From bc3a6535b642243eafa95f28302763c273f0ca38 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Thu, 12 Sep 2024 17:43:04 +1000 Subject: [PATCH] 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 | 6 ++-- classdesc | 2 +- include/Makefile | 4 +-- include/graph.h | 4 +-- include/nauty.h | 2 +- include/pack.h | 4 +-- include/random_basic.h | 3 +- include/random_gsl.h | 8 ++--- include/random_unuran.h | 15 ++++------ models/ecolab_model.cc | 64 +++++++++++++++++++++++++++++++++++----- models/ecolab_model.h | 8 ++--- models/migration_test.py | 33 +++++++++++++++++++++ models/spatial_ecolab.py | 5 ++++ src/automorph.cc | 2 +- src/ecolab.cc | 2 ++ 15 files changed, 125 insertions(+), 37 deletions(-) create mode 100644 models/migration_test.py diff --git a/Makefile b/Makefile index 702492a..fa13540 100644 --- a/Makefile +++ b/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 diff --git a/classdesc b/classdesc index aa059c0..8ea0a32 160000 --- a/classdesc +++ b/classdesc @@ -1 +1 @@ -Subproject commit aa059c01ba1461114b517bdd16fcf332214e529b +Subproject commit 8ea0a324bb0ecef26d50925d1f11b45aa820a8ae diff --git a/include/Makefile b/include/Makefile index dc14a2e..f2e552e 100644 --- a/include/Makefile +++ b/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, diff --git a/include/graph.h b/include/graph.h index a8ef6cc..ce37087 100644 --- a/include/graph.h +++ b/include/graph.h @@ -107,8 +107,8 @@ namespace ecolab }; /// iterator over edges - class const_iterator: public classdesc::poly/*, - public std::iterator*/ + class const_iterator: public classdesc::poly /*, + public std::iterator*/ { typedef classdesc::poly super; Graph::const_iterator_base& base() {return super::operator*();} diff --git a/include/nauty.h b/include/nauty.h index 77122ca..b393efb 100644 --- a/include/nauty.h +++ b/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 /* auto generate these macros */ +//#include /* auto generate these macros */ #define HAVE_CONST 1 /* compiler properly supports const */ diff --git a/include/pack.h b/include/pack.h index 11bed7d..142e813 100644 --- a/include/pack.h +++ b/include/pack.h @@ -11,8 +11,8 @@ #ifndef PACK_H #define PACK_H #include "pack_base.h" -#include "ref.h" -#include +//#include "ref.h" +//#include #include "pack_stl.h" #include "pack_graph.h" #endif diff --git a/include/random_basic.h b/include/random_basic.h index 84a376b..2714f7a 100644 --- a/include/random_basic.h +++ b/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 diff --git a/include/random_gsl.h b/include/random_gsl.h index bda40b8..15766c6 100644 --- a/include/random_gsl.h +++ b/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; diff --git a/include/random_unuran.h b/include/random_unuran.h index 805c4f8..85334cb 100644 --- a/include/random_unuran.h +++ b/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); diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 6179f58..e017e22 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -141,12 +141,12 @@ void PanmicticModel::condense() void SpatialModel::condense() { array total_density(species.size()); - for (auto& i: *this) total_density+=i->as()->density; + for (auto& i: *this) total_density+=i->as()->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()->condense(mask, mask_true); + for (auto& i: *this) i->as()->condense(mask, mask_true); } @@ -174,7 +174,7 @@ void SpatialModel::mutate() array num_new_sp; for (auto& i: *this) { - auto new_species_in_cell=i->as()->mutate(mut_scale); + auto new_species_in_cell=i->as()->mutate(mut_scale); num_new_sp<<=new_species_in_cell.size(); new_sp <<= new_species_in_cell; } @@ -182,7 +182,7 @@ void SpatialModel::mutate() // assign 1 for all new species created in this cell, 0 for the others for (auto& i: *this) { - auto& density=i->as()->density; + auto& density=i->as()->density; density<<=array(new_sp.size(),0); for (size_t j=0; jas()->generate(niter,*this); + for (auto& i: *this) i->as()->generate(niter,*this); tstep+=niter; } +void SpatialModel::migrate() +{ + /* each cell gets a distinct random salt value */ + for (auto& i: *this) + (*this)[i]->as()->salt=array_urand.rand(); + + // prepareNeighbours + + vector > delta(size(), array(species.size(),0)); + + for (size_t i=0; ias(); + /* loop over neighbours */ + for (auto& n: *(*this)[i]) + { + auto& nbr=*n->as(); + array m( double(tstep-last_mig_tstep) * migration * + (nbr.density - cell.density) ); + double salt=(*this)[i].id()(m + array(m!=0.0)*(2*(m>0.0)-1)) * salt; + } + } + last_mig_tstep=tstep; + for (size_t i=0; ias()->density+=delta[i]; + + /* assertion testing that population numbers are conserved */ +#ifndef NDEBUG + array ssum(species.size()), s(species.size()); + unsigned mig=0, i; + for (ssum=0, i=0; ias()->density.size()); + for (auto& i: *this) nsp=max(nsp,i->as()->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()->density<<=array(nsp-i->as()->density.size(),0); + i->as()->density<<=array(nsp-i->as()->density.size(),0); ModelData::makeConsistent(nsp); } diff --git a/models/ecolab_model.h b/models/ecolab_model.h index b39bd01..1b5d546 100644 --- a/models/ecolab_model.h +++ b/models/ecolab_model.h @@ -79,15 +79,15 @@ struct PanmicticModel: public ModelData, public EcolabPoint array lifetimes(); }; -struct EcoLabCell: public EcolabPoint, public graphcode::Object {}; +struct EcolabCell: public EcolabPoint, public graphcode::Object {}; -class SpatialModel: public ModelData, public graphcode::Graph +class SpatialModel: public ModelData, public graphcode::Graph { 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 nsp() const; @@ -96,6 +96,6 @@ class SpatialModel: public ModelData, public graphcode::Graph void generate() {generate(1);} void condense(); void mutate(); -// void migrate(); + void migrate(); }; diff --git a/models/migration_test.py b/models/migration_test.py new file mode 100644 index 0000000..ce58a27 --- /dev/null +++ b/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) + + + diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py index 6d5ccc7..b751dc8 100644 --- a/models/spatial_ecolab.py +++ b/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}') diff --git a/src/automorph.cc b/src/automorph.cc index 1edceb9..6512b75 100644 --- a/src/automorph.cc +++ b/src/automorph.cc @@ -8,7 +8,7 @@ #include #include -#include +//#include #include "ecolab_epilogue.h" #include diff --git a/src/ecolab.cc b/src/ecolab.cc index 31f544e..0275334 100644 --- a/src/ecolab.cc +++ b/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); }