From 5b2c22287377868276c3b85dfa84fc3429c56898 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Tue, 10 Dec 2024 12:00:32 +1100 Subject: [PATCH] generate step now running, although all cells behave identitically... --- include/arrays.h | 19 ++++++++++-- include/ecolab.h | 30 ++++++++++++------ include/graph.h | 4 +-- include/sparse_mat.h | 16 ++++++++-- models/ecolab_model.cc | 67 +++++++++++++++++++++++----------------- models/ecolab_model.h | 35 +++++++++------------ models/spatial_ecolab.py | 8 ++--- src/ecolab.cc | 2 +- 8 files changed, 109 insertions(+), 72 deletions(-) diff --git a/include/arrays.h b/include/arrays.h index d66492c..70de12f 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1629,6 +1629,12 @@ namespace ecolab ~array() {release();} const Allocator& allocator() const {return m_allocator;} + const Allocator& allocator(const Allocator& alloc) { + array tmp(size(),alloc); + asg_v(tmp.data(),size(),data()); + swap(tmp); + return m_allocator; + } /// resize array to \a s elements void resize(size_t s) { @@ -1648,6 +1654,7 @@ namespace ecolab void swap(array& x) { std::swap(dt, x.dt); + std::swap(m_allocator,x.m_allocator); } T& operator[](size_t i) {assert(i typename @@ -2453,10 +2461,17 @@ namespace ecolab /// fill array with uniformly random numbers from [0,1) template void fillrand(array& x); + template void fillrand(array& x) + {array tmp(x.size()); fillrand(tmp); x=tmp;} /// fill array with gaussian random numbers from \f$N(0,1)\propto\exp(-x^2/2)\f$ template void fillgrand(array& x); + template void fillgrand(array& x) + {array tmp(x.size()); fillgrand(tmp); x=tmp;} /// fill array with exponential random numbers \f$x\leftarrow-\ln\xi,\,\xi\in [0,1)\f$ template void fillprand(array& x); + template void fillprand(array& x) + {array tmp(x.size()); fillprand(tmp); x=tmp;} + /// fill with uniform numbers drawn from [0...\a max] without replacement void fill_unique_rand(array& x, unsigned max); diff --git a/include/ecolab.h b/include/ecolab.h index 8412a50..f3f4be2 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -20,6 +20,14 @@ #include "device/Ouroboros_impl.dp.hpp" #include "InstanceDefinitions.dp.hpp" #include "device/MemoryInitialization.dp.hpp" +#if 1 //def __INTEL_LLVM_COMPLER +#warning "Engaging experimental printf" +template +void syclPrintf(Args... args) {sycl::ext::oneapi::experimental::printf(args...);} +sycl::nd_item<1> syclItem() {return sycl::ext::oneapi::this_work_item::get_nd_item<1>();} +#else +#define syclPrintf(...) +#endif #endif #include @@ -137,22 +145,19 @@ namespace ecolab { public: Ouro::SyclDesc<1>*const* desc=nullptr; - MemAllocator*const* memAlloc; - const sycl::stream* out=nullptr; + MemAllocator*const* memAlloc=nullptr; template friend class Allocator; CellAllocator()=default; - CellAllocator(Ouro::SyclDesc<1>* const& desc, MemAllocator*const& memAlloc, const sycl::stream* out): - desc(&desc), memAlloc(&memAlloc), out(out) { + CellAllocator(Ouro::SyclDesc<1>* const& desc, MemAllocator*const& memAlloc): + desc(&desc), memAlloc(&memAlloc) { } template CellAllocator(const CellAllocator& x): - desc(x.desc) {} + desc(x.desc), memAlloc(x.memAlloc) {} T* allocate(size_t sz) { if (memAlloc && *memAlloc && desc && *desc) { auto r=reinterpret_cast((*memAlloc)->malloc(**desc,sz*sizeof(T))); - if (out) (*out)<<"allocated "< CellAllocator allocator() const { - return CellAllocator(desc,memAlloc,out); + return CellAllocator(desc,memAlloc); } #else template using CellAllocator=std::allocator; @@ -170,6 +175,14 @@ namespace ecolab #endif }; + template + void printAllocator(const char* prefix, const T&) {} + + template + void printAllocator(const char* prefix, const CellBase::CellAllocator& x) + {syclPrintf("%s: allocator.desc=%p memAlloc=%p\n",prefix,x.desc,x.memAlloc);} + + class SyclGraphBase { protected: @@ -215,7 +228,6 @@ namespace ecolab h.parallel_for(sycl::nd_range<1>(range*workGroupSize, workGroupSize), [=,this](auto i) { auto idx=i.get_global_linear_id(); if (idxsize()) { - out<<"idx="<get()<<" as "<<(*this)[idx]->template as()<template as(); Ouro::SyclDesc<> desc(i,{}); cell.desc=&desc; diff --git a/include/graph.h b/include/graph.h index 0adc213..ddda759 100644 --- a/include/graph.h +++ b/include/graph.h @@ -408,8 +408,8 @@ bool GraphAdaptor::directed() const {return false;} class sparse_mat_graph: public ConcreteGraph { public: - template - const sparse_mat_graph& operator=(const sparse_mat& mat) { + template class A> + const sparse_mat_graph& operator=(const sparse_mat& mat) { clear(); std::map, double > weights; for (size_t i=0; i +using sycl::ext::oneapi::experimental::printf; + namespace ecolab { /// sparse matrix class @@ -37,6 +39,13 @@ namespace ecolab diag=x.diag; val=x.val; row=x.row; col=x.col; } + void setAllocators(const A& ialloc, const A& falloc) { + diag.allocator(falloc); + val.allocator(falloc); + row.allocator(ialloc); + col.allocator(ialloc); + } + /*matrix multiplication*/ template typename array_ns::enable_if @@ -52,15 +61,18 @@ namespace ecolab { auto alloc=array_ns::makeAllocator(x.allocator()); array_ns::array r(alloc); - assert(row.size()==col.size() && row.size()==val.size()); + assert(row.size()==col.size() && row.size()==val.size()); if (diag.size()>0) { assert(diag.size()==x.size()); - r = diag*x; + r = diag*x; + assert(r.size()==x.size()); } else r.resize(rowsz,0); + assert(r.size()>max(row) && x.size()>max(col)); r[row]+=val*x[col]; + return r; } /* initialise sparse mat in a random configuration */ diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 55b2155..641a6c4 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -64,21 +64,11 @@ int EcolabPoint::ROUND(Float x) Float dum; if (x<0) x=0; if (x>std::numeric_limits::max()-1) x=std::numeric_limits::max()-1; + //syclPrintf("ROUND inner x=%g, modf=%g, rand()=%g, rand.max=%g, rand.min=%g\n",x,std::fabs(std::modf(x,&dum)),rand(),rand.max(),rand.min()); return std::fabs(std::modf(x,&dum))*(rand.max()-rand.min()) > (rand()-rand.min()) ? (int)x+1 : (int)x; } -//inline int ROUND(float x) {return ROUND(double(x));} - -// template -// inline array EcolabPoint::RoundArray(const E& x) -// { -// array r(x.size()); -// for (size_t i=0; i struct RoundArray { @@ -87,21 +77,27 @@ struct RoundArray RoundArray(P& point, const E& expr): expr(expr), point(point) {} using value_type=int; size_t size() const {return expr.size();} - int operator[](size_t i) const {return point.ROUND(expr[i]);} + int operator[](size_t i) const //{return point.ROUND(expr[i]);} + { + auto r=point.ROUND(expr[i]); + //syclPrintf("ROUND: %d, %g=%d\n",i,expr[i],r); + return r; + } }; namespace ecolab::array_ns {template struct is_expression>: public true_type {};} -template -RoundArray roundArray(P& point, const E& expr) -{return RoundArray(point,expr);} +template +template +RoundArray> EcolabPoint::roundArray(const E& expr) +{return RoundArray>(*this,expr);} template void EcolabPoint::generate(unsigned niter, const ModelData& model) -{ +{ for (unsigned i=0; i @@ -224,7 +220,7 @@ array EcolabPoint::mutate(const array& mut_scale) { array speciations; /* calculate the number of mutants each species produces */ - speciations = roundArray(*this, mut_scale * density); + speciations = roundArray(mut_scale * density); /* generate index list of old species that mutate to the new */ array new_sp = gen_index(speciations); @@ -485,7 +481,9 @@ void SpatialModel::setGrid(size_t nx, size_t ny) void SpatialModel::generate(unsigned niter) { if (tstep==0) makeConsistent(); - forAll([=,this](EcolabCell& c) {c.generate(niter,*this);}); + forAll([=,this](EcolabCell& c) { + c.generate(niter,*this); + }); tstep+=niter; } @@ -540,6 +538,27 @@ void SpatialModel::migrate() } +void ModelData::makeConsistent(size_t nsp) +{ +#ifdef SYCL_LANGUAGE_VERSION + FAlloc falloc(syclQ(),sycl::usm::alloc::shared); + species.allocator(graphcode::Allocator(syclQ(),sycl::usm::alloc::shared)); + create.allocator(falloc); + repro_rate.allocator(falloc); + mutation.allocator(falloc); + migration.allocator(falloc); + interaction.setAllocators + (graphcode::Allocator(syclQ(),sycl::usm::alloc::shared),falloc); +#endif + + if (!species.size()) + species=pcoord(nsp); + + if (!create.size()) create.resize(species.size(),0); + if (!mutation.size()) mutation.resize(species.size(),0); + if (!migration.size()) migration.resize(species.size(),0); +} + void SpatialModel::makeConsistent() { @@ -549,23 +568,13 @@ void SpatialModel::makeConsistent() #ifdef MPI_SUPPORT MPI_Allreduce(MPI_IN_PLACE,&nsp,1,MPI_UNSIGNED_LONG,MPI_MAX,MPI_COMM_WORLD); #endif - cout<c.density.size()) { array> tmp(nsp,0,c.allocator()); asg_v(tmp.data(),c.density.size(),c.density); - //(*c.out)<<"tmp.size:"<as(); - //cout< { - array density; - sparse_mat connections; + array> density; + sparse_mat connections; bool redraw(int x0, int y0, int width, int height) override; void requestRedraw() const {if (surface) surface->requestRedraw();} - void update(const array& d, const sparse_mat& c) { + void update(const array>& d, const sparse_mat& c) { density=d; connections=c; requestRedraw(); @@ -36,33 +36,25 @@ struct ConnectionPlot: public Object struct ModelData { - array species; - array create; - array repro_rate, mutation, migration; - sparse_mat interaction; + using FAlloc=graphcode::Allocator; + array> species; + array create; + array repro_rate, mutation, migration; + sparse_mat interaction; sparse_mat_graph foodweb; unsigned long long tstep=0, last_mut_tstep=0, last_mig_tstep=0; //mutation parameters float sp_sep=0.1, mut_max=0, repro_min=0, repro_max=1, odiag_min=0, odiag_max=1; - void makeConsistent(size_t nsp) - { - if (!species.size()) - species=pcoord(nsp); - - if (!create.size()) create.resize(species.size(),0); - if (!mutation.size()) mutation.resize(species.size(),0); - if (!migration.size()) migration.resize(species.size(),0); - } - + void makeConsistent(size_t nsp); void random_interaction(unsigned conn, double sigma); - void condense(const array& mask, size_t mask_true); void mutate(const array&); - double complexity() {return ::complexity(foodweb);} }; +template struct RoundArray; + /* ecolab cell */ template struct EcolabPoint: public Exclude @@ -74,8 +66,9 @@ struct EcolabPoint: public Exclude array mutate(const array&); unsigned nsp() const; ///< number of living species in this cell /// Rounding function, randomly round up or down, in the range 0..INT_MAX - int ROUND(Float x); - template array RoundArray(const E& x); + int ROUND(Float x); + template RoundArray roundArray(const E& expr); + //template array RoundArray(const E& x); Exclude rand; // random number generator }; diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py index 9660786..a5bbab5 100644 --- a/models/spatial_ecolab.py +++ b/models/spatial_ecolab.py @@ -1,6 +1,4 @@ from ecolab_model import spatial_ecolab as ecolab -# spatial_ecolab is a smart ptr to possible Device accessible memery, so must be dereferenced -#ecolab=spatial_ecolab() from random import random, seed as randomSeed @@ -28,8 +26,8 @@ def randomList(num, min, max): ecolab.species(range(nsp)) -numX=2 -numY=2 +numX=8 +numY=8 ecolab.setGrid(numX,numY) ecolab.partitionObjects() @@ -38,9 +36,7 @@ def randomList(num, min, max): for i in range(numX): for j in range(numY): - print(i,j,ecolab.cell(i,j).density()) ecolab.cell(i,j).density(nsp*[100]) -# print(i,j,ecolab.cell(i,j).density()) ecolab.repro_rate(randomList(nsp, ecolab.repro_min(), ecolab.repro_max())) ecolab.interaction.diag(randomList(nsp, -1e-3, -1e-3)) diff --git a/src/ecolab.cc b/src/ecolab.cc index 9e3444f..ccc9311 100644 --- a/src/ecolab.cc +++ b/src/ecolab.cc @@ -64,7 +64,7 @@ namespace ecolab { bool interpExiting=false; void interpExitProc(ClientData cd) {} - + #ifdef MPI_SUPPORT unsigned myid() {return graphcode::myid();} unsigned nprocs() {return graphcode::nprocs();}