Return number of extinctions from condense.

Renable migrate code.
Make makeId() work for some negative arguments.
Code to run epochs of high and low migration.

Author: highperformancecoder

models/ecolab_model.cc

@@ -200,16 +200,17 @@ void ModelData::condense(const array<bool>& mask, size_t mask_true)   /* remove
   foodweb = interaction;
 }
 
-void PanmicticModel::condense()
+unsigned PanmicticModel::condense()
 {
   auto mask=density != 0;
   size_t mask_true=sum(mask);
-  if (mask.size()==mask_true) return; /* no change ! */
+  if (mask.size()==mask_true) return 0; /* no change ! */
   ModelData::condense(mask,mask_true);
   EcolabPoint::condense(mask, mask_true);
+  return mask.size()-mask_true;
 }
 
-void SpatialModel::condense()
+unsigned SpatialModel::condense()
 {
   array<int> total_density(species.size());
   for (auto& i: *this) total_density+=i->as<EcolabCell>()->density; // TODO
@@ -221,9 +222,10 @@ void SpatialModel::condense()
   auto mask=total_density != 0;
 
   size_t mask_true=sum(mask);
-  if (mask.size()==mask_true) return; /* no change ! */
+  if (mask.size()==mask_true) return 0; /* no change ! */
   ModelData::condense(mask,mask_true);
   for (auto& i: *this) i->as<EcolabCell>()->condense(mask, mask_true);
+  return mask.size()-mask_true;
 }
 
 
@@ -633,59 +635,71 @@ void SpatialModel::generate(unsigned niter)
 
 void SpatialModel::migrate()
 {
-//  /* each cell gets a distinct random salt value */
-//  forAll([=,this](EcolabCell& c) {c.salt=c.rand();});
-//  
-//  prepareNeighbours();
-//
-////#ifdef SYCL_LANGUAGE_VERSION
-////  using ArrayAlloc=graphcode::Allocator<int>;
-////  using Array=vector<int,ArrayAlloc>;
-////  using DeltaAlloc=graphcode::Allocator<Array>;
-////  Array init(species.size(),0,ArrayAlloc(syclQ(),sycl::usm::alloc::device));
-////  vector<Array,DeltaAlloc> delta(size(), init, DeltaAlloc(syclQ(),sycl::usm::alloc::device));
-////#else
-////  vector<array<int>> delta(size(), array<int>(species.size(),0));
-////#endif
-//
+  /* each cell gets a distinct random salt value */
+  hostForAll([=,this](EcolabCell& c) {c.salt=c.rand();});
+  
+  prepareNeighbours();
+
+//#ifdef SYCL_LANGUAGE_VERSION
+//  using ArrayAlloc=graphcode::Allocator<int>;
+//  using Array=vector<int,ArrayAlloc>;
+//  using DeltaAlloc=graphcode::Allocator<Array>;
+//  Array init(species.size(),0,ArrayAlloc(syclQ(),sycl::usm::alloc::device));
+//  vector<Array,DeltaAlloc> delta(size(), init, DeltaAlloc(syclQ(),sycl::usm::alloc::device));
+//#else
+  vector<array<int>> delta(size(), array<int>(species.size(),0));
+//#endif
+
 //  for (auto& c: *this)
 //    c->as<EcolabCell>()->delta.resize(species.size());
-//
-//  groupedForAll([=,this](EcolabCell& c) {
-//    /* loop over neighbours */
-//    for (auto& n: c) 
-//      {
-//        auto& nbr=*n->as<EcolabCell>();
-//        Float salt=c.id<nbr.id? c.salt: nbr.salt;
-//        array_ns::groupAsg(species.size(), [&](size_t i) {
-//            Float m=(tstep-last_mig_tstep) * migration[i] * (nbr.density[i] - c.density[i]);
-//            c.delta[i]+=m + (m!=0.0)*(2*(m>0.0)-1) * salt;
-//        });
-//      }
-//  });
-//  last_mig_tstep=tstep;
-//  syncThreads();
-//  groupedForAll([=,this](EcolabCell& c) {
-//    c.density+=c.delta;
-//  });
+
+  hostForAll([&,this](EcolabCell& c) {
+    /* loop over neighbours */
+    for (auto& n: c) 
+      {
+        auto& nbr=*n->as<EcolabCell>();
+        Float salt=c.id<nbr.id? c.salt: nbr.salt;
+        array<Float> m=(tstep-last_mig_tstep) * migration * (nbr.density - c.density);
+        delta[c.idx()]+=m + array<Float>((m!=0.0)*(2*(m>0.0)-1)) * salt;
+      }
+  });
+
+  array<int> ssum(species.size(),0);
+  hostForAll([&,this](EcolabCell& c) {
+    c.density+=delta[c.idx()];
+#if !defined(NDEBUG)
+#pragma omp critical
+    ssum+=delta[c.idx()];
+#endif
+    });
+  last_mig_tstep=tstep;
 
   /* assertion testing that population numbers are conserved */
 #if !defined(NDEBUG) && !defined(SYCL_LANGUAGE_VERSION)
-  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
+  array<int> s(species.size()); 
   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 (sum(ssum==0)!=int(ssum.size()))
+    {
+      for (size_t i=0; i<ssum.size(); ++i)
+        if (ssum[i])
+          {
+            cout<<"species "<<i<<":"<<endl;
+            for (size_t idx=0; idx<size(); ++idx)
+              if (delta[idx][i])
+                {
+                  auto& c=*(*this)[idx]->as<EcolabCell>();
+                  cout<<"  delta "<<c.id<<"="<<delta[idx][i]<<" n="<<c.density[i]<<endl;
+                  for (auto& n: c)
+                    {
+                      auto& nbr=*n->as<EcolabCell>();
+                      cout<<"    nbr:"<<nbr.id<<"="<<nbr.density[i]<<endl;
+                    }
+                }                    
+          }
+    }
   if (myid()==0) assert(sum(ssum==0)==int(ssum.size()));
 #endif
 }

