---

kernel.c

Go to the documentation of this file.

#include <math.h>
#include <vector>
//#include <unistd.h> // Temporary; for sync()

#include "ind_types.h"
#include "ipc.h"
#include "cmdparam.h"
#include "binarysave.h"
#include "shuffledrand.h"
#include "globals.h"
#include "kernel.h"
#include "kernelwrapper.h"
#include "stringutils.h"
#include "neuralregion.h"

/* Currently multiprocessor support is unimplemented; the calls to
   ipc_put and ipc_get use dynamically-allocated data, whose address
   is not known on remote PEs.  They would need to be changed to
   broadcast addresses into known locations before broadcasting
   the data. */
#if (NPES>1)
#error "Multiprocessor support is not implemented for kernel.c"
#endif



/******************************************************************************/
/* Defines and typedefs                                                       */
/******************************************************************************/

/*  For speed, define NUM_EYES_IS_CONSTANT=<num_eyes> to declare that
    num_eyes is constant and equal to the given value.  Skip this if
    flexibility is more important than speed.  However, if you can
    make this parameter constant, it speeds up certain inner loops
    considerably, because automatic loop unrolling can then eliminate
    the entire loop in many cases.  Historically, the speedup has been
    about 10-15%.  It is your responsibility to set num_eyes to match,
    although warnings should be generated if it does not.
*/
/* #define NUM_EYES_IS_CONSTANT=1 */



#define ACTLIST_START(actlist) 1
#define ACTLIST_COUNT(actlist) ((actlist)[0])

/* 12-09-95 Added by J.A.B. -- Reload from file immediately after save.
   Allows ASCII saving in middle of training by replacing the weights
   that were destroyed, or can be used for debugging.
*/
/* #define RELOAD_AFTER_WEIGHT_SAVE */

#define VERIFY_WEIGHT_SAVES_ERROR_LIMIT 10




/* Macros to assist in handling circular receptive fields                      */

#define WITHIN_RADIUS(xdistance,ydistance,radius_sq) \
(((xdistance)*(xdistance) + (ydistance)*(ydistance)) <= (radius_sq))

#define WITHIN_RADIUS_SQ(xdistance_sq,ydistance_sq,radius_sq) \
((xdistance_sq) + (ydistance_sq) <= (radius_sq))



/******************************************************************************/
/* LissomMap class                                                            */
/******************************************************************************/

class LissomMap : public BasicLissomMap, public InternalNeuralRegion {
  typedef std::vector<int> actlist;
  
public:
  LissomMap(const string& name_i, Dimensions dims, 
            bool init_weights_, ParamMap* params_);

  virtual ~LissomMap() {}


  virtual void activate(bool learn=false, bool settle=true, bool activatefn=true);

  /* For OBJ-LISSOM, but also useful for the backproject command */
  virtual void backproject();
  virtual void backproject(const string& name, const double gammaaff=1.0);
#ifdef OBJ_LISSOM  
  virtual void resp_to_residual();
  virtual void resp_to_residual(const string& name, const double gammaaff=1.0);
  virtual void train_obj_weights();
#endif /* OBJ_LISSOM */
  


  virtual const WeightMatrix get_weights(const string& name, int ui=0, int uj=0) const;
  virtual NeuralRegion* get_input(const string& name);
  virtual void prune();
  virtual void save_state(const string&);
  virtual bool restore_state(const int, const string&);
  virtual double size_connection_bytes()   const;
  virtual double size_unique_connections() const;
  virtual void   add_input(const string& name_, NeuralRegion& upstream_region,
                           WeightFunction& fn, Length size_scale);
  virtual Dimensions input_dimensions(WeightFunction& fn, Length size_scale=1.0);
  virtual WeightBounds* get_weights_bounds(const string& name, int ui=0, int uj=0) const;
  virtual int    num_aff_inputs() const
    {  return aff_inputs.size();  }

protected:
  void   present_inputs(bool learn, bool settle=true, bool activate=true);
  void   init_weights(bool init_weights);
  void   collect_activation_data(int dest_pe=Uninitialized);
  double dist_lat_wt(double value, double amt_of_randomness);
  double lat_resp(int ui, int uj, int radius, int array_width, l_weight *weights, int weight_idx_bound);
  double lat_resp_rad1(int ui, int uj, int radius, int array_width, int array_width_2, l_weight *weights, int weight_idx_bound);
  double lat_resp_named(int ui, int uj, const string& name);
  double sigmoid(double activity, double sigdelta, double sigbeta);
  double truesigmoid(double activity, double beta);
  double truesigmoidprime(double sigmoidactivity, double beta);
  double positivesigmoid(double activity, double beta);
  double logsigmoid(double activity, double beta);
  double logsigmoidprime(double activity, double beta);
  int    binary_search(const actlist& list, int lower, int upper, int value);
  double change_lateral_exc_radius(double old_radius, double new_radius);
  void   clear_weights(void);
  void   randomize_lat_wts(void);
  void   compute_responses(int ts, bool activate, bool continuousdynamics);
  inline  void   crop_actlist(const actlist&, int lowval, int highval, int *lowidx, int *highidx);
  void   initialize_actlists(void);
  void   initialize_markers(void);
  void   modify_all_weights(void);
  void   modify_lat_wts(double alpha_exc, double alpha_inh);
  void   modify_latwt_loop(int i, int j, int radius, int array_width, l_weight *weights, int weight_idx_bound, double alpha);
  void   modify_latwt_loop_killed(int i, int j, int radius, int array_width, l_weight *weights, int weight_idx_bound, double alpha);
  void   modify_latwt_loop_rad1(int ui, int uj, int radius, int array_width, int array_width_2, l_weight *weights, int weight_idx_bound, double alpha);
  void   modify_weights(double alpha, bool normalize);
  void   normalize_afferent_wts(int lowk, int highk, int lowl, int highl, int row, int j, double length);
  void   normalize_lat_wts(int i, int j, int radius, int array_width, l_weight *weights);
  double prune_circular_aff_weights(int ui, int uj, int radius, AfferentWeightsList& weights,
                                    bool smooth, double radius_trim, double smooth_trim,
                                    double length);
  void   prune_circular_lat_weights(int ui, int uj, int radius_sq, int array_width, l_weight *weights,
                                    bool smooth, double radius_trim, double smooth_trim);
  void   receptor_normalize(void); 
  void   reduce_lateral_radius(int ui, int uj, double old_radius, double new_radius, int array_width, l_weight *weights,
                               bool circular, bool smooth, double smooth_trim);
  void   response_to_input(bool continuousdynamics);
  void   set_markers(int radius); 
  void   settle_responses(bool learn, int tsettle, bool activate, bool continuousdynamics);
  void   setup_afferent_weights(Neuron& neuron, AfferentWeightsList& awts,
                                int ui, int uj, int radius, double radius_sq,
                                int preset, double sigma, double amt_of_randomness,
                                bool circular, bool smooth, 
                                double radius_trim, double smooth_trim,
                                double force_cx=-1, double force_cy=-1);
  void   setup_lateral_weights(int ui, int uj, int radius,
                               int array_width, l_weight *weights,
                               int preset, double sigma, double amt_of_randomness,
                               bool circular, bool smooth, double radius_trim,
                               double smooth_trim, int* con_killed_flag);
//void   verify_actlist(const actlist& list);
//void   verify_actlist_array(const actlist& lists[],int num);
  void   verify_saved_lateral_weights( int i, int j, int radius, int array_width, l_weight *old_weights, l_weight *new_weights, double tolerance, const char *description);
  void   verify_saved_weights( double tolerance );


  SETFNOBJ_DECLARE2(LissomMap,exc_rad_setfn2);

  OwningPointer<ParamMap> params;

#ifndef NUM_EYES_IS_CONSTANT
  const int num_aff;
#endif
  int    exc_connections_killed;          
  int    inh_connections_killed;          
  const int exc_array_width_2;            
  const int RN_sq;                        
  ActivityMatrix gather_activity;         
  ActivityMatrix init_activity;           
  ActivityMatrix map_activity;            
  ActivityMatrix resp_to_inp;             
  ActivityMatrix prev_map_activity;       
  /* Using a vector of vectors results in a slight performance
     decrease compared to an array of vectors (approximately 5%
     overall for small (30MB) networks), but is significantly easier
     to manage.  If the speedup is desired, arrays could be allocated
     for these variables instead.
  */
  typedef std::vector<actlist> ActivityListTable;
  ActivityListTable activity_list;        
  ActivityListTable prev_activity_list;   
  // Should be <bool> but for some reason that doesn't appear to work well with MTL (using mtl-2.1.2-15)
  typedef MatrixType<int>::rectangular  MarkerType;
  MarkerType activity_marker;             
  typedef std::vector<double> SumType;
  SumType aff_norm_sum;                   
  SumType afferent_sum;                   
  SumType temp_sum;                       
  typedef NeuralRegion::ActivityMatrix        aff_input_type;
  typedef std::vector<const aff_input_type*>  aff_inputs_type;
  aff_inputs_type aff_inputs;
  std::vector<string> aff_names;

  typedef Boundary::BoundingCircle<int,int> BoundingCircle;
  typedef Boundary::AARBoundingBox<int,int> BoundingRectangle;

  int array_radius_int(const ParamMap& parammap, const string& radius_param) const
    {  return array_radius_int(parammap.get_with_default(radius_param.c_str(),0.0));  }

  int array_radius_int(const double radius_float) const
    //Probably would be better, but would require extra work to support:
    //{  return Generic::round(radius_float);  }
    {  return int(radius_float);  }



/* yschoe: the following is an (almost) literal copy from fixedwtregion */
protected: 

  struct Input {
    Input(const string& name_i, const WeightMatrix& kernel_i,
          NeuralRegion* input_i, const Length size_scale_i)
      : name_str(name_i), kernel(kernel_i), input(input_i), size_scale(size_scale_i) { }
    const string& name() const {  return name_str;  }

    const double size_bytes() const {  return mat::size(kernel)*sizeof(WeightMatrix::value_type);  }
    const double size_conns() const {  return mat::size(kernel);  }

    string name_str;              
    WeightMatrix kernel;          
    NeuralRegion* input;    
    Length size_scale;            
  };

  inline Subscript input_row(const Input& input, Subscript row) const
    {  return Subscript(row*(input.size_scale)+(input.kernel.nrows()/2.0));  }

  inline Subscript input_col(const Input& input, Subscript col) const
    {  return Subscript(col*(input.size_scale)+(input.kernel.ncols()/2.0));  }

  typedef std::vector<Input*> inputs_type;

  const Activity xo,yo; 
  inputs_type            inputs;

};



LissomMap::LissomMap(const string& name_i, Dimensions dims,
                     bool init_weights_, ParamMap* params_) :
  BasicLissomMap(dims.height,dims.width,
                 params_->get_with_default("num_aff_inputs",0),
                 params_->get_with_default("RN",0),
                 array_radius_int(*params_,"rf_radius"),
                 array_radius_int(*params_,"exc_rad"  ),
                 array_radius_int(*params_,"inh_rad"  ),
                 NPEs),
  InternalNeuralRegion(name_i,dims),
  params(params_,true),
#ifndef NUM_EYES_IS_CONSTANT
  num_aff(params_->get_with_default("num_aff_inputs",0)),
#else
#define num_aff NUM_EYES_IS_CONSTANT
#endif
  exc_connections_killed(False), inh_connections_killed(False),
    
  /* The array widths set below don't change even if exc_radius is decreased. */
  exc_array_width_2(2*exc_array_width),
  RN_sq(RN*RN),
  
  // Should these arguments be swapped?  This order is correct
  // for a matrix, but does not match the swapped conventions that
  // this file is using.
  gather_activity(dims.height+1,dims.width+1), // Was [NMAX][NMAX+1]; just being safe here
  init_activity(dims.height,dims.width),
  map_activity(dims.height,dims.width),
  resp_to_inp(dims.height,dims.width),
  prev_map_activity(dims.height,dims.width),
  
  activity_list(perows,actlist(N+1)),
  prev_activity_list(N,actlist(N+1)),
    
  activity_marker(dims.height,dims.width),
  
  aff_norm_sum(num_aff*RN*RN),
  afferent_sum(num_aff*RN*RN),
  temp_sum(num_aff*RN*RN)
  
{
  /* Make sure that our exc_rad will be changed when someone sets it in
     our parameter set */
  ParameterMap<>* pm = dynamic_cast< ParameterMap<>* >(params.valueptr());
  pm->set_local("exc_rad",params->get_with_default("exc_rad",0.0));
  pm->lookup("exc_rad").add_setfn(new setfnobj2_exc_rad_setfn2(this));
  
  /* Initialize activity in the map */
  for(int i=0; i< N; i++)   
    for(int j=0; j< N; j++)         
        prev_map_activity[i][j]=map_activity[i][j]=0.0;

  ipc_notify(IPC_ONE,IPC_SUMMARY,
             "%s %d-connection (%dMB%s) network",
             (init_weights_? "Initializing" : "Declaring"),
             int(size_unique_connections()),(int)(size_connection_bytes()/1024.0/1024.0+1),
             (NPEs>1 ? "/PE" : ""));
  
  init_weights(init_weights_);

  /* Register map for InputResponse */
  typedef mat::CArrayMatrixAdapter<double,ActivityMatrix>                       DoubleMatWrapper;
  typedef mat::MatrixInterfaceAdapter<DoubleMatWrapper,NeuralRegion::Magnitude> DoubleMatAdapter;
  DoubleMatWrapper* inpresp = new DoubleMatWrapper(init_activity, (unsigned int)N, (unsigned int)N);
  table_set("InputResponse",  new DoubleMatAdapter(inpresp)); // Swaps coord order for plotting
}



NeuralRegion* LissomMap_create(const string& name_i, NeuralRegion::Dimensions dims,
                               bool init_weights_, ParamMap* params)
{  return new LissomMap(name_i,dims,init_weights_,params);  }



cmdstat LissomMap::exc_rad_setfn2(const string& name, Typeless& existing, const Typeless& requested)
{
#ifdef NO_WEIGHTS
  (void)name;
  existing.set(requested);
  
#else 
  double oldvalue       = (double)existing.numericvalue();
  double requestedvalue = (double)requested.numericvalue();

  if (!network_initialized)
    existing.set(requested);
  else if (oldvalue != requestedvalue) {
    ipc_notify(IPC_ONE,IPC_STD,"Renormalizing %s's weights to change %s from %g to %g",
               LissomMap::name().c_str(),name.c_str(),double(oldvalue),double(requestedvalue));
    double actualvalue = change_lateral_exc_radius(oldvalue, requestedvalue);
    existing.set(Polymorph<double>(actualvalue));
  }
  exc_rad=array_radius_int(*params,"exc_rad");
#endif
  
  return CMD_NO_ERROR;
}



