Refactored so that classdesc and graphcode are not pushed into ecolab…

… source heirarchy.

Random number generator API refactored.
Spatial Ecolab model migration operator added

Makefile

@@ -205,7 +205,8 @@ endif
 endif
 
 $(CLASSDESC):
-	cd classdesc; $(MAKE) PREFIX=$(ECOLAB_HOME) XDR=$(XDR) install
+	cd classdesc; $(MAKE) PREFIX=$(ECOLAB_HOME) XDR=$(XDR)
+	ln -sf `pwd`/classdesc/classdesc $(CLASSDESC)
 
 src/xdr_pack.cc: classdesc/xdr_pack.cc
 	-cp $< $@
@@ -264,9 +265,10 @@ install: all-without-models
 	$(MAKE) bin/ecolab$(ECOLIBS_EXT)
 # if installing on top of a different version, remove directory completely
 	if [ -f $(PREFIX)/include/version.h ] && ! diff -q $(PREFIX)/include/version.h include/version.h; then rm -rf $(PREFIX); fi
-	mkdir -p $(PREFIX)
+	mkdir -p $(PREFIX)/bin
 	cp -r include $(PREFIX)
 	cp classdesc/*.{h,cd} $(PREFIX)/include
+	cp classdesc/classdesc $(PREFIX)/bin
 	cp graphcode/*.{h,cd} $(PREFIX)/include
 # ensure cd files are more up-to-date than their sources.
 	touch $(PREFIX)/include/*.cd

include/Makefile

@@ -47,7 +47,7 @@ include $(ECOLAB_HOME)/include/Makefile.config
 endif
 
 .SUFFIXES: .cd .d .rc $(SUFFIXES)
-VPATH=$(ECOLAB_HOME)/include
+VPATH=$(ECOLAB_HOME)/classdesc:$(ECOLAB_HOME)/classdesc/json5_parser/json5_parser:$(ECOLAB_HOME)/graphcode:$(ECOLAB_HOME)/include
 PATH:=$(ECOLAB_HOME)/bin:$(PATH):$(HOME)/usr/bin
 CLASSDESC=$(ECOLAB_HOME)/bin/classdesc
 PKG_CONFIG_PATH:=$(HOME)/usr/lib/pkgconfig:/usr/local/lib/pkgconfig:/opt/local/lib/pkgconfig:$(PKG_CONFIG_PATH)
@@ -182,7 +182,7 @@ ECOLIBS=$(ECOLAB_HOME)/lib/libecolab$(ECOLIBS_EXT).a
 
 # why is boost_thread required here?
 LIBS+=-L$(ECOLAB_HOME)/lib -lecolab$(ECOLIBS_EXT) -lboost_thread
-FLAGS+=-I. -I$(ECOLAB_HOME)/classdesc -I$(ECOLAB_HOME)/include -DHASH_TCL_hash
+FLAGS+=-I. -I$(ECOLAB_HOME)/classdesc -I$(ECOLAB_HOME)/classdesc/json5_parser/json5_parser -I$(ECOLAB_HOME)/graphcode -I$(ECOLAB_HOME)/include -DHASH_TCL_hash
 
 # The following section uses GNU Make specific syntax. If not using
 # GNU Make, edit the  FLAGS, LIBS, CC,

include/graph.h

@@ -107,8 +107,8 @@ namespace ecolab
     };
 
     /// iterator over edges
-    class const_iterator: public classdesc::poly<const_iterator_base>/*,
-                                                                       public std::iterator<std::forward_iterator_tag,Edge>*/
+    class const_iterator: public classdesc::poly<const_iterator_base> /*,
+                          public std::iterator<std::forward_iterator_tag,Edge>*/
     {
       typedef classdesc::poly<Graph::const_iterator_base> super;
       Graph::const_iterator_base& base() {return super::operator*();}

include/nauty.h

@@ -35,7 +35,7 @@ it is necessary to check they are correct.
 #define SIZEOF_LONG 4
 #define SIZEOF_LONG_LONG 8   /* 0 if nonexistent */
 #endif
-#include <nauty_sizes.h>  /* auto generate these macros */
+//#include <nauty_sizes.h>  /* auto generate these macros */
 
 #define HAVE_CONST 1    /* compiler properly supports const */
 

include/pack.h

@@ -11,8 +11,8 @@
 #ifndef PACK_H
 #define PACK_H
 #include "pack_base.h"
-#include "ref.h"
-#include <ref.cd>
+//#include "ref.h"
+//#include <ref.cd>
 #include "pack_stl.h"
 #include "pack_graph.h"
 #endif

include/random_basic.h

@@ -12,8 +12,7 @@ namespace ecolab
   {
   public:
     double rand();
-    void seed(TCL_args args) {Seed(args);}
-    void Seed(int s) {srand(s);}
+    void seed(int s) {srand(s);}
   };
 
   class gaussrand: public random_gen

include/random_gsl.h

@@ -19,17 +19,18 @@ namespace ecolab
 
   class urand: public random_gen
   {
-    void operator=(urand&);
+    void operator=(const urand&)=delete;
+    urand(const urand&)=delete;
     CLASSDESC_ACCESS(urand);
   public:
     gsl_rng *gen;
     urand(const gsl_rng_type *descr=gsl_rng_mt19937) {gen=gsl_rng_alloc(descr);}
     ~urand() {gsl_rng_free(gen);}
-    void Seed(int s) {gsl_rng_set(gen,s);}
+    void seed(int s) {gsl_rng_set(gen,s);}
     double rand();
     /* select a different uniform random generator according the GSL's rng
        string interface */
-    void Set_gen(const std::string& descr)
+    void setGen(const std::string& descr)
     {
       static const gsl_rng_type ** rngTypes=gsl_rng_types_setup();
       const gsl_rng_type **g=rngTypes;
@@ -46,7 +47,6 @@ namespace ecolab
 
   class gaussrand: public random_gen
   {
-    void operator=(gaussrand&);
     CLASSDESC_ACCESS(gaussrand);
   public:
     urand uni;

include/random_unuran.h

@@ -24,34 +24,31 @@ namespace ecolab
     UNUR_URNG *gen;
     friend class gaussrand;
     friend class unuran;
-    void operator=(const urand&);
-    urand(const urand&);
+    operator=(const urand&)=delete;
+    urand(const urand&)=delete;
     CLASSDESC_ACCESS(urand);
   public:
     urand();
     ~urand();
     void seed(int s) {unur_urng_seed(gen,s);}
-    void Seed(int s) {seed(s);} // for backwards compatibility
     double rand() {return unur_urng_sample(gen);};
-    // for backwards compatibility
-    void Set_gen(const char* descr) {set_gen(descr);} 
     /** select a different uniform random generator according the prng's
         string interface. If PRNG not avail, this routine is a nop. */
-    void set_gen(const char* descr);
+    void setGen(const char* descr);
   };
 
   /// universal non-uniform random generator
   class unuran: public random_gen
   {
     UNUR_GEN   *gen;
-    void operator=(unuran&);
+    void operator=(const unuran&)=delete;
+    unuran(const unuran&)=delete;
     CLASSDESC_ACCESS(unuran);
   public:
     urand uni;
     /* specify a random generator according to unuran's string interface */
-    void set_gen(TCL_args args) {Set_gen(args);}
     UNUR_GEN *get_gen() {return gen;}
-    void Set_gen(const char *descr)
+    void setGen(const char *descr)
     {
       if (gen) unur_free(gen);
       gen=unur_str2gen(descr); 

models/ecolab_model.cc

@@ -141,12 +141,12 @@ void PanmicticModel::condense()
 void SpatialModel::condense()
 {
   array<int> total_density(species.size());
-  for (auto& i: *this) total_density+=i->as<EcoLabCell>()->density;
+  for (auto& i: *this) total_density+=i->as<EcolabCell>()->density;
   auto mask=total_density != 0;
   size_t mask_true=sum(mask);
   if (mask.size()==mask_true) return; /* no change ! */
   ModelData::condense(mask,mask_true);
-  for (auto& i: *this) i->as<EcoLabCell>()->condense(mask, mask_true);
+  for (auto& i: *this) i->as<EcolabCell>()->condense(mask, mask_true);
 }
 
 
@@ -174,15 +174,15 @@ void SpatialModel::mutate()
   array<unsigned> num_new_sp;
   for (auto& i: *this)
     {
-      auto new_species_in_cell=i->as<EcoLabCell>()->mutate(mut_scale);
+      auto new_species_in_cell=i->as<EcolabCell>()->mutate(mut_scale);
       num_new_sp<<=new_species_in_cell.size();
       new_sp <<= new_species_in_cell;
     }
   unsigned offset=species.size(), offi=0;
   // assign 1 for all new species created in this cell, 0 for the others
   for (auto& i: *this)
     {
-      auto& density=i->as<EcoLabCell>()->density;
+      auto& density=i->as<EcolabCell>()->density;
       density<<=array<int>(new_sp.size(),0);
       for (size_t j=0; j<num_new_sp[offi]; ++j)
         {
@@ -450,19 +450,69 @@ void SpatialModel::setGrid(size_t nx, size_t ny)
 void SpatialModel::generate(unsigned niter)
 {
   if (tstep==0) makeConsistent();
-  for (auto& i: *this) i->as<EcoLabCell>()->generate(niter,*this);
+  for (auto& i: *this) i->as<EcolabCell>()->generate(niter,*this);
   tstep+=niter;
 }
 
+void SpatialModel::migrate()
+{
+  /* each cell gets a distinct random salt value */
+  for (auto& i: *this)
+    (*this)[i]->as<EcolabCell>()->salt=array_urand.rand();
+
+  // prepareNeighbours
+
+  vector<array<int> > delta(size(), array<int>(species.size(),0));
+
+  for (size_t i=0; i<size(); i++)
+    { 
+      auto& cell=*(*this)[i]->as<EcolabCell>();
+      /* loop over neighbours */ 
+      for (auto& n: *(*this)[i]) 
+	{
+          auto& nbr=*n->as<EcolabCell>();
+	  array<double> m( double(tstep-last_mig_tstep) * migration * 
+                           (nbr.density - cell.density) );
+          double salt=(*this)[i].id()<n.id()? cell.salt: nbr.salt;
+          delta[i] += array<int>(m + array<double>(m!=0.0)*(2*(m>0.0)-1)) * salt;
+	}
+    }
+  last_mig_tstep=tstep;
+  for (size_t i=0; i<size(); i++)
+    (*this)[i]->as<EcolabCell>()->density+=delta[i];
+
+  /* assertion testing that population numbers are conserved */
+#ifndef NDEBUG
+  array<int> ssum(species.size()), s(species.size()); 
+  unsigned mig=0, i;
+  for (ssum=0, i=0; i<size(); i++)
+    {
+      ssum+=delta[i];
+      for (size_t j=0; j<delta[i].size(); j++)
+	mig+=abs(delta[i][j]);
+    }
+#ifdef MPI_SUPPORT
+  MPI_Reduce(ssum.data(),s.data(),s.size(),MPI_INT,MPI_SUM,0,MPI_COMM_WORLD);
+  ssum=s;
+  int m;
+  MPI_Reduce(&mig,&m,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD);
+  mig=m;
+#endif
+  if (myid==0) assert(sum(ssum==0)==int(ssum.size()));
+#endif
+
+}
+
+
 void SpatialModel::makeConsistent()
 {
   // all cells must have same number of species. Pack out with zero density if necessary
   unsigned long nsp=0;
-  for (auto& i: *this) nsp=max(nsp,i->as<EcoLabCell>()->density.size());
+  for (auto& i: *this) nsp=max(nsp,i->as<EcolabCell>()->density.size());
 #ifdef MPI_SUPPORT
   MPI_AllReduce(MPI_IN_PLACE,&nsp,1,MPI_UNSIGNED_LONG,MPI_MAX,MPI_COMM_WORLD);
 #endif
   for (auto& i: *this)
-    i->as<EcoLabCell>()->density<<=array<int>(nsp-i->as<EcoLabCell>()->density.size(),0);
+    i->as<EcolabCell>()->density<<=array<int>(nsp-i->as<EcolabCell>()->density.size(),0);
   ModelData::makeConsistent(nsp);
 }

models/ecolab_model.h

@@ -79,15 +79,15 @@ struct PanmicticModel: public ModelData, public EcolabPoint
   array<double> lifetimes();
 };
 
-struct EcoLabCell: public EcolabPoint, public graphcode::Object<EcoLabCell> {};
+struct EcolabCell: public EcolabPoint, public graphcode::Object<EcolabCell> {};
 
-class SpatialModel: public ModelData, public graphcode::Graph<EcoLabCell>
+class SpatialModel: public ModelData, public graphcode::Graph<EcolabCell>
 {
   size_t numX=1, numY=1;
 public:
   size_t makeId(size_t x, size_t y) const {return x%numX + numX*(y%numY);}
   void setGrid(size_t nx, size_t ny);
-  EcoLabCell& cell(size_t x, size_t y) {
+  EcolabCell& cell(size_t x, size_t y) {
     return *objects[makeId(x,y)];
   }
   array<unsigned> nsp() const;
@@ -96,6 +96,6 @@ class SpatialModel: public ModelData, public graphcode::Graph<EcoLabCell>
   void generate() {generate(1);}
   void condense();
   void mutate();
-//  void migrate();
+  void migrate();
 };
 

models/migration_test.py

@@ -0,0 +1,33 @@
+from ecolab_model import spatial_ecolab as ecolab
+from random import random
+
+from ecolab import array_urand
+array_urand.seed(10)
+
+# initial number of species
+nsp=100
+
+ecolab.species(range(nsp))
+
+numX=10
+numY=10
+ecolab.setGrid(numX,numY)
+ecolab.partitionObjects()
+ecolab.cell(4,4).density(nsp*[100])
+
+ecolab.migration(nsp*[1e-5])
+
+from plot import plot
+from GUI import gui, statusBar, windows
+
+def step():
+    ecolab.migrate()
+    ecolab.tstep(ecolab.tstep()+1)
+    plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())
+    plot(f'Density(4,4)',ecolab.tstep(),ecolab.cell(4,4).density(), pens=ecolab.species())
+    plot(f'Density(2,2)',ecolab.tstep(),ecolab.cell(2,2).density(), pens=ecolab.species())
+
+gui(step)
+
+
+

models/spatial_ecolab.py

@@ -1,5 +1,9 @@
 from ecolab_model import spatial_ecolab as ecolab
 from random import random
+
+from ecolab import array_urand
+array_urand.seed(10)
+
 # initial number of species
 nsp=100
 
@@ -39,6 +43,7 @@ def randomList(num, min, max):
 def step():
     ecolab.generate()
     ecolab.mutate()
+    ecolab.migrate()
     ecolab.condense()
     nsp=len(ecolab.species)
     statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')

src/automorph.cc

@@ -8,7 +8,7 @@
 
 #include <graph.h>
 #include <netcomplexity.h>
-#include <ref.h>
+//#include <ref.h>
 #include "ecolab_epilogue.h"
 
 #include <vector>

src/ecolab.cc

@@ -31,6 +31,8 @@ namespace ecolab
 
   std::string ecolabHome=ECOLAB_HOME;
   CLASSDESC_ADD_GLOBAL(ecolabHome);
+  CLASSDESC_ADD_GLOBAL(array_urand);
+  CLASSDESC_ADD_GLOBAL(array_grand);
   CLASSDESC_DECLARE_TYPE(Plot);
   CLASSDESC_PYTHON_MODULE(ecolab);
 }