models/ecolab_model.h

@@ -95,7 +95,8 @@ struct PanmicticModel: public ModelData, public EcolabPoint<AllocatorBase>, publ
   void seed(unsigned x) {rand.seed(x);}
   void generate(unsigned niter);
   void generate() {generate(1);}
-  void condense();
+  /// returns number of extinctions
+  unsigned condense();
   void mutate();
   array<double> lifetimes();
 };
@@ -113,7 +114,8 @@ class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>,
   CLASSDESC_ACCESS(SpatialModel);
 public:
   static constexpr size_t log2MaxNsp=10;
-  size_t makeId(size_t x, size_t y) const {return x%numX + numX*(y%numY);}
+  // function valid for x∈(-numX,∞], y∈(-numY,∞]
+  size_t makeId(size_t x, size_t y) const {return (x+numX)%numX + numX*((y+numY)%numY);}
   void setGrid(size_t nx, size_t ny);
   EcolabCell& cell(size_t x, size_t y) {
     return *objects[makeId(x,y)];
@@ -126,7 +128,8 @@ class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>,
   void seed(unsigned x) {forAll([=](EcolabCell& cell){cell.rand.seed(x);});}
   void generate(unsigned niter);
   void generate() {generate(1);}
-  void condense();
+  /// returns number of extinctions
+  unsigned condense();
   void mutate();
   void migrate();
 };

models/spatial_ecolab.py

@@ -7,7 +7,7 @@
 print(device())
 
 # we want initialisation to be identical across all processes
-randomSeed(1)
+randomSeed(10)
 
 # we want the array operations to have a different seed across processes
 array_urand.seed(10+myid())
@@ -28,14 +28,15 @@ def randomList(num, min, max):
 
 ecolab.species(range(nsp))
 
-numX=2
-numY=2
+numX=3
+numY=3
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 
 print("initialising density")
 for i in range(numX):
     for j in range(numY):
+        #ecolab.cell(i,j).density([int(100*random()) for i in range(nsp)])
         ecolab.cell(i,j).density(nsp*[100])
 print("density initialised")
 
@@ -46,20 +47,30 @@ def randomList(num, min, max):
 ecolab.interaction.val(randomList(len(ecolab.interaction.val), ecolab.odiag_min(), ecolab.odiag_max()))
 
 ecolab.mutation(nsp*[ecolab.mut_max()])
-ecolab.migration(nsp*[1e-5])
+ecolab.migration(nsp*[1e-4])
 
 from plot import plot
 from GUI import gui, statusBar, windows
 
+epoch=2000000
+mut_factor=1000
 
+extinctions=0
 def stepImpl():
     #ecolab.setDensitiesDevice()
     ecolab.generate(100)
     ecolab.mutate()
-    #    ecolab.migrate()
-    #    ecolab.condense()
+
+    epochTs=ecolab.tstep()%epoch
+    if (epochTs==0):
+        ecolab.migration([x*mut_factor for x in ecolab.migration()])
+    if (epochTs==epoch//2):
+        ecolab.migration([x/mut_factor for x in ecolab.migration()])
+
+    #ecolab.migrate()
+    global extinctions
+    extinctions+=ecolab.condense()
     #ecolab.syncThreads()
-    print("tstep=",ecolab.tstep())
     #print(ecolab.nsp()())
     #ecolab.setDensitiesShared()
     #ecolab.gather()
@@ -73,15 +84,20 @@ def stepImpl():
 print(timeit('stepImpl()', globals=globals(), number=10))
 
 def step():
-    stepImpl()
+    global extinctions
+    extinctions=0
+    for i in range(epoch//10000):
+        stepImpl()
     if myid()==0:
         nsp=len(ecolab.species)
         statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
-        plot('No. species',ecolab.tstep(),nsp)
+        plot('No. species',ecolab.tstep(),nsp,200*(ecolab.tstep()%epoch<0.5*epoch))
+        #plot('No. species',ecolab.tstep(),nsp)
         plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())
-        for i in range(numX):
-            for j in range(numY):
-                plot(f'Density({i},{j})',ecolab.tstep(),ecolab.cell(i,j).density(), pens=ecolab.species())
+        plot('Extinctions',ecolab.tstep(),extinctions)
+#        for i in range(numX):
+#            for j in range(numY):
+#                plot(f'Density({i},{j})',ecolab.tstep(),ecolab.cell(i,j).density(), pens=ecolab.species())
 
 gui(step)