void LissomMap::activate(bool learn, bool settle, bool activatefn)
{
  collect_activation_data();
  present_inputs(learn,settle,activatefn);

  /* Copy the activity to our local matrix to allow clean plotters to be developed */
  /* Swaps the index order to match that assumed by ppm_draw.c */
  for (Subscript i=0; i<output.nrows(); i++)
    for (Subscript j=0; j<output.ncols(); j++)
      output[i][j] = bp[i][j] = prev_map_activity[j][i];
  //mat::gnuplot(output, name());

}
  

#define assert_equals(x,y) \
if (x!=y) {  failed=true; ipc_notify(IPC_ALL,IPC_ERROR,"Assertion failed: " #x "!=" #y " (%d!=%d)",x,y);  }
#define assert_ge(x,y) \
if (!(x>=y)) {  failed=true; ipc_notify(IPC_ALL,IPC_ERROR,"Assertion failed: !(" #x "=>" #y ") !(%d=>%d)",x,y);  }

/* This implementation is essentially just a stub, but it does check that
   the values match what we already assume anyway */
void LissomMap::add_input(const string& name_, NeuralRegion& source, WeightFunction& fn, Length size_scale)
{


  /* yschoe: fixedwtregion style add_input to retain vector of upstream 
     regions (source): copied from fixedwtregion.h */
  WeightMatrix kernel = (fn)(size_scale);
  /* Kernel height and width must be odd */
  assert(kernel.nrows()%2==1);
  assert(kernel.ncols()%2==1);
    
  Generic::insert_named(inputs,name_,new Input(name_,kernel,&source,size_scale));
  /* end of fixedwtregion style add_input */



  const int      afferent_edge_buffer     = params->get_with_default("retina_edge_buffer",0);

  bool failed=false;
  if (String::non_numeric_basename_matches<string>(name_,"Afferent")) {
    WeightMatrix sampleweights = (fn)(size_scale);
    assert_equals(int(sampleweights.nrows()-1)/2,rf_radius);
    assert_equals(int(sampleweights.ncols()-1)/2,rf_radius);
    assert_equals(int(source.const_activity().nrows()),RN);
    assert_equals(int(source.const_activity().ncols()),RN);
    const int calculated_N = (int(size_scale*(RN-2*afferent_edge_buffer)+0.0001));
    if (calculated_N!=N)
      ipc_notify(IPC_ALL,IPC_WARNING,"Mismatched size_scale: %d!=%d",calculated_N,N);
    
    /* Does one useful thing: adds this input to the list of afferent regions. */
    /* Temporarily uses the number being zero as the clue that we
       need to replace the list of inputs with a new one. */
    if (!failed) {
      const int aff   = String::numeric_extension<int>(name_);
      if (aff==0) {  aff_inputs.resize(0);  aff_names.resize(0);  }
      aff_inputs.push_back(&(source.const_activity()));
      aff_names.push_back(name_);
    }
  }
  else if (name_ == "LateralExcitatory") {
    WeightMatrix sampleweights = (fn)(size_scale);
    assert_equals(int(sampleweights.nrows()-1)/2,exc_rad);
    assert_equals(int(sampleweights.ncols()-1)/2,exc_rad);
    assert_equals(int(source.const_activity().nrows()),N);
    assert_equals(int(source.const_activity().ncols()),N);
    assert_equals(int(100*size_scale+0.0001),100); // size_scale==1.0
  }
  else if (name_ == "LateralInhibitory") {
    WeightMatrix sampleweights = (fn)(size_scale);
    assert_equals(int(sampleweights.nrows()-1)/2,inh_rad);
    assert_equals(int(sampleweights.ncols()-1)/2,inh_rad);
    assert_equals(int(source.const_activity().nrows()),N);
    assert_equals(int(source.const_activity().ncols()),N);
    assert_equals(int(100*size_scale+0.0001),100); // size_scale==1.0
  }
  else {
    ipc_notify(IPC_ALL,IPC_ERROR,"Could not add unsupported weight type (%s)",name_.c_str());    
  }

  if (failed)
    ipc_notify(IPC_ALL,IPC_ERROR,"Unable to add connection %s",name_.c_str());
}


NeuralRegion::Dimensions LissomMap::input_dimensions(WeightFunction&, Length size_scale)
{
  const int      afferent_edge_buffer     = params->get_with_default("retina_edge_buffer",0);
  const int calculated_N = (int(size_scale*(RN-2*afferent_edge_buffer)+0.0001));
  if (calculated_N!=N)
    ipc_notify(IPC_ALL,IPC_WARNING,"Mismatched size_scale: %d!=%d",calculated_N,N);
  
  return NeuralRegion::Dimensions(Subscript(RN),Subscript(RN));
}


const LissomMap::WeightMatrix LissomMap::get_weights(const string& name, int ui, int uj) const
{
  bool failed=false;

  const bool   allow_negative_wts = params->get_with_default("allow_negative_wts",False);

  /* Check assumptions */
  assert_ge(ui,0);
  assert_ge(uj,0);
  assert_ge(N,ui);
  assert_ge(N,uj);

  // Multiprocessor note:
  /*
    At present plotting weights on multiprocessor machines requires that
    the plot_pe correctly specify the particular PE on which the weights
    reside, which is not something that the user is ordinarily aware of.
    Instead we should just fetch the weights from the PE where they
    are stored and return them.
  */
  assert_equals(ROWISLOCAL(ui),true);

  if (failed) return WeightMatrix();
  
  /* The array indeces are swapped to translate between the correct and the legacy formats. */
  const Neuron&              neuron = cortex_map[uj][ui];
  const AfferentWeightsList& awts   = *(neuron.weights);
  //const AfferentWeightsList& awts   = reference_wts.weights;

  // pre-allocate these matrices so that they do not get allocated
  // subsequently.
  /*
  if (aff_w == NULL) 
        aff_w = new WeightMatrix(2*rf_radius+1,2*rf_radius+1);
  
  if (lat_w == NULL) 
        lat_w = new WeightMatrix(2*MAX(inh_rad,exc_rad)+1,2*MAX(inh_rad,exc_rad)+1);
  */

  /* Select appropriate set of weights to copy */
  if (String::non_numeric_basename_matches<string>(name,"Afferent")) {
    const int aff   = std::find(ISEQ(aff_names),name) - aff_names.begin();
    assert (aff<num_aff);

    const int lowk  = MAX(neuron.centery-rf_radius,0   );
    const int highk = MIN(neuron.centery+rf_radius,RN-1);
    const int lowl  = MAX(neuron.centerx-rf_radius,0   );
    const int highl = MIN(neuron.centerx+rf_radius,RN-1);
    
    WeightMatrix w(highk-lowk+1,highl-lowl+1);
    // WeightMatrix& w=(*aff_w);
    
    for (int k=lowk; k <=highk; k++)
      for (int l=lowl; l<=highl; l++) {
        a_weight aw = awts[aff][l-lowl][k-lowk];
        w[k-lowk][l-lowl] = ((aw>(-DEATH_FLAG) || allow_negative_wts)?aw:0.0);
      }
    
    return w;
  }
  else if (name == "LateralExcitatory") {
    const int lowi  = MAX(uj-exc_rad,0);
    const int highi = MIN(uj+exc_rad,N-1);
    const int lowj  = MAX(ui-exc_rad,0);
    const int highj = MIN(ui+exc_rad,N-1);
    
    WeightMatrix w(highj-lowj+1,highi-lowi+1);
    // WeightMatrix& w=(*lat_w);
    
    for (int i=lowi; i<=highi; i++)
      for (int j=lowj; j<=highj; j++) {
        l_weight lw = neuron.lat_exc_wts[LAT_INDEX(uj,ui,i,j,
                                        exc_rad,exc_array_width)]; 
        w[j-lowj][i-lowi] = ((lw>(-DEATH_FLAG) || allow_negative_wts)?lw:0.0);
      }
    
    return w;
  }
  else if (name == "LateralInhibitory") {
    const int lowi  = MAX(uj-inh_rad,0);
    const int highi = MIN(uj+inh_rad,N-1);
    const int lowj  = MAX(ui-inh_rad,0);
    const int highj = MIN(ui+inh_rad,N-1);
    
    WeightMatrix w(highj-lowj+1,highi-lowi+1);
    // WeightMatrix& w=(*lat_w);
    
    for (int i=lowi; i<=highi; i++)
      for (int j=lowj; j<=highj; j++) {
        l_weight lw = neuron.lat_inh_wts[LAT_INDEX(uj,ui,i,j,
                                        inh_rad,inh_array_width)]; 
        w[j-lowj][i-lowi] = ((lw>(-DEATH_FLAG) || allow_negative_wts)?lw:0.0);
      }
      
    return w;
  }

  /* Unknown set of weights */
  return WeightMatrix();
}



NeuralRegion* LissomMap::get_input(const string& name)
{
  Input* i = Generic::find_named<Input>(inputs,name);

  return (i? i->input : 0);

}



LissomMap::WeightBounds* LissomMap::get_weights_bounds(const string& name, int ui, int uj) const
{

  /* The array indeces are swapped to translate between the correct and the legacy formats. */
  int ci=ui;
  int cj=uj;
  string radparam = "rf_radius";
  NeuralRegion::Bounds b=bounds();
  const bool     circular_aff_wts    = params->get_with_default("circular_aff_wts",True);
  const bool     circular_lat_wts    = params->get_with_default("circular_lat_wts",True);
  bool circular=circular_lat_wts;
  bool failed=false;

  /* Check assumptions */
  assert_ge(ui,0);
  assert_ge(uj,0);
  assert_ge(N,ui);
  assert_ge(N,uj);
  assert_equals(ROWISLOCAL(ui),true);
  if (failed) return new NeuralRegion::Bounds();
  
  if      (name=="LateralExcitatory") radparam = "exc_rad";
  else if (name=="LateralInhibitory") radparam = "inh_rad";
  else if (String::non_numeric_basename_matches(name,string("Afferent"))) {
    circular=circular_aff_wts;
    radparam = "rf_radius";
    const Neuron& neuron = cortex_map[MAPROW(uj)][ui];    
    ci = neuron.centery;
    cj = neuron.centerx;
    b  = Bounds(0,0,RN-1,RN-1);
  }
  else ipc_notify(IPC_ALL,IPC_ERROR,"Don't know weight type: %s",name.c_str());

  if (circular) {
    const double radius = params->get_with_default(radparam.c_str(),0.0);
    const BoundingCircle e(ci,cj,radius);
    return new Boundary::Intersection<int,int,double,BoundingCircle,Bounds>(e,b);
  }
  else {
    const int rad = array_radius_int(*params,radparam);
    const BoundingRectangle e(ci-rad,cj-rad,ci+rad,cj+rad);
    return new Boundary::Intersection<int,int,double,BoundingRectangle,Bounds>(e,b);
  }
}


/******************************************************************************/
/* Core routines                                                              */
/******************************************************************************/

void LissomMap::present_inputs(bool learn, bool settle, bool activatefn)
{
  int i,pe,j,count;
  
  const bool continuousdynamics = params->get_with_default("dynamics",False);
  const int  tsettle            = params->get_with_default("tsettle",1);
  
  /* Present and settle activity into bubbles */
  settle_responses(learn,(settle? tsettle : 1),activatefn,continuousdynamics);
  
  for(pe=0; pe< NPEs; pe++)       /* Collect map_activity from each processor */
    for(i=0; i<perows; i++)
      ipc_get_to(&(prev_map_activity[ARBITRARY_MAPROW(i,pe)][0]),
                 &(map_activity[ARBITRARY_MAPROW(i,pe)][0]),IPC_DOUBLE,N,pe);
  
  /* gather_activity is indexed from [1..count] like the activity list */
  for (i=0; i<N; i++)
    if ((count=ACTLIST_COUNT(prev_activity_list[i])) > 0)
      for (j=ACTLIST_START(prev_activity_list); j<=count; j++)
        gather_activity[i][j]=prev_map_activity[i][prev_activity_list[i][j]];
  
  if (learn && !continuousdynamics)      /* Otherwise, modify at every step */
    modify_all_weights();
}



void LissomMap::settle_responses(bool learn, int tsettle, bool activatefn, bool continuousdynamics)
{
  int ts;
  const bool   allow_negative_wts = params->get_with_default("allow_negative_wts",False);

  response_to_input(continuousdynamics);
  initialize_markers();
  initialize_actlists();

  ipc_barrier(); /* To make sure dynamics = 0 init has taken place */

  for(ts=1; ts<=tsettle; ts++){ /* Let network settle into an activity bubble */
    compute_responses(ts,activatefn,continuousdynamics);  /* Add up input and lateral activations */
    ipc_barrier();          /* To ensure all activations have been computed  */
    
    if (learn && continuousdynamics) { /* Modify at each step   */
      (void)allow_negative_wts;
      assert(!allow_negative_wts);
      modify_all_weights();
    }
  }
}



void LissomMap::response_to_input(bool continuousdynamics)
{
  const bool   divide_input_sum      = params->get_with_default("divide_input_sum",False);
  const bool   divide_reference_sum  = params->get_with_default("divide_reference_sum",False);
  const double inp_sum_scale         = params->get_with_default("divide_input_sum_scale",1.0);
  const double inp_sum_offset        = params->get_with_default("divide_input_sum_offset",0.0);
  const double gammasur              = params->get_with_default("gammasur",1.0);
#ifdef OBJ_LISSOM
  const bool   l2_norm  = params->get_with_default("l2_norm",False);
#endif /* OBJ_LISSOM */

  if (int(aff_inputs.size()) != num_aff) {
    ipc_notify(IPC_ALL,IPC_ERROR,"No response; the number of afferent inputs was declared as %d but is %d (see num_aff_inputs)",
               num_aff, int(aff_inputs.size()));
    return;
  }
      
  for(int i=0; i<perows; i++){
    const int row = MAPROW(i);                      /* calculate row index  */
    for(int j=0; j<N; j++)                              /* Init all activations */
      { 
        const Neuron&              neuron = cortex_map[row][j];
        const AfferentWeightsList& awts   = *(neuron.weights);
        const AfferentWeightsList& ref_awts = reference_wts.weights;

        const int lowk   = MAX(neuron.centerx-rf_radius,0   );
        const int highk  = MIN(neuron.centerx+rf_radius,RN-1);
        const int lowl   = MAX(neuron.centery-rf_radius,0   );
        const int highl  = MIN(neuron.centery+rf_radius,RN-1);

        const int rlowk  = neuron.centerx-rf_radius;
        const int rlowl  = neuron.centery-rf_radius;

        double resp =0;
        double resp1=0;
        double inp  =0;
        double inp1 =0;
        double ref  =0;
        double ref1 =0;
        double num_inp =0;
        
        /* Blank out previous activity in network if appropriate. */
        if (!continuousdynamics) prev_map_activity[row][j] = map_activity[row][j] =0.0;

#ifdef OBJ_LISSOM
        /* calculate l2_norm if necessary */
        if (l2_norm) {
          for (inputs_type::iterator i=inputs.begin();i != inputs.end(); ++i) {
           if (String::non_numeric_basename_matches<string>((*i)->name(),"Afferent")){
            NeuralRegion* inp = dynamic_cast<NeuralRegion*>(get_input((*i)->name()));
            ActivityMatrix& in = inp->output;

            int nr = inp->output.nrows(), nc = inp->output.ncols();
            double l2_sum=0.0;
            /* gather sum sqr */
            for (int k=0; k<nr ; ++k) 
              for (int l=0; l<nc ; ++l) {
                  l2_sum += in[k][l]*in[k][l];
              }

            l2_sum = 1.0/sqrt(l2_sum);

            /* normalize input */
            for (int k=0; k<nr ; ++k) 
              for (int l=0; l<nc ; ++l) {
                  in[k][l] *= l2_sum;
              }

           }
          }
        }
#endif /* OBJ_LISSOM */ 
        /* This is a time-critical loop; make sure that the compiler unrolls it */

        

#ifndef NUM_EYES_IS_CONSTANT
        /* Use canonical but slower version of algorithm */
        
#ifdef SPARSE_MATRICES
        /* Basic version */
        for (int k=lowk; k <=highk; k++) {
          for (int l=lowl; l<=highl; l++)
            for (int aff=0; aff<num_aff; aff++) {
              const aff_input_type::value_type& in = (*aff_inputs[aff])[l][k];
              resp  += in * awts[aff][k-lowk][l-lowl];
              ref   += in * ref_awts[aff][k-rlowk][l-rlowl];
              inp   += in;
            }
        }
        
#else     
        /* Optimized assuming dense rectangular layout */
        const aff_inputs_type::value_type*     ap = &(aff_inputs[0]);
        const AfferentWeightsList::value_type* wp = &(awts[0]);
        const AfferentWeightsList::value_type* rwp= &(ref_awts[0]);
        for (int aff=0; aff<num_aff; aff++,ap++,wp++,rwp++) {
          for (int k=lowk; k <=highk; k++) {
            /* This could be further optimized if both arrays had the same index order... */
            //const aff_input_type::value_type*  al=&((**ap)[l][lowk]);
            const AfferentWeights::value_type* wl=&((*wp)[k-lowk][0]);
            const AfferentWeights::value_type* rwl=&((*rwp)[k-rlowk][lowl-rlowl]);
            for (int l=lowl; l<=highl; l++,wl++,rwl++) {
              const aff_input_type::value_type& in = (**ap)[l][k];
              //resp  += (*al) * (*wl);
              resp  += in * (*wl);
              ref   += in * (*rwl);
              inp   += in;
              num_inp++;
            }
          }
        }
#endif

        
#else /* num_aff is constant */
#error divide_input_sum and divide_reference_sum have not yet been implemented when NUM_EYES_IS_CONSTANT
        // (Implementing them is trivial (see above) but has not been done.)
        
        for (int k=lowk; k <=highk; k++){

          /*
            This is an extremely time-critical loop, so it has been optimized
            by unrolling the inner loop by hand for the typical case of
            num_aff==1 and num_aff==2.
          */
          for (int l=lowl; l<=highl; l++){
            // Would need to update inp and ref here...
            resp  += (*aff_inputs[0])[l][k] * awts[0][k-lowk][l-lowl];
#if   (num_aff>1)
            resp1 += (*aff_inputs[1])[l][k] * awts[1][k-lowk][l-lowl];
#endif
#if   (num_aff>2)
            {
              for (int aff=2; aff<num_aff; aff++)
                resp  += (*aff_inputs[aff])[l][k] * awts[aff][k-lowk][l-lowl];
            }
#endif    
          }
        }
#endif /* num_aff is constant */

        const double response  = resp + resp1;
        const double reference = ref  + ref1;
        /* Effectively normalizes input_sum weights to 1.0, all equal */
        const double inp_avg  = (num_inp? (inp + inp1)/num_inp : 0.0);
        const double inp_resp = 
          (divide_input_sum?      (inp_avg==0   ? response : response/(inp_sum_scale*inp_avg+inp_sum_offset)) : 
           (divide_reference_sum? (reference==0 ? response : response/(inp_sum_scale*ref    +inp_sum_offset)) :
            response));

        resp_to_inp[row][j] = inp_resp - gammasur * inp_avg;
      }
  }
}



void LissomMap::backproject()
{
  // NeuralRegion* inp = get_input("Afferent00");
  // cout << "- retrieved upstream region: " << inp->name() << endl;

  const double gammaaff = params->get_with_default("gammaaff",1.0);
  const double beta = params->get_with_default("beta",1.0);
  // Hack; should pass verbose through backproject() instead
  const bool   verbose = params->get_with_default("::cmd::backproject::verbose",False);

  // set up current bp matrix by passing it through activation function.
  // we are assuming that bp[][] in its final stage did not go through
  // activation function: this applies to both present_input() and 
  // also to resp_to_residual().
  int nr = bp.nrows(), nc=bp.ncols();
  for (int x=0; x<nr; ++x) 
     for (int y=0; y<nc; ++y) {
        // bp[x][y] = logsigmoid(bp[x][y],beta);
        bp[x][y] = truesigmoid(bp[x][y],beta);
        // These functions are not available on all platforms
        //assert(!isnan(bp[x][y]));
        //assert(!isinf(bp[x][y]));
     }

  // clean upstream bp and residual matrices
  for (inputs_type::iterator i=inputs.begin();i != inputs.end(); ++i) {
    if (String::non_numeric_basename_matches<string>((*i)->name(),"Afferent")){
      NeuralRegion* inp = dynamic_cast<NeuralRegion*>(get_input((*i)->name()));
      int nr = inp->bp.nrows(), nc = inp->bp.ncols();
      ActivityMatrix& r = inp->residual;
      ActivityMatrix& b = inp->bp;
      for (int x=0; x<nr; ++x) 
        for (int y=0; y<nc; ++y) {
          // do not touch inp->output[][] !
          b[x][y] = 0.0;
          r[x][y] = 0.0;
        }
    }
  }

  // backproject current bp to the upstream bp matrix (using addition)
  // and update residual
  for (inputs_type::iterator i=inputs.begin();i != inputs.end(); ++i) {
    if (String::non_numeric_basename_matches<string>((*i)->name(),"Afferent")){

      backproject((*i)->name(),gammaaff);

      NeuralRegion* inp = dynamic_cast<NeuralRegion*>(get_input((*i)->name()));

      double err = 0.0, err_tmp=0.0, bp_tmp=0.0, output_tmp=0.0;
      double v1_len = 0.0, v2_len = 0.0, dot_p = 0.0, vector_err = 0.0;
      int nr = inp->bp.nrows(), nc = inp->bp.ncols();

      ActivityMatrix& r = inp->residual;
      ActivityMatrix& b = inp->bp;
      ActivityMatrix& o = inp->output;
      for (int x=0; x<nr; ++x) 
        for (int y=0; y<nc; ++y) {
            bp_tmp      = b[x][y];
            output_tmp  = o[x][y];

            // sqrt(sum sqr err)
            err_tmp     = output_tmp - bp_tmp;
            err                += err_tmp * err_tmp;
            r[x][y]     += err_tmp;

            // vector error
            dot_p += bp_tmp * output_tmp;
            v1_len += bp_tmp*bp_tmp;
            v2_len += output_tmp*output_tmp;
        }

        vector_err = acos(dot_p/(sqrt(v1_len * v2_len)));

        if (verbose)
          ipc_notify(IPC_ALL,IPC_STD,"%s: sqrt(sum err sqr) = %f : vector error = %f\n",
                     ((*i)->name()).c_str(),sqrt(err), vector_err);
        else 
          ipc_log(IPC_ONE,"%s: sqrt(sum err sqr) = %f : vector error = %f\n",
                  ((*i)->name()).c_str(),sqrt(err),vector_err);
        //sync();
    }
  }

}


#ifdef OBJ_LISSOM

void LissomMap::resp_to_residual()
{

  // How the various matrices are used: the usage is pretty
  //    brittle and depends on a lot of conditions outside
  //    of this function: e.g., backproject() being called immediately
  //    before the call to the current function.
  //
  // precondition:
  //    output:                 a
  //    bp:                     f(a)    (backproj() should have applied f())
  //    residual (upstream):    r       (calculated by backproject)
  //    * backproject() should have been called prior to this function call!
  //    * repeated calls of backproject() and resp_to_residual() will
  //      update the activity 'a' (inner loop optimization).
  // 
  // intermediate:
  //    output:                 r*W
  //    resp_to_inp:            a
  //    prev_map_activity:      f(a)
  //    map_activity:           f'(a)
  //    residual (current):     -dE/da = alpha_settle * f'(a)[r*W+exc-inh]
  //
  // postcondition:
  //    output = prev_map_activity:     new a
  //    bp:                             new a (f will be applied by backproj)
  //
  // As a result, output[][] and bp[][] should contain the new a_i.
  // 
  //    * weight update need to be based on the new 'f(a)' values
  //      and the to-be-calculated residual 'r' based on the current 'a'. 
  //      All these necessary steps can be taken care of by backproject() 
  //      which will put 'f(a)' in bp[][] and 'r' in the upstream residual[][].

  const double gammaaff = params->get_with_default("gammaaff",1.0);
  const double gammaexc = params->get_with_default("gammaexc",1.0);
  const double gammainh = params->get_with_default("gammainh",1.0);
  const double alpha_settle = params->get_with_default("alpha_settle",1.0);
  const double beta = params->get_with_default("beta",1.0);
  const bool obj_lateral_on = params->get_with_default("obj_lateral_on",True);

  // Initialize.
  int nr = output.nrows(), nc = output.ncols();
  for (int x=0; x<nr; ++x) 
    for (int y=0; y<nc; ++y) {
      resp_to_inp[y][x] = output[x][y];
      // prev_map_activity[y][x] = logsigmoid(resp_to_inp[y][x],beta);
      // map_activity[y][x] = logsigmoidprime(resp_to_inp[y][x],beta);
      prev_map_activity[y][x] = truesigmoid(resp_to_inp[y][x],beta);
      map_activity[y][x] = truesigmoidprime(prev_map_activity[y][x],beta);
      output[x][y] = 0.0;
    } // a, f(a), and f'(a) prepared. output (for r*W) initialized to 0.

  // calculate response to residual: r*W
  for (inputs_type::iterator i=inputs.begin();i != inputs.end(); ++i) {
    if (String::non_numeric_basename_matches<string>((*i)->name(),"Afferent")){

      // Calculate r*W, the projection of residual through afferent weights
      // - assume backproject() was called and all "bp" matrices hold
      //   updated values.
      // - write r*W to "output": each call is cumulative.
      resp_to_residual((*i)->name(),gammaaff);
    }
  } // r*W ready in output
 
  // calculate lateral response and update all matrices
  for (int x=0; x<nr; ++x) {
    for (int y=0; y<nc; ++y) {

      // lateral response should use prev_map_activity[][] as f(a)

      // update residual, which holds -dE/da
      residual[x][y] = alpha_settle 
                     * map_activity[y][x] 
                     * ( output[x][y] +
                         ((obj_lateral_on)
                          ?( gammaexc*lat_resp_named(x,y,"LateralExcitatory")
                            -gammainh*lat_resp_named(x,y,"LateralInhibitory") )
                          :0.0
                         ) 
                       );

      // update output(for plotting) and bp(as new a, for mor back projection) 
      output[x][y] = bp[x][y] = resp_to_inp[y][x] + residual[x][y];

    }
  }

}


void LissomMap::train_obj_weights() {

  const double alpha_aff     = params->get_with_default("alpha_input",0.0);
  const double alpha_exc = params->get_with_default("alpha_exc",0.0);
  const double alpha_inh = params->get_with_default("alpha_inh",0.0);
  const bool   crop_negative_wts = params->get_with_default("crop_negative_wts",False);
  const bool   allow_negative_wts = params->get_with_default("allow_negative_wts",False);
  const bool   l2_norm = params->get_with_default("l2_norm",False);
  
  // 0. Set up 'f(a)' and 'r' properly in bp[][] and upstream 
  //    residual[][] respectively.
  backproject(); 

  int nr = output.nrows();
  int nc = output.ncols();

  // 1. For all afferent regions, adapt weights
  for (inputs_type::iterator i=inputs.begin();i != inputs.end(); ++i) {
    if (String::non_numeric_basename_matches<string>((*i)->name(),"Afferent")){
  
      // calculate the afferent region number and get the input region
      string name = (*i)->name();
      int aff   = std::find(ISEQ(aff_names),name) - aff_names.begin();
      NeuralRegion* inp = dynamic_cast<NeuralRegion*>(get_input((*i)->name()));
      ActivityMatrix& r = inp->residual;
  
      // loop through target units
      for (int ui=0; ui<nr; ++ui) {
        for (int uj=0; uj<nc; ++uj) { 
  
          /* Array idx swapped to translt b/w correct and legacy formats. */
          Neuron&              neuron = cortex_map[uj][ui];
          AfferentWeightsList& awts   = *(neuron.weights);
      
          int lowk  = MAX(neuron.centery-rf_radius,0   );
          int highk = MIN(neuron.centery+rf_radius,RN-1);
          int lowl  = MAX(neuron.centerx-rf_radius,0   );
          int highl = MIN(neuron.centerx+rf_radius,RN-1);
          
          if (l2_norm) {
                  // adapt weight
                  double sum=0.0; 
                  for (int k=lowk; k <=highk; k++)
                    for (int l=lowl; l<=highl; l++) {
                        a_weight* const wptr = &(awts[aff][l-lowl][k-lowk]);
                        if (*wptr > -DEATH_FLAG || allow_negative_wts) {
                          *wptr += alpha_aff * r[k][l] * bp[ui][uj];    
                          sum += (*wptr)*(*wptr);
                        }
                    }

                  // normalize
                  sum=1.0/sqrt(sum);
                  for (int k=lowk; k <=highk; k++)
                    for (int l=lowl; l<=highl; l++) {
                        a_weight* const wptr = &(awts[aff][l-lowl][k-lowk]);
                          *wptr *= sum;
                    }

          } else { // L1-norm 

                  // adapt weight
                  double sum=0.0; 
                  for (int k=lowk; k <=highk; k++)
                    for (int l=lowl; l<=highl; l++) {
                        a_weight* const wptr = &(awts[aff][l-lowl][k-lowk]);
                        if (((*wptr) > -DEATH_FLAG) || allow_negative_wts) {
                          (*wptr) += alpha_aff * r[k][l] * bp[ui][uj];  
                          if (crop_negative_wts && ((*wptr)<0.0)) *wptr=0.0;

                          // jbednar: Why fabs?  To avoid divide by zero when doing 1.0/sum,
                          // when crop_negative_wts==False?
                          sum += fabs((*wptr));
                        }
                    }

                  // normalize
                  assert(sum!=0.0); 
                  sum = 1.0/sum;
                  for (int k=lowk; k <=highk; k++)
                    for (int l=lowl; l<=highl; l++) {
                        a_weight* const wptr = &(awts[aff][l-lowl][k-lowk]);
                        if (((*wptr) > -DEATH_FLAG) || allow_negative_wts) {
                                (*wptr) *= sum;
                        }
                    }

          } // end else

        } // uj loop
      } // ui loop

    } // endf afferent check if
  } // end input regions loop

  // 2. Adapt excitatory lateral connections
  for (int ui=0; ui<nr; ++ui) {
    for (int uj=0; uj<nc; ++uj) {

        /* Array idx swapped to translt b/w correct and legacy formats. */
        Neuron&              neuron = cortex_map[uj][ui];
      
        const int lowi  = MAX(uj-exc_rad,0);
        const int highi = MIN(uj+exc_rad,N-1);
        const int lowj  = MAX(ui-exc_rad,0);
        const int highj = MIN(ui+exc_rad,N-1);
        
        double sum=0.0;
        for (int i=lowi; i<=highi; i++)
          for (int j=lowj; j<=highj; j++) {
            l_weight* const wptr = &(neuron.lat_exc_wts[LAT_INDEX(uj,ui,i,j,exc_rad,exc_array_width)]);
            (*wptr) += alpha_exc * bp[j][i] * bp[ui][uj];
            sum += fabs((*wptr));
          }

        assert(sum!=0.0);
        sum = 1.0/sum;
        for (int i=lowi; i<=highi; i++)
          for (int j=lowj; j<=highj; j++) {
            l_weight* const wptr = &(neuron.lat_exc_wts[LAT_INDEX(uj,ui,i,j,exc_rad,exc_array_width)]);
            (*wptr) *= sum; 
          }
        
    }
  }

  // 3. Adapt inhibitory lateral connections
  for (int ui=0; ui<nr; ++ui) {
    for (int uj=0; uj<nc; ++uj) {

        /* Array idx swapped to translt b/w correct and legacy formats. */
        Neuron&              neuron = cortex_map[uj][ui];
      
        const int lowi  = MAX(uj-inh_rad,0);
        const int highi = MIN(uj+inh_rad,N-1);
        const int lowj  = MAX(ui-inh_rad,0);
        const int highj = MIN(ui+inh_rad,N-1);
        
        double sum=0.0;
        for (int i=lowi; i<=highi; i++)
          for (int j=lowj; j<=highj; j++) {
            l_weight* const wptr = &(neuron.lat_inh_wts[LAT_INDEX(uj,ui,i,j,inh_rad,inh_array_width)]);
            (*wptr) -= alpha_inh * bp[j][i] * bp[ui][uj];
            sum += fabs((*wptr));
           
          }
        
        assert(sum!=0.0);
        sum=1.0/sum;
        for (int i=lowi; i<=highi; i++)
          for (int j=lowj; j<=highj; j++) {
            l_weight* const wptr = &(neuron.lat_inh_wts[LAT_INDEX(uj,ui,i,j,inh_rad,inh_array_width)]);
            (*wptr) *= sum;
            // this fails! : assert((*wptr)>0.0);
          }
    }
  }

}
#endif /* OBJ_LISSOM */


void LissomMap::compute_responses(int ts, bool activatefn, bool continuousdynamics)
{
  const double gammaaff = params->get_with_default("gammaaff",1.0);
  const double gammaexc = params->get_with_default("gammaexc",1.0);
  const double gammainh = params->get_with_default("gammainh",1.0);
  const Tristate ls_type= params->get_with_default("lateral_sign_type",Uninitialized);
  const double ls_scale = params->get_with_default("lateral_sign_scale",0.0);
  const double ls_offset= params->get_with_default("lateral_sign_offset",-1.0);
  const double delta    = params->get_with_default("delta",0.0);
  const double beta     = params->get_with_default("beta",1.0);
  const double noise    = params->get_with_default("noise",0.0);
  const bool   allow_negative_wts = params->get_with_default("allow_negative_wts",False);

  int i,pe,j,row,count; double resp;
 
  if ((ts==1) && (!continuousdynamics)){ /* Don't add lateral interaction for first */
    for(i=0; i<perows; i++){  /* Add up input and lateral activations. */
      row = MAPROW(i);
      for(j=0; j<N; j++)
        {
          resp = gammaaff * resp_to_inp[row][j];
          if (noise > 0.0)   /* If neuron is noisy, add some zero mean noise */
            resp += noise * (shuffled_rand()-0.5);
          map_activity[row][j] = init_activity[row][j] = activatefn ? sigmoid(resp,delta,beta) : resp;
          if (map_activity[row][j] > 0.0)         
            activity_list[i][++ACTLIST_COUNT(activity_list[i])] = j;
        }
      assert(ACTLIST_COUNT(activity_list[i]) >= 0); assert(ACTLIST_COUNT(activity_list[i]) <= N);
    }
    ipc_barrier(); /* a barrier before ipc_put is essential to avoid races */
    for(i=0; i<perows; i++){  /* Store all the activity lists */
      row = MAPROW(i);
      for (pe=0; pe < NPEs; pe++)    /* Store activity list to all processors */
        ipc_put_to(&(prev_activity_list[row][0]), &(activity_list[i][0]),
                   IPC_INT,ACTLIST_COUNT(activity_list[i])+1,pe);
      ACTLIST_COUNT(activity_list[i]) = 0;
    }
    return;
  }

  (void)allow_negative_wts;
  assert(!allow_negative_wts);
  /*============================= OTHERWISE ============================= */
  /* Copy the the activity from the previous settling iteration.
   * This is used to compute the responses in the next iteration.
   */
  for(pe=0; pe< NPEs; pe++)          /* From each processor     */
    for (i=0; i<perows; i++){          /* Fetch each row to my PE */
      row = ARBITRARY_MAPROW(i,pe);
      ipc_get_to(&(prev_map_activity[row][0]),&(map_activity[row][0]),IPC_DOUBLE,N,pe);
    }

  /* Set markers after two (heuristic) iterations are complete */
  if (ts==3){
    if (exc_rad > 1) set_markers(exc_rad);
    else set_markers(2);
  }
  /* Gather the activity from each row to the gather_activity array. This
   * increases cache_locality and helps increase performance in the unrolled 
   * loop. gather_activity is indexed from [1..count] like the activity list *
   */
  for (i=0; i<N; i++)
    if ((count=ACTLIST_COUNT(prev_activity_list[i])) > 0)
      for (j=ACTLIST_START(prev_activity_list); j<=count; j++)
        gather_activity[i][j]=prev_map_activity[i][prev_activity_list[i][j]];

  for(i=0; i<perows; i++){  /* Add up input and lateral activations. */
    row = MAPROW(i);
    for(j=0; j<N; j++){
      
      if (activity_marker[row][j]){ /* Activity exists in neighborhood*/
        resp = gammaaff * resp_to_inp[row][j];
        
        if (noise > 0.0)   /* If neuron is noisy, add some zero mean noise */
          resp += noise * (shuffled_rand()-0.5);

        const Neuron& neuron = cortex_map[row][j];
        if (exc_rad <= 1)/*   Use optimized version of the lat_resp function*/
          resp += gammaexc * lat_resp_rad1(row, j, exc_rad, exc_array_width, exc_array_width_2,
                                           neuron.lat_exc_wts, lat_exc_dimension);

        else
          resp += gammaexc * lat_resp(row, j, exc_rad, exc_array_width, 
                                      neuron.lat_exc_wts, lat_exc_dimension);
    
        /* Add inhibition only if fixed-sign and activity is already
           greater than zero, or if variable sign is being used (and
           thus actual sign is not predictable) */
        const double act = activatefn ? sigmoid(resp,delta,beta) : resp;
        if ( act>0.0 || ls_type>Uninitialized) {
          const double inh_resp = lat_resp(row,j,inh_rad,inh_array_width,neuron.lat_inh_wts, lat_inh_dimension);
          const double scaled_inh_resp =
            /* Standard: fixed-sign inhibitory connections */
            (ls_type==Uninitialized?
             gammainh*inh_resp :
             
             /* Untested: scales sign of long-range connections by
                strength of incoming inh activations */
             // Consider removing the first inh_resp; right now
             // responses are squared for high inh or low ls_offset
             (ls_type==False?
              gammainh*inh_resp*(ls_scale*inh_resp-ls_offset) : 
              
              /* Untested: scales sign of long-range connections by resp_to_inp */
              gammainh*inh_resp*(ls_scale*resp_to_inp[row][j]-ls_offset) ));
          
          const double inp      = resp-scaled_inh_resp;
          map_activity[row][j] = activatefn ? sigmoid(inp, delta, beta) : inp;
        }
        
        if (map_activity[row][j] > 0.0)   
          activity_list[i][++ACTLIST_COUNT(activity_list[i])] = j;
        
      }
      else map_activity[row][j]=0.0;
    }
  }

  ipc_barrier(); /* a barrier before ipc_put is essential to avoid races */
  for(i=0; i<perows; i++){  /* Store all the activity lists */
    row = MAPROW(i);
    for (pe=0; pe < NPEs; pe++)   /* Store activity list to all processors */
      ipc_put_to(&(prev_activity_list[row][0]), &(activity_list[i][0]),
                 IPC_INT,ACTLIST_COUNT(activity_list[i])+1,pe);
    ACTLIST_COUNT(activity_list[i]) = 0;
  }
}



int LissomMap::binary_search(const actlist& list, int lower, int upper, int value)
{
  /* Don't bother to search if entirely outside the range */
  if (value <= list[lower]) return lower;    /* value lower than entire range */
  if (value >  list[upper]) return upper+1;  /* value higher than entire range */

  /* Binary search within range */
  while(lower < upper) {
    const int middle = (lower+upper) >> 1;  /* same as divide by 2 */
    const int compare_to = list[middle];
    if      (value > compare_to) lower = middle + 1;
    else if (value < compare_to) upper = middle;
    else return(middle);   /* value has been found at index middle */
  }
  return(lower); /* value not found, but is in range */
}



void LissomMap::crop_actlist(const actlist& list, int lowval, int highval, int *lowidx, int *highidx)
{
  (*lowidx )=ACTLIST_START(list);
  (*highidx)=ACTLIST_COUNT(list);
  
  if (lowval == 0)
    (*lowidx) = 1;
  else
    (*lowidx) = binary_search(list,(*lowidx),(*highidx),lowval);

  if (highval==(N-1))
    (*highidx) = ACTLIST_COUNT(list); 
  else {
    (*highidx) = binary_search(list,(*lowidx),(*highidx),highval);

    if ((*highidx) > ACTLIST_COUNT(list) || (list[(*highidx)]> highval))
      (*highidx)--;
  }
}



/*
void LissomMap::crop_actlist_orig(const actlist& list, int lowval, int highval, int *lowidx, int *highidx)
{
  (*lowidx )=ACTLIST_START(list);
  (*highidx)=ACTLIST_COUNT(list);

  while ((*highidx) >= (*lowidx)  &&  list[(*lowidx )] < lowval ) (*lowidx )++;
  while ((*highidx) >= (*lowidx)  &&  list[(*highidx)] > highval) (*highidx)--;
} 
*/



/*
void LissomMap::verify_actlist(const actlist& list)
{
  int i;
  int count=ACTLIST_COUNT(list);

  if (count <0 || count > N) {
    ipc_notify(IPC_ALL,IPC_ERROR,"List has invalid count (%d)",count);
    return;
  }
  for (i=ACTLIST_START(list); i<count; i++)
    if (list[i]>= list[i+1]) {
      ipc_notify(IPC_ALL,IPC_ERROR,"List is not in proper order -- for %d elements, list[%d]=%d > list[%d]=%d",
                count,i,list[i],i+1,list[i+1]);
      return;
    }
  return;
}
*/

/*
void LissomMap::verify_actlist_array(const actlist& lists[],int num)
{
  int i;
  for (i=0; i<num; i++)
    verify_actlist(lists[i]);
}
*/


void LissomMap::initialize_actlists( void )
{
  int i;
  for(i=0; i<perows; i++) {
    ACTLIST_COUNT(prev_activity_list[MAPROW(i)]) = 0;
    ACTLIST_COUNT(activity_list[i]) = 0;
  }
}



void LissomMap::initialize_markers( void )
{
  int i,j;

  for(i=0; i<perows; i++){
    const int row = MAPROW(i);
    for(j=0; j<N; j++)  
      activity_marker[row][j] = True;
  }
}



void LissomMap::set_markers(int radius)
{
  int i,j,k,l;
  const double act_threshold = params->get_with_default("activity_threshold",0.0);
  const bool   has_threshold = act_threshold>=0;
  
  for(i=0; i<perows; i++){
    const int row   = MAPROW(i);                      
    const int lowk  = MAX(row-radius,0);
    const int highk = MIN(row+radius,N-1);

    for(j=0; j<N; j++){
      const int lowl  = MAX(j-radius,0);
      const int highl = MIN(j+radius,N-1);
      double marker_resp = 0.0;

      for(k=lowk; k<= highk; k++)
        for(l=lowl; l<=highl; l++)
          marker_resp += prev_map_activity[k][l];

      activity_marker[row][j] = (!has_threshold || (marker_resp > act_threshold));
    }
  }
}



double LissomMap::lat_resp(int ui, int uj, int radius, int array_width, l_weight *weights, int weight_idx_bound)
{
  double resp=0;
  
  const int lowi  = MAX(ui-radius,0); 
  const int highi = MIN(ui+radius,N-1);
  
  (void)weight_idx_bound; /* Parameter used only with assertions */
  
  for (int i=lowi; i<=highi; i++)
    if (ACTLIST_COUNT(prev_activity_list[i]) > 0){
        const int lowj  = MAX(uj-radius,0); 
        const int highj = MIN(uj+radius,N-1); 
        const int partial_idx = PARTIAL_LAT_INDEX(ui,uj,i,radius,array_width); 
        int high_idx, low_idx;
        
        crop_actlist(prev_activity_list[i],lowj,highj,&low_idx,&high_idx);
        assert (high_idx <= ACTLIST_COUNT(prev_activity_list[i]));
 
        /* This is a time-critical loop; make sure that the compiler unrolls it */
        for (int j=low_idx; j <= high_idx; j++) {
          const int idx = FULL_LAT_INDEX(partial_idx,prev_activity_list[i][j]);
          assert (prev_activity_list[i][j] >= lowj);
          assert (prev_activity_list[i][j] <= highj);
          assert (idx >= 0); assert (idx <= weight_idx_bound);
          
          resp += gather_activity[i][j] * weights[idx];
        }
    }
  return(resp);
}



double LissomMap::lat_resp_rad1(int ui, int uj, int radius, int array_width, int array_width_2, l_weight *weights, int weight_idx_bound)
{
  double resp, resp1, resp2;

  /* This parameter is for error checking and consistency with
     modify_latwt_loop; it is currently unused */
  (void)weight_idx_bound;
  
  if (radius==0){
    resp = prev_map_activity[ui][uj]; /* Weights are 1.0 */
    return(resp);
  }
  resp = resp1 = resp2=0.0;
  if ((ui > 0) && (ui < N-1))    /* Case where 'ui' is away from border      */
    {
      if ((uj >0) && (uj < N-1)) /* Case where 'uj' is also away from border */
        {
          /* FULL_LAT_INDEX need not be called when either argument is 0 */
          resp   = prev_map_activity[ui-1][uj-1] * weights[                             0 ];
          resp1  = prev_map_activity[ui-1][uj  ] * weights[                             1 ];
          resp2  = prev_map_activity[ui-1][uj+1] * weights[                             2 ];

          resp  += prev_map_activity[ui  ][uj-1] * weights[               array_width     ];
          resp1 += prev_map_activity[ui  ][uj  ] * weights[FULL_LAT_INDEX(array_width,  1)];
          resp2 += prev_map_activity[ui  ][uj+1] * weights[FULL_LAT_INDEX(array_width,  2)];

          resp  += prev_map_activity[ui+1][uj-1] * weights[               array_width_2   ];
          resp1 += prev_map_activity[ui+1][uj  ] * weights[FULL_LAT_INDEX(array_width_2,1)];
          resp2 += prev_map_activity[ui+1][uj+1] * weights[FULL_LAT_INDEX(array_width_2,2)];
        }
      else{                     /*  'uj' is at some border */
        const int lowj  = MAX(uj-1,0);
        const int highj = MIN(uj+1,N-1);
        const int partial_index0 =               - uj + 1;
        const int partial_index1 = array_width   - uj + 1; 
        const int partial_index2 = array_width_2 - uj + 1;
        int j;
        
        for(j=lowj; j<=highj; j++){
          resp  += prev_map_activity[ui-1][j] * weights[FULL_LAT_INDEX(partial_index0,j)];
          resp1 += prev_map_activity[ui  ][j] * weights[FULL_LAT_INDEX(partial_index1,j)];
          resp2 += prev_map_activity[ui+1][j] * weights[FULL_LAT_INDEX(partial_index2,j)];
        }
      }
      resp += resp1 + resp2;
    }
  else if (ui > 0){            /* Case where 'ui' is at top border (N-1)     */
    const int lowj  = MAX(uj-1,0);
    const int highj = MIN(uj+1,N-1);
    const int partial_index0 = -uj + 1;
    const int partial_index1 = array_width - uj + 1;
    int j;
      
    for(j=lowj; j<=highj; j++){
      resp += prev_map_activity[ui-1][j] * weights[FULL_LAT_INDEX(partial_index0,j)];
      resp1+= prev_map_activity[ui  ][j] * weights[FULL_LAT_INDEX(partial_index1,j)];
    }
    resp += resp1;
  }
  else{                        /* Case where 'ui' is at bottom border (0)    */
    const int lowj  = MAX(uj-1,0);
    const int highj = MIN(uj+1,N-1);
    const int partial_index0 = array_width   - uj + 1;
    const int partial_index1 = array_width_2 - uj + 1;
    int j;
      
    for(j=lowj; j<=highj; j++){
      resp += prev_map_activity[ui  ][j] * weights[FULL_LAT_INDEX(partial_index0,j)];
      resp1+= prev_map_activity[ui+1][j] * weights[FULL_LAT_INDEX(partial_index1,j)];
    }
    resp += resp1;
  }
  return(resp);
}


double LissomMap::lat_resp_named(int ui, int uj, const string& name)
{
 
  double sum = 0.0;

  // Assumptions: prev_map_activity[][] is used to calculate
  // the lateral contributions.

  Neuron&              neuron = cortex_map[uj][ui];

  if (name == "LateralExcitatory") {
    /* Array idx swapped to translt b/w correct and legacy formats. */
    const int lowi  = MAX(uj-exc_rad,0);
    const int highi = MIN(uj+exc_rad,N-1);
    const int lowj  = MAX(ui-exc_rad,0);
    const int highj = MIN(ui+exc_rad,N-1);
        
    for (int i=lowi; i<=highi; i++)
      for (int j=lowj; j<=highj; j++)
         sum += prev_map_activity[i][j] 
             * neuron.lat_exc_wts[LAT_INDEX(uj,ui,i,j,exc_rad,exc_array_width)];

  } else if (name == "LateralInhibitory") {
    /* Array idx swapped to translt b/w correct and legacy formats. */
    const int lowi  = MAX(uj-inh_rad,0);
    const int highi = MIN(uj+inh_rad,N-1);
    const int lowj  = MAX(ui-inh_rad,0);
    const int highj = MIN(ui+inh_rad,N-1);
        
    for (int i=lowi; i<=highi; i++)
      for (int j=lowj; j<=highj; j++)
         sum += prev_map_activity[i][j] 
             * neuron.lat_inh_wts[LAT_INDEX(uj,ui,i,j,inh_rad,inh_array_width)];
  }

  return sum;
}



void LissomMap::modify_all_weights(void)
{
  const double alpha_aff     = params->get_with_default("alpha_input",0.0);
  const bool   normalize_aff = params->get_with_default("normalize_aff",False);
  if (alpha_aff!=0) modify_weights(alpha_aff,normalize_aff);

  const double alpha_exc = params->get_with_default("alpha_exc",0.0);
  const double alpha_inh = params->get_with_default("alpha_inh",0.0);
  if (alpha_exc!=0 || alpha_inh!=0) modify_lat_wts(alpha_exc,alpha_inh);
  
  ipc_barrier();      /* To ensure all wts have been modified */
}



void LissomMap::modify_weights(double alpha, bool normalize)
{                         
  const bool   allow_negative_wts = params->get_with_default("allow_negative_wts",False);

  (void)allow_negative_wts;
  assert(!allow_negative_wts);

  if (int(aff_inputs.size()) != num_aff) {
    ipc_notify(IPC_ALL,IPC_ERROR,"No response; the number of afferent inputs is required to be %d but is %d",
               num_aff, int(aff_inputs.size()));
    return;
  }
      
  int i,j,k,l,aff;
  
  if (normalize)    /* If afferents from each receptor are normalized */
    for(i=0; i<RN; i++)
      for(j=0; j<RN; j++)
        for (aff=0; aff<num_aff; aff++)
          afferent_sum[INP_INDEX(aff,i,j)] = 0.0;
  
  for(i=0; i<perows; i++){
    const int row = MAPROW(i);
    
    for(j=0; j<N; j++)    /* Only wts of active neurons need to be modified */
      if (normalize || (prev_map_activity[row][j] > 0.0)){
        Neuron&              neuron = cortex_map[row][j];
        AfferentWeightsList& awts = *(neuron.weights);
        
        /* The weights must be modified only within the receptive field. */
        const int lowk   = MAX(neuron.centerx-rf_radius,0   );
        const int highk  = MIN(neuron.centerx+rf_radius,RN-1);
        const int lowl   = MAX(neuron.centery-rf_radius,0   );
        const int highl  = MIN(neuron.centery+rf_radius,RN-1);
        
        if (prev_map_activity[row][j] > 0.0){
          
          double deltaw = alpha * prev_map_activity[row][j]; 
          double length = 0.0;

          for (aff=0; aff<num_aff; aff++)         
            for (k=lowk; k<=highk; k++){
              for (l=lowl; l<=highl; l++) {
                a_weight* const wtptr = &(awts[aff][k-lowk][l-lowl]);

                if (*wtptr > -DEATH_FLAG) {
                  const double load = (*aff_inputs[aff])[l][k];
                
                  if ( load > 0.0){
                    const double delt = load * deltaw;
                    *wtptr += delt;
                    length += delt;
                  }
                }
              }
            }
          
          if (!normalize) /* Do only if receptor_normalize is not called */
            normalize_afferent_wts(lowk,highk,lowl,highl,row,j,1.0/(1.0+length));
        }
        
        if (normalize) /* Do this whether map activity is >0 or not */
          for (k=lowk; k<=highk; k++)
            for (l=lowl; l<=highl; l++)
              for (aff=0; aff<num_aff; aff++) {
                assert (INP_INDEX(aff,k,l) >= INP_INDEX(aff,  0,0));
                assert (INP_INDEX(aff,k,l) <  INP_INDEX(aff+1,0,0));
                
                afferent_sum[INP_INDEX(aff,k,l)] += awts[aff][k-lowk][l-lowl];
              }
      }
  }
  
  if (normalize) receptor_normalize();
}



void LissomMap::receptor_normalize(void)
{
  double result;
  int i,j,k,l,pe,aff;

  for(k=0; k<RN; k++)  /* Store local data in temp array */
    for(l=0; l<RN; l++)
      for (aff=0; aff<num_aff; aff++)
        temp_sum[INP_INDEX(aff,k,l)] = afferent_sum[INP_INDEX(aff,k,l)];

  ipc_barrier();

  for(pe=0; pe<NPEs; pe++)
    if (!PEISME(pe)){
      ipc_get_to(&aff_norm_sum[0], &afferent_sum[0], IPC_DOUBLE,num_aff*RN*RN,pe);
      
      for(k=0; k<RN; k++)
        for(l=0; l<RN; l++)
          for (aff=0; aff<num_aff; aff++)
            temp_sum[INP_INDEX(aff,k,l)] += aff_norm_sum[INP_INDEX(aff,k,l)];
    }
  
  ipc_barrier(); /* To ensure that afferent_sum is not changed before fetched */

  for(k=0; k<RN; k++)
    for (l=0; l<RN; l++)
      for (aff=0; aff<num_aff; aff++)
        afferent_sum[INP_INDEX(aff,k,l)] = 1.0/temp_sum[INP_INDEX(aff,k,l)];
  
  for(i=0; i<perows; i++){
    const int row = MAPROW(i);
    for(j=0; j<N; j++){
      Neuron&              neuron = cortex_map[row][j];
      AfferentWeightsList& awts = *(neuron.weights);

      const int lowk   = MAX(neuron.centerx-rf_radius,0   );
      const int highk  = MIN(neuron.centerx+rf_radius,RN-1);
      const int lowl   = MAX(neuron.centery-rf_radius,0   );
      const int highl  = MIN(neuron.centery+rf_radius,RN-1);

      double length = 0.0;

      /* This is a time-critical loop; make sure that the compiler unrolls it */
      for (k=lowk; k<=highk; k++)
        for (l=lowl; l<=highl; l++)
          for (aff=0; aff<num_aff; aff++){
            a_weight* const wtptr = &(awts[aff][k-lowk][l-lowl]);
            assert (INP_INDEX(aff,k,l) >= INP_INDEX(aff,  0,0));
            assert (INP_INDEX(aff,k,l) <  INP_INDEX(aff+1,0,0));
            
            result  = (*wtptr) * afferent_sum[INP_INDEX(aff,k,l)];
            length += result;
            
            *wtptr = result;
          }
      
      normalize_afferent_wts(lowk,highk,lowl,highl,row,j,1.0/length);
    }
  }
}



void LissomMap::normalize_afferent_wts(int lowk, int highk, int lowl, int highl, int row, int j, double length)
{
  int k,l,aff;
  Neuron&              neuron = cortex_map[row][j];
  AfferentWeightsList& awts   = *(neuron.weights);
  
  /* This is a time-critical loop; make sure that the compiler unrolls it */
  for (k=lowk; k<=highk; k++)
    for (l=lowl; l<=highl; l++)
      for (aff=0; aff<num_aff; aff++) {
        a_weight* const wtptr = &(awts[aff][k-lowk][l-lowl]);
        if (*wtptr > -DEATH_FLAG)
          *wtptr *= length;
      }
}



void LissomMap::modify_lat_wts(double alpha_exc, double alpha_inh)
{
  int i,j;

  for(i=0; i<perows; i++){
    const int row = MAPROW(i);

    for(j=0; j<N; j++)     /* Modify only neurons with nonzero activity */
      if (prev_map_activity[row][j] > 0.0){
        Neuron& neuron = cortex_map[row][j];
        
        if (exc_rad==1 && !exc_connections_killed)
          modify_latwt_loop_rad1(row,j,exc_rad, exc_array_width, exc_array_width_2,
                                 neuron.lat_exc_wts, lat_exc_dimension, alpha_exc);
        else 
          modify_latwt_loop(     row,j, exc_rad, exc_array_width,
                                 neuron.lat_exc_wts, lat_exc_dimension, alpha_exc);
        
        modify_latwt_loop(       row,j, inh_rad, inh_array_width,
                                 neuron.lat_inh_wts, lat_inh_dimension, alpha_inh);
      }
  }
}



void LissomMap::modify_latwt_loop(int i, int j, int radius, int array_width, l_weight *weights, int weight_idx_bound, double alpha)
{
  int k,m;
  const int lowk  = MAX(i-radius,0  );
  const int highk = MIN(i+radius,N-1);
  const int lowl  = MAX(j-radius,0  );
  const int highl = MIN(j+radius,N-1);
  const double loop_const=alpha * prev_map_activity[i][j];
  double new_length=0;

  (void)weight_idx_bound; /* Parameter used only with assertions */

  for (k=lowk; k<=highk; k++)
    if ( ACTLIST_COUNT(prev_activity_list[k]) > 0) {
      const int partial_idx = PARTIAL_LAT_INDEX(i,j,k,radius,array_width);
      int low_idx, high_idx;
      crop_actlist(prev_activity_list[k],lowl,highl,&low_idx,&high_idx);

      /* This is a time-critical loop; make sure that the compiler unrolls it */
      for (m=low_idx; m <= high_idx; m++) {
        const int l   = prev_activity_list[k][m];
        const int idx = FULL_LAT_INDEX(partial_idx,l);
        assert (l >= lowl); assert (l <= highl);
        assert (idx >= 0);  assert (idx <= weight_idx_bound);
        l_weight* const wt = &(weights[idx]);
        if (*wt>=0) {
          const double delt = gather_activity[k][m] * loop_const;
          new_length   += delt;
          *wt += delt;
        }
      }
    }
  new_length = 1.0/(1.0+new_length);

  for (k=lowk; k<=highk; k++) {
    int idx = LAT_INDEX(i,j,k,lowl,radius,array_width);
    int l;
    assert (idx >= 0);

    /* This is a time-critical loop; make sure that the compiler unrolls it */
    for (l=lowl; l<=highl; l++,idx++) {
      assert (idx <= weight_idx_bound);
      l_weight* const wt = &(weights[idx]);
      if (*wt>=0)
        *wt *= new_length;
    }
  }
}



void LissomMap::modify_latwt_loop_rad1(int ui, int uj, int radius, int array_width, int array_width_2, l_weight *weights, int weight_idx_bound, double alpha)
{
  int j;
  double length=0;
  double delta0=0;
  double delta1=0;
  double delta2=0;
  double delta3=0;
  const double loop_constant = alpha * prev_map_activity[ui][uj];
  
  /* These parameters are for error checking and consistency with
     modify_latwt_loop; they are currently unused */
  (void)radius;
  (void)weight_idx_bound;
  
  if ((ui > 0) && (ui < N-1))    /* Case where 'ui' is away from border      */
    {
      if ((uj >0) && (uj < N-1)) /* Case where 'uj' is also away from border */
        {
          /* FULL_LAT_INDEX need not be called when either argument is 0 */
          delta0 = prev_map_activity[ui-1][uj-1] * loop_constant; weights[                           0   ] += delta0; 
          delta1 = prev_map_activity[ui-1][uj  ] * loop_constant; weights[                           1   ] += delta1;
          delta2 = prev_map_activity[ui-1][uj+1] * loop_constant; weights[                           2   ] += delta2;
          length = delta0 + delta1 + delta2;

          delta3 = prev_map_activity[ui  ][uj-1] * loop_constant; weights[               array_width     ]+= delta3; 
          delta0 = prev_map_activity[ui  ][uj  ] * loop_constant; weights[FULL_LAT_INDEX(array_width,1)  ]+= delta0; 
          delta1 = prev_map_activity[ui  ][uj+1] * loop_constant; weights[FULL_LAT_INDEX(array_width,2)  ]+= delta1;
          length += delta3 + delta0 + delta1;

          delta2 = prev_map_activity[ui+1][uj-1] * loop_constant; weights[               array_width_2   ]+= delta2;
          delta3 = prev_map_activity[ui+1][uj  ] * loop_constant; weights[FULL_LAT_INDEX(array_width_2,1)]+= delta3;
          delta0 = prev_map_activity[ui+1][uj+1] * loop_constant; weights[FULL_LAT_INDEX(array_width_2,2)]+= delta0;
          length = 1.0/(1.0+length + delta2 + delta3 + delta0);

          for(j=0; j<3; j++){
            weights[                             j ] *= length;
            weights[FULL_LAT_INDEX(array_width  ,j)] *= length;
            weights[FULL_LAT_INDEX(array_width_2,j)] *= length;
          }
        }
      else{                     /*  'uj' is at some border */
        const int lowj  = MAX(uj-1,0); 
        const int highj = MIN(uj+1,N-1);
        const int index0 =               -uj+1;
        const int index1 = array_width   -uj+1; 
        const int index2 = array_width_2 -uj+1;

        for(j=lowj; j<=highj; j++){
          delta0 = prev_map_activity[ui-1][j] * loop_constant; weights[FULL_LAT_INDEX(index0,j)] += delta0; 
          delta1 = prev_map_activity[ui  ][j] * loop_constant; weights[FULL_LAT_INDEX(index1,j)] += delta1; 
          delta2 = prev_map_activity[ui+1][j] * loop_constant; weights[FULL_LAT_INDEX(index2,j)] += delta2; 
          length += delta0 + delta1 + delta2; 
        }
        length = 1.0/(1.0+length);
        for(j=lowj; j<=highj; j++){
          weights[FULL_LAT_INDEX(index0,j)] *= length;
          weights[FULL_LAT_INDEX(index1,j)] *= length;
          weights[FULL_LAT_INDEX(index2,j)] *= length;
        }
      }
    }
  else if (ui > 0){            /* Case where 'ui' is at top border (N-1)     */
    const int lowj  = MAX(uj-1,0); 
    const int highj = MIN(uj+1,N-1);
    const int index0 = -uj + 1;
    const int index1 = array_width - uj + 1;

    for(j=lowj; j<=highj; j++){
      delta0 = prev_map_activity[ui-1][j] * loop_constant; weights[FULL_LAT_INDEX(index0,j)] += delta0; length += delta0; 
      delta1 = prev_map_activity[ui  ][j] * loop_constant; weights[FULL_LAT_INDEX(index1,j)] += delta1; length += delta1;
    }
    length = 1.0/(1.0+length);
    for(j=lowj; j<=highj; j++){
      weights[FULL_LAT_INDEX(index0,j)] *= length;
      weights[FULL_LAT_INDEX(index1,j)] *= length;
    }
  }
  else{                        /* Case where 'ui' is at bottom border (0)    */
    const int lowj  = MAX(uj-1,0); 
    const int highj = MIN(uj+1,N-1);
    const int index1 = array_width   -uj+1;
    const int index2 = array_width_2 -uj+1;

    for(j=lowj; j<=highj; j++){
      delta0 = prev_map_activity[ui  ][j] * loop_constant; weights[FULL_LAT_INDEX(index1,j)] += delta0; length += delta0;
      delta1+= prev_map_activity[ui+1][j] * loop_constant; weights[FULL_LAT_INDEX(index2,j)] += delta1; length += delta1;
    }
    length = 1.0/(1.0+length);
    for(j=lowj; j<=highj; j++){
      weights[FULL_LAT_INDEX(index1,j)] *= length;
      weights[FULL_LAT_INDEX(index2,j)] *= length;
    }
  }  
}



void LissomMap::normalize_lat_wts(int i, int j, int radius, int array_width, l_weight *weights)
{
  int k,l;
  const int lowk   = MAX(i-radius,0  );
  const int highk  = MIN(i+radius,N-1);
  const int lowl   = MAX(j-radius,0  );
  const int highl  = MIN(j+radius,N-1);
  double length = 0;

  for(k=lowk; k<= highk; k++){
    int idx = LAT_INDEX(i,j,k,lowl,radius,array_width);    
    assert (idx >= 0);
    
    /* This is a time-critical loop; make sure that the compiler unrolls it */
    for (l=lowl; l <=highl; l++,idx++) {
      /* assert (idx <= weight_idx_bound); */
      length += weights[idx];
    }
  }
  length = 1.0/length;
  
  for(k=lowk; k<= highk; k++){
    int idx = LAT_INDEX(i,j,k,lowl,radius,array_width);    
    assert (idx >= 0);
    
    /* This is a time-critical loop; make sure that the compiler unrolls it */
    for (l=lowl; l <=highl; l++,idx++){
      /* assert (idx <= weight_idx_bound); */
      weights[idx] *= length;
    }
  }
}



double LissomMap::sigmoid(double activity, double sigdelta, double sigbeta) /* Linear approx of sigmoid */
                                      /* Lower & upper thresholds */    
{
#ifdef DEBUG_SIGMOID
  ipc_notify(IPC_ALL,IPC_VERBOSE,"activity=%.6f; sigdelta=%.6f, sigbeta=%.6f",
             (double)activity,(double)sigdelta,(double)sigbeta);
#endif
  if      (activity < sigdelta) return(0.0);
  else if (activity > sigbeta)  return(1.0);
  else return( (activity-sigdelta)/(sigbeta-sigdelta) );
}

double LissomMap::truesigmoid(double activity, double beta)
{
        double ex = exp(-2*beta*activity);
        return ((1-ex)/(1+ex));
}

double LissomMap::truesigmoidprime(double sigmoidactivity, double beta)
{
        return beta*(1-sigmoidactivity*sigmoidactivity);        
}

double LissomMap::positivesigmoid(double activity, double beta)
{
        return 1.0/(1.0+exp(-(activity-0.5)*beta));
}

double LissomMap::logsigmoid(double activity, double beta)
{
        return (activity>0.0)?(log(1+activity*beta)):0.0;
}

double LissomMap::logsigmoidprime(double activity, double beta)
{
        return (activity>0.0)?(beta/(1+activity*beta)):0.0;
}

/******************************************************************************/
/* Other routines                                                             */
/******************************************************************************/

double LissomMap::size_connection_bytes() const
{
  int memsize = (N*N * ( (2*rf_radius +1)*(2*rf_radius +1)*sizeof(a_weight)*num_aff +
                         (2*exc_rad   +1)*(2*exc_rad   +1)*sizeof(l_weight) +
                         (2*inh_rad   +1)*(2*inh_rad   +1)*sizeof(l_weight) ))/NPEs;
  
  return memsize;
}


double LissomMap::size_unique_connections() const
{
  return        N*N * ( (2*rf_radius +1)*(2*rf_radius +1)*num_aff +
                        (2*exc_rad   +1)*(2*exc_rad   +1) +
                        (2*inh_rad   +1)*(2*inh_rad   +1) );
}


void LissomMap::clear_weights(void)
{
  int local_row, j, k, aff, x, y;

  for(local_row=0; local_row<perows; local_row++){
    const int map_row = MAPROW(local_row);

    for(j=0; j<N; j++) {
      Neuron&              neuron = cortex_map[map_row][j];
      AfferentWeightsList& awts = *(neuron.weights);

      neuron.centerx = 0;
      neuron.centery = 0;
      
      for(k=0; k<lat_exc_dimension; k++)
        neuron.lat_exc_wts[k]   = -DEATH_FLAG;
      for(k=0; k<lat_inh_dimension; k++)
        neuron.lat_inh_wts[k]   = -DEATH_FLAG; 
      for(aff=0; aff<num_aff; aff++)
        for(x=0; x<aff_array_width; x++)
          for(y=0; y<aff_array_width; y++)
            awts[aff][x][y] = -DEATH_FLAG;  
    }
  }
}



void LissomMap::init_weights(bool init_weights)
{
  const double   afferent_randomness = params->get_with_default("randomness",0.0);
  const double   lat_exc_randomness  = params->get_with_default("lat_exc_randomness",0.0);
  const double   lat_inh_randomness  = params->get_with_default("lat_inh_randomness",0.0);
  const Tristate preset_lat_wts      = params->get_with_default("preset_lat_wts",False);
  const Tristate preset_aff_wts      = params->get_with_default("preset_aff_wts",False);
  const double   preset_sigma_aff    = params->get_with_default("preset_sigma_aff",0.0);
  const double   preset_sigma_lat    = params->get_with_default("preset_sigma_lat",0.0);
  const double   preset_sigma_exc    = params->get_with_default("preset_sigma_exc",0.0);
  const double   preset_sigma_inh    = params->get_with_default("preset_sigma_inh",0.0);
  const bool     circular_aff_wts    = params->get_with_default("circular_aff_wts",True);
  const bool     circular_lat_wts    = params->get_with_default("circular_lat_wts",True);
  const bool     smooth              = params->get_with_default("smooth_circular_outlines",False);
  const double   smooth_trim         = params->get_with_default("smooth_circular_radius_trim",0.0);
  const double   rf_radius_float     = params->get_with_default("rf_radius",0.0);
  const double   exc_rad_float       = params->get_with_default("exc_rad",0.0);
  const double   inh_rad_float       = params->get_with_default("inh_rad",0.0);

  const double preset_sigma_e =
   (preset_sigma_exc  != Uninitialized ? preset_sigma_exc : 
    (preset_sigma_lat != Uninitialized ? preset_sigma_lat :
     0.9*exc_rad));

  const double preset_sigma_i =
   (preset_sigma_inh  != Uninitialized ? preset_sigma_inh : 
    (preset_sigma_lat != Uninitialized ? preset_sigma_lat*MIN(1,sqrt(inh_rad)) :
     0.9*inh_rad));

  /* Ensure that all weights are known. */
  if (init_weights) {
    clear_weights();
    
    for(int i=0; i<perows; i++){
      const int row = MAPROW(i);
      /* To make weights generated independent of NPEs */
      shuffled_rand_reset(ShuffledRandom<double>::default_seed+123456+row);
      for(int j=0; j<N; j++) {
        /* Set up each set of weights */
        Neuron& neuron = cortex_map[row][j];
        
        setup_afferent_weights(neuron, *(neuron.weights),
                               row,j, rf_radius, rf_radius_float*rf_radius_float,
                               preset_aff_wts, preset_sigma_aff, afferent_randomness,
                               circular_aff_wts, smooth, rf_radius_float-rf_radius, smooth_trim);

        setup_lateral_weights(row,j, exc_rad, exc_array_width,neuron.lat_exc_wts,
                              preset_lat_wts, preset_sigma_e, lat_exc_randomness,
                              circular_lat_wts, smooth, exc_rad_float-exc_rad, smooth_trim, &exc_connections_killed);
        setup_lateral_weights(row,j, inh_rad, inh_array_width,neuron.lat_inh_wts,
                              preset_lat_wts, preset_sigma_i, lat_inh_randomness,
                              circular_lat_wts, smooth, inh_rad_float-inh_rad, smooth_trim, &inh_connections_killed);
      }
    }
  }
  
  /* Store one neuron with pre-set weights for comparison later. */
  {
    Neuron neuron; /* Dummy; unused */
    const int pos=int(N/2)-1; /* Makes sure the RF is centered on a neuron */

    setup_afferent_weights(neuron, reference_wts.weights,
                           pos,pos, rf_radius, rf_radius_float*rf_radius_float,
                           true, preset_sigma_aff, afferent_randomness,
                           circular_aff_wts, smooth, rf_radius_float-rf_radius, smooth_trim, RN/2+0.5, RN/2+0.5);

    setup_lateral_weights(pos,pos, exc_rad, exc_array_width,&(reference_wts.lat_exc_wts[0]),
                          true, preset_sigma_e, lat_exc_randomness,
                          circular_lat_wts, smooth, exc_rad_float-exc_rad, smooth_trim, &exc_connections_killed);
    setup_lateral_weights(pos,pos, inh_rad, inh_array_width,&(reference_wts.lat_inh_wts[0]),
                          true, preset_sigma_i, lat_inh_randomness,
                          circular_lat_wts, smooth, inh_rad_float-inh_rad, smooth_trim, &inh_connections_killed);
  }
}


void LissomMap::setup_afferent_weights(Neuron& neuron, AfferentWeightsList& awts,
                                       int ui, int uj, int radius, double radius_sq,
                                       int preset, double sigma, double amt_of_randomness,
                                       bool circular, bool smooth,
                                       double radius_trim, double smooth_trim,
                                       double force_cx, double force_cy)
{
  const int    afferent_edge_buffer = params->get_with_default("retina_edge_buffer",0);
  const bool   match_weights        = params->get_with_default("force_matching_initial_weights",True);
  const double cortical_image_width = RN-2*afferent_edge_buffer;
  
  /* Choose receptive field center with a random degree of scatter. */
  const double cx = (force_cx>=0? force_cx : afferent_edge_buffer + cortical_image_width *
    ((ui+0.5)/double(N) + amt_of_randomness*(2.0*shuffled_rand()-1.0)));
  
  const double cy = (force_cy>=0? force_cy : afferent_edge_buffer + cortical_image_width *
    ((uj+0.5)/double(N) + amt_of_randomness*(2.0*shuffled_rand()-1.0)));

  /* The RF center is constrained to be an integer and to be located 
     within the bounds of the input surface, because of how the lowk,
     highk, etc. are calculated and used in the rest of the program. */

  neuron.centerx = MAX(0,MIN(RN-1,int(cx)));
  neuron.centery = MAX(0,MIN(RN-1,int(cy)));

  const int lowk  = (force_cx>=0? neuron.centerx-radius : MAX(neuron.centerx-radius,0   ));
  const int highk = (force_cx>=0? neuron.centerx+radius : MIN(neuron.centerx+radius,RN-1));
  const int lowl  = (force_cy>=0? neuron.centery-radius : MAX(neuron.centery-radius,0   ));
  const int highl = (force_cy>=0? neuron.centery+radius : MIN(neuron.centery+radius,RN-1));
  
  /* Distribute weights in the receptive field randomly, or (if
     preset_weights is on), preset the weights to gaussians of given
     sigma.  Circular receptive fields are accomplished by killing
     all the connections that fall outside a circular radius. */
  double length=0.0;
  if (preset){                                /* Initialize to Gaussian */
    for (int k=lowk; k <=highk; k++) 
      for (int l=lowl; l<=highl; l++){
        double value;
        double xsq = std::abs(cx-k-0.5); /* Uses full floating-point precision */
        double ysq = std::abs(cy-l-0.5); /*  of the center point */
        xsq *= xsq;
        ysq *= ysq;
        value = exp(-((xsq+ysq)/(sigma*sigma)));
        for (int aff=0; aff<num_aff; aff++)
          awts[aff][k-lowk][l-lowl] = value;
        length += num_aff*value;
      }
  }
  
  else{                                       /* Initialize randomly    */
    for (int k=lowk; k <=highk; k++) 
      for (int l=lowl; l<=highl; l++){
        if (match_weights) {
          double value = shuffled_rand();
          for (int aff=0; aff<num_aff; aff++)
            awts[aff][k-lowk][l-lowl] = value;
          length += num_aff*value;
        }
        else {
          for (int aff=0; aff<num_aff; aff++) {
            double value = shuffled_rand();
            awts[aff][k-lowk][l-lowl] = value;
            length += value;
          }
        }
      }
  }
  
  /* Initialize all connections outside the circle to negative, if appropriate */
  if (circular)
    length = prune_circular_aff_weights(neuron.centerx,neuron.centery,
                                        radius, awts,
                                        smooth, radius_trim, smooth_trim, length);
  
  /* Normalize */
  length = 1.0/length;
  for (int k=lowk; k <=highk; k++) 
    for (int l=lowl; l<=highl; l++)
      if (!circular || WITHIN_RADIUS(k-neuron.centerx,l-neuron.centery,radius_sq))
        for (int aff=0; aff<num_aff; aff++) 
          awts[aff][k-lowk][l-lowl] *= length;
}



void LissomMap::setup_lateral_weights(int ui, int uj, int radius,
                                      int array_width, l_weight *weights,
                                      int preset, double sigma, double amt_of_randomness,
                                      bool circular, bool smooth, 
                                      double radius_trim, double smooth_trim,
                                      int* con_killed_flag)
{
  int k,l;
  
  const int lowk  = MAX(ui-radius,0  );
  const int highk = MIN(ui+radius,N-1);
  const int lowl  = MAX(uj-radius,0  );
  const int highl = MIN(uj+radius,N-1);
  
  const int max_num_connections = (highk-lowk+1)*(highl-lowl+1);
  const double mean_wt = 1.0/max_num_connections;
  
  if (preset)
    for(k=lowk; k<=highk; k++)        
      for(l=lowl; l<=highl; l++){
        double xsq = (double)(k-ui); 
        double ysq = (double)(l-uj);
        xsq *= xsq;
        ysq *= ysq;
        weights[LAT_INDEX(ui,uj,k,l,radius,array_width)] = 
          mean_wt * exp( -((xsq+ysq)/(sigma*sigma)));
      }
  else /* Init weights randomly */
    for(k=lowk; k<=highk; k++)        
      for(l=lowl; l<=highl; l++)
        weights[LAT_INDEX(ui,uj,k,l,radius,array_width)] = 
          dist_lat_wt(mean_wt,amt_of_randomness);
  
  if (circular) {
    prune_circular_lat_weights(ui,uj,radius, array_width, weights,
                               smooth, radius_trim, smooth_trim);
    *con_killed_flag=True;
  }
  
  normalize_lat_wts(ui,uj,radius,array_width,weights);
}



void LissomMap::prune_circular_lat_weights(int ui, int uj, int radius, 
                                           int array_width, l_weight *weights,
                                           bool smooth, double radius_trim, double smooth_trim )
{
  const double radius_float    = radius+radius_trim;
  const double radius_sq       = radius_float*radius_float;
  const double smooth_start    = radius+smooth_trim;
  const double smooth_start_sq = smooth_start*smooth_start;
  const double smooth_rad      = (radius_float-smooth_start)/2;
  const double smooth_rad_sq   = smooth_rad*smooth_rad;
    
  const int lowk  = MAX(ui-radius,0  );
  const int highk = MIN(ui+radius,N-1);
  const int lowl  = MAX(uj-radius,0  );
  const int highl = MIN(uj+radius,N-1);

  int k,l;
  
  for (k=lowk; k <=highk; k++)
    for (l=lowl; l<=highl; l++) {
      const double xdist    = k-ui;
      const double ydist    = l-uj;
      const double xdist_sq = xdist*xdist;
      const double ydist_sq = ydist*ydist;
      const double rdist_sq = xdist_sq+ydist_sq;
      
      if (radius>1 && !WITHIN_RADIUS_SQ(xdist_sq,ydist_sq,radius_sq))
        weights[LAT_INDEX(ui,uj,k,l,radius,array_width)] = -DEATH_FLAG;
      
      else if (smooth && (rdist_sq >= smooth_start_sq)) {
        double dist_sq = sqrt(rdist_sq)-smooth_start;
        dist_sq*=dist_sq;
        
        weights[LAT_INDEX(ui,uj,k,l,radius,array_width)] = 
          weights[LAT_INDEX(ui,uj,k,l,radius,array_width)] *
          exp( -(dist_sq/smooth_rad_sq) );
      }
    }
}



double LissomMap::prune_circular_aff_weights(int ui, int uj, int radius, 
                                             AfferentWeightsList& weights,
                                             bool smooth, double radius_trim, double smooth_trim,
                                             double length)
{
  const double radius_float    = radius+radius_trim;
  const double radius_sq       = radius_float*radius_float;
  const double smooth_start    = radius+smooth_trim;
  const double smooth_start_sq = smooth_start*smooth_start;
  const double smooth_rad      = (radius_float-smooth_start)/2;
  const double smooth_rad_sq   = smooth_rad*smooth_rad;
    
  const int lowk  = MAX(ui-radius,0   );
  const int highk = MIN(ui+radius,RN-1);
  const int lowl  = MAX(uj-radius,0   );
  const int highl = MIN(uj+radius,RN-1);

  int k,l;
  
  for (k=lowk; k <=highk; k++)
    for (l=lowl; l<=highl; l++) {
      const double xdist    = k-ui;
      const double ydist    = l-uj;
      const double xdist_sq = xdist*xdist;
      const double ydist_sq = ydist*ydist;
      const double rdist_sq = xdist_sq+ydist_sq;
      
      int aff;
      
      if (radius>1 && !WITHIN_RADIUS_SQ(xdist_sq,ydist_sq,radius_sq))
        for (aff=0; aff<num_aff; aff++) {
          /* Kill connection */
          length -= weights[aff][k-lowk][l-lowl];
          weights[aff][k-lowk][l-lowl] = -DEATH_FLAG;
        }
      
      else if (smooth && (rdist_sq >= smooth_start_sq)) {
        for (aff=0; aff<num_aff; aff++) {
          double dist_sq = sqrt(rdist_sq)-smooth_start;
          a_weight init_weight = weights[aff][k-lowk][l-lowl];
          dist_sq*=dist_sq;
          
          weights[aff][k-lowk][l-lowl] = 
            init_weight * exp( -(dist_sq/smooth_rad_sq) );
          length += weights[aff][k-lowk][l-lowl]-init_weight;
        }
      }
    }

  return length;
}



double LissomMap::dist_lat_wt(double value, double amt_of_randomness)
/* Randomness is in range 0-1 */
{
  double low, add;

  low = value * (1.0 - amt_of_randomness)  ;
  add = 2*value*amt_of_randomness*(ShuffledRandom<double>::plain_rand());

  return(low+add);
}



void LissomMap::prune()
{
  int i,j,k;
  int living_inh_connections = 0;
  int dead_inh_connections   = 0;
  int killed_inh_connections = 0;
  //int theoretical_total_inh_connections  = perows*N*lat_inh_dimension;

  const double i_death = params->get_with_default("inh_death",1.0);
  const double e_death = params->get_with_default("exc_death",1.0);
  const bool   e_kill  = (e_death > 0);
  const bool   i_kill  = (i_death > 0);
  
  if (e_kill)  exc_connections_killed = True;
  if (i_kill)  inh_connections_killed = True;

  const int max_num_inh = (2*inh_rad+1)*(2*inh_rad+1);
  const int max_num_exc = (2*exc_rad+1)*(2*exc_rad+1);
    
  for (i=0; i<perows; i++){                       /* For each row in the PE   */
    const int row = MAPROW(i);
    for(j=0; j<N; j++) {                          /* For each column          */
      
      Neuron& neuron = cortex_map[row][j];
      
      if (e_kill) {
        /* Adjust killing to account for actual number of connections */
        const int radius= exc_rad;
        const int lowk  = MAX(i-radius,0  );
        const int highk = MIN(i+radius,N-1);
        const int lowl  = MAX(j-radius,0  );
        const int highl = MIN(j+radius,N-1);
        const int num   = (highk-lowk+1)*(highl-lowl+1);
        const double death = e_death*double(max_num_exc)/double(num);
          
        for(k=0; k<lat_exc_dimension; k++) {
          l_weight& wt = neuron.lat_exc_wts[k];
          if (wt < death) wt = -DEATH_FLAG;
        }
        normalize_lat_wts(i,j,exc_rad,exc_array_width,neuron.lat_exc_wts);
      }
        
      if (i_kill) {
        /* Adjust killing to account for actual number of connections */
        const int radius= inh_rad;
        const int lowk  = MAX(i-radius,0  );
        const int highk = MIN(i+radius,N-1);
        const int lowl  = MAX(j-radius,0  );
        const int highl = MIN(j+radius,N-1);
        const int num   = (highk-lowk+1)*(highl-lowl+1);
        const double death = i_death*double(max_num_inh)/double(num);

        for(k=0; k<lat_inh_dimension; k++) {
          l_weight& wt = neuron.lat_inh_wts[k];
          if (wt < death) {
            dead_inh_connections++;
            /* count this connection as killed if it was live before */
            if (wt > 0) killed_inh_connections++;
            wt = -DEATH_FLAG;
          }
          else living_inh_connections++;
        }
        normalize_lat_wts(i,j,inh_rad,inh_array_width,neuron.lat_inh_wts);
      }
    }
  }
  
  ipc_pe_msg_delay(1.0);
  int total_inh_connections = living_inh_connections + dead_inh_connections;
  if (i_kill)
    ipc_notify(IPC_ALL,IPC_STD,"%s: Killed %g%% (%d/%d) of the inhibitory connections, leaving %d",
               name().c_str(),(double)(total_inh_connections>0 ?
                                       (100*((float)killed_inh_connections)/total_inh_connections) : 0.0),
               killed_inh_connections, total_inh_connections, living_inh_connections);
}



double LissomMap::change_lateral_exc_radius(double old_radius, double new_radius)
{
  int i,j;
  const bool   smooth       = params->get_with_default("smooth_circular_outlines",False);
  const bool   circular     = params->get_with_default("circular_lat_wts",True);
  const double smooth_trim  = params->get_with_default("smooth_circular_radius_trim",0.0);
  
  /* Rearrange and renormalize lateral weights, if needed */
  if (new_radius > old_radius) {
    ipc_notify(IPC_ONE,IPC_ERROR,"Increasing excitatory radius is not supported");
    new_radius=old_radius;
  }
  
  else if (new_radius < old_radius)
    for (i=0; i<perows; i++) {
      const int row = MAPROW(i);
      for (j=0; j<N; j++)
        reduce_lateral_radius(row,j,old_radius,new_radius,exc_array_width,
                              cortex_map[row][j].lat_exc_wts,circular,
                              smooth, smooth_trim);
    }
  
  return new_radius;
}


void LissomMap::reduce_lateral_radius(int ui, int uj,
                                      double old_radius_float, double new_radius_float,
                                      int array_width, l_weight *weights, bool circular,
                                      bool smooth, double smooth_trim)
{
  int k,l;
  const int old_radius=array_radius_int(old_radius_float);
  const int new_radius=array_radius_int(new_radius_float);
  const int lowk  = MAX(ui-new_radius,0  );
  const int highk = MIN(ui+new_radius,N-1);
  const int lowl  = MAX(uj-new_radius,0  );
  const int highl = MIN(uj+new_radius,N-1);
        
  for(k=lowk; k<= highk; k++){
    const int partial_index_new = PARTIAL_LAT_INDEX(ui,uj,k,new_radius,array_width);
    const int partial_index_old = PARTIAL_LAT_INDEX(ui,uj,k,old_radius,array_width);
    
    for(l=lowl; l<=highl; l++)
      weights[FULL_LAT_INDEX(partial_index_new,l)] = 
        weights[FULL_LAT_INDEX(partial_index_old,l)];
  }

  if (circular)
    prune_circular_lat_weights(ui,uj,new_radius, array_width, weights,
                               smooth, new_radius_float-new_radius, smooth_trim);
      
  /* Renormalize the copied weights */
  normalize_lat_wts(ui,uj,new_radius,array_width,weights);
}



void LissomMap::collect_activation_data(int dest_pe)
{
  (void)dest_pe;
#if (NPES>1)
  int i,pe;
  ipc_barrier(); /* Ensure that data is ready to copy */

  /* Currently collects initial activity from each processor */

  if (dest_pe==Uninitialized)
    for(pe=0; pe<NPEs; pe++) {
      if (!PEISME(pe))
        for(i=0; i<perows; i++)  {
          ipc_put(init_activity[MAPROW(i)], IPC_DOUBLE, N, pe);
          /* ipc_pe_msg_delay(1.0);
          ipc_notify(IPC_ALL,IPC_STD,"Copying data from row %d to PE%d",MAPROW(i),pe); */
        }
    }
  
  else if (!PEISME(dest_pe))
    for(i=0; i<perows; i++)  
      ipc_put(init_activity[MAPROW(i)], IPC_DOUBLE, N, dest_pe);
  
  ipc_barrier(); /* Ensure that data has arrived */
#endif
} 



void LissomMap::save_state(const string& filename)
{
#ifdef NO_IO
#else
  FILE *fp=NULL;

#ifdef VERIFY_WEIGHT_SAVES 
  {
    int row, i, j, k;
    for(i=0; i<perows; i++){
      row = MAPROW(i);
      for(j=0; j<N; j++) {
        Neuron&              neuron   = cortex_map[row][j];
        AfferentWeightsList& awts     = *(neuron.weights);
        Neuron&              neuron_o = Orig_map[row][j];
        AfferentWeightsList& awts_o   = *(neuron_o.weights);

        /* Copy weights from original arrays */
        neuron_o.centerx = neuron.centerx;
        neuron_o.centery = neuron.centery;
        for(k=0; k<lat_exc_dimension; k++)
          neuron_o.lat_exc_wts[k] = neuron.lat_exc_wts[k];
        for(k=0; k<lat_inh_dimension; k++)
          neuron_o.lat_inh_wts[k] = neuron.lat_inh_wts[k];

        for(int aff=0; aff<num_aff; aff++)
          for(int x=0; x<aff_array_width; x++)
            for(int y=0; y<aff_array_width; y++)
              awts_o[aff][x][y] = awts[aff][x][y];
      }
    }
  }
  verify_saved_weights(0.0);
#endif

  /* The last PE is in charge of writing to the file */
  if (AMYOUNGESTPE) {
    if ((fp=fopen(filename.c_str(), "wb+"))==NULL){
      ipc_abort(IPC_EXIT_FILE_PROBLEM,"Weights file %s couldn't be opened for saving", filename.c_str());
    }

    ipc_notify(IPC_ALL,IPC_SUMMARY,"Saving weights to file %s",filename.c_str());
  }

  
  
  /* All PEs process the data */
  LISSOMBinaryStateSaver lss(*this); // When this file becomes a class, lss should be a member variable.
  StateSaver& ss=lss;
  ss.write(fp);
  

  if (fp!=NULL)
    fclose(fp);
  
  /* Allows saving in middle of training by reloading the weights destroyed by saving */
#ifdef RELOAD_AFTER_WEIGHT_SAVE
  ipc_barrier();
  clear_weights(); /* Just to be sure it's working */
  const bool interactive = params->get_with_default("interactive",False);
  if ( !restore_state(t,filename) && !interactive)
    ipc_exit(IPC_EXIT_FILE_PROBLEM,"Problem with weights file");
#endif

#endif
}  



bool LissomMap::restore_state(const int, const string& filename)
{
  if (filename=="")
    ipc_abort(IPC_EXIT_FILE_PROBLEM,"No weight file specified");
  
  bool status=true;
#ifdef NO_IO
#else
  FILE *fp=NULL;

  if (AMPARENTPE){ /* Parent does the file reading and input generation */

    if ((fp=fopen(filename.c_str(), "rb" ))==NULL) {
      /* Could add synchronization mechanism for NPES>1 instead of aborting */
      const bool interactive = params->get_with_default("interactive",False);
      if (!interactive || NPEs>1) 
        ipc_abort(IPC_EXIT_FILE_PROBLEM,"Wt file %s couldn't be opened", filename.c_str());
      else {
        ipc_notify(IPC_ONE,IPC_ERROR,"Wt file %s couldn't be opened", filename.c_str());
        return false;
      }
    }
  }
  
  /* Ensure that all weights are re-initialized before loading because
     some, e.g. killed inhibitory weights, are not initialized during
     the reload */
  // Temporarily skips clearing if loading only the afferent weights
  const bool aff_only = params->get_with_default("load_afferent_weights_only",False);
  if (!aff_only)
    clear_weights();
  
  /* Read & distribute the network to all PEs     */
  ipc_notify(IPC_ALL,IPC_SUMMARY,"Reading weights from file %s",filename.c_str());
  LISSOMBinaryStateSaver lss(*this); // When this file becomes a class, lss should be a member variable.
  StateSaver& ss=lss;
  if(ss.read(fp)) status=false;

  if (AMPARENTPE)  fclose(fp);/* Only the parent has to close the file        */

  /*
    Set these flags for safety, assuming that there might be some
    connections killed in the file (presumably should add checking for it,
    but it's almost always true that connections were killed, and it
    only affects the performance, and even that only slightly.)
  */
  exc_connections_killed=True;
  inh_connections_killed=True;
  

  /* Test every weight */
#ifdef VERIFY_WEIGHT_SAVES
  verify_saved_weights(0.0000001);
#endif

#endif

  return status;
}



#ifdef VERIFY_WEIGHT_SAVES 
void LissomMap::verify_saved_weights( double tolerance )
{
  int row, i, j, aff, x, y;
  int aff_errors_found=0;

  for(i=0; i<perows; i++){
    row = MAPROW(i);
    for(j=0; j<N; j++) {
      Neuron&              neuron   = cortex_map[row][j];
      AfferentWeightsList& awts     = *(neuron.weights);
      Neuron&              neuron_o = Orig_map[row][j];
      AfferentWeightsList& awts_o   = *(neuron_o.weights);
      
      if (neuron_o.centerx - neuron.centerx > tolerance)
        ipc_notify(IPC_ALL,IPC_ERROR,"row:%d j:%d       -- Centerx changed: %d -> %d",
                  row,j,(int)neuron_o.centerx,(int)neuron.centerx);
      if (neuron_o.centery - neuron.centery > tolerance)
        ipc_notify(IPC_ALL,IPC_ERROR,"row:%d j:%d       -- Centery changed: %d -> %d",
                  row,j,(int)neuron_o.centery,(int)neuron.centery);
      
      for(aff=0; aff<num_aff; aff++)
        for(x=0; x<aff_array_width; x++)
          for(y=0; y<aff_array_width; y++)
            if (awts_o[aff][x][y] - awts[aff][x][y] >tolerance) {
              if (++aff_errors_found < VERIFY_WEIGHT_SAVES_ERROR_LIMIT)
                ipc_notify(IPC_ALL,IPC_ERROR,"row:%d j:%d aff:%d x:%d y:%d -- Afferent weights differ: %e -> %e",
                          row,j,aff,x,y,
                          (double)(awts_o[aff][x][y]),
                          (double)(awts[aff][x][y]));
              else if (aff_errors_found == VERIFY_WEIGHT_SAVES_ERROR_LIMIT)
                ipc_notify(IPC_ALL,IPC_ERROR,"Reached limit for VERIFY_WEIGHT_SAVES error reporting; rest will be discarded");
            }

      verify_saved_lateral_weights(row,j, exc_rad, exc_array_width,
                                   neuron_o.lat_exc_wts,
                                   neuron.lat_exc_wts,
                                   tolerance,"lat_exc_wts");

      verify_saved_lateral_weights(row,j, inh_rad, inh_array_width,
                                   neuron_o.lat_inh_wts,
                                   neuron.lat_inh_wts,
                                   tolerance,"lat_inh_wts");
      
    }
  }
}
#endif



#ifdef VERIFY_WEIGHT_SAVES 
void LissomMap::verify_saved_lateral_weights( int i, int j, int radius, int array_width,
                                   l_weight *old_weights, l_weight *new_weights,
                                   double tolerance, const char *description)
{
  const int lowk  = MAX(i-radius,0  );
  const int highk = MIN(i+radius,N-1);
  const int lowl  = MAX(j-radius,0  );
  const int highl = MIN(j+radius,N-1);

  static int errors_found=0;
  int k,l;
  
  for(k=lowk; k<= highk; k++) {
    const int partial_idx = PARTIAL_LAT_INDEX(i,j,k,radius,array_width);
    for (l=lowl; l<=highl; l++) {
      const int idx = FULL_LAT_INDEX(partial_idx,l);
      
      if (old_weights[idx] - new_weights[idx] > tolerance) {
        if (++errors_found < VERIFY_WEIGHT_SAVES_ERROR_LIMIT)
          ipc_notify(IPC_ALL,IPC_ERROR,"i:%d j:%d k:%d l:%d -- %s differ: %e -> %e",
                    i,j,k,l,description,
                    (double)old_weights[idx],
                    (double)new_weights[idx]);
        else if (errors_found == VERIFY_WEIGHT_SAVES_ERROR_LIMIT)
          ipc_notify(IPC_ALL,IPC_ERROR,"Reached limit for VERIFY_WEIGHT_SAVES error reporting; rest will be discarded");
      }
      
    }
  }
}
#endif

void LissomMap::backproject(const string& name, const double gammaaff) {

#ifdef OBJ_LISSOM
  const bool   l2_norm = params->get_with_default("l2_norm",False);
#endif /* OBJ_LISSOM */

  // This function should be called by derived classes' 
  // backproject(noarg), which will know which map to call.
  // --- cout << " - " << name << endl;

  // 1. get named input (upstream's output). Cleaning up the 
  //    backprojection region should be done prior to calling this
  //    function, so it is omitted here.
  NeuralRegion* inp = dynamic_cast<NeuralRegion*>(get_input(name));
  if (!inp) {
    ipc_notify(IPC_ALL,IPC_ERROR,"Error: input region nonexistent");        
    return;
  }
  ActivityMatrix& inp_bp = inp->bp;
  int aff   = std::find(ISEQ(aff_names),name) - aff_names.begin();

  // 2. Loop thorough all neurons in the current NeuralRegion bp 

  int nr = output.nrows(), nc = output.ncols();
  
#ifdef OBJ_LISSOM
  double l2_sum=0.0;

  // pre-calculate l2_norm if necessary
  if (l2_norm) {
    for (int i=0; i<nr; i++)
      for (int j=0; j<nc; j++) 
        l2_sum += output[i][j]*output[i][j];
    l2_sum = 1.0/sqrt(l2_sum);
    if (l2_sum==0.0) {
      ipc_notify(IPC_ALL,IPC_ERROR,"Panic: l2_norm value is 0.0");        
    }
  }
#endif /* OBJ_LISSOM */

         
  for (int i=0; i<nr; i++)
    for (int j=0; j<nc; j++) {
  
      double cur_output_x_gammaaff = bp[i][j]/gammaaff;
#ifdef OBJ_LISSOM
      if (l2_norm) 
        cur_output_x_gammaaff *= l2_sum;
#endif /* OBJ_LISSOM */
  
      /* Array idx swapped to translt b/w correct and legacy formats. */
      Neuron&              neuron = cortex_map[j][i];
      AfferentWeightsList& awts   = *(neuron.weights);
      
      int lowk  = MAX(neuron.centery-rf_radius,0   );
      int highk = MIN(neuron.centery+rf_radius,RN-1);
      int lowl  = MAX(neuron.centerx-rf_radius,0   );
      int highl = MIN(neuron.centerx+rf_radius,RN-1);

      // backproject
      for (int k=lowk; k <=highk; k++)
        for (int l=lowl; l<=highl; l++)
          inp_bp[k][l] += cur_output_x_gammaaff*awts[aff][l-lowl][k-lowk];
        }
  // 2. end

};


#ifdef OBJ_LISSOM

void LissomMap::resp_to_residual(const string& name, const double gammaaff) {

        // This function should be called by derived classes' 
        // resp_to_residual(noarg), which will know which map to call.
        //
        // precondition: upstream "residual" must have been properly calculated
        // assumption: result will be ADDED to the current "output"
        //      so "output" should have been cleaned if this is the first time.

        // 1. propagation starts at the current "residual"

        // 2. Loop thorough all neurons in the current NeuralRegion residual
        NeuralRegion* inp = dynamic_cast<NeuralRegion*>(get_input(name));

        ActivityMatrix& inp_residual = inp->residual;
        int aff   = std::find(ISEQ(aff_names),name) - aff_names.begin();

        int nr = output.nrows(), nc = output.ncols();
        for (int i=0; i<nr; i++)
          for (int j=0; j<nc; j++) {

            // 2.1 get named weights and bounds (upstream->downstream)
            Neuron&              neuron = cortex_map[j][i];
            AfferentWeightsList& awts   = *(neuron.weights);

            int lowk  = MAX(neuron.centery-rf_radius,0   );
            int highk = MIN(neuron.centery+rf_radius,RN-1);
            int lowl  = MAX(neuron.centerx-rf_radius,0   );
            int highl = MIN(neuron.centerx+rf_radius,RN-1);

            // 2.2 ffwd residual
            double sum = 0.0;
            for (int k=lowk; k <=highk; k++)
              for (int l=lowl; l<=highl; l++)
                sum += inp_residual[k][l] * awts[aff][l-lowl][k-lowk];

            output[i][j] += sum / gammaaff;

          }
        // 2. end

}

void LissomMap::randomize_lat_wts(void)
{
  int map_row, j, k;
  const bool   allow_negative_wts = params->get_with_default("allow_negative_wts",False);

  for(map_row=0; map_row<N; map_row++){

    for(j=0; j<N; j++) {
      Neuron&              neuron = cortex_map[map_row][j];

      double sum=0.0;
      for(k=0; k<lat_exc_dimension; k++) {
        if (allow_negative_wts || neuron.lat_exc_wts[k] > -DEATH_FLAG || 1) {
                neuron.lat_exc_wts[k]   = shuffled_rand();
                sum +=  neuron.lat_exc_wts[k];
        } 
      }
      // normalize
      for(k=0; k<lat_exc_dimension; k++) {
        if (allow_negative_wts || neuron.lat_exc_wts[k] > -DEATH_FLAG || 1) {
                neuron.lat_exc_wts[k] /= sum;
        } 
      }

      sum=0.0;
      for(k=0; k<lat_inh_dimension; k++) {
        if (allow_negative_wts || neuron.lat_inh_wts[k] > -DEATH_FLAG || 1) {
                neuron.lat_inh_wts[k]   = shuffled_rand();
                sum +=  neuron.lat_inh_wts[k];
        } 
      }
      // normalize
      for(k=0; k<lat_inh_dimension; k++) {
        if (allow_negative_wts || neuron.lat_inh_wts[k] > -DEATH_FLAG || 1) {
                neuron.lat_inh_wts[k] /= sum;
        } 
      }


    }
  }
}

#endif /* OBJ_LISSOM */

---