From c1ac797dfc4d366fc4ce95f32d57c0827cd35b3c Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 8 Nov 2023 17:19:12 +0100 Subject: [PATCH 01/37] Add ltrim/rtrim hoc functions to avoid some regex (#2603) --- azure-pipelines.yml | 2 +- src/ivoc/strfun.cpp | 28 +++++++++++++++++++ .../hoc_python/test_StringFunctions.py | 20 +++++++++++++ 3 files changed, 49 insertions(+), 1 deletion(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index f98807a7f2..56de36855f 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -81,7 +81,7 @@ stages: # Jobs to build OSX wheels natively - job: 'MacOSWheels' - timeoutInMinutes: 45 + timeoutInMinutes: 60 pool: vmImage: 'macOS-11' strategy: diff --git a/src/ivoc/strfun.cpp b/src/ivoc/strfun.cpp index 5e8fd7192c..f21820b2d8 100644 --- a/src/ivoc/strfun.cpp +++ b/src/ivoc/strfun.cpp @@ -68,6 +68,32 @@ static double l_tail(void*) { return double(i); } +static double l_ltrim(void*) { + std::string s(gargstr(1)); + std::string chars = " \r\n\t\f\v"; + if (ifarg(3)) { + chars = gargstr(3); + } + s.erase(0, s.find_first_not_of(chars)); + + char** ret = hoc_pgargstr(2); + hoc_assign_str(ret, s.c_str()); + return 0.; +} + +static double l_rtrim(void*) { + std::string s(gargstr(1)); + std::string chars = " \r\n\t\f\v"; + if (ifarg(3)) { + chars = gargstr(3); + } + s.erase(s.find_last_not_of(chars) + 1); + + char** ret = hoc_pgargstr(2); + hoc_assign_str(ret, s.c_str()); + return 0.; +} + static double l_left(void*) { std::string text(gargstr(1)); std::string newtext = text.substr(0, int(chkarg(2, 0, strlen(gargstr(1))))); @@ -311,6 +337,8 @@ static Member_func l_members[] = {{"substr", l_substr}, {"len", l_len}, {"head", l_head}, {"tail", l_tail}, + {"ltrim", l_ltrim}, + {"rtrim", l_rtrim}, {"right", l_right}, {"left", l_left}, {"is_name", l_is_name}, diff --git a/test/unit_tests/hoc_python/test_StringFunctions.py b/test/unit_tests/hoc_python/test_StringFunctions.py index 0cfcb44dc4..ecea51955d 100644 --- a/test/unit_tests/hoc_python/test_StringFunctions.py +++ b/test/unit_tests/hoc_python/test_StringFunctions.py @@ -55,6 +55,26 @@ def test_tail(): assert tail[0] == "" +def text_rtrim(): + text = "bar\t; \t\n" + out = h.ref("") + sf.rtrim(text, out) + assert out[0] == "bar\t;" + + sf.rtrim(text, out, " \t\n\f\v\r;") + assert out[0] == "bar" + + +def test_ltrim(): + text = " \t \n# foo" + out = h.ref("") + sf.ltrim(text, out) + assert out[0] == "# foo" + + sf.ltrim(text, out, " \t\n\f\r\v#") + assert out[0] == "foo" + + def test_right(): s = h.ref("foobarshi") sf.right(s, 6) From 9700e91d1a2efe4d49ba21fa386fc4ed316d0761 Mon Sep 17 00:00:00 2001 From: WeinaJi Date: Thu, 9 Nov 2023 14:02:31 +0100 Subject: [PATCH 02/37] Raise error when dividing vector by zero (#2605) --- src/ivoc/ivocvect.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index 31b71e5944..e662334649 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -2369,7 +2369,11 @@ static Object** v_mul(void* v1) { static Object** v_div(void* v1) { Vect* x = (Vect*) v1; if (hoc_argtype(1) == NUMBER) { - std::for_each(x->begin(), x->end(), [](double& d) { d /= *getarg(1); }); + if (*getarg(1) == 0.0) { + hoc_execerror("Vector", "Division by zero"); + } else { + std::for_each(x->begin(), x->end(), [](double& d) { d /= *getarg(1); }); + } } if (hoc_is_object_arg(1)) { Vect* y = vector_arg(1); From ae9b284134771e52a3e12066a9ef05810cf7a059 Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Sat, 11 Nov 2023 22:53:22 +0100 Subject: [PATCH 03/37] Initialize `prop_id_` in `NPyMechObj`. (#2604) * Default initialize and destroy `NPyMechObj`. When creating a `NPyMechObj` the Python C API will perform the equivalent of a `malloc(sizeof(NPyMechObj))`. This will return an invalid object (unitialized prop_id_), that must be fixed using placement new. Details: * `PyObject_New(NPyMechObj, pmech_generic_type)` does not call `tp_new`. * `tp_free` is the mental equivalent of `free`. * Creating a valid NPyMechObj object should be done via `new_mechobj()`. --------- Co-authored-by: nrnhines --- src/nrnpython/nrnpy_nrn.cpp | 33 +++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/src/nrnpython/nrnpy_nrn.cpp b/src/nrnpython/nrnpy_nrn.cpp index dfd7da6949..e1f82848bf 100644 --- a/src/nrnpython/nrnpy_nrn.cpp +++ b/src/nrnpython/nrnpy_nrn.cpp @@ -55,19 +55,21 @@ typedef struct { PyObject_HEAD NPySegObj* pyseg_; Prop* prop_; - neuron::container::non_owning_identifier_without_container prop_id_{}; + // Following cannot be initialized when NPyMechObj allocated by Python. See new_pymechobj + // wrapper. + neuron::container::non_owning_identifier_without_container prop_id_; int type_; } NPyMechObj; typedef struct { PyObject_HEAD - NPyMechObj* pymech_{}; + NPyMechObj* pymech_; } NPyMechOfSegIter; typedef struct { PyObject_HEAD - NPyMechObj* pymech_{}; - NPyDirectMechFunc* f_{}; + NPyMechObj* pymech_; + NPyDirectMechFunc* f_; } NPyMechFunc; typedef struct { @@ -257,12 +259,26 @@ static void NPyRangeVar_dealloc(NPyRangeVar* self) { static void NPyMechObj_dealloc(NPyMechObj* self) { // printf("NPyMechObj_dealloc %p %s\n", self, self->ob_type->tp_name); Py_XDECREF(self->pyseg_); + // Must manually call destructor since it was manually constructed in new_pymechobj wrapper + self->prop_id_.~non_owning_identifier_without_container(); ((PyObject*) self)->ob_type->tp_free((PyObject*) self); } +static NPyMechObj* new_pymechobj() { + NPyMechObj* m = PyObject_New(NPyMechObj, pmech_generic_type); + if (m) { + // Use "placement new" idiom since Python C allocation cannot call the initializer to start + // it as "null". So later `a = b` might segfault because copy constructor decrements the + // refcount of `a`s nonsense memory. + new (&m->prop_id_) neuron::container::non_owning_identifier_without_container; + } + + return m; +} + // Only call if p is valid static NPyMechObj* new_pymechobj(NPySegObj* pyseg, Prop* p) { - NPyMechObj* m = PyObject_New(NPyMechObj, pmech_generic_type); + NPyMechObj* m = new_pymechobj(); if (!m) { return NULL; } @@ -436,6 +452,7 @@ static PyObject* NPyMechObj_new(PyTypeObject* type, PyObject* args, PyObject* kw // printf("NPyMechObj_new %p %s\n", self, // ((PyObject*)self)->ob_type->tp_name); if (self != NULL) { + new (self) NPyMechObj; self->pyseg_ = pyseg; Py_INCREF(self->pyseg_); } @@ -1669,7 +1686,7 @@ static NPyRangeVar* rvnew(Symbol* sym, NPySecObj* sec, double x) { if (!r) { return NULL; } - r->pymech_ = PyObject_New(NPyMechObj, pmech_generic_type); + r->pymech_ = new_pymechobj(); r->pymech_->pyseg_ = PyObject_New(NPySegObj, psegment_type); r->pymech_->pyseg_->pysec_ = sec; Py_INCREF(sec); @@ -1916,7 +1933,7 @@ static PyObject* segment_getattro(NPySegObj* self, PyObject* pyname) { sym = ((NPyRangeVar*) rv)->sym_; if (ISARRAY(sym)) { NPyRangeVar* r = PyObject_New(NPyRangeVar, range_type); - r->pymech_ = PyObject_New(NPyMechObj, pmech_generic_type); + r->pymech_ = new_pymechobj(); r->pymech_->pyseg_ = self; Py_INCREF(r->pymech_->pyseg_); r->sym_ = sym; @@ -1944,7 +1961,7 @@ static PyObject* segment_getattro(NPySegObj* self, PyObject* pyname) { sym->type == RANGEVAR) { if (ISARRAY(sym)) { NPyRangeVar* r = PyObject_New(NPyRangeVar, range_type); - r->pymech_ = PyObject_New(NPyMechObj, pmech_generic_type); + r->pymech_ = new_pymechobj(); r->pymech_->pyseg_ = self; Py_INCREF(self); r->sym_ = sym; From 6f2e0871e36779cb4beede52d63790ce10e0f923 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Mon, 13 Nov 2023 06:31:53 -0500 Subject: [PATCH 04/37] Cvode version 2 cannot use IDASolve for interpolation. (#2606) --- src/nrncvode/nrndaspk.cpp | 5 +-- test/hoctests/tests/test_cvinterp.py | 60 ++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 3 deletions(-) create mode 100644 test/hoctests/tests/test_cvinterp.py diff --git a/src/nrncvode/nrndaspk.cpp b/src/nrncvode/nrndaspk.cpp index 4dbde077d3..5434cdc031 100644 --- a/src/nrncvode/nrndaspk.cpp +++ b/src/nrncvode/nrndaspk.cpp @@ -348,13 +348,12 @@ int Daspk::advance_tn(double tstop) { int Daspk::interpolate(double tt) { // printf("Daspk::interpolate %.15g\n", tt); assert(tt >= cv_->t0_ && tt <= cv_->tn_); - IDASetStopTime(mem_, tt); - int ier = IDASolve(mem_, tt, &cv_->t_, cv_->y_, yp_, IDA_NORMAL); + int ier = IDAGetSolution(mem_, tt, cv_->y_, yp_); if (ier < 0) { Printf("DASPK interpolate error\n"); return ier; } - assert(MyMath::eq(tt, cv_->t_, NetCvode::eps(cv_->t_))); + cv_->t_ = tt; // interpolation does not call res. So we have to. res_gvardt(cv_->t_, cv_->y_, yp_, delta_, cv_); // if(MyMath::eq(t, cv_->t_, NetCvode::eps(cv_->t_))) { diff --git a/test/hoctests/tests/test_cvinterp.py b/test/hoctests/tests/test_cvinterp.py new file mode 100644 index 0000000000..afcbdc2cbb --- /dev/null +++ b/test/hoctests/tests/test_cvinterp.py @@ -0,0 +1,60 @@ +# cvode and ida should interpolate correctly +# Consider a simulation with linearly increasing solution (charging capacitor) +# If during a free running large time step, one requests the value of +# the voltage in the middle of that time step, the voltage should be +# the average of the voltages at the beginning and end of the time step. +# Prior to this PR, IDA failed this test + +from neuron import h + +h.load_file("stdrun.hoc") +import math + + +def model(): + s = h.Section(name="s") + s.L = s.diam = h.sqrt(100.0 / h.PI) + ic = h.IClamp(s(0.5)) + ic.dur = 10 + ic.amp = 0.01 + return s, ic + + +def run(ida): + h.cvode_active(1) + h.cvode.use_daspk(ida) + h.v_init = 0.0 + h.run() + + +def points(vecs): + return [x for x in zip(vecs[0], vecs[1])] + + +def test(): + m = model() + vref = m[0](0.5)._ref_v + freerun = [h.Vector().record(h._ref_t), h.Vector().record(vref)] + for ida in [0, 1]: + run(ida) + n = len(freerun[0]) + std = [v.c() for v in freerun] + pts = points(freerun)[-3:-1] + midpnt = [(pts[0][i] + pts[1][i]) / 2 for i in range(2)] + trec = h.Vector([midpnt[0]]) + rec = [h.Vector().record(h._ref_t, trec), h.Vector().record(vref, trec)] + run(ida) + print(midpnt) + print(points(rec)) + + assert len(freerun[0]) == n # rec does not add to original freerun + for i, s in enumerate(std): # and rec did not affect freerun. + for j, v in enumerate(s): + assert math.isclose(v, freerun[i][j]) + + for i in range(2): + assert math.isclose(midpnt[i], rec[i][0]) + + +if __name__ == "__main__": + test() From 15bbeb25f299fe53b4b8602f690bee4825f5f41a Mon Sep 17 00:00:00 2001 From: WeinaJi Date: Wed, 15 Nov 2023 15:34:11 +0100 Subject: [PATCH 05/37] Raise errors in import3d_gui for invalide soma (#2609) --- share/lib/hoc/import3d/import3d_gui.hoc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/share/lib/hoc/import3d/import3d_gui.hoc b/share/lib/hoc/import3d/import3d_gui.hoc index acde58f828..ba3b8a0104 100755 --- a/share/lib/hoc/import3d/import3d_gui.hoc +++ b/share/lib/hoc/import3d/import3d_gui.hoc @@ -1237,6 +1237,9 @@ proc contour2centroid() {local i, j, imax, imin, ok localobj mean, pts, d, max, // minor is normal and in xy plane minor = m.getcol(3-tobj.min_ind-tobj.max_ind) minor.x[2] = 0 + if (minor.mag == 0.) { + execerror("Failed to compute soma centroid from contour.") + } minor.div(minor.mag) if (g != nil) { g.beginline(4, 3) g.line(mean.x[0], mean.x[1]) From bd9426d9647ffc8e684c48d320a5d8cb37dde3bd Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Wed, 15 Nov 2023 18:07:10 +0100 Subject: [PATCH 06/37] Revert moving copy to build time. (#2610) Some of these files might be needed at the new location during the build (not just after), due to missing dependencies, it might be better to create the files during the configure phase. --- bin/CMakeLists.txt | 3 +-- cmake/SanitizerHelper.cmake | 4 ++-- src/coreneuron/CMakeLists.txt | 12 ++++++------ 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/bin/CMakeLists.txt b/bin/CMakeLists.txt index 561c594b81..a60c5888f4 100644 --- a/bin/CMakeLists.txt +++ b/bin/CMakeLists.txt @@ -37,8 +37,7 @@ include(CMakeListsNrnMech) # nrnmech_makefile (based on coreneuron Configure templates) # ============================================================================= nrn_configure_file(nrngui bin) -cpp_cc_build_time_copy(INPUT ${CMAKE_CURRENT_SOURCE_DIR}/sortspike - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/sortspike) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/sortspike ${CMAKE_CURRENT_BINARY_DIR}/sortspike COPYONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/nrnivmodl_makefile_cmake.in ${PROJECT_BINARY_DIR}/bin/nrnmech_makefile @ONLY) string(JOIN " " NRN_PYTHON_VERSIONS_STRING ${NRN_PYTHON_VERSIONS}) diff --git a/cmake/SanitizerHelper.cmake b/cmake/SanitizerHelper.cmake index 28d3b8af65..e88a1db7ea 100644 --- a/cmake/SanitizerHelper.cmake +++ b/cmake/SanitizerHelper.cmake @@ -48,8 +48,8 @@ if(NRN_SANITIZERS) # directory foreach(sanitizer ${nrn_sanitizers}) if(EXISTS "${PROJECT_SOURCE_DIR}/.sanitizers/${sanitizer}.supp") - cpp_cc_build_time_copy(INPUT "${PROJECT_SOURCE_DIR}/.sanitizers/${sanitizer}.supp" - OUTPUT "${PROJECT_BINARY_DIR}/share/nrn/sanitizers/${sanitizer}.supp") + configure_file("${PROJECT_SOURCE_DIR}/.sanitizers/${sanitizer}.supp" + "${PROJECT_BINARY_DIR}/share/nrn/sanitizers/${sanitizer}.supp" COPYONLY) install(FILES "${PROJECT_BINARY_DIR}/share/nrn/sanitizers/${sanitizer}.supp" DESTINATION "${CMAKE_INSTALL_PREFIX}/share/nrn/sanitizers") endif() diff --git a/src/coreneuron/CMakeLists.txt b/src/coreneuron/CMakeLists.txt index 55bee86717..52b0de72d9 100644 --- a/src/coreneuron/CMakeLists.txt +++ b/src/coreneuron/CMakeLists.txt @@ -652,15 +652,15 @@ file( RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" *.h *.hpp) -cpp_cc_build_time_copy(INPUT ${CMAKE_BINARY_DIR}/generated/coreneuron/config/neuron_version.hpp - OUTPUT ${CMAKE_BINARY_DIR}/include/coreneuron/config/neuron_version.hpp) +configure_file(${CMAKE_BINARY_DIR}/generated/coreneuron/config/neuron_version.hpp + ${CMAKE_BINARY_DIR}/include/coreneuron/config/neuron_version.hpp COPYONLY) foreach(header ${main_headers}) - cpp_cc_build_time_copy(INPUT ${CMAKE_CURRENT_SOURCE_DIR}/${header} - OUTPUT ${CMAKE_BINARY_DIR}/include/coreneuron/${header}) + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${header} + ${CMAKE_BINARY_DIR}/include/coreneuron/${header} COPYONLY) endforeach() -cpp_cc_build_time_copy(INPUT ${CMAKE_CURRENT_SOURCE_DIR}/utils/profile/profiler_interface.h - OUTPUT ${CMAKE_BINARY_DIR}/include/coreneuron/nrniv/profiler_interface.h) +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/utils/profile/profiler_interface.h + ${CMAKE_BINARY_DIR}/include/coreneuron/nrniv/profiler_interface.h COPYONLY) # main program required for building special-core cpp_cc_build_time_copy(INPUT ${CMAKE_CURRENT_SOURCE_DIR}/apps/coreneuron.cpp From 7cfdc1f8efe4dd36ee6b170cb59f4ec7b48ba21f Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Thu, 16 Nov 2023 16:47:14 +0100 Subject: [PATCH 07/37] Fixup for #2609. (#2613) Improves the floating point comparison with zero. The strategy is to use a relative tolerance for the ratio `minor.mag/major.mag`. The softening constant `1e-100` is chosen such that: 1. If `minor.mag == 0.0` then the ratio is always `0.0` and no division by zero occurs. 2. The value of `major.mag + 1e-100 == major.mag` for all biophysically sensible values of `major.mag`. If both `minor.mag` and `major.mag` are zero as floating point number, then check may fail to raise a exception. --- share/lib/hoc/import3d/import3d_gui.hoc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/share/lib/hoc/import3d/import3d_gui.hoc b/share/lib/hoc/import3d/import3d_gui.hoc index ba3b8a0104..205d24a1a8 100755 --- a/share/lib/hoc/import3d/import3d_gui.hoc +++ b/share/lib/hoc/import3d/import3d_gui.hoc @@ -1237,7 +1237,7 @@ proc contour2centroid() {local i, j, imax, imin, ok localobj mean, pts, d, max, // minor is normal and in xy plane minor = m.getcol(3-tobj.min_ind-tobj.max_ind) minor.x[2] = 0 - if (minor.mag == 0.) { + if (minor.mag / (major.mag + 1e-100) < 1e-6) { execerror("Failed to compute soma centroid from contour.") } minor.div(minor.mag) From bc5ed0576598fa30c30607620473ce1effb7f2e1 Mon Sep 17 00:00:00 2001 From: Omar Awile Date: Mon, 20 Nov 2023 15:13:07 +0100 Subject: [PATCH 08/37] Simplifying mech register function (#2600) This is just a small commit to somewhat simplify the code around mechanism registration; mostly as a means for me to better understand what is happening there. --- src/nrnoc/init.cpp | 433 +++++++++++++++++++++++++-------------------- 1 file changed, 243 insertions(+), 190 deletions(-) diff --git a/src/nrnoc/init.cpp b/src/nrnoc/init.cpp index 6a67507a8e..4feddaa3df 100644 --- a/src/nrnoc/init.cpp +++ b/src/nrnoc/init.cpp @@ -167,28 +167,28 @@ void hoc_reg_watch_allocate(int type, NrnWatchAllocateFunc_t waf) { bbcore_write_t* nrn_bbcore_write_; bbcore_write_t* nrn_bbcore_read_; -void hoc_reg_bbcore_write(int type, bbcore_write_t f) { - nrn_bbcore_write_[type] = f; +void hoc_reg_bbcore_write(int mechtype, bbcore_write_t f) { + nrn_bbcore_write_[mechtype] = f; } -void hoc_reg_bbcore_read(int type, bbcore_write_t f) { - nrn_bbcore_read_[type] = f; +void hoc_reg_bbcore_read(int mechtype, bbcore_write_t f) { + nrn_bbcore_read_[mechtype] = f; } const char** nrn_nmodl_text_; -void hoc_reg_nmodl_text(int type, const char* txt) { - nrn_nmodl_text_[type] = txt; +void hoc_reg_nmodl_text(int mechtype, const char* txt) { + nrn_nmodl_text_[mechtype] = txt; } const char** nrn_nmodl_filename_; -void hoc_reg_nmodl_filename(int type, const char* filename) { - nrn_nmodl_filename_[type] = filename; +void hoc_reg_nmodl_filename(int mechtype, const char* filename) { + nrn_nmodl_filename_[mechtype] = filename; } -void add_nrn_has_net_event(int type) { +void add_nrn_has_net_event(int mechtype) { ++nrn_has_net_event_cnt_; nrn_has_net_event_ = (int*) erealloc(nrn_has_net_event_, nrn_has_net_event_cnt_ * sizeof(int)); - nrn_has_net_event_[nrn_has_net_event_cnt_ - 1] = type; + nrn_has_net_event_[nrn_has_net_event_cnt_ - 1] = mechtype; } /* values are type numbers of mechanisms which have FOR_NETCONS statement */ @@ -208,9 +208,9 @@ void add_nrn_fornetcons(int type, int indx) { short* nrn_is_artificial_; short* nrn_artcell_qindex_; -void add_nrn_artcell(int type, int qi) { - nrn_is_artificial_[type] = 1; - nrn_artcell_qindex_[type] = qi; +void add_nrn_artcell(int mechtype, int qi) { + nrn_is_artificial_[mechtype] = 1; + nrn_artcell_qindex_[mechtype] = qi; } int nrn_is_artificial(int pnttype) { @@ -453,6 +453,23 @@ void initnrn(void) { static int pointtype = 1; /* starts at 1 since 0 means not point in pnt_map*/ int n_memb_func; + +void reallocate_mech_data(int mechtype); +void initialize_memb_func(int mechtype, + nrn_cur_t cur, + nrn_jacob_t jacob, + Pvmp alloc, + nrn_state_t stat, + nrn_init_t initialize, + int vectorized); +void check_mech_version(const char** m); +std::pair count_variables_in_mechanism(const char** m2, int modltypemax); +void register_mech_vars(const char** var_buffers, + int modltypemax, + Symbol* mech_symbol, + int mechtype, + int nrnpointerindex); + /* if vectorized then thread_data_size added to it */ void nrn_register_mech_common(const char** m, Pvmp alloc, @@ -463,15 +480,135 @@ void nrn_register_mech_common(const char** m, int nrnpointerindex, /* if -1 then there are none */ int vectorized) { // initialize at first entry, it will be incremented at exit of the function - static int type = 2; /* 0 unused, 1 for cable section */ - int j, k, modltype, pindx, modltypemax; - Symbol* s; + static int mechtype = 2; /* 0 unused, 1 for cable section */ + int modltype; + int modltypemax; + Symbol* mech_symbol; const char** m2; nrn_load_name_check(m[1]); - if (type >= memb_func_size_) { - // we exhausted the allocated space in the tables for the mechanism type data - // so reallocate + + reallocate_mech_data(mechtype); + + initialize_memb_func(mechtype, cur, jacob, alloc, stat, initialize, vectorized); + + check_mech_version(m); + + mech_symbol = hoc_install(m[1], MECHANISM, 0.0, &hoc_symlist); + mech_symbol->subtype = mechtype; + memb_func[mechtype].sym = mech_symbol; + m2 = m + 2; + if (nrnpointerindex == -1) { + modltypemax = STATE; + } else { + modltypemax = NRNPOINTER; + } + auto [nvartypes, nvars] = count_variables_in_mechanism(m2, modltypemax); + mech_symbol->s_varn = nvars; + mech_symbol->u.ppsym = (Symbol**) emalloc((unsigned) (nvartypes * sizeof(Symbol*))); + + register_mech_vars(m2, modltypemax, mech_symbol, mechtype, nrnpointerindex); + ++mechtype; + n_memb_func = mechtype; + // n_memb_func has changed, so any existing NrnThread do not know about the + // new mechanism + v_structure_change = 1; +} + +void register_mech_vars(const char** var_buffers, + int modltypemax, + Symbol* mech_symbol, + int mechtype, + int nrnpointerindex) { + /* this is set up for the possiblility of overloading range variables. + We are currently not allowing this. Hence the #if. + If never overloaded then no reason for list of symbols for each mechanism. + */ + /* the indexing is confusing because k refers to index in the range indx list + and j refers to index in mechanism list which has 0 elements to separate + nrnocCONST, DEPENDENT, and STATE */ + /* variable pointers added on at end, if they exist */ + /* allowing range variable arrays. Must extract dimension info from name[%d]*/ + /* pindx refers to index into the p-array */ + int pindx = 0; + int modltype; + int j, k; + for (j = 0, k = 0, modltype = nrnocCONST; modltype <= modltypemax; modltype++, j++) { + for (; var_buffers[j]; j++, k++) { + Symbol* var_symbol; + std::string varname(var_buffers[j]); // copy out the varname to allow modifying it + int index = 1; + unsigned nsub = 0; + auto subscript = varname.find('['); + if (subscript != varname.npos) { +#if EXTRACELLULAR + if (varname[subscript + 1] == 'N') { + index = nlayer; + } else +#endif // EXTRACELLULAR + { + index = std::stoi(varname.substr(subscript + 1)); + } + nsub = 1; + varname.erase(subscript); + } + /*SUPPRESS 624*/ + if ((var_symbol = hoc_lookup(varname.c_str()))) { +#if 0 + if (var_symbol->subtype != RANGEVAR) { + IGNORE(fprintf(stderr, CHKmes, + varname.c_str())); + } +#else // not 0 + IGNORE(fprintf(stderr, CHKmes, varname.c_str())); +#endif // not 0 + } else { + var_symbol = hoc_install(varname.c_str(), RANGEVAR, 0.0, &hoc_symlist); + var_symbol->subtype = modltype; + var_symbol->u.rng.type = mechtype; + var_symbol->cpublic = 1; + if (modltype == NRNPOINTER) { /* not in p array */ + var_symbol->u.rng.index = nrnpointerindex; + } else { + var_symbol->u.rng.index = pindx; + } + if (nsub) { + var_symbol->arayinfo = (Arrayinfo*) emalloc(sizeof(Arrayinfo) + + nsub * sizeof(int)); + var_symbol->arayinfo->a_varn = nullptr; + var_symbol->arayinfo->refcount = 1; + var_symbol->arayinfo->nsub = nsub; + var_symbol->arayinfo->sub[0] = index; + } + if (modltype == NRNPOINTER) { + if (nrn_dparam_ptr_end_[mechtype] == 0) { + nrn_dparam_ptr_start_[mechtype] = nrnpointerindex; + } + nrnpointerindex += index; + nrn_dparam_ptr_end_[mechtype] = nrnpointerindex; + } else { + pindx += index; + } + } + mech_symbol->u.ppsym[k] = var_symbol; + } + } +} + +std::pair count_variables_in_mechanism(const char** m2, int modltypemax) { + int j, k, modltype; + // count the number of variables registered in this mechanism + for (j = 0, k = 0, modltype = nrnocCONST; modltype <= modltypemax; modltype++) { + // while we have not encountered a 0 (sentinel for variable type) + while (m2[j++]) { + k++; + } + } + return std::make_pair(j, k); +} + +void reallocate_mech_data(int mechtype) { + if (mechtype >= memb_func_size_) { memb_func_size_ += 20; pointsym = (Symbol**) erealloc(pointsym, memb_func_size_ * sizeof(Symbol*)); point_process = (Point_process**) erealloc(point_process, @@ -506,7 +643,7 @@ void nrn_register_mech_common(const char** m, nrn_watch_allocate_ = (NrnWatchAllocateFunc_t*) erealloc(nrn_watch_allocate_, memb_func_size_ * sizeof(NrnWatchAllocateFunc_t)); - for (j = memb_func_size_ - 20; j < memb_func_size_; ++j) { + for (int j = memb_func_size_ - 20; j < memb_func_size_; ++j) { pnt_map[j] = 0; point_process[j] = (Point_process*) 0; pointsym[j] = (Symbol*) 0; @@ -525,41 +662,52 @@ void nrn_register_mech_common(const char** m, } nrn_mk_prop_pools(memb_func_size_); } +} + +void initialize_memb_func(int mechtype, + nrn_cur_t cur, + nrn_jacob_t jacob, + Pvmp alloc, + nrn_state_t stat, + nrn_init_t initialize, + int vectorized) { + assert(mechtype >= memb_list.size()); + memb_list.resize(mechtype + 1); + memb_func.resize(mechtype + 1); + nrn_prop_param_size_[mechtype] = 0; /* fill in later */ + nrn_prop_dparam_size_[mechtype] = 0; /* fill in later */ + nrn_dparam_ptr_start_[mechtype] = 0; /* fill in later */ + nrn_dparam_ptr_end_[mechtype] = 0; /* fill in later */ + memb_func[mechtype].current = cur; + memb_func[mechtype].jacob = jacob; + memb_func[mechtype].alloc = alloc; + memb_func[mechtype].state = stat; + memb_func[mechtype].set_initialize(initialize); + memb_func[mechtype].destructor = nullptr; + memb_func[mechtype].vectorized = vectorized ? 1 : 0; + memb_func[mechtype].thread_size_ = vectorized ? (vectorized - 1) : 0; + memb_func[mechtype].thread_mem_init_ = nullptr; + memb_func[mechtype].thread_cleanup_ = nullptr; + memb_func[mechtype].thread_table_check_ = nullptr; + memb_func[mechtype].is_point = 0; + memb_func[mechtype].hoc_mech = nullptr; + memb_func[mechtype].setdata_ = nullptr; + memb_func[mechtype].dparam_semantics = nullptr; + memb_order_[mechtype] = mechtype; + memb_func[mechtype].ode_count = nullptr; + memb_func[mechtype].ode_map = nullptr; + memb_func[mechtype].ode_spec = nullptr; + memb_func[mechtype].ode_matsol = nullptr; + memb_func[mechtype].ode_synonym = nullptr; + memb_func[mechtype].singchan_ = nullptr; +} - assert(type >= memb_list.size()); - memb_list.resize(type + 1); - memb_func.resize(type + 1); - nrn_prop_param_size_[type] = 0; /* fill in later */ - nrn_prop_dparam_size_[type] = 0; /* fill in later */ - nrn_dparam_ptr_start_[type] = 0; /* fill in later */ - nrn_dparam_ptr_end_[type] = 0; /* fill in later */ - memb_func[type].current = cur; - memb_func[type].jacob = jacob; - memb_func[type].alloc = alloc; - memb_func[type].state = stat; - memb_func[type].set_initialize(initialize); - memb_func[type].destructor = nullptr; - memb_func[type].vectorized = vectorized ? 1 : 0; - memb_func[type].thread_size_ = vectorized ? (vectorized - 1) : 0; - memb_func[type].thread_mem_init_ = nullptr; - memb_func[type].thread_cleanup_ = nullptr; - memb_func[type].thread_table_check_ = nullptr; - memb_func[type].is_point = 0; - memb_func[type].hoc_mech = nullptr; - memb_func[type].setdata_ = nullptr; - memb_func[type].dparam_semantics = nullptr; - memb_order_[type] = type; - memb_func[type].ode_count = nullptr; - memb_func[type].ode_map = nullptr; - memb_func[type].ode_spec = nullptr; - memb_func[type].ode_matsol = nullptr; - memb_func[type].ode_synonym = nullptr; - memb_func[type].singchan_ = nullptr; +void check_mech_version(const char** m) { /* as of 5.2 nmodl translates so that the version string - is the first string in m. This allows the neuron application - to determine if nmodl c files are compatible with this version - Note that internal mechanisms have a version of "0" and are - by nature consistent. + is the first string in m. This allows the neuron application + to determine if nmodl c files are compatible with this version + Note that internal mechanisms have a version of "0" and are + by nature consistent. */ /*printf("%s %s\n", m[0], m[1]);*/ @@ -586,102 +734,6 @@ It's version %s \"c\" code is incompatible with this neuron version.\n", nrn_exit(1); } } - - s = hoc_install(m[1], MECHANISM, 0.0, &hoc_symlist); - s->subtype = type; - memb_func[type].sym = s; - /* printf("%s type=%d\n", s->name, type);*/ - m2 = m + 2; - if (nrnpointerindex == -1) { - modltypemax = STATE; - } else { - modltypemax = NRNPOINTER; - } - for (k = 0, j = 0, modltype = nrnocCONST; modltype <= modltypemax; modltype++, j++) { - /*EMPTY*/ - for (; m2[j]; j++, k++) { - ; - } - } - s->s_varn = k; - s->u.ppsym = (Symbol**) emalloc((unsigned) (j * sizeof(Symbol*))); - /* this is set up for the possiblility of overloading range variables. - We are currently not allowing this. Hence the #if. - If never overloaded then no reason for list of symbols for each mechanism. - */ - /* the indexing is confusing because k refers to index in the range indx list - and j refers to index in mechanism list which has 0 elements to separate - nrnocCONST, DEPENDENT, and STATE */ - /* variable pointers added on at end, if they exist */ - /* allowing range variable arrays. Must extract dimension info from name[%d]*/ - /* pindx refers to index into the p-array */ - pindx = 0; - for (j = 0, k = 0, modltype = nrnocCONST; modltype <= modltypemax; modltype++, j++) { - for (; m2[j]; j++, k++) { - Symbol* s2; - char buf[200], *cp; - int indx; - unsigned nsub = 0; - strcpy(buf, m2[j]); /* not allowed to change constant string */ - indx = 1; - cp = strchr(buf, '['); - if (cp) { -#if EXTRACELLULAR - if (cp[1] == 'N') { - indx = nlayer; - } else -#endif // EXTRACELLULAR - { - sscanf(cp + 1, "%d", &indx); - } - nsub = 1; - *cp = '\0'; - } - /*SUPPRESS 624*/ - if ((s2 = hoc_lookup(buf))) { -#if 0 - if (s2->subtype != RANGEVAR) { - IGNORE(fprintf(stderr, CHKmes, - buf)); - } -#else // not 0 - IGNORE(fprintf(stderr, CHKmes, buf)); -#endif // not 0 - } else { - s2 = hoc_install(buf, RANGEVAR, 0.0, &hoc_symlist); - s2->subtype = modltype; - s2->u.rng.type = type; - s2->cpublic = 1; - if (modltype == NRNPOINTER) { /* not in p array */ - s2->u.rng.index = nrnpointerindex; - } else { - s2->u.rng.index = pindx; - } - if (nsub) { - s2->arayinfo = (Arrayinfo*) emalloc(sizeof(Arrayinfo) + nsub * sizeof(int)); - s2->arayinfo->a_varn = (unsigned*) 0; - s2->arayinfo->refcount = 1; - s2->arayinfo->nsub = nsub; - s2->arayinfo->sub[0] = indx; - } - if (modltype == NRNPOINTER) { - if (nrn_dparam_ptr_end_[type] == 0) { - nrn_dparam_ptr_start_[type] = nrnpointerindex; - } - nrnpointerindex += indx; - nrn_dparam_ptr_end_[type] = nrnpointerindex; - } else { - pindx += indx; - } - } - s->u.ppsym[k] = s2; - } - } - ++type; - n_memb_func = type; - // n_memb_func has changed, so any existing NrnThread do not know about the - // new mechanism - v_structure_change = 1; } void register_mech(const char** m, @@ -692,24 +744,24 @@ void register_mech(const char** m, nrn_init_t initialize, int nrnpointerindex, /* if -1 then there are none */ int vectorized) { - int type = n_memb_func; + int mechtype = n_memb_func; nrn_register_mech_common(m, alloc, cur, jacob, stat, initialize, nrnpointerindex, vectorized); if (nrnpy_reg_mech_p_) { - (*nrnpy_reg_mech_p_)(type); + (*nrnpy_reg_mech_p_)(mechtype); } } -void nrn_writes_conc(int type, int unused) { +void nrn_writes_conc(int mechtype, int unused) { static int lastion = EXTRACELL + 1; int i; for (i = n_memb_func - 2; i >= lastion; --i) { memb_order_[i + 1] = memb_order_[i]; } - memb_order_[lastion] = type; + memb_order_[lastion] = mechtype; #if 0 - printf("%s reordered from %d to %d\n", memb_func[type].sym->name, type, lastion); + printf("%s reordered from %d to %d\n", memb_func[mechtype].sym->name, mechtype, lastion); #endif // 0 - if (nrn_is_ion(type)) { + if (nrn_is_ion(mechtype)) { ++lastion; } } @@ -756,17 +808,18 @@ int dparam_semantics_to_int(std::string_view name) { } // namespace namespace neuron::mechanism::detail { -void register_data_fields(int type, +void register_data_fields(int mechtype, std::vector> const& param_info, std::vector> const& dparam_info) { - nrn_prop_param_size_[type] = param_info.size(); - nrn_prop_dparam_size_[type] = dparam_info.size(); - delete[] std::exchange(memb_func[type].dparam_semantics, nullptr); + nrn_prop_param_size_[mechtype] = param_info.size(); + nrn_prop_dparam_size_[mechtype] = dparam_info.size(); + delete[] std::exchange(memb_func[mechtype].dparam_semantics, nullptr); if (!dparam_info.empty()) { - memb_func[type].dparam_semantics = new int[dparam_info.size()]; + memb_func[mechtype].dparam_semantics = new int[dparam_info.size()]; for (auto i = 0; i < dparam_info.size(); ++i) { // dparam_info[i].first is the name of the variable, currently unused... - memb_func[type].dparam_semantics[i] = dparam_semantics_to_int(dparam_info[i].second); + memb_func[mechtype].dparam_semantics[i] = dparam_semantics_to_int( + dparam_info[i].second); } } // Translate param_info into the type we want to use internally now we're fully inside NEURON @@ -781,14 +834,14 @@ void register_data_fields(int type, // Create a per-mechanism data structure as part of the top-level // neuron::model() structure. auto& model = neuron::model(); - model.delete_mechanism(type); // e.g. extracellular can call hoc_register_prop_size multiple - // times - auto& mech_data = model.add_mechanism(type, - memb_func[type].sym->name, // the mechanism name + model.delete_mechanism(mechtype); // e.g. extracellular can call hoc_register_prop_size + // multiple times + auto& mech_data = model.add_mechanism(mechtype, + memb_func[mechtype].sym->name, // the mechanism name std::move(param_info_new)); // names and array dimensions // of double-valued // per-instance variables - memb_list[type].set_storage_pointer(&mech_data); + memb_list[mechtype].set_storage_pointer(&mech_data); } } // namespace neuron::mechanism::detail namespace neuron::mechanism { @@ -829,8 +882,8 @@ int get_field_count(int mech_type) { * create a per mechanism map of f members that can be called directly * without prior call to setmech. */ -void hoc_register_npy_direct(int type, NPyDirectMechFunc* f) { - auto& fmap = nrn_mech2funcs_map[type] = {}; +void hoc_register_npy_direct(int mechtype, NPyDirectMechFunc* f) { + auto& fmap = nrn_mech2funcs_map[mechtype] = {}; for (int i = 0; f[i].name; ++i) { fmap[f[i].name] = &f[i]; } @@ -842,9 +895,9 @@ std::unordered_map nrn_mech2funcs_map; * * Superseded by neuron::mechanism::register_data_fields. */ -void hoc_register_prop_size(int type, int psize, int dpsize) { - assert(nrn_prop_param_size_[type] == psize); - assert(nrn_prop_dparam_size_[type] == dpsize); +void hoc_register_prop_size(int mechtype, int psize, int dpsize) { + assert(nrn_prop_param_size_[mechtype] == psize); + assert(nrn_prop_dparam_size_[mechtype] == dpsize); } /** @@ -852,8 +905,8 @@ void hoc_register_prop_size(int type, int psize, int dpsize) { * * Superseded by neuron::mechanism::register_data_fields. */ -void hoc_register_dparam_semantics(int type, int ix, const char* name) { - assert(memb_func[type].dparam_semantics[ix] == dparam_semantics_to_int(name)); +void hoc_register_dparam_semantics(int mechtype, int ix, const char* name) { + assert(memb_func[mechtype].dparam_semantics[ix] == dparam_semantics_to_int(name)); } void hoc_register_cvode(int i, @@ -954,14 +1007,14 @@ void state_discontinuity(int i, double* pd, double d) { } } -void hoc_register_limits(int type, HocParmLimits* limits) { +void hoc_register_limits(int mechtype, HocParmLimits* limits) { int i; Symbol* sym; for (i = 0; limits[i].name; ++i) { sym = (Symbol*) 0; - if (type && memb_func[type].is_point) { + if (mechtype && memb_func[mechtype].is_point) { Symbol* t; - t = hoc_lookup(memb_func[type].sym->name); + t = hoc_lookup(memb_func[mechtype].sym->name); sym = hoc_table_lookup(limits[i].name, t->u.ctemplate->symtable); } if (!sym) { @@ -971,14 +1024,14 @@ void hoc_register_limits(int type, HocParmLimits* limits) { } } -void hoc_register_units(int type, HocParmUnits* units) { +void hoc_register_units(int mechtype, HocParmUnits* units) { int i; Symbol* sym; for (i = 0; units[i].name; ++i) { sym = (Symbol*) 0; - if (type && memb_func[type].is_point) { + if (mechtype && memb_func[mechtype].is_point) { Symbol* t; - t = hoc_lookup(memb_func[type].sym->name); + t = hoc_lookup(memb_func[mechtype].sym->name); sym = hoc_table_lookup(units[i].name, t->u.ctemplate->symtable); } if (!sym) { @@ -1037,14 +1090,14 @@ void _cvode_abstol(Symbol** s, double* tol, int i) { } } -void hoc_register_tolerance(int type, HocStateTolerance* tol, Symbol*** stol) { +void hoc_register_tolerance(int mechtype, HocStateTolerance* tol, Symbol*** stol) { int i; Symbol* sym; - /*printf("register tolerance for %s\n", memb_func[type].sym->name);*/ + /*printf("register tolerance for %s\n", memb_func[mechtype].sym->name);*/ for (i = 0; tol[i].name; ++i) { - if (memb_func[type].is_point) { + if (memb_func[mechtype].is_point) { Symbol* t; - t = hoc_lookup(memb_func[type].sym->name); + t = hoc_lookup(memb_func[mechtype].sym->name); sym = hoc_table_lookup(tol[i].name, t->u.ctemplate->symtable); } else { sym = hoc_lookup(tol[i].name); @@ -1052,16 +1105,16 @@ void hoc_register_tolerance(int type, HocStateTolerance* tol, Symbol*** stol) { hoc_symbol_tolerance(sym, tol[i].tolerance); } - if (memb_func[type].ode_count) { - if (auto const n = memb_func[type].ode_count(type); n > 0) { + if (memb_func[mechtype].ode_count) { + if (auto const n = memb_func[mechtype].ode_count(mechtype); n > 0) { auto* const psym = new Symbol* [n] {}; Node node{}; // dummy node node.sec_node_index_ = 0; - prop_alloc(&(node.prop), MORPHOLOGY, &node); /* in case we need diam */ - auto* p = prop_alloc(&(node.prop), type, &node); /* this and any ions */ + prop_alloc(&(node.prop), MORPHOLOGY, &node); /* in case we need diam */ + auto* p = prop_alloc(&(node.prop), mechtype, &node); /* this and any ions */ // Fill `pv` with pointers to `2*n` parameters inside `p` std::vector> pv(2 * n); - memb_func[type].ode_map(p, 0, pv.data(), pv.data() + n, nullptr, type); + memb_func[mechtype].ode_map(p, 0, pv.data(), pv.data() + n, nullptr, mechtype); // The first n elements of `pv` are "pv", the second n are "pvdot" for (int i = 0; i < n; ++i) { // `index` is the legacy index of `pv[i]` inside mechanism instance `p` @@ -1121,8 +1174,8 @@ void _nrn_setdata_reg(int i, void (*call)(Prop*)) { memb_func[i].setdata_ = call; } /* there is some question about the _extcall_thread variables, if any. */ -double nrn_call_mech_func(Symbol* s, int narg, Prop* p, int type) { - void (*call)(Prop*) = memb_func[type].setdata_; +double nrn_call_mech_func(Symbol* s, int narg, Prop* p, int mechtype) { + void (*call)(Prop*) = memb_func[mechtype].setdata_; if (call) { (*call)(p); } From 27bec3eefef22b01c69c24ac1449f6404a095d54 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 22 Nov 2023 12:20:29 +0100 Subject: [PATCH 09/37] Rewrite resources to use std::vector and update NMODL submodule (#2616) * Rewrite resources to use std::vector * Update nmodl submodule to latest master --- external/nmodl | 2 +- src/ivos/InterViews/resource.h | 6 ++-- src/ivos/resource.cpp | 55 ++++++++++++++-------------------- 3 files changed, 26 insertions(+), 37 deletions(-) diff --git a/external/nmodl b/external/nmodl index 6f6db3b6f2..6c89dc639c 160000 --- a/external/nmodl +++ b/external/nmodl @@ -1 +1 @@ -Subproject commit 6f6db3b6f25e066db46d8a5fc85a1697363b995c +Subproject commit 6c89dc639c1ee59e8cfdd69719780e0185f1ec08 diff --git a/src/ivos/InterViews/resource.h b/src/ivos/InterViews/resource.h index ffa147deca..bf487f2dfb 100644 --- a/src/ivos/InterViews/resource.h +++ b/src/ivos/InterViews/resource.h @@ -33,8 +33,8 @@ class Resource { public: - Resource(); - virtual ~Resource(); + Resource() = default; + virtual ~Resource() = default; virtual void ref() const; virtual void unref() const; @@ -54,7 +54,7 @@ class Resource { virtual void Reference() const { ref(); } virtual void Unreference() const { unref(); } private: - unsigned refcount_; + unsigned refcount_{}; private: /* prohibit default assignment */ Resource& operator =(const Resource&); diff --git a/src/ivos/resource.cpp b/src/ivos/resource.cpp index 29fe48360b..d3718e7ce0 100644 --- a/src/ivos/resource.cpp +++ b/src/ivos/resource.cpp @@ -25,24 +25,19 @@ * OF THIS SOFTWARE. */ -#include -#include +#include -declarePtrList(ResourceList,Resource) -implementPtrList(ResourceList,Resource) +#include class ResourceImpl { friend class Resource; static bool deferred_; - static ResourceList* deletes_; + static std::vector deletes_; }; bool ResourceImpl::deferred_ = false; -ResourceList* ResourceImpl::deletes_; - -Resource::Resource() { refcount_ = 0; } -Resource::~Resource() { } +std::vector ResourceImpl::deletes_; void Resource::ref() const { Resource* r = (Resource*)this; @@ -52,29 +47,26 @@ void Resource::ref() const { void Resource::unref() const { Resource* r = (Resource*)this; if (r->refcount_ != 0) { - r->refcount_ -= 1; + r->refcount_ -= 1; } if (r->refcount_ == 0) { - r->cleanup(); - delete r; + r->cleanup(); + delete r; } } void Resource::unref_deferred() const { Resource* r = (Resource*)this; if (r->refcount_ != 0) { - r->refcount_ -= 1; + r->refcount_ -= 1; } if (r->refcount_ == 0) { - r->cleanup(); - if (ResourceImpl::deferred_) { - if (ResourceImpl::deletes_ == nil) { - ResourceImpl::deletes_ = new ResourceList; - } - ResourceImpl::deletes_->append(r); - } else { - delete r; - } + r->cleanup(); + if (ResourceImpl::deferred_) { + ResourceImpl::deletes_.push_back(r); + } else { + delete r; + } } } @@ -82,41 +74,38 @@ void Resource::cleanup() { } void Resource::ref(const Resource* r) { if (r != nil) { - r->ref(); + r->ref(); } } void Resource::unref(const Resource* r) { if (r != nil) { - r->unref(); + r->unref(); } } void Resource::unref_deferred(const Resource* r) { if (r != nil) { - r->unref_deferred(); + r->unref_deferred(); } } bool Resource::defer(bool b) { bool previous = ResourceImpl::deferred_; if (b != previous) { - flush(); - ResourceImpl::deferred_ = b; + flush(); + ResourceImpl::deferred_ = b; } return previous; } void Resource::flush() { - ResourceList* list = ResourceImpl::deletes_; - if (list != nil) { bool previous = ResourceImpl::deferred_; ResourceImpl::deferred_ = false; - for (ListItr(ResourceList) i(*list); i.more(); i.next()) { - Resource* r = i.cur(); + for (auto& r: ResourceImpl::deletes_) { delete r; } - list->remove_all(); + ResourceImpl::deletes_.clear(); + ResourceImpl::deletes_.shrink_to_fit(); ResourceImpl::deferred_ = previous; - } } From 5721cc9119172d0f35782729113cde94baf7b500 Mon Sep 17 00:00:00 2001 From: Ioannis Magkanaris Date: Wed, 22 Nov 2023 16:23:31 +0100 Subject: [PATCH 10/37] NEURON 9.0+ random API changes and Random123 documentation (#2618) * Added a paragraph about NEURON random API * Added an example with netstim.mod and mentioned CoreNEURON --- docs/guide/porting_mechanisms_to_cpp.rst | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/docs/guide/porting_mechanisms_to_cpp.rst b/docs/guide/porting_mechanisms_to_cpp.rst index 5b2df166ec..5835725a31 100644 --- a/docs/guide/porting_mechanisms_to_cpp.rst +++ b/docs/guide/porting_mechanisms_to_cpp.rst @@ -324,3 +324,27 @@ If your MOD files produce these deprecation warnings, make sure that the relevant method (``vector_capacity`` in this example) is being called with an argument of the correct type (``IvocVect*``), and not a type that is implicitly converted to ``void*``. + +Legacy random number generators and API +--------------------------------------- + +Various changes have also been done in the API of NEURON functions related to random +number generators. + +First, in |neuron_with_cpp_mechanisms| parameters passed to the functions need to be +of the correct type as it was already mentioned in :ref:`function-decls-with-incorrect-types`. +The most usual consequence of that is that NEURON random API functions that were taking as +an argument a ``void*`` now need to called with a ``Rand*``. An example of the changes needed +to fix this issue is given in `182129 `_. + +Another related change is with ``scop_rand()`` function that is usually used for defining a +``URAND`` ``FUNCTION`` in mod files to return a number based on a uniform distribution from 0 to 1. +This function now takes no argument anymore. An example of this change can also be found in +`182129 `_. + +Finally, the preferred random number generator is ``Random123``. You can find more information +about that in :meth:`Random.Random123` and :ref:`Randomness in NEURON models`. An example of the +usage of ``Random123`` can be seen in `netstim.mod `_ +and its `corresponding test `_.` +Another important aspect of ``Random123`` is that it's supported in CoreNEURON as well. For more +information about this see :ref:`Random Number Generators: Random123 vs MCellRan4`. \ No newline at end of file From 1009732c7849ce160c968721a6153cad6cc8c212 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 27 Nov 2023 11:48:41 +0100 Subject: [PATCH 11/37] Rewrite observers to use std::vector (#2615) * Rewrite observers to use std::vector --- src/ivos/InterViews/observe.h | 14 +++++----- src/ivos/observe.cpp | 52 ++++++++--------------------------- src/nrnpython/CMakeLists.txt | 9 ++++-- 3 files changed, 24 insertions(+), 51 deletions(-) diff --git a/src/ivos/InterViews/observe.h b/src/ivos/InterViews/observe.h index 8a733226df..48490e717f 100755 --- a/src/ivos/InterViews/observe.h +++ b/src/ivos/InterViews/observe.h @@ -32,30 +32,30 @@ #include #include +#include class Observer; -class ObserverList; class Observable { public: - Observable(); + Observable() = default; virtual ~Observable(); virtual void attach(Observer*); virtual void detach(Observer*); virtual void notify(); private: - ObserverList* observers_; + std::vector observers_; }; class Observer { protected: - Observer(); + Observer() = default; public: - virtual ~Observer(); + virtual ~Observer() = default; - virtual void update(Observable*); - virtual void disconnect(Observable*); + virtual void update(Observable*) {}; + virtual void disconnect(Observable*) {}; }; #include diff --git a/src/ivos/observe.cpp b/src/ivos/observe.cpp index c6b09e53bb..e5e991233f 100755 --- a/src/ivos/observe.cpp +++ b/src/ivos/observe.cpp @@ -30,58 +30,28 @@ */ #include -#include - -declarePtrList(ObserverList,Observer) -implementPtrList(ObserverList,Observer) - -Observable::Observable() { - observers_ = nil; -} +#include "utils/enumerate.h" Observable::~Observable() { - ObserverList* list = observers_; - if (list != nil) { - // in case a disconnect removes items from the ObserverList - for (long i = list->count() - 1; i >= 0; --i) { - list->item(i)->disconnect(this); - if (i > list->count()) { i = list->count(); } + // in case a disconnect removes items from the observers + for (long long i = static_cast(observers_.size()) - 1; i >= 0; --i) { + observers_[i]->disconnect(this); + if (i > observers_.size()) { + i = observers_.size(); + } } - delete list; - } } void Observable::attach(Observer* o) { - ObserverList* list = observers_; - if (list == nil) { - list = new ObserverList(5); - observers_ = list; - } - list->append(o); + observers_.push_back(o); } void Observable::detach(Observer* o) { - ObserverList* list = observers_; - if (list != nil) { - for (ListUpdater(ObserverList) i(*list); i.more(); i.next()) { - if (i.cur() == o) { - i.remove_cur(); - break; - } - } - } + erase_first(observers_, o); } void Observable::notify() { - ObserverList* list = observers_; - if (list != nil) { - for (ListItr(ObserverList) i(*list); i.more(); i.next()) { - i.cur()->update(this); - } + for (auto& obs: observers_) { + obs->update(this); } } - -Observer::Observer() { } -Observer::~Observer() { } -void Observer::update(Observable*) { } -void Observer::disconnect(Observable*) { } diff --git a/src/nrnpython/CMakeLists.txt b/src/nrnpython/CMakeLists.txt index 1a2c045b72..6f114673f6 100644 --- a/src/nrnpython/CMakeLists.txt +++ b/src/nrnpython/CMakeLists.txt @@ -24,14 +24,17 @@ if(NRN_ENABLE_PYTHON_DYNAMIC) ../nrnoc ../ivoc ../nrniv - ../ivos ../gnu ../mesch ../nrnmpi ${PROJECT_BINARY_DIR}/src/nrnpython ${PROJECT_BINARY_DIR}/src/ivos - ${PROJECT_BINARY_DIR}/src/oc - ${IV_INCLUDE_DIR}) + ${PROJECT_BINARY_DIR}/src/oc) + if(NRN_ENABLE_INTERVIEWS) + list(APPEND INCLUDE_DIRS ${IV_INCLUDE_DIR}) + else() + list(APPEND INCLUDE_DIRS ../ivos) + endif() foreach(val RANGE ${NRN_PYTHON_ITERATION_LIMIT}) list(GET NRN_PYTHON_VERSIONS ${val} pyver) list(GET NRN_PYTHON_INCLUDES ${val} pyinc) From 170a3acfd90ff61eb93adad3bc84433a2cbade95 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 27 Nov 2023 18:48:40 +0100 Subject: [PATCH 12/37] Remove Lists from ivos (#2577) - After #2514 , #2615 and #2616 Lists are not used anymore in the code and can be removed --- cmake/NeuronFileLists.cmake | 2 +- src/ivoc/datapath.cpp | 1 - src/ivoc/ivocvect.cpp | 2 - src/ivoc/nrnsymdiritem.h | 2 - src/ivoc/symdir.cpp | 1 - src/ivos/OS/_defines.h | 2 - src/ivos/OS/_undefs.h | 2 - src/ivos/OS/list.h | 312 ------------------------------------ src/ivos/listimpl.cpp | 61 ------- src/nrncvode/cvodeobj.cpp | 1 - src/nrncvode/vrecitem.h | 1 - src/nrniv/finithnd.cpp | 1 - src/nrniv/glinerec.cpp | 1 - src/nrniv/kschan.cpp | 1 - src/nrniv/kssingle.cpp | 1 - src/nrniv/linmod.h | 1 - src/nrniv/matrixmap.h | 3 - src/nrniv/nrndae.h | 3 - src/nrniv/prcellstate.cpp | 1 - src/nrniv/vrecord.cpp | 1 - 20 files changed, 1 insertion(+), 399 deletions(-) delete mode 100644 src/ivos/OS/list.h delete mode 100755 src/ivos/listimpl.cpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 2b6ee3d1ee..8cd9a06553 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -427,7 +427,7 @@ set(NMODL_FILES_LIST units.cpp version.cpp) -set(IVOS_FILES_LIST listimpl.cpp observe.cpp regexp.cpp resource.cpp) +set(IVOS_FILES_LIST observe.cpp regexp.cpp resource.cpp) set(MPI_DYNAMIC_INCLUDE nrnmpi_dynam.h nrnmpi_dynam_cinc nrnmpi_dynam_wrappers.inc) diff --git a/src/ivoc/datapath.cpp b/src/ivoc/datapath.cpp index b7c65cdcdf..cee0f49d5b 100644 --- a/src/ivoc/datapath.cpp +++ b/src/ivoc/datapath.cpp @@ -2,7 +2,6 @@ #include #include #include -#include #include "hoclist.h" #if HAVE_IV #include "graph.h" diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index e662334649..cc6565bea3 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -26,8 +26,6 @@ //#include #include -#else -#include #endif #if defined(SVR4) diff --git a/src/ivoc/nrnsymdiritem.h b/src/ivoc/nrnsymdiritem.h index 7a88be6c56..bb7ab96871 100644 --- a/src/ivoc/nrnsymdiritem.h +++ b/src/ivoc/nrnsymdiritem.h @@ -3,8 +3,6 @@ // allow communication between src/ivoc/symdir.cpp and src/nrniv/pysecname.cpp -#include - class SymbolItem { public: SymbolItem(const char*, int whole_array = 0); diff --git a/src/ivoc/symdir.cpp b/src/ivoc/symdir.cpp index 9cc3d82535..5fb09510a9 100644 --- a/src/ivoc/symdir.cpp +++ b/src/ivoc/symdir.cpp @@ -1,7 +1,6 @@ #include <../../nrnconf.h> #include #include -#include #include #include "ocobserv.h" #include "utils/enumerate.h" diff --git a/src/ivos/OS/_defines.h b/src/ivos/OS/_defines.h index 4bbe5cc7c5..ca24eed4cb 100755 --- a/src/ivos/OS/_defines.h +++ b/src/ivos/OS/_defines.h @@ -1,3 +1 @@ #define u_char _lib_os(u_char) -#define List _lib_os(List) -#define PtrList _lib_os(PtrList) diff --git a/src/ivos/OS/_undefs.h b/src/ivos/OS/_undefs.h index 8998a61956..0c2e5c7b88 100755 --- a/src/ivos/OS/_undefs.h +++ b/src/ivos/OS/_undefs.h @@ -1,3 +1 @@ #undef u_char -#undef List -#undef PtrList diff --git a/src/ivos/OS/list.h b/src/ivos/OS/list.h deleted file mode 100644 index f5f09ca845..0000000000 --- a/src/ivos/OS/list.h +++ /dev/null @@ -1,312 +0,0 @@ -/* - * Copyright (c) 1987, 1988, 1989, 1990, 1991 Stanford University - * Copyright (c) 1991 Silicon Graphics, Inc. - * - * Permission to use, copy, modify, distribute, and sell this software and - * its documentation for any purpose is hereby granted without fee, provided - * that (i) the above copyright notices and this permission notice appear in - * all copies of the software and related documentation, and (ii) the names of - * Stanford and Silicon Graphics may not be used in any advertising or - * publicity relating to the software without the specific, prior written - * permission of Stanford and Silicon Graphics. - * - * THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, - * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY - * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. - * - * IN NO EVENT SHALL STANFORD OR SILICON GRAPHICS BE LIABLE FOR - * ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, - * OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, - * WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF - * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE - * OF THIS SOFTWARE. - */ - -/* - * Generic list implemented as dynamic array - */ - -#ifndef os_list_h -#define os_list_h - -#include - -extern void ListImpl_range_error(long index); -extern long ListImpl_best_new_count(long count, unsigned int size, unsigned int m = 1); - -#if 1 || defined(__STDC__) || defined(__ANSI_CPP__) -#define __ListItr(List) List##_Iterator -#define ListItr(List) __ListItr(List) -#define __ListUpdater(List) List##_Updater -#define ListUpdater(List) __ListUpdater(List) -#else -#define __ListItr(List) List/**/_Iterator -#define ListItr(List) __ListItr(List) -#define __ListUpdater(List) List/**/_Updater -#define ListUpdater(List) __ListUpdater(List) -#endif - -#define declareList(List,T) \ -class List { \ -public: \ - List(long size = 0); \ - ~List(); \ -\ - long count() const; \ - T item(long index) const; \ - T& item_ref(long index) const; \ -\ - void prepend(const T&); \ - void append(const T&); \ - void insert(long index, const T&); \ - void remove(long index); \ - void remove_all(); \ -private: \ - T* items_; \ - long size_; \ - long count_; \ - long free_; \ -}; \ -\ -inline long List::count() const { return count_; } \ -\ -inline T List::item(long index) const { \ - if (index < 0 || index >= count_) { \ - ListImpl_range_error(index); \ - } \ - long i = index < free_ ? index : index + size_ - count_; \ - return items_[i]; \ -} \ -inline T& List::item_ref(long index) const { \ - if (index < 0 || index >= count_) { \ - ListImpl_range_error(index); \ - } \ - long i = index < free_ ? index : index + size_ - count_; \ - return items_[i]; \ -} \ -\ -inline void List::append(const T& item) { insert(count_, item); } \ -inline void List::prepend(const T& item) { insert(0, item); } \ -\ -class ListItr(List) { \ -public: \ - ListItr(List)(const List&); \ -\ - bool more() const; \ - T cur() const; \ - T& cur_ref() const; \ - void next(); \ -private: \ - const List* list_; \ - long cur_; \ -}; \ -\ -inline bool ListItr(List)::more() const { return cur_ < list_->count(); } \ -inline T ListItr(List)::cur() const { return list_->item(cur_); } \ -inline T& ListItr(List)::cur_ref() const { \ - return list_->item_ref(cur_); \ -} \ -inline void ListItr(List)::next() { ++cur_; } \ -\ -class ListUpdater(List) { \ -public: \ - ListUpdater(List)(List&); \ -\ - bool more() const; \ - T cur() const; \ - T& cur_ref() const; \ - void remove_cur(); \ - void next(); \ -private: \ - List* list_; \ - long cur_; \ -}; \ -\ -inline bool ListUpdater(List)::more() const { \ - return cur_ < list_->count(); \ -} \ -inline T ListUpdater(List)::cur() const { return list_->item(cur_); } \ -inline T& ListUpdater(List)::cur_ref() const { \ - return list_->item_ref(cur_); \ -} \ -inline void ListUpdater(List)::remove_cur() { list_->remove(cur_); } \ -inline void ListUpdater(List)::next() { ++cur_; } - -/* - * Lists of pointers - * - * Don't ask me to explain the AnyPtr nonsense. C++ compilers - * have a hard time deciding between (const void*)& and const (void*&). - * Typedefs help, though still keep me guessing. - */ - -typedef void* __AnyPtr; - -declareList(__AnyPtrList,__AnyPtr) - -#define declarePtrList(PtrList,T) \ -class PtrList { \ -public: \ - PtrList(long size = 0); \ -\ - long count() const; \ - T* item(long index) const; \ -\ - void prepend(T*); \ - void append(T*); \ - void insert(long index, T*); \ - void remove(long index); \ - void remove_all(); \ -private: \ - __AnyPtrList impl_; \ -}; \ -\ -inline PtrList::PtrList(long size) : impl_(size) { } \ -inline long PtrList::count() const { return impl_.count(); } \ -inline T* PtrList::item(long index) const { return (T*)impl_.item(index); } \ -inline void PtrList::append(T* item) { insert(impl_.count(), item); } \ -inline void PtrList::prepend(T* item) { insert(0, item); } \ -inline void PtrList::remove(long index) { impl_.remove(index); } \ -inline void PtrList::remove_all() { impl_.remove_all(); } \ -\ -class ListItr(PtrList) { \ -public: \ - ListItr(PtrList)(const PtrList&); \ -\ - bool more() const; \ - T* cur() const; \ - void next(); \ -private: \ - const PtrList* list_; \ - long cur_; \ -}; \ -\ -inline bool ListItr(PtrList)::more() const { \ - return cur_ < list_->count(); \ -} \ -inline T* ListItr(PtrList)::cur() const { return list_->item(cur_); } \ -inline void ListItr(PtrList)::next() { ++cur_; } \ -\ -class ListUpdater(PtrList) { \ -public: \ - ListUpdater(PtrList)(PtrList&); \ -\ - bool more() const; \ - T* cur() const; \ - void remove_cur(); \ - void next(); \ -private: \ - PtrList* list_; \ - long cur_; \ -}; \ -\ -inline bool ListUpdater(PtrList)::more() const { \ - return cur_ < list_->count(); \ -} \ -inline T* ListUpdater(PtrList)::cur() const { return list_->item(cur_); } \ -inline void ListUpdater(PtrList)::remove_cur() { list_->remove(cur_); } \ -inline void ListUpdater(PtrList)::next() { ++cur_; } - -/* - * List implementation - */ - -#define implementList(List,T) \ -List::List(long size) { \ - if (size > 0) { \ - size_ = ListImpl_best_new_count(size, sizeof(T)); \ - items_ = new T[size_]; \ - } else { \ - size_ = 0; \ - items_ = 0; \ - } \ - count_ = 0; \ - free_ = 0; \ -} \ -\ -List::~List() { \ - delete [] items_; \ -} \ -\ -void List::insert(long index, const T& item) { \ - if (count_ == size_) { \ - long size = ListImpl_best_new_count(size_ + 1, sizeof(T), 2); \ - T* items = new T[size]; \ - if (items_ != 0) { \ - long i; \ - for (i = 0; i < free_; ++i) { \ - items[i] = items_[i]; \ - } \ - for (i = 0; i < count_ - free_; ++i) { \ - items[free_ + size - count_ + i] = \ - items_[free_ + size_ - count_ + i]; \ - } \ - delete [] items_; \ - } \ - items_ = items; \ - size_ = size; \ - } \ - if (index >= 0 && index <= count_) { \ - if (index < free_) { \ - for (long i = free_ - index - 1; i >= 0; --i) { \ - items_[index + size_ - count_ + i] = items_[index + i]; \ - } \ - } else if (index > free_) { \ - for (long i = 0; i < index - free_; ++i) { \ - items_[free_ + i] = items_[free_ + size_ - count_ + i]; \ - } \ - } \ - free_ = index + 1; \ - count_ += 1; \ - items_[index] = item; \ - } \ -} \ -\ -void List::remove(long index) { \ - if (index >= 0 && index <= count_) { \ - if (index < free_) { \ - for (long i = free_ - index - 2; i >= 0; --i) { \ - items_[size_ - count_ + index + 1 + i] = \ - items_[index + 1 + i]; \ - } \ - } else if (index > free_) { \ - for (long i = 0; i < index - free_; ++i) { \ - items_[free_ + i] = items_[free_ + size_ - count_ + i]; \ - } \ - } \ - free_ = index; \ - count_ -= 1; \ - } \ -} \ -\ -void List::remove_all() { \ - count_ = 0; \ - free_ = 0; \ -} \ -\ -ListItr(List)::ListItr(List)(const List& list) { \ - list_ = &list; \ - cur_ = 0; \ -} \ -\ -ListUpdater(List)::ListUpdater(List)(List& list) { \ - list_ = &list; \ - cur_ = 0; \ -} - -#define implementPtrList(PtrList,T) \ -void PtrList::insert(long index, T* item) { \ - const __AnyPtr p = item; \ - impl_.insert(index, p); \ -} \ -ListItr(PtrList)::ListItr(PtrList)(const PtrList& list) { \ - list_ = &list; \ - cur_ = 0; \ -} \ -\ -ListUpdater(PtrList)::ListUpdater(PtrList)(PtrList& list) { \ - list_ = &list; \ - cur_ = 0; \ -} - -#endif diff --git a/src/ivos/listimpl.cpp b/src/ivos/listimpl.cpp deleted file mode 100755 index 814cb733ea..0000000000 --- a/src/ivos/listimpl.cpp +++ /dev/null @@ -1,61 +0,0 @@ -#ifdef HAVE_CONFIG_H -#include <../../config.h> -#endif -/* - * Copyright (c) 1991 Stanford University - * Copyright (c) 1991 Silicon Graphics, Inc. - * - * Permission to use, copy, modify, distribute, and sell this software and - * its documentation for any purpose is hereby granted without fee, provided - * that (i) the above copyright notices and this permission notice appear in - * all copies of the software and related documentation, and (ii) the names of - * Stanford and Silicon Graphics may not be used in any advertising or - * publicity relating to the software without the specific, prior written - * permission of Stanford and Silicon Graphics. - * - * THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, - * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY - * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. - * - * IN NO EVENT SHALL STANFORD OR SILICON GRAPHICS BE LIABLE FOR - * ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, - * OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, - * WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF - * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE - * OF THIS SOFTWARE. - */ - -/* - * Support routines for lists. - */ - -#include -#include -#include - -implementList(__AnyPtrList,__AnyPtr) - -static long ListImpl_best_new_sizes[] = { - 48, 112, 240, 496, 1008, 2032, 4080, 8176, - 16368, 32752, 65520, 131056, 262128, 524272, 1048560, - 2097136, 4194288, 8388592, 16777200, 33554416, 67108848, - 134217712, 268435440, 536870896, 1073741808, 2147483632 -}; - -long ListImpl_best_new_count(long count, unsigned int size, unsigned int m) { - for (int i = 0; i < sizeof(ListImpl_best_new_sizes)/sizeof(long); i++) { - if (count * size < ListImpl_best_new_sizes[i]) { - return ListImpl_best_new_sizes[i] / size; - } - } - return count*m; -} - -void ListImpl_range_error(long i) { -#if defined(WIN32) - printf("internal error: list index %ld out of range\n", i); -#else - fprintf(stderr, "internal error: list index %ld out of range\n", i); -#endif - abort(); -} diff --git a/src/nrncvode/cvodeobj.cpp b/src/nrncvode/cvodeobj.cpp index e9f879fad2..0579d71964 100644 --- a/src/nrncvode/cvodeobj.cpp +++ b/src/nrncvode/cvodeobj.cpp @@ -24,7 +24,6 @@ extern int hoc_return_type_code; #include "tqueue.h" #include "mymath.h" #include "htlist.h" -#include #include #if NRN_ENABLE_THREADS diff --git a/src/nrncvode/vrecitem.h b/src/nrncvode/vrecitem.h index 6d30a0ec17..17b8b208a0 100644 --- a/src/nrncvode/vrecitem.h +++ b/src/nrncvode/vrecitem.h @@ -1,7 +1,6 @@ #ifndef vrecitem_h #define vrecitem_h -#include #include #include #include diff --git a/src/nrniv/finithnd.cpp b/src/nrniv/finithnd.cpp index c9d237f935..0bf58a2e32 100644 --- a/src/nrniv/finithnd.cpp +++ b/src/nrniv/finithnd.cpp @@ -15,7 +15,6 @@ Type 3 are at the very beginning of finitialize. ie structure changes #include #include -#include #include #include #include diff --git a/src/nrniv/glinerec.cpp b/src/nrniv/glinerec.cpp index 58cf63099e..f224b658b2 100644 --- a/src/nrniv/glinerec.cpp +++ b/src/nrniv/glinerec.cpp @@ -13,7 +13,6 @@ #undef begin #undef add -#include #include "nrnoc2iv.h" #include "vrecitem.h" #include "netcvode.h" diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index 655d4b376a..f5e2fe0ac4 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -1,7 +1,6 @@ #include <../../nrnconf.h> #include #include -#include #include #include "nrnoc2iv.h" #include "classreg.h" diff --git a/src/nrniv/kssingle.cpp b/src/nrniv/kssingle.cpp index 2883fb480b..f23267ee33 100644 --- a/src/nrniv/kssingle.cpp +++ b/src/nrniv/kssingle.cpp @@ -1,7 +1,6 @@ #include <../../nrnconf.h> #include #include -#include #include "nrnoc2iv.h" #include "kschan.h" #include "kssingle.h" diff --git a/src/nrniv/linmod.h b/src/nrniv/linmod.h index fada6c69a5..090d054456 100644 --- a/src/nrniv/linmod.h +++ b/src/nrniv/linmod.h @@ -1,7 +1,6 @@ #ifndef linmod_h #define linmod_h -#include #include "ocmatrix.h" #include "ivocvect.h" #include "nrnoc2iv.h" diff --git a/src/nrniv/matrixmap.h b/src/nrniv/matrixmap.h index 202ba1081d..4838d24af4 100644 --- a/src/nrniv/matrixmap.h +++ b/src/nrniv/matrixmap.h @@ -1,9 +1,6 @@ #ifndef matrixmap_h #define matrixmap_h -// this defines things needed by ocmatrix -#include - #include "ocmatrix.h" #include "nrnoc2iv.h" diff --git a/src/nrniv/nrndae.h b/src/nrniv/nrndae.h index 06a55cde5d..bf0817e1f8 100644 --- a/src/nrniv/nrndae.h +++ b/src/nrniv/nrndae.h @@ -11,9 +11,6 @@ #ifndef nrndae_h #define nrndae_h -// this defines things needed by ocmatrix -#include - #include "ivocvect.h" #include "matrixmap.h" diff --git a/src/nrniv/prcellstate.cpp b/src/nrniv/prcellstate.cpp index 46e4a0952b..840d22f306 100644 --- a/src/nrniv/prcellstate.cpp +++ b/src/nrniv/prcellstate.cpp @@ -4,7 +4,6 @@ #include "nrniv_mf.h" #include "netcon.h" #include -#include "OS/list.h" #include "neuron.h" #include "utils/enumerate.h" diff --git a/src/nrniv/vrecord.cpp b/src/nrniv/vrecord.cpp index 3100b1ffba..131595070c 100644 --- a/src/nrniv/vrecord.cpp +++ b/src/nrniv/vrecord.cpp @@ -1,6 +1,5 @@ #include <../../nrnconf.h> -#include #if HAVE_IV #include "ivoc.h" #endif From 2e7c56101a8553da2e873c6da08231fc604cd5ac Mon Sep 17 00:00:00 2001 From: Ioannis Magkanaris Date: Tue, 5 Dec 2023 15:29:22 +0100 Subject: [PATCH 13/37] Add CLI options in nrnivmodl for selecting NMODL binary and flags (#2626) Note: this is for development purposes and will be removed/refactored in near future. --- bin/nrnivmodl.in | 17 +++++++++++++++-- bin/nrnivmodl_makefile_cmake.in | 7 +++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/bin/nrnivmodl.in b/bin/nrnivmodl.in index 0c5379b7c4..100c173cc9 100755 --- a/bin/nrnivmodl.in +++ b/bin/nrnivmodl.in @@ -46,6 +46,8 @@ LinkCoreNEURON=false UserINCFLAGS="" UserLDFLAGS="" UserCOREFLAGS=() +UserNMODLBIN="" +UserNMODLFLAGS="" # - options come first but can be in any order. while [ "$1" ] ; do @@ -68,6 +70,17 @@ while [ "$1" ] ; do UserCOREFLAGS+=(-l "${2}") shift shift;; + -nmodl) + echo "[NMODL][warning] Code generation with NMODL is pre-alpha, lacks features and is intended only for development use" + UserNMODLBIN="$2" + UserNMODLFLAGS="--neuron $UserNMODLFLAGS" + shift + shift;; + -nmodlflags) + echo "[NMODL][warning] If sympy is enabled, NMODL needs to be found in PYTHONPATH" + UserNMODLFLAGS="$UserNMODLFLAGS $2" + shift + shift;; -*) echo "$1 unrecognized" exit 1;; @@ -159,7 +172,7 @@ for i in "${files[@]}" ; do echo "\ ./${base_name// /\\ }.cpp: ${f}.mod @printf \" -> \$(C_GREEN)NMODL\$(C_RESET) \$<\\\n\" - (cd \"$dir_name\"; @NRN_NOCMODL_SANITIZER_ENVIRONMENT_STRING@ MODLUNIT=\$(NRNUNITS) \$(NOCMODL) \"$base_name.mod\" -o \"$mdir\") + (cd \"$dir_name\"; @NRN_NOCMODL_SANITIZER_ENVIRONMENT_STRING@ MODLUNIT=\$(NRNUNITS) \$(NOCMODL) \"$base_name.mod\" -o \"$mdir\" $UserNMODLFLAGS) ./${base_name// /\\ }.o: ./${base_name// /\\ }.cpp @printf \" -> \$(C_GREEN)Compiling\$(C_RESET) \$<\\\n\" @@ -246,5 +259,5 @@ if [ "$LinkCoreNEURON" = true ] ; then fi fi -make -j 4 -f "${bindir}/nrnmech_makefile" "ROOT=${prefix}" "MODOBJFILES=$MODOBJS" "UserLDFLAGS=$UserLDFLAGS" "UserINCFLAGS=$UserINCFLAGS" "LinkCoreNEURON=$LinkCoreNEURON" special && +make -j 4 -f "${bindir}/nrnmech_makefile" "ROOT=${prefix}" "MODOBJFILES=$MODOBJS" "UserLDFLAGS=$UserLDFLAGS" "UserINCFLAGS=$UserINCFLAGS" "LinkCoreNEURON=$LinkCoreNEURON" "UserNMODLBIN=$UserNMODLBIN" "UserNMODLFLAGS=$UserNMODLFLAGS" special && echo "Successfully created $MODSUBDIR/special" diff --git a/bin/nrnivmodl_makefile_cmake.in b/bin/nrnivmodl_makefile_cmake.in index 4f6e27b42a..0bbf3d9375 100644 --- a/bin/nrnivmodl_makefile_cmake.in +++ b/bin/nrnivmodl_makefile_cmake.in @@ -4,6 +4,8 @@ # UserLDFLAGS # UserINCFLAGS # LinkCoreNEURON +# UserNMODLBIN +# UserNMODLFLAGS # Rules to build MODOBJFILES from mod files are found in makemod2c_inc # Mechanisms version are by default 0.0, but should be overriden @@ -15,6 +17,7 @@ OUTPUT = . DESTDIR = UserINCFLAGS = UserLDFLAGS = +UserNMODLBIN = # install dirs bindir := ${ROOT}/@CMAKE_INSTALL_BINDIR@ @@ -67,7 +70,11 @@ CCOMPILE = $(CC) $(CFLAGS) @NRN_COMPILE_DEFS_STRING@ @NRN_COMPILE_FLAGS_STRING@ CXX_LINK_EXE = $(CXX) $(CXXFLAGS) @CMAKE_EXE_LINKER_FLAGS@ @NRN_LINK_FLAGS_STRING@ CXX_LINK_SHARED = $(CXX) $(CXXFLAGS) @CMAKE_SHARED_LIBRARY_CREATE_CXX_FLAGS@ @CMAKE_SHARED_LIBRARY_CXX_FLAGS@ @CMAKE_SHARED_LINKER_FLAGS@ @NRN_LINK_FLAGS_STRING@ +ifeq ($(UserNMODLBIN), ) NOCMODL = $(bindir)/nocmodl +else +NOCMODL = $(UserNMODLBIN) +endif NRNUNITS = $(datadir_lib)/nrnunits.lib # File path config (internal) From f94a4652dc5c7c957e403f6e487e58ce75206fa7 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 11 Dec 2023 19:20:50 +0100 Subject: [PATCH 14/37] Fix a bug with generating argc/argv for MPI_Init and move from macos11, 12 to 12, 13 (#2634) * Move from macos to 12, 13 * Unlink brew python formulas before installing * Fix a bug with generating argc/argv for MPI_Init Co-authored-by: Pramod S Kumbhar --- .github/workflows/neuron-ci.yml | 8 ++++++-- azure-pipelines.yml | 10 +++++----- src/coreneuron/apps/main1.cpp | 5 ++++- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/.github/workflows/neuron-ci.yml b/.github/workflows/neuron-ci.yml index ecc62ed386..9472e60686 100644 --- a/.github/workflows/neuron-ci.yml +++ b/.github/workflows/neuron-ci.yml @@ -44,7 +44,7 @@ jobs: strategy: matrix: - os: [ macOS-11, ubuntu-20.04] + os: [ macOS-12, ubuntu-20.04] config: - { matrix_eval : "CC=gcc-9 CXX=g++-9", build_mode: "setuptools"} - { matrix_eval : "CC=gcc-10 CXX=g++-10", build_mode: "cmake", music: ON} @@ -78,7 +78,7 @@ jobs: cmake_option: -DNRN_ENABLE_CORENEURON=ON -DNRN_ENABLE_MPI=OFF -DCORENRN_ENABLE_OPENMP=OFF -DNRN_ENABLE_RX3D=OFF sanitizer: thread - - os: macOS-12 + - os: macOS-13 config: build_mode: cmake # TODO: investigate rxd test timeouts in this build and re-enable them @@ -96,7 +96,11 @@ jobs: - name: Install homebrew packages if: startsWith(matrix.os, 'macOS') run: | + # Unlink and re-link to prevent errors when GitHub macOS runner images + # install Python outside of brew; See actions/setup-python#577 and BlueBrain/libsonata/pull/317 + brew list -1 | grep python | while read formula; do brew unlink $formula; brew link --overwrite $formula; done brew install ccache coreutils doxygen flex bison mpich ninja xz autoconf autoconf automake libtool + # We use both for dynamic mpi in nrn brew unlink mpich brew install openmpi brew install --cask xquartz diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 56de36855f..d65673ed84 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -83,24 +83,24 @@ stages: - job: 'MacOSWheels' timeoutInMinutes: 60 pool: - vmImage: 'macOS-11' + vmImage: 'macOS-12' strategy: matrix: Python38: python.version: '3.8' python.org.version: '3.8.10' - python.installer.name: 'macosx10.9.pkg' + python.installer.name: 'macos11.pkg' Python39: python.version: '3.9' python.org.version: '3.9.13' - python.installer.name: 'macosx10.9.pkg' + python.installer.name: 'macos11.pkg' Python310: python.version: '3.10' - python.org.version: '3.10.5' + python.org.version: '3.10.11' python.installer.name: 'macos11.pkg' Python311: python.version: '3.11' - python.org.version: '3.11.1' + python.org.version: '3.11.7' python.installer.name: 'macos11.pkg' steps: diff --git a/src/coreneuron/apps/main1.cpp b/src/coreneuron/apps/main1.cpp index 6bea87dc5f..c5d5b3d2db 100644 --- a/src/coreneuron/apps/main1.cpp +++ b/src/coreneuron/apps/main1.cpp @@ -89,7 +89,7 @@ char* prepare_args(int& argc, char**& argv, std::string& args) { free(first); // now build char*argv - argv = new char*[argc]; + argv = new char*[argc + 1]; first = strdup(args.c_str()); token = strtok(first, sep); for (int i = 0; token; i++) { @@ -97,6 +97,9 @@ char* prepare_args(int& argc, char**& argv, std::string& args) { token = strtok(nullptr, sep); } + // make sure argv is terminated by NULL! + argv[argc] = nullptr; + // return actual data to be freed return first; } From b3ebf693cd973496843c1ece77448cf0324f2844 Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Mon, 11 Dec 2023 20:42:38 +0100 Subject: [PATCH 15/37] Change release tarball name from full-src-package.tar.gz to nrn-full-src-package.tar.gz (#2632) --- .github/workflows/release.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index b52bb98149..994c3c0028 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -50,7 +50,7 @@ jobs: name: Release ${{ env.REL_TAG }} prerelease: true - full-src-package: + nrn-full-src-package: runs-on: ubuntu-latest needs: tag-n-release steps: @@ -72,13 +72,13 @@ jobs: cmake -DNRN_ENABLE_PYTHON=OFF -DNRN_ENABLE_RX3D=OFF -DNRN_ENABLE_MPI=OFF -DNRN_ENABLE_INTERVIEWS=OFF ../nrn make nrnversion_h VERBOSE=1 - - name: Create full-src-package + - name: Create nrn-full-src-package id: tar run: | - tar -czvf full-src-package-${REL_TAG}.tar.gz nrn - echo "asset_file=full-src-package-${REL_TAG}.tar.gz" >> $GITHUB_OUTPUT + tar -czvf nrn-full-src-package-${REL_TAG}.tar.gz nrn + echo "asset_file=nrn-full-src-package-${REL_TAG}.tar.gz" >> $GITHUB_OUTPUT - - name: Upload full-src-package to release + - name: Upload nrn-full-src-package to release run: | gh release upload ${{ needs.tag-n-release.outputs.rel_tag }} ${{ steps.tar.outputs.asset_file }} env: From 8ba72ff15e149cb0a7dc79ac3577c40d93337167 Mon Sep 17 00:00:00 2001 From: Luc Grosheintz Date: Tue, 12 Dec 2023 06:21:49 +0100 Subject: [PATCH 16/37] Remove false comment about being generated. (#2587) --- src/coreneuron/mechanism/capac.cpp | 3 --- src/coreneuron/mechanism/eion.cpp | 2 -- 2 files changed, 5 deletions(-) diff --git a/src/coreneuron/mechanism/capac.cpp b/src/coreneuron/mechanism/capac.cpp index b6454377b6..4342803ccb 100644 --- a/src/coreneuron/mechanism/capac.cpp +++ b/src/coreneuron/mechanism/capac.cpp @@ -1,6 +1,3 @@ -/*** - THIS FILE IS AUTO GENERATED DONT MODIFY IT. - ***/ /* # ============================================================================= # Copyright (c) 2016 - 2021 Blue Brain Project/EPFL diff --git a/src/coreneuron/mechanism/eion.cpp b/src/coreneuron/mechanism/eion.cpp index d0f5e01a0c..e468ee3e01 100644 --- a/src/coreneuron/mechanism/eion.cpp +++ b/src/coreneuron/mechanism/eion.cpp @@ -6,8 +6,6 @@ # =============================================================================. */ -/// THIS FILE IS AUTO GENERATED DONT MODIFY IT. - #include #include From b087a926b943248c4774e493ed2e3d84f390fb1d Mon Sep 17 00:00:00 2001 From: nrnhines Date: Tue, 12 Dec 2023 01:25:46 -0500 Subject: [PATCH 17/37] Migrate macosx pkg build to use notarytool (instead of altool) (#2584) * Migrate to notarytool * Update notarization documentation. * Test basic functionality before notarization. --- bldnrnmacpkgcmake.sh | 49 ++++++++++++---------- docs/install/mac_pkg.md | 66 +++++++++++++++++------------ src/mac/nrn_notarize.sh | 92 ++++++++++++++++++++++------------------- 3 files changed, 116 insertions(+), 91 deletions(-) diff --git a/bldnrnmacpkgcmake.sh b/bldnrnmacpkgcmake.sh index e29c3f7738..2f0a56af8d 100644 --- a/bldnrnmacpkgcmake.sh +++ b/bldnrnmacpkgcmake.sh @@ -1,7 +1,7 @@ #!/usr/bin/env bash set -ex -default_pythons="python3.8 python3.9 python3.10" +default_pythons="python3.8 python3.9 python3.10 python3.11" # distribution built with # bash bldnrnmacpkgcmake.sh # without args, default are the pythons above. @@ -11,6 +11,10 @@ default_pythons="python3.8 python3.9 python3.10" # All pythons must have the same macos version and that will become # the MACOSX_DEPLOYMENT_TARGET +# On my machine, to build nrn-x.x.x-macosx-10.9-universal2-py-38-39-310-311.pkg +# I built my own versions of 3.8 in $HOME/soft/python3.8, and +export PATH=$HOME/soft/python3.8/bin:$PATH + CPU=`uname -m` universal="yes" # changes to "no" if any python not universal @@ -45,11 +49,15 @@ if test "$archs" != "universal2" ; then universal=no fi +# Arrgh. Recent changes to nrn source require at least 10.15! +macosver=10.15 +mac_platform=macosx-$macosver-$archs + export MACOSX_DEPLOYMENT_TARGET=$macosver echo "MACOSX_DEPLOYMENT_TARGET=$MACOSX_DEPLOYMENT_TARGET" if test "$NRN_SRC" == "" ; then - NRN_SRC=$HOME/neuron/nrn + NRN_SRC=`pwd` fi NRN_BLD=$NRN_SRC/build NSRC=$NRN_SRC @@ -64,7 +72,7 @@ mkdir -p $NRN_BLD rm -r -f $NRN_BLD/* cd $NRN_BLD -PYVS="py" # will be part of package file name, eg. py-37-38-39-310 +PYVS="py" # will be part of package file name, eg. py-38-39-310-311 pythons="" # will be arg value of NRN_PYTHON_DYNAMIC for i in $args ; do PYVER=`$i -c 'from sys import version_info as v ; print (str(v.major) + str(v.minor)); quit()'` @@ -83,7 +91,7 @@ fi # from brew install gedit). User installations are expected to have the # former and would only accidentally have the latter. -cmake .. -DCMAKE_INSTALL_PREFIX=$NRN_INSTALL \ +cmake .. -G Ninja -DCMAKE_INSTALL_PREFIX=$NRN_INSTALL \ -DNRN_ENABLE_MPI_DYNAMIC=ON \ -DPYTHON_EXECUTABLE=`which python3` -DNRN_ENABLE_PYTHON_DYNAMIC=ON \ -DNRN_PYTHON_DYNAMIC="$pythons" \ @@ -94,7 +102,7 @@ cmake .. -DCMAKE_INSTALL_PREFIX=$NRN_INSTALL \ -DCMAKE_PREFIX_PATH=/usr/X11 \ -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -make -j install +ninja install if test "$universal" = "yes" ; then _temp="`lipo -archs $NRN_INSTALL/share/nrn/demo/release/$CPU/special`" @@ -121,32 +129,29 @@ chk () { ) } +# test basic functionality +for i in $args ; do + chk $i +done + #/Applications/Packages.app from # http://s.sudre.free.fr/Software/Packages/about.html # For mac to do a productsign, need my developerID_installer.cer # and Neurondev.p12 file. To add to the keychain, double click each # of those files. By default, I added my certificates to the login keychain. -make macpkg # will sign the binaries, construct below - # mentioned PACKAGE_FILE_NAME, and request notarization from - # Apple. At the end it will print a stapling request that you - # should run manually after receiving a "success" email from - # Apple. +ninja macpkg # will sign the binaries, construct below + # mentioned PACKAGE_FILE_NAME, request notarization from + # Apple, and staple the package. -# test basic functionality -for i in $args ; do - chk $i -done - -# upload package to neuron.yale.edu -ALPHADIR='hines@neuron.yale.edu:/home/htdocs/ftp/neuron/versions/alpha' +# Copy the package to $HOME/$PACKAGE_FULL_NAME +# You should then manually upload that to github. describe="`sh $NRN_SRC/nrnversion.sh describe`" macos=macos${MACOSX_DEPLOYMENT_TARGET} PACKAGE_FULL_NAME=nrn-${describe}-${mac_platform}-${PYVS}.pkg -PACKATE_DOWNLOAD_NAME=$ALPHADIR/$PACKAGE_FULL_NAME PACKAGE_FILE_NAME=$NRN_BLD/src/mac/build/NEURON.pkg + +cp $PACKAGE_FILE_NAME $HOME/$PACKAGE_FULL_NAME + echo " - Until we figure out how to automatically staple the notarization - the following two commands must be executed manually. - xcrun stapler staple $PACKAGE_FILE_NAME - cp $PACKAGE_FILE_NAME $HOME/$PACKAGE_FULL_NAME + Manually upload $HOME/$PACKAGE_FULL_NAME to github " diff --git a/docs/install/mac_pkg.md b/docs/install/mac_pkg.md index 17c30ee70d..8f69c73c64 100644 --- a/docs/install/mac_pkg.md +++ b/docs/install/mac_pkg.md @@ -17,7 +17,7 @@ On an Apple M1 or x86_64, the script, by default, creates, e.g., ```nrn-8.0a-726-gb9a811a32-macosx-11-universal2-py-38-39-310.pkg``` where the information between nrn and macosx comes from ```git describe```, -the number after macos refers to the ```MACOSX_DEPLOYMENT_TARGET=11``` +the number after macosx refers to the ```MACOSX_DEPLOYMENT_TARGET=11``` the next item(s) before py indicate the architectures on which the program can run (i.e. ```arm64```, ```x86_64```, or ```universal2``` for both) @@ -32,7 +32,7 @@ the same MACOSX_DEPLOYMENT_TARGET. You can check both of these with A space separated list of python executable arguments can be used in place of the internal default lists. ```$NRN_SRC``` is the location of the -repository, default ```$HOME/neuron/nrn```. The script makes sure ```$NRN_SRC/build``` +repository, default is the current working directory (hopefully you launch at the location of this script, e.g. ```$HOME/neuron/nrn```). The script makes sure ```$NRN_SRC/build``` exists and uses that folder to configure and build the software. The software is installed in ```/Applications/NEURON```. @@ -50,9 +50,12 @@ cmake .. -DCMAKE_INSTALL_PREFIX=$NRN_INSTALL \ -DCMAKE_PREFIX_PATH=/usr/X11 \ -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ ``` +[^1] +[^1]: NRN_RX3D_OPT_LEVEL=2 can build VERY slowly (cython translated cpp file can take a half hour or more). So for testing, we generally copy the script to temp.sh and modify to NRN_RX3D_OPT_LEVEL=0 + The default variables above will be ``` -pythons="python3.8;python3.9;python3.10" +pythons="python3.8;python3.9;python3.10;python3.11" archs_cmake='-DCMAKE_OSX_ARCHITECTURES=arm64;x86_64' ``` @@ -62,7 +65,11 @@ The MPI library does not have to be universal if configuring with ```-DNRN_ENABLE_MPI_DYNAMIC=ON``` as the library is not linked against during the build. -```make -j install``` is used to build and install in ```/Applications/NEURON``` +```make -j install``` [^2] is used to build and install in ```/Applications/NEURON``` + +[^2]: Instead of make, the script uses Ninja. ```ninja install``` and the cmake line begins with ```cmake .. -G Ninja ...``` + +Runs some basic tests to verify a successful build ```make macpkg``` ( see ```src/mac/CMakeLists.txt``` ) is used to: @@ -90,42 +97,47 @@ during the build. - request Apple to notarize NEURON.pkg ```src/macnrn_notarize.sh``` + If notarization succeeds it will finish with a few lines of the form + ``` + Current status: Accepted........Processing complete + id: bda82fa9-e0ab-4f86-aad8-3fee1d8c2369 + status: Accepted + ``` + and staple the package. + If notarizaton fails, it is occasionally due to Apple changing the contracts and demanding that "You must first sign the relevant contracts online. (1048)". In that case, go to [appstoreconnect.apple.com](https://appstoreconnect.apple.com) - to accept the legal docs. For other notarization failures, one must consult - the LogFileURL which can be obtained with + to accept the legal docs an try again. For other notarization failures, one must consult the LogFileURL which can be obtained with [^3] + [^3]: altool has been replaced by notarytool for purposes of notarization. See + https://developer.apple.com/documentation/technotes/tn3147-migrating-to-the-latest-notarization-tool + and ```nrn/src/mac/nrn_notarize.sh``` + + ``` + % xcrun notarytool log \ + --apple-id "$apple_id" \ + --team-id "$team_id" \ + --password "$app_specific_password" \ + "$id" ``` - % xcrun altool --notarization-info $RequestIdentifier \ - --username "michael.hines@yale.edu" \ - --password "`cat ~/.ssh/notarization-password`" - No errors getting notarization info. - - Date: 2022-01-02 23:38:12 +0000 - Hash: 7254157952d4f3573c2804879cf6da8d... - LogFileURL: https://osxapps-ssl.itunes.apple.com/itunes-assets... - RequestUUID: 152f0f0e-af58-4d22-b291-6a441825dd20 - Status: invalid - Status Code: 2 - Status Message: Package Invalid + where $id was printed by the notarization request. The apple_id, team_id, and + app_specfic_password are set by the script to the contents of the files ``` - where RequestIdentifer (the RequestUUID) appears in the email sent - back in response to the notarization request. + $HOME/.ssh/apple-id + $HOME/.ssh/apple-team-id + $HOME/.ssh/apple-notarization-password + ``` + The script ends by printing: ``` - Until we figure out how to automatically staple the notarization - the following two commands must be executed manually. - xcrun stapler staple $PACKAGE_FILE_NAME cp $PACKAGE_FILE_NAME $HOME/$PACKAGE_FULL_NAME + Manually upload $HOME/nrn-9.0a-88-gc6d8b9af6-macosx-10.15-universal2-py-39-310-311.pkg to github ``` -where the variables are filled out and can be copy/pasted to your -terminal after Apple sends an email declaring that notarization was successful. -The email from Apple usually takes just a few minutes but can be hours. I've been uploading the ```$PACKAGE_FULL_NAME``` as an artifact for a -Release version from https://github.com/neuronsimulatior/nrn by choosing +Release version from https://github.com/neuronsimulator/nrn by choosing Releases, choosing the proper Release tag (e.g. Release 8.0a), Edit release, and clicking in the "Attach binaries ..." near the bottom of the page. diff --git a/src/mac/nrn_notarize.sh b/src/mac/nrn_notarize.sh index 7e58fc9649..d99c8ca1d6 100755 --- a/src/mac/nrn_notarize.sh +++ b/src/mac/nrn_notarize.sh @@ -1,64 +1,72 @@ #!/usr/bin/env bash set -e -# App specific password used to request notarization -password_path="$HOME/.ssh/notarization-password" + +# See https://developer.apple.com/documentation/technotes/tn3147-migrating-to-the-latest-notarization-tool + +# App specific password and credentials used to request notarization +password_path="$HOME/.ssh/apple-notarization-password" +apple_id_path="$HOME/.ssh/apple-id" +team_id_path="$HOME/.ssh/apple-team-id" + if test -f "$password_path" ; then app_specific_password=`cat "$password_path"` else echo "\"$password_path\" does not exist" exit 1 fi +if test -f "$apple_id_path" ; then + apple_id=`cat "$apple_id_path"` +else + echo "\"$apple_id_path\" does not exist" + exit 1 +fi +if test -f "$team_id_path" ; then + team_id=`cat "$team_id_path"` +else + echo "\"$team_id_path\" does not exist" + exit 1 +fi pkg="$1" pkgname="$2" echo "Notarize request" -xcrun altool --notarize-app \ - --primary-bundle-id "edu.yale.neuron.pkg.$pkgname" \ - --username "michael.hines@yale.edu" \ +xcrun notarytool submit \ + --wait \ + --apple-id "$apple_id" \ + --team-id "$team_id" \ --password "$app_specific_password" \ - --file "$pkg" + "$pkg" -# you should get -# 2021-01-31 17:43:58.537 altool[3021:3231297] CFURLRequestSetHTTPCookieStorageAcceptPolicy_block_invoke: no longer implemented and should not be called -# No errors uploading '/Users/michaelhines/neuron/notarize/build/src/mac/build/NEURON.pkg'. -# RequestUUID = dc9dd4dd-a942-475a-ab59-4996f938d606 - -# and eventually get an email in just a few minutes saying that -# Your Mac software has been notarized. +# this should end with something like +#Processing complete +# id: 87c2e5e9-f37f-4a95-9941-8c2205fc90bd +# status: Accepted # At this point the software can be installed. However it is recommended -# xcrun stapler staple "$pkg" +# to staple the pkg file # so that Gatekeeper will be able to find the whitelist in the file itself # without the need to perform an online check. -# However it is unclear how long the wait will be if you do a wait loop over -# xcrun altool --notarization-info $RequestUUID \ -# --username "michael.hines@yale.edu" \ -# --password "$app_specific_password" -# until -# Status: in progress -# changes to -# Status: success -# at which point it is ok to run the stapler. - -# For now, echo a suggestion about waiting for an email and what stapler -# command to execute manually. +xcrun stapler staple "$pkg" -echo " -After getting an email from Apple stating that: - Your Mac software has been notarized. -manually execute +# Lastly, copy the pkg file to its proper name in $HOME - xcrun stapler staple \"$pkg\" - -so that Gatekeeper will be able to find the whitelist in the file itself -without the need to perform an online check. -" +# I've read that it is a good idea to check for warnings by +# fetching the notary log (using the id generated by the above submit) +# And if there are errors, then you certainly need to fetch the log, +# address the issues it shows and try again. +if false ; then + xcrun notarytool log \ + --apple-id "$apple_id" \ + --team-id "$team_id" \ + --password "$app_specific_password" \ + "$id" +fi -# If the notarization fails. E.g. the message from Apple is -# "Your Mac software was not notarized" -# then review the notarization LogFileURL by obtaining the log url with -# xcrun altool --notarization-info $RequestUUID \ -# --username "michael.hines@yale.edu" \ -# --password "$app_specific_password" -# and address the issues it shows and try again. +# To get a history of notarizations... +if false ; then + xcrun notarytool history \ + --apple-id "$apple_id" \ + --team-id "$team_id" \ + --password "$app_specific_password" +fi From 6fe551a57e871e50ff7ae8deb0106622e2590418 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 11:48:06 +0100 Subject: [PATCH 18/37] Fully remove Meschach (#2619) Mostly the code/cmake/docs that were still referring to meschach --- cmake/CompilerHelper.cmake | 5 +---- docs/hoc/programming/math/matrix.rst | 13 +++++-------- docs/python/programming/math/matrix.rst | 13 +++++-------- src/ivoc/matrix.cpp | 4 ---- src/nrniv/CMakeLists.txt | 7 ------- src/nrnpython/CMakeLists.txt | 1 - 6 files changed, 11 insertions(+), 32 deletions(-) diff --git a/cmake/CompilerHelper.cmake b/cmake/CompilerHelper.cmake index 74bc939c36..54adf1ce31 100644 --- a/cmake/CompilerHelper.cmake +++ b/cmake/CompilerHelper.cmake @@ -37,16 +37,13 @@ if(CMAKE_C_COMPILER_ID MATCHES "PGI" OR CMAKE_C_COMPILER_ID MATCHES "NVHPC") # "src/modlunit/units.cpp", warning #170-D: pointer points outside of underlying object # "src/nrnpython/grids.cpp", warning #174-D: expression has no effect # "src/nmodl/nocpout.cpp", warning #177-D: variable "j" was declared but never referenced - # "src/mesch/conjgrad.c", warning #180-D: argument is incompatible with formal parameter # "src/nrniv/partrans.cpp", warning #186-D: pointless comparison of unsigned integer with zero - # "src/mesch/machine.h", warning #301-D: typedef name has already been declared (with same type) # "src/nrnpython/rxdmath.cpp", warning #541-D: allowing all exceptions is incompatible with previous function # "src/nmodl/nocpout.cpp", warning #550-D: variable "sion" was set but never used # "src/gnu/neuron_gnu_builtin.h", warning #816-D: type qualifier on return type is meaningless" # "src/modlunit/consist.cpp", warning #2465-D: conversion from a string literal to "char *" is deprecated # ~~~ - list(APPEND NRN_COMPILE_FLAGS - --diag_suppress=1,47,111,128,170,174,177,180,186,301,541,550,816,2465) + list(APPEND NRN_COMPILE_FLAGS --diag_suppress=1,47,111,128,170,174,177,186,541,550,816,2465) endif() list(APPEND NRN_COMPILE_FLAGS -noswitcherror) list(APPEND NRN_LINK_FLAGS -noswitcherror) diff --git a/docs/hoc/programming/math/matrix.rst b/docs/hoc/programming/math/matrix.rst index 15a31632b2..617ad13fa4 100644 --- a/docs/hoc/programming/math/matrix.rst +++ b/docs/hoc/programming/math/matrix.rst @@ -53,20 +53,17 @@ Matrix By default, a new Matrix is of type MFULL (= 1) and allocates storage for all nrow*ncol elements. Scaffolding is in place for matrices of storage type MSPARSE (=2) and MBAND (=3) but not many methods have been interfaced - to the meschach library at this time. If a method is called on a matrix type + to the eigen library at this time. If a method is called on a matrix type whose method has not been implemented, an error message will be printed. It is intended that implemented methods will be transparent to the user, eg m*x=b (``x = m.solv(b)`` ) will solve the linear system regardless of the type of m and v1 = m*v2 (``v1 = m.mulv(v2)`` ) will perform the vector multiplication. - Matrix is implemented using the - `meschach c library by David E. Stewart `_ - (discovered at http://www.netlib.org/c/index.html\ ) which contains a large collection - of routines for sparse, banded, and full matrices. Many of the useful - routines have not - been interfaced with the hoc interpreter but can be easily added on request - or you can add it yourself + Matrix is implemented using the `eigen3 library `_ + which contains a large collection of routines for sparse, banded, and full matrices. + Many of the useful routines have not been interfaced with the hoc + interpreter but can be easily added on request or you can add it yourself by analogy with the code in ``nrn/src/ivoc/(matrix.c ocmatrix.[ch])`` At this time the MFULL matrix type is complete enough to do useful work and MSPARSE can be used to multiply a matrix by a vector and solve diff --git a/docs/python/programming/math/matrix.rst b/docs/python/programming/math/matrix.rst index 62e7f1b044..f494cb50b6 100755 --- a/docs/python/programming/math/matrix.rst +++ b/docs/python/programming/math/matrix.rst @@ -53,20 +53,17 @@ Matrix By default, a new Matrix is of type MFULL (= 1) and allocates storage for all nrow*ncol elements. Scaffolding is in place for matrices of storage type MSPARSE (=2) and MBAND (=3) but not many methods have been interfaced - to the meschach library at this time. If a method is called on a matrix type + to the eigen library at this time. If a method is called on a matrix type whose method has not been implemented, an error message will be printed. It is intended that implemented methods will be transparent to the user, eg m*x=b (``x = m.solv(b)`` ) will solve the linear system regardless of the type of m and v1 = m*v2 (``v1 = m.mulv(v2)`` ) will perform the vector multiplication. - Matrix is implemented using the - `meschach c library by David E. Stewart `_ - (discovered at http://www.netlib.org/c/index.html\ ) which contains a large collection - of routines for sparse, banded, and full matrices. Many of the useful - routines have not - been interfaced with the hoc interpreter but can be easily added on request - or you can add it yourself + Matrix is implemented using the `eigen3 library `_ + which contains a large collection of routines for sparse, banded, and full matrices. + Many of the useful routines have not been interfaced with the hoc + interpreter but can be easily added on request or you can add it yourself by analogy with the code in ``nrn/src/ivoc/(matrix.c ocmatrix.[ch])`` At this time the MFULL matrix type is complete enough to do useful work and MSPARSE can be used to multiply a matrix by a vector and solve diff --git a/src/ivoc/matrix.cpp b/src/ivoc/matrix.cpp index 67e2367578..3c66406403 100644 --- a/src/ivoc/matrix.cpp +++ b/src/ivoc/matrix.cpp @@ -358,10 +358,6 @@ static Object** m_muls(void* v) { if (ifarg(2)) { out = matrix_arg(2); } - // I believe meschach does this for us - // if (out != m) { - // out->resize(... - // } m->muls(*getarg(1), out); return temp_objvar(out); } diff --git a/src/nrniv/CMakeLists.txt b/src/nrniv/CMakeLists.txt index 2af3fc458d..babfe56cc4 100644 --- a/src/nrniv/CMakeLists.txt +++ b/src/nrniv/CMakeLists.txt @@ -214,7 +214,6 @@ set(NRN_INCLUDE_DIRS ${PROJECT_SOURCE_DIR}/external/Random123/include ${PROJECT_SOURCE_DIR}/src ${PROJECT_SOURCE_DIR}/src/gnu - ${PROJECT_SOURCE_DIR}/src/mesch ${PROJECT_SOURCE_DIR}/src/nrncvode ${PROJECT_SOURCE_DIR}/src/nrnmpi ${PROJECT_SOURCE_DIR}/src/nrnpython @@ -327,12 +326,6 @@ if(NRN_USE_BACKWARD) endif() if(NRN_HAVE_NVHPC_COMPILER) - # NVHPC/21.7 cannot compile znorm.c with -O2 or above. See also: - # https://forums.developer.nvidia.com/t/nvc-21-7-regression-internal-compiler-error-can-only-coerce-indirect-args/184847 - if(${CMAKE_CXX_COMPILER_VERSION} VERSION_EQUAL 21.7) - set_source_files_properties(${PROJECT_SOURCE_DIR}/src/mesch/znorm.c PROPERTIES COMPILE_OPTIONS - -Mnovect) - endif() # For NVHPC we will rely on FE exceptions as opposed to errno in order to make use of faster # builtins. One caveat is that if we use an optimization level greater than -O1, the FE exception # is not raised. See https://github.com/neuronsimulator/nrn/pull/1930 diff --git a/src/nrnpython/CMakeLists.txt b/src/nrnpython/CMakeLists.txt index 6f114673f6..79a24c67fb 100644 --- a/src/nrnpython/CMakeLists.txt +++ b/src/nrnpython/CMakeLists.txt @@ -25,7 +25,6 @@ if(NRN_ENABLE_PYTHON_DYNAMIC) ../ivoc ../nrniv ../gnu - ../mesch ../nrnmpi ${PROJECT_BINARY_DIR}/src/nrnpython ${PROJECT_BINARY_DIR}/src/ivos From f71c64e2a0e56aa628a6196bc9110937ad91c47e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 14:38:59 +0100 Subject: [PATCH 19/37] Simplify solving of sparse matrix (#2639) If we give a matrix to the SparseLU ctor it already analyze and factorize the matrix. Don't do it twice! --- src/ivoc/ocmatrix.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index 2823c58b50..2fd427fe11 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -282,8 +282,6 @@ void OcSparseMatrix::solv(Vect* in, Vect* out, bool use_lu) { if (!lu_ || !use_lu || lu_->rows() != m_.rows()) { m_.makeCompressed(); lu_ = std::make_unique>(m_); - lu_->analyzePattern(m_); - lu_->factorize(m_); } auto v1 = Vect2VEC(in); auto v2 = Vect2VEC(out); From 46be3b77f19fc0819aa7003bef52dc1ace5ab86a Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 12 Dec 2023 17:55:45 +0100 Subject: [PATCH 20/37] Support Python 3.12 wheels (#2630) * Update azure macOS and python for wheels * Add a numpy version * python 3.12 needs setuptools --- azure-pipelines.yml | 6 ++++++ nrn_requirements.txt | 3 +++ packaging/python/build_wheels.bash | 1 + packaging/python/test_wheels.sh | 4 +++- 4 files changed, 13 insertions(+), 1 deletion(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index d65673ed84..4a237d6b43 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -42,6 +42,8 @@ stages: python.version: '3.10' Python311: python.version: '3.11' + Python312: + python.version: '3.12' steps: # Secure files documentation: @@ -102,6 +104,10 @@ stages: python.version: '3.11' python.org.version: '3.11.7' python.installer.name: 'macos11.pkg' + Python312: + python.version: '3.12' + python.org.version: '3.12.0' + python.installer.name: 'macos11.pkg' steps: diff --git a/nrn_requirements.txt b/nrn_requirements.txt index 77ccb32ddc..06c1f1920a 100644 --- a/nrn_requirements.txt +++ b/nrn_requirements.txt @@ -13,3 +13,6 @@ pytest-cov mpi4py numpy find_libpython + +# For python 3.12 +setuptools diff --git a/packaging/python/build_wheels.bash b/packaging/python/build_wheels.bash index ed42822fbf..a2c9bec5f2 100755 --- a/packaging/python/build_wheels.bash +++ b/packaging/python/build_wheels.bash @@ -62,6 +62,7 @@ pip_numpy_install() { 39) numpy_ver="numpy==1.19.3" ;; 310) numpy_ver="numpy==1.21.3" ;; 311) numpy_ver="numpy==1.23.5" ;; + 312) numpy_ver="numpy==1.26.0" ;; *) echo "Error: numpy version not specified for this python!" && exit 1;; esac diff --git a/packaging/python/test_wheels.sh b/packaging/python/test_wheels.sh index dd015076d7..3e458d1f2a 100755 --- a/packaging/python/test_wheels.sh +++ b/packaging/python/test_wheels.sh @@ -256,7 +256,9 @@ $python_exe -m pip install --upgrade pip # install numpy, pytest and neuron -$python_exe -m pip install numpy pytest +# we install setuptools because since python 3.12 it is no more installed +# by default +$python_exe -m pip install numpy pytest setuptools $python_exe -m pip install $python_wheel $python_exe -m pip show neuron || $python_exe -m pip show neuron-nightly From c483af882ab62e93c416ceb37d52811033f3b8fe Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 14 Dec 2023 13:29:18 +0100 Subject: [PATCH 21/37] Bump Eigen to stable version 3.4.0 (#2642) The actual version was not an official one and got one bug. --- external/eigen | 2 +- src/ivoc/ocmatrix.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/external/eigen b/external/eigen index 328b5f9085..3147391d94 160000 --- a/external/eigen +++ b/external/eigen @@ -1 +1 @@ -Subproject commit 328b5f90858f93344ebc1484df8cadfd2a5da6dd +Subproject commit 3147391d946bb4b6c68edd901f2add6ac1f31f8c diff --git a/src/ivoc/ocmatrix.cpp b/src/ivoc/ocmatrix.cpp index 2fd427fe11..7afa87498a 100644 --- a/src/ivoc/ocmatrix.cpp +++ b/src/ivoc/ocmatrix.cpp @@ -128,7 +128,7 @@ void OcFullMatrix::symmeigen(Matrix* mout, Vect* vout) { void OcFullMatrix::svd1(Matrix* u, Matrix* v, Vect* d) { auto v1 = Vect2VEC(d); - Eigen::JacobiSVD svd(m_); + Eigen::JacobiSVD svd(m_, Eigen::ComputeFullU | Eigen::ComputeFullV); v1 = svd.singularValues(); if (u) { u->full()->m_ = svd.matrixU().transpose(); From 2cf724869ddac97efe9db76f1980fbd2ce86a0f9 Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Thu, 21 Dec 2023 16:27:26 +0100 Subject: [PATCH 22/37] Bring SIMD relared flags from Spack recipe to CMake (#2636) See flags added in BlueBrain/spack/pull/2023 and discussion in BlueBrain/spack/pull/2136/files#r1340030026 * I am bringing these flags at the NEURON level as we will need them as we want to soon enable SIMD/new-code generation for NEURON. * Introduce NRN_ENABLE_MATH_OPT cmake option to hide/disable extra math optimisations. * Rename NRN_WHEEL_BUILD CMake option to NRN_BINARY_DIST_BUILD In order to use the flag for mac or windows installer, rename it. * Add FastDebug as CMAKE_BUILD_TYPE Co-authored-by: Ioannis Magkanaris --- CMakeLists.txt | 10 +++--- ci/win_build_cmake.sh | 1 + cmake/BuildOptionDefaults.cmake | 3 +- cmake/CompilerFlagsHelpers.cmake | 22 ++++++++++++- cmake/ReleaseDebugAutoFlags.cmake | 50 +++++++++++++++++++++++++++--- docs/cmake_doc/options.rst | 16 ++++++++-- packaging/python/build_wheels.bash | 4 +-- src/coreneuron/CMakeLists.txt | 2 +- 8 files changed, 93 insertions(+), 15 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2630543e6c..44983a771b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -74,6 +74,8 @@ option(NRN_ENABLE_RX3D "Enable rx3d support" ${NRN_ENABLE_RX3D_DEFAULT}) option(NRN_ENABLE_CORENEURON "Enable CoreNEURON support" ${NRN_ENABLE_CORENEURON_DEFAULT}) option(NRN_ENABLE_BACKTRACE "Enable pretty-printed backtraces" ${NRN_ENABLE_BACKTRACE_DEFAULT}) option(NRN_ENABLE_TESTS "Enable unit tests" ${NRN_ENABLE_TESTS_DEFAULT}) +option(NRN_ENABLE_MATH_OPT "Enable extra math optimisations (to enable SIMD)" + ${NRN_ENABLE_MATH_OPT_DEFAULT}) set(NRN_ENABLE_MODEL_TESTS "${NRN_ENABLE_MODEL_TESTS_DEFAULT}" CACHE STRING "Comma-separated list of detailed models to enable tests of.") @@ -94,11 +96,11 @@ option(NRN_ENABLE_MOD_COMPATIBILITY "Enable CoreNEURON compatibility for MOD fil ${NRN_ENABLE_MOD_COMPATIBILITY_DEFAULT}) option(NRN_ENABLE_REL_RPATH "Use relative RPATH in binaries. for relocatable installs/Python" ${NRN_ENABLE_REL_RPATH_DEFAULT}) -option(NRN_WHEEL_BUILD ${NRN_WHEEL_BUILD_DEFAULT}) +option(NRN_BINARY_DIST_BUILD ${NRN_BINARY_DIST_BUILD_DEFAULT}) option(NRN_WHEEL_STATIC_READLINE "Use static readline libraries for the wheels." ${NRN_WHEEL_STATIC_READLINE_DEFAULT}) mark_as_advanced(NRN_ENABLE_REL_RPATH) -mark_as_advanced(NRN_WHEEL_BUILD) +mark_as_advanced(NRN_BINARY_DIST_BUILD) # ============================================================================= # Build options (string) @@ -214,7 +216,7 @@ include(cmake/modules/FindPythonModule.cmake) include(cmake/Coverage.cmake) # set CMAKE_BUILD_TYPE and associated flags using allowableBuildTypes and CMAKE_BUILD_TYPE_DEFAULT -set(allowableBuildTypes Custom Debug Release RelWithDebInfo Fast) +set(allowableBuildTypes Custom Debug Release RelWithDebInfo Fast FastDebug) include(ReleaseDebugAutoFlags) # Try and emit an intelligent warning if the version number currently set in the CMake project(...) @@ -929,7 +931,7 @@ if(BUILD_TYPE_UPPER MATCHES "CUSTOM") else() set(COMPILER_FLAGS "${CMAKE_CXX_FLAGS_${BUILD_TYPE_UPPER}}") endif() -string(JOIN " " COMPILER_FLAGS "${COMPILER_FLAGS}" ${NRN_COMPILE_FLAGS}) +string(JOIN " " COMPILER_FLAGS "${COMPILER_FLAGS}" ${NRN_COMPILE_FLAGS} ${CMAKE_CXX_FLAGS}) message(STATUS "") message(STATUS "Configured NEURON ${PROJECT_VERSION}") diff --git a/ci/win_build_cmake.sh b/ci/win_build_cmake.sh index 942e0fd71a..3bf55bc842 100755 --- a/ci/win_build_cmake.sh +++ b/ci/win_build_cmake.sh @@ -29,6 +29,7 @@ cd $BUILD_SOURCESDIRECTORY/build -DNRN_ENABLE_PYTHON=ON \ -DNRN_ENABLE_RX3D=ON \ -DNRN_RX3D_OPT_LEVEL=2 \ + -DNRN_BINARY_DIST_BUILD=ON \ -DPYTHON_EXECUTABLE=/c/Python38/python.exe \ -DNRN_ENABLE_PYTHON_DYNAMIC=ON \ -DNRN_PYTHON_DYNAMIC='c:/Python38/python.exe;c:/Python39/python.exe;c:/Python310/python.exe;c:/Python311/python.exe' \ diff --git a/cmake/BuildOptionDefaults.cmake b/cmake/BuildOptionDefaults.cmake index d52b65acd8..e3e1f3dd63 100644 --- a/cmake/BuildOptionDefaults.cmake +++ b/cmake/BuildOptionDefaults.cmake @@ -28,6 +28,7 @@ set(NRN_ENABLE_REL_RPATH_DEFAULT OFF) set(NRN_AVOID_ABSOLUTE_PATHS_DEFAULT OFF) set(NRN_NMODL_CXX_FLAGS_DEFAULT "-O0") set(NRN_SANITIZERS_DEFAULT "") +set(NRN_ENABLE_MATH_OPT_DEFAULT OFF) # Some distributions may set the prefix. To avoid errors, unset it set(NRN_PYTHON_DYNAMIC_DEFAULT "") @@ -43,7 +44,7 @@ set(PYTHON_EXECUTABLE_DEFAULT "") set(IV_LIB_DEFAULT "") # For wheel deployment -set(NRN_WHEEL_BUILD_DEFAULT OFF) +set(NRN_BINARY_DIST_BUILD_DEFAULT OFF) set(NRN_WHEEL_STATIC_READLINE_DEFAULT OFF) # we add some coreneuron options in order to check support like GPU diff --git a/cmake/CompilerFlagsHelpers.cmake b/cmake/CompilerFlagsHelpers.cmake index 942c614ff4..7272084057 100644 --- a/cmake/CompilerFlagsHelpers.cmake +++ b/cmake/CompilerFlagsHelpers.cmake @@ -14,7 +14,8 @@ set(SUPPORTED_COMPILER_LANGUAGE_LIST "C;CXX") foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) if(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "XL") set(CMAKE_${COMPILER_LANGUAGE}_COMPILER_IS_XLC ON) - elseif(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "Intel") + elseif(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "Intel" + OR CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "IntelLLVM") set(CMAKE_${COMPILER_LANGUAGE}_COMPILER_IS_ICC ON) elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") set(CMAKE_${COMPILER_LANGUAGE}_COMPILER_IS_MSVC) @@ -27,6 +28,10 @@ foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) endif() endforeach() +set(UNSAFE_MATH_FLAG + "-ffinite-math-only -fno-math-errno -funsafe-math-optimizations -fno-associative-math") +set(FASTDEBUG_FLAG "-g -O1") + # Set optimization flags for each compiler foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) @@ -39,6 +44,7 @@ foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) set(CMAKE_${COMPILER_LANGUAGE}_OPT_NONE "-O0") set(CMAKE_${COMPILER_LANGUAGE}_OPT_NORMAL "-O2") set(CMAKE_${COMPILER_LANGUAGE}_OPT_FAST "-O3") + set(CMAKE_${COMPILER_LANGUAGE}_OPT_FASTDEBUG ${FASTDEBUG_FLAG}) set(CMAKE_${COMPILER_LANGUAGE}_STACK_PROTECTION "-qstackprotect") set(CMAKE_${COMPILER_LANGUAGE}_POSITION_INDEPENDENT "-qpic=small") set(CMAKE_${COMPILER_LANGUAGE}_VECTORIZE "-qhot") @@ -59,9 +65,11 @@ foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) set(CMAKE_${COMPILER_LANGUAGE}_OPT_NONE "-O0") set(CMAKE_${COMPILER_LANGUAGE}_OPT_NORMAL "-O2") set(CMAKE_${COMPILER_LANGUAGE}_OPT_FAST "-O3") + set(CMAKE_${COMPILER_LANGUAGE}_OPT_FASTDEBUG ${FASTDEBUG_FLAG}) set(CMAKE_${COMPILER_LANGUAGE}_STACK_PROTECTION "-fstack-protector") set(CMAKE_${COMPILER_LANGUAGE}_POSITION_INDEPENDENT "-fPIC") set(CMAKE_${COMPILER_LANGUAGE}_VECTORIZE "-ftree-vectorize") + set(CMAKE_${COMPILER_LANGUAGE}_UNSAFE_MATH ${UNSAFE_MATH_FLAG}) set(IGNORE_UNKNOWN_PRAGMA_FLAGS "-Wno-unknown-pragmas") if(CMAKE_${COMPILER_LANGUAGE}_COMPILER_VERSION VERSION_GREATER "4.7.0") @@ -82,6 +90,11 @@ foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) set(CMAKE_${COMPILER_LANGUAGE}_OPT_NONE "-O0") set(CMAKE_${COMPILER_LANGUAGE}_OPT_NORMAL "-O2") set(CMAKE_${COMPILER_LANGUAGE}_OPT_FAST "-O3") + if(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "IntelLLVM") + set(CMAKE_${COMPILER_LANGUAGE}_OPT_FASTDEBUG "${FASTDEBUG_FLAG} -fp-model precise") + else() + set(CMAKE_${COMPILER_LANGUAGE}_OPT_FASTDEBUG "${FASTDEBUG_FLAG} -fp-model consistent") + endif() set(CMAKE_${COMPILER_LANGUAGE}_STACK_PROTECTION "-fstack-protector") set(CMAKE_${COMPILER_LANGUAGE}_POSITION_INDEPENDENT "-fpic") set(CMAKE_${COMPILER_LANGUAGE}_VECTORIZE "") @@ -94,11 +107,18 @@ foreach(COMPILER_LANGUAGE ${SUPPORTED_COMPILER_LANGUAGE_LIST}) set(CMAKE_${COMPILER_LANGUAGE}_OPT_NONE "-O0") set(CMAKE_${COMPILER_LANGUAGE}_OPT_NORMAL "-O2") set(CMAKE_${COMPILER_LANGUAGE}_OPT_FAST "-O3") + set(CMAKE_${COMPILER_LANGUAGE}_OPT_FASTDEBUG "-g -O1") set(CMAKE_${COMPILER_LANGUAGE}_STACK_PROTECTION "") set(CMAKE_${COMPILER_LANGUAGE}_POSITION_INDEPENDENT "-fPIC") set(CMAKE_${COMPILER_LANGUAGE}_VECTORIZE "") if(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "PGI") set(CMAKE_${COMPILER_LANGUAGE}_WARNING_ALL "") endif() + if(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "Clang") + set(CMAKE_${COMPILER_LANGUAGE}_UNSAFE_MATH ${UNSAFE_MATH_FLAG}) + endif() + if(CMAKE_${COMPILER_LANGUAGE}_COMPILER_ID STREQUAL "NVHPC") + set(CMAKE_${COMPILER_LANGUAGE}_OPT_FASTDEBUG "${FASTDEBUG_FLAG} -fno-omit-frame-pointer") + endif() endif() endforeach() diff --git a/cmake/ReleaseDebugAutoFlags.cmake b/cmake/ReleaseDebugAutoFlags.cmake index 61565e60f4..747b8801e6 100644 --- a/cmake/ReleaseDebugAutoFlags.cmake +++ b/cmake/ReleaseDebugAutoFlags.cmake @@ -21,6 +21,9 @@ endif() # Release : Release mode, no debuginfo # RelWithDebInfo : Distribution mode, basic optimizations for portable code with debuginfos # Fast : Maximum level of optimization. Target native architecture, not portable code +# FastDebug: Similar to Debug with a bit higher level optimisations (-O1) and other compiler +# flags so that it's faster than -O0 but still produces consistent results for +# testing and debugging purposes. # ~~~ include(CompilerFlagsHelpers) @@ -33,17 +36,56 @@ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_DEBUGINFO_FLAGS} ${CMAKE_CXX_OPT_NONE} ${CMAKE_CXX_STACK_PROTECTION} ${CMAKE_CXX_IGNORE_WARNINGS}" ) -set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_OPT_NORMAL} ${CMAKE_C_IGNORE_WARNINGS}") -set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_OPT_NORMAL} ${CMAKE_CXX_IGNORE_WARNINGS}") +set(C_UNSAFE_MATH_FLAGS "") +set(CXX_UNSAFE_MATH_FLAGS "") +if(NRN_ENABLE_MATH_OPT) + set(C_UNSAFE_MATH_FLAGS ${CMAKE_C_UNSAFE_MATH}) + set(CXX_UNSAFE_MATH_FLAGS ${CMAKE_CXX_UNSAFE_MATH}) +endif() + +set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_OPT_NORMAL} ${CMAKE_C_IGNORE_WARNINGS} ${C_UNSAFE_MATH_FLAGS}") +set(CMAKE_CXX_FLAGS_RELEASE + "${CMAKE_CXX_OPT_NORMAL} ${CMAKE_CXX_IGNORE_WARNINGS} ${CXX_UNSAFE_MATH_FLAGS}") set(CMAKE_C_FLAGS_RELWITHDEBINFO - "${CMAKE_C_DEBUGINFO_FLAGS} ${CMAKE_C_OPT_NORMAL} ${CMAKE_C_IGNORE_WARNINGS}") + "${CMAKE_C_DEBUGINFO_FLAGS} ${CMAKE_C_OPT_NORMAL} ${CMAKE_C_IGNORE_WARNINGS} ${C_UNSAFE_MATH_FLAGS}" +) set(CMAKE_CXX_FLAGS_RELWITHDEBINFO - "${CMAKE_CXX_DEBUGINFO_FLAGS} ${CMAKE_CXX_OPT_NORMAL} ${CMAKE_CXX_IGNORE_WARNINGS}") + "${CMAKE_CXX_DEBUGINFO_FLAGS} ${CMAKE_CXX_OPT_NORMAL} ${CMAKE_CXX_IGNORE_WARNINGS} ${CXX_UNSAFE_MATH_FLAGS}" +) set(CMAKE_C_FLAGS_FAST "${CMAKE_C_OPT_FAST} ${CMAKE_C_LINK_TIME_OPT} ${CMAKE_C_GEN_NATIVE} ${CMAKE_C_IGNORE_WARNINGS}") set(CMAKE_CXX_FLAGS_FAST "${CMAKE_CXX_OPT_FAST} ${CMAKE_CXX_LINK_TIME_OPT} ${CMAKE_CXX_GEN_NATIVE} ${CMAKE_CXX_IGNORE_WARNINGS}" ) + +set(CMAKE_C_FLAGS_FASTDEBUG "${CMAKE_C_OPT_FASTDEBUG} ${CMAKE_C_IGNORE_WARNINGS}") +set(CMAKE_CXX_FLAGS_FASTDEBUG "${CMAKE_CXX_OPT_FASTDEBUG} ${CMAKE_CXX_IGNORE_WARNINGS}") # ~~~ + +# for binary distributions, avoid addition of OpenMP specific flag as compiler on end-user machine +# may not support it. +if(NOT DEFINED NRN_BINARY_DIST_BUILD OR NOT NRN_BINARY_DIST_BUILD) + include(CheckCXXCompilerFlag) + # Check support for OpenMP SIMD constructs + set(SIMD_FLAGS "") + + if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") + set(SIMD_FLAGS "-qopenmp-simd") + elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") + set(SIMD_FLAGS "-openmp:experimental") + elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang") + set(SIMD_FLAGS "-openmp-simd") + else() # not ICC, MSVC, or Clang => GCC and others + set(SIMD_FLAGS "-fopenmp-simd") + endif() + + check_cxx_compiler_flag("${SIMD_FLAGS}" COMPILER_SUPPORT_SIMD) + if(COMPILER_SUPPORT_SIMD) + set(CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${SIMD_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${SIMD_FLAGS}") + else() + message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no support for OpenMP SIMD construct") + endif() +endif() diff --git a/docs/cmake_doc/options.rst b/docs/cmake_doc/options.rst index 7dca20f1e9..83497d5c0b 100644 --- a/docs/cmake_doc/options.rst +++ b/docs/cmake_doc/options.rst @@ -12,14 +12,14 @@ together using the following instructions: .. code-block:: shell git clone https://github.com/neuronsimulator/nrn # latest development branch - git clone https://github.com/neuronsimulator/nrn -b 8.0.0 # specific release version 8.0.0 + git clone https://github.com/neuronsimulator/nrn -b 8.2.3 # specific release version 8.2.3 cd nrn .. .. warning:: To build NEURON from source you either need to clone the NEURON Git repository or download a source code archive that includes - Git submodules, such as the ``full-src-package-X.Y.Z.tar.gz`` file in + Git submodules, such as the ``nrn-full-src-package-X.Y.Z.tar.gz`` file in the `NEURON releases `__ on GitHub. The tarballs like ``Source code (tar.gz)`` or @@ -648,3 +648,15 @@ NRN_PYTHON_EXTRA_FOR_TESTS:STRING= built with support for, for use in tests of error messages and reporting. For these purposes, minor versions (3.X and 3.Y) are considered different and patch versions (3.8.X and 3.8.Y) are considered to be the same. + +NRN_ENABLE_MATH_OPT:BOOL=OFF +------------------------------------- + Enable extra math optimisations. + + When using compilers like GCC and Clang, one needs to explicitly use compiler + flags like `-funsafe-math-optimizations` in order to generate SIMD/vectorised + code using vector math library. This flag adds these extra compiler flags + to enable SIMD code. + + Note: Compilers like Intel, NVHPC, Cray etc enable such optimisations + by default. diff --git a/packaging/python/build_wheels.bash b/packaging/python/build_wheels.bash index a2c9bec5f2..007fa581ce 100755 --- a/packaging/python/build_wheels.bash +++ b/packaging/python/build_wheels.bash @@ -91,7 +91,7 @@ build_wheel_linux() { CMAKE_DEFS="NRN_MPI_DYNAMIC=$3" if [ "$USE_STATIC_READLINE" == "1" ]; then - CMAKE_DEFS="$CMAKE_DEFS,NRN_WHEEL_BUILD=ON,NRN_WHEEL_STATIC_READLINE=ON" + CMAKE_DEFS="$CMAKE_DEFS,NRN_BINARY_DIST_BUILD=ON,NRN_WHEEL_STATIC_READLINE=ON" fi if [ "$2" == "coreneuron" ]; then @@ -148,7 +148,7 @@ build_wheel_osx() { CMAKE_DEFS="NRN_MPI_DYNAMIC=$3" if [ "$USE_STATIC_READLINE" == "1" ]; then - CMAKE_DEFS="$CMAKE_DEFS,NRN_WHEEL_BUILD=ON,NRN_WHEEL_STATIC_READLINE=ON" + CMAKE_DEFS="$CMAKE_DEFS,NRN_BINARY_DIST_BUILD=ON,NRN_WHEEL_STATIC_READLINE=ON" fi # We need to "fix" the platform tag if the Python installer is universal2 diff --git a/src/coreneuron/CMakeLists.txt b/src/coreneuron/CMakeLists.txt index 52b0de72d9..a18bb38d7d 100644 --- a/src/coreneuron/CMakeLists.txt +++ b/src/coreneuron/CMakeLists.txt @@ -318,7 +318,7 @@ endif() # Do not set this when building wheels. The nrnivmodl workflow means that we do not know what # compiler will be invoked with these flags, so we have to use flags that are as generic as # possible. -if(NOT DEFINED NRN_WHEEL_BUILD OR NOT NRN_WHEEL_BUILD) +if(NOT DEFINED NRN_BINARY_DIST_BUILD OR NOT NRN_BINARY_DIST_BUILD) list(APPEND CORENRN_EXTRA_CXX_FLAGS "${IGNORE_UNKNOWN_PRAGMA_FLAGS}") endif() From b045e504e044c6084759ada9d9981c85fcb5e1d7 Mon Sep 17 00:00:00 2001 From: nrnhines Date: Sun, 7 Jan 2024 04:03:07 -0500 Subject: [PATCH 23/37] cell builder session files involving extracellular could not be saved. (#2651) --- src/ivoc/datapath.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/ivoc/datapath.cpp b/src/ivoc/datapath.cpp index cee0f49d5b..c4da677e77 100644 --- a/src/ivoc/datapath.cpp +++ b/src/ivoc/datapath.cpp @@ -462,7 +462,12 @@ void HocDataPathImpl::search(Prop* prop, double x) { if (memb_func[type].hoc_mech) { pd = prop->ob->u.dataspace[ir].pval; } else { - pd = static_cast(prop->param_handle_legacy(ir)); + if (type == EXTRACELL && ir == neuron::extracellular::vext_pseudoindex()) { + // skip as it was handled by caller + continue; + } else { + pd = static_cast(prop->param_handle_legacy(ir)); + } } imax = hoc_total_array_data(psym, 0); for (i = 0; i < imax; ++i) { From 5a91ef533009b4de757364fcb2466ee0fe8c6b2c Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso <41900536+jorblancoa@users.noreply.github.com> Date: Sun, 7 Jan 2024 16:45:35 +0100 Subject: [PATCH 24/37] Soma summation fix for SONATA simulations (#2647) - Both Cell/Soma SectionType should work - BBP internal ticket: https://bbpteam.epfl.ch/project/issues/browse/BBPBGLIB-1107 --- src/coreneuron/io/reports/report_handler.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/coreneuron/io/reports/report_handler.cpp b/src/coreneuron/io/reports/report_handler.cpp index 8089b803ab..0fd48d7a0a 100644 --- a/src/coreneuron/io/reports/report_handler.cpp +++ b/src/coreneuron/io/reports/report_handler.cpp @@ -297,7 +297,8 @@ VarsToReport ReportHandler::get_summation_vars_to_report( if (report.section_type == SectionType::All) { double* variable = report_variable + segment_id; to_report.emplace_back(VarWithMapping(section_id, variable)); - } else if (report.section_type == SectionType::Cell) { + } else if (report.section_type == SectionType::Cell || + report.section_type == SectionType::Soma) { summation_report.gid_segments_[gid].push_back(segment_id); } } From d2c8c4b2ea0c55bc522917cd224a2f43bb038f77 Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Sun, 7 Jan 2024 16:47:37 +0100 Subject: [PATCH 25/37] Gitlab CI: allocate GPU jobs exclusively for CI user (#2648) With this change we will see if CI failures on BB5 phase 2 GPU nodes will be mitigated. --- .gitlab-ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 9ff3ce0594..d6baa64fb0 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -187,6 +187,7 @@ simulation_stack: bb5_constraint: volta bb5_cpus_per_task: 2 bb5_partition: prod # assume this is a good source of GPU nodes + bb5_exclusive: user # allocate gpu node exclusively for the CI user (to avoid errors from oversubscription) .test_neuron: extends: [.ctest] variables: From 1e0d8be3b56c142973c8f82e42340d093fdf464e Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 8 Jan 2024 12:59:49 +0100 Subject: [PATCH 26/37] Make have_to_want with templates (#2645) * Having a hpp file instead of cpp file, * Use template so that the algorithm can be included twice, * Remove usage of 3 configurations variables by preprocessor (before the algorithm was paramretrized by preprocessor variable leading to some weird inclusion), * Use std::vector<> instead of new T*, * Create a new structure Data to group data, count and displ. --- src/nrniv/have2want.cpp | 259 ---------------------------------------- src/nrniv/have2want.hpp | 202 +++++++++++++++++++++++++++++++ src/nrniv/partrans.cpp | 133 +++++++++++---------- 3 files changed, 272 insertions(+), 322 deletions(-) delete mode 100644 src/nrniv/have2want.cpp create mode 100644 src/nrniv/have2want.hpp diff --git a/src/nrniv/have2want.cpp b/src/nrniv/have2want.cpp deleted file mode 100644 index 5ac55709aa..0000000000 --- a/src/nrniv/have2want.cpp +++ /dev/null @@ -1,259 +0,0 @@ -/* -To be included by a file that desires rendezvous rank exchange functionality. -Need to define HAVEWANT_t, HAVEWANT_alltoallv, and HAVEWANT2Int -The latter is a map or unordered_map. -E.g. std::unordered_map -*/ - -#ifdef have2want_cpp -#error "This implementation can only be included once" -// The static function names used to involve a macro name (NrnHash) but now, -// with the use of std::..., it may be the case this could be included -// multiple times or even transformed into a template. -#endif - -#define have2want_cpp - -/* - -A rank owns a set of HAVEWANT_t keys and wants information associated with -a set of HAVEWANT_t keys owned by unknown ranks. Owners do not know which -ranks want their information. Ranks that want info do not know which ranks -own that info. - -The have_to_want function returns two new vectors of keys along with -associated count and displacement vectors of length nhost and nhost+1 -respectively. Note that a send_to_want_displ[i+1] = - send_to_want_cnt[i] + send_to_want_displ[i] . - -send_to_want[send_to_want_displ[i] to send_to_want_displ[i+1]] contains -the keys from this rank for which rank i wants information. - -recv_from_have[recv_from_have_displ[i] to recv_from_have_displ[i+1] contains -the keys from which rank i is sending information to this rank. - -Note that on rank i, the order of keys in the rank j area of send_to_want -is the same order of keys on rank j in the ith area in recv_from_have. - -The rendezvous_rank function is used to parallelize this computation -and minimize memory usage so that no single rank ever needs to know all keys. -*/ - -#ifndef HAVEWANT_t -#define HAVEWANT_t int -#endif - -// round robin default rendezvous rank function -static int default_rendezvous(HAVEWANT_t key) { - return key % nrnmpi_numprocs; -} - -static int* cnt2displ(int* cnt) { - int* displ = new int[nrnmpi_numprocs + 1]; - displ[0] = 0; - for (int i = 0; i < nrnmpi_numprocs; ++i) { - displ[i + 1] = displ[i] + cnt[i]; - } - return displ; -} - -static int* srccnt2destcnt(int* srccnt) { - int* destcnt = new int[nrnmpi_numprocs]; - nrnmpi_int_alltoall(srccnt, destcnt, 1); - return destcnt; -} - -static void rendezvous_rank_get(HAVEWANT_t* data, - int size, - HAVEWANT_t*& sdata, - int*& scnt, - int*& sdispl, - HAVEWANT_t*& rdata, - int*& rcnt, - int*& rdispl, - int (*rendezvous_rank)(HAVEWANT_t)) { - int nhost = nrnmpi_numprocs; - - // count what gets sent - scnt = new int[nhost]; - for (int i = 0; i < nhost; ++i) { - scnt[i] = 0; - } - for (int i = 0; i < size; ++i) { - int r = (*rendezvous_rank)(data[i]); - ++scnt[r]; - } - - sdispl = cnt2displ(scnt); - rcnt = srccnt2destcnt(scnt); - rdispl = cnt2displ(rcnt); - sdata = new HAVEWANT_t[sdispl[nhost] + 1]; // ensure not 0 size - rdata = new HAVEWANT_t[rdispl[nhost] + 1]; // ensure not 0 size - // scatter data into sdata by recalculating scnt. - for (int i = 0; i < nhost; ++i) { - scnt[i] = 0; - } - for (int i = 0; i < size; ++i) { - int r = (*rendezvous_rank)(data[i]); - sdata[sdispl[r] + scnt[r]] = data[i]; - ++scnt[r]; - } - HAVEWANT_alltoallv(sdata, scnt, sdispl, rdata, rcnt, rdispl); -} - -static void have_to_want(HAVEWANT_t* have, - int have_size, - HAVEWANT_t* want, - int want_size, - HAVEWANT_t*& send_to_want, - int*& send_to_want_cnt, - int*& send_to_want_displ, - HAVEWANT_t*& recv_from_have, - int*& recv_from_have_cnt, - int*& recv_from_have_displ, - int (*rendezvous_rank)(HAVEWANT_t)) { - // 1) Send have and want to the rendezvous ranks. - // 2) Rendezvous rank matches have and want. - // 3) Rendezvous ranks tell the want ranks which ranks own the keys - // 4) Ranks that want tell owner ranks where to send. - - int nhost = nrnmpi_numprocs; - - // 1) Send have and want to the rendezvous ranks. - HAVEWANT_t *have_s_data, *have_r_data; - int *have_s_cnt, *have_s_displ, *have_r_cnt, *have_r_displ; - rendezvous_rank_get(have, - have_size, - have_s_data, - have_s_cnt, - have_s_displ, - have_r_data, - have_r_cnt, - have_r_displ, - rendezvous_rank); - delete[] have_s_cnt; - delete[] have_s_displ; - delete[] have_s_data; - // assume it is an error if two ranks have the same key so create - // hash table of key2rank. Will also need it for matching have and want - HAVEWANT2Int havekey2rank = HAVEWANT2Int(have_r_displ[nhost] + 1); // ensure not empty. - for (int r = 0; r < nhost; ++r) { - for (int i = 0; i < have_r_cnt[r]; ++i) { - HAVEWANT_t key = have_r_data[have_r_displ[r] + i]; - if (havekey2rank.find(key) != havekey2rank.end()) { - hoc_execerr_ext( - "internal error in have_to_want: key %lld owned by multiple ranks\n", - (long long) key); - } - havekey2rank[key] = r; - } - } - delete[] have_r_data; - delete[] have_r_cnt; - delete[] have_r_displ; - - HAVEWANT_t *want_s_data, *want_r_data; - int *want_s_cnt, *want_s_displ, *want_r_cnt, *want_r_displ; - rendezvous_rank_get(want, - want_size, - want_s_data, - want_s_cnt, - want_s_displ, - want_r_data, - want_r_cnt, - want_r_displ, - rendezvous_rank); - - // 2) Rendezvous rank matches have and want. - // we already have made the havekey2rank map. - // Create an array parallel to want_r_data which contains the ranks that - // have that data. - int n = want_r_displ[nhost]; - int* want_r_ownerranks = new int[n]; - for (int r = 0; r < nhost; ++r) { - for (int i = 0; i < want_r_cnt[r]; ++i) { - int ix = want_r_displ[r] + i; - HAVEWANT_t key = want_r_data[ix]; - auto search = havekey2rank.find(key); - if (search == havekey2rank.end()) { - hoc_execerr_ext( - "internal error in have_to_want: key = %lld is wanted but does not exist\n", - (long long) key); - } - want_r_ownerranks[ix] = search->second; - } - } - delete[] want_r_data; - - // 3) Rendezvous ranks tell the want ranks which ranks own the keys - // The ranks that want keys need to know the ranks that own those keys. - // The want_s_ownerranks will be parallel to the want_s_data. - // That is, each item defines the rank from which information associated - // with that key is coming from - int* want_s_ownerranks = new int[want_s_displ[nhost]]; - if (nrn_sparse_partrans > 0) { - nrnmpi_int_alltoallv_sparse(want_r_ownerranks, - want_r_cnt, - want_r_displ, - want_s_ownerranks, - want_s_cnt, - want_s_displ); - } else { - nrnmpi_int_alltoallv(want_r_ownerranks, - want_r_cnt, - want_r_displ, - want_s_ownerranks, - want_s_cnt, - want_s_displ); - } - - delete[] want_r_ownerranks; - delete[] want_r_cnt; - delete[] want_r_displ; - - // 4) Ranks that want tell owner ranks where to send. - // Finished with the rendezvous ranks. The ranks that want keys know the - // owner ranks for those keys. The next step is for the want ranks to - // tell the owner ranks where to send. - // The parallel want_s_ownerranks and want_s_data are now uselessly ordered - // by rendezvous rank. Reorganize so that want ranks can tell owner ranks - // what they want. - n = want_s_displ[nhost]; - delete[] want_s_displ; - for (int i = 0; i < nhost; ++i) { - want_s_cnt[i] = 0; - } - HAVEWANT_t* old_want_s_data = want_s_data; - want_s_data = new HAVEWANT_t[n]; - // compute the counts - for (int i = 0; i < n; ++i) { - int r = want_s_ownerranks[i]; - ++want_s_cnt[r]; - } - want_s_displ = cnt2displ(want_s_cnt); - for (int i = 0; i < nhost; ++i) { - want_s_cnt[i] = 0; - } // recount while filling - for (int i = 0; i < n; ++i) { - int r = want_s_ownerranks[i]; - HAVEWANT_t key = old_want_s_data[i]; - want_s_data[want_s_displ[r] + want_s_cnt[r]] = key; - ++want_s_cnt[r]; - } - delete[] want_s_ownerranks; - delete[] old_want_s_data; - want_r_cnt = srccnt2destcnt(want_s_cnt); - want_r_displ = cnt2displ(want_r_cnt); - want_r_data = new HAVEWANT_t[want_r_displ[nhost]]; - HAVEWANT_alltoallv( - want_s_data, want_s_cnt, want_s_displ, want_r_data, want_r_cnt, want_r_displ); - // now the want_r_data on the have_ranks are grouped according to the ranks - // that want those keys. - - send_to_want = want_r_data; - send_to_want_cnt = want_r_cnt; - send_to_want_displ = want_r_displ; - recv_from_have = want_s_data; - recv_from_have_cnt = want_s_cnt; - recv_from_have_displ = want_s_displ; -} diff --git a/src/nrniv/have2want.hpp b/src/nrniv/have2want.hpp new file mode 100644 index 0000000000..11b22d28ed --- /dev/null +++ b/src/nrniv/have2want.hpp @@ -0,0 +1,202 @@ +#pragma once + +#include +#include +#include + +/* +A rank owns a set of T keys and wants information associated with +a set of T keys owned by unknown ranks. Owners do not know which +ranks want their information. Ranks that want info do not know which ranks +own that info. + +The have_to_want function returns two new vectors of keys along with +associated count and displacement vectors of length nhost and nhost+1 +respectively. Note that a send_to_want_displ[i+1] = + send_to_want_cnt[i] + send_to_want_displ[i] . + +send_to_want[send_to_want_displ[i] to send_to_want_displ[i+1]] contains +the keys from this rank for which rank i wants information. + +recv_from_have[recv_from_have_displ[i] to recv_from_have_displ[i+1] contains +the keys from which rank i is sending information to this rank. + +Note that on rank i, the order of keys in the rank j area of send_to_want +is the same order of keys on rank j in the ith area in recv_from_have. + +The rendezvous_rank function is used to parallelize this computation +and minimize memory usage so that no single rank ever needs to know all keys. +*/ + +// round robin rendezvous rank function +template +int rendezvous_rank(const T& key) { + return key % nrnmpi_numprocs; +} + +template +struct Data { + std::vector data{}; + std::vector cnt{}; + std::vector displ{}; +}; + +static std::vector cnt2displ(const std::vector& cnt) { + std::vector displ(nrnmpi_numprocs + 1); + std::partial_sum(cnt.cbegin(), cnt.cend(), displ.begin() + 1); + return displ; +} + +static std::vector srccnt2destcnt(std::vector srccnt) { + std::vector destcnt(nrnmpi_numprocs); + nrnmpi_int_alltoall(srccnt.data(), destcnt.data(), 1); + return destcnt; +} + +template +static std::tuple, Data> rendezvous_rank_get(const std::vector& data, + F alltoall_function) { + int nhost = nrnmpi_numprocs; + + Data s; + // count what gets sent + s.cnt = std::vector(nhost); + + for (const auto& e: data) { + int r = rendezvous_rank(e); + s.cnt[r] += 1; + } + + s.displ = cnt2displ(s.cnt); + s.data.resize(s.displ[nhost] + 1); + + Data r; + r.cnt = srccnt2destcnt(s.cnt); + r.displ = cnt2displ(r.cnt); + r.data.resize(r.displ[nhost]); + // scatter data into sdata by recalculating s.cnt. + std::fill(s.cnt.begin(), s.cnt.end(), 0); + for (const auto& e: data) { + int rank = rendezvous_rank(e); + s.data[s.displ[rank] + s.cnt[rank]] = e; + s.cnt[rank] += 1; + } + alltoall_function(s, r); + return {s, r}; +} + +template +std::pair, Data> have_to_want(const std::vector& have, + const std::vector& want, + F alltoall_function) { + // 1) Send have and want to the rendezvous ranks. + // 2) Rendezvous rank matches have and want. + // 3) Rendezvous ranks tell the want ranks which ranks own the keys + // 4) Ranks that want tell owner ranks where to send. + + int nhost = nrnmpi_numprocs; + + // 1) Send have and want to the rendezvous ranks. + + // hash table of key2rank. Will also need it for matching have and want + std::unordered_map havekey2rank{}; + { + auto [_, have_r] = rendezvous_rank_get(have, alltoall_function); + // assume it is an error if two ranks have the same key so create + havekey2rank.reserve(have_r.displ[nhost] + 1); + for (int r = 0; r < nhost; ++r) { + for (int i = 0; i < have_r.cnt[r]; ++i) { + T key = have_r.data[have_r.displ[r] + i]; + if (havekey2rank.find(key) != havekey2rank.end()) { + hoc_execerr_ext( + "internal error in have_to_want: key %lld owned by multiple ranks\n", + (long long) key); + } + havekey2rank[key] = r; + } + } + } + + auto [want_s, want_r] = rendezvous_rank_get(want, alltoall_function); + + // 2) Rendezvous rank matches have and want. + // we already have made the havekey2rank map. + // Create an array parallel to want_r_data which contains the ranks that + // have that data. + int n = want_r.displ[nhost]; + std::vector want_r_ownerranks(n); + for (int r = 0; r < nhost; ++r) { + for (int i = 0; i < want_r.cnt[r]; ++i) { + int ix = want_r.displ[r] + i; + T key = want_r.data[ix]; + auto search = havekey2rank.find(key); + if (search == havekey2rank.end()) { + hoc_execerr_ext( + "internal error in have_to_want: key = %lld is wanted but does not exist\n", + (long long) key); + } + want_r_ownerranks[ix] = search->second; + } + } + + // 3) Rendezvous ranks tell the want ranks which ranks own the keys + // The ranks that want keys need to know the ranks that own those keys. + // The want_s_ownerranks will be parallel to the want_s_data. + // That is, each item defines the rank from which information associated + // with that key is coming from + std::vector want_s_ownerranks(want_s.displ[nhost]); + if (nrn_sparse_partrans > 0) { + nrnmpi_int_alltoallv_sparse(want_r_ownerranks.data(), + want_r.cnt.data(), + want_r.displ.data(), + want_s_ownerranks.data(), + want_s.cnt.data(), + want_s.displ.data()); + } else { + nrnmpi_int_alltoallv(want_r_ownerranks.data(), + want_r.cnt.data(), + want_r.displ.data(), + want_s_ownerranks.data(), + want_s.cnt.data(), + want_s.displ.data()); + } + + // 4) Ranks that want tell owner ranks where to send. + // Finished with the rendezvous ranks. The ranks that want keys know the + // owner ranks for those keys. The next step is for the want ranks to + // tell the owner ranks where to send. + // The parallel want_s_ownerranks and want_s_data are now uselessly ordered + // by rendezvous rank. Reorganize so that want ranks can tell owner ranks + // what they want. + n = want_s.displ[nhost]; + for (int i = 0; i < nhost; ++i) { + want_s.cnt[i] = 0; + } + std::vector old_want_s_data(n); + std::swap(old_want_s_data, want_s.data); + // compute the counts + for (int i = 0; i < n; ++i) { + int r = want_s_ownerranks[i]; + ++want_s.cnt[r]; + } + want_s.displ = cnt2displ(want_s.cnt); + for (int i = 0; i < nhost; ++i) { + want_s.cnt[i] = 0; + } // recount while filling + for (int i = 0; i < n; ++i) { + int r = want_s_ownerranks[i]; + T key = old_want_s_data[i]; + want_s.data[want_s.displ[r] + want_s.cnt[r]] = key; + ++want_s.cnt[r]; + } + + Data new_want_r{}; + new_want_r.cnt = srccnt2destcnt(want_s.cnt); + new_want_r.displ = cnt2displ(new_want_r.cnt); + new_want_r.data.resize(new_want_r.displ[nhost]); + alltoall_function(want_s, new_want_r); + // now the want_r_data on the have_ranks are grouped according to the ranks + // that want those keys. + + return {new_want_r, want_s}; +} diff --git a/src/nrniv/partrans.cpp b/src/nrniv/partrans.cpp index 1f34f9bbf2..1ebfdf335d 100644 --- a/src/nrniv/partrans.cpp +++ b/src/nrniv/partrans.cpp @@ -17,23 +17,48 @@ #include #include +#if NRNMPI +#include "have2want.hpp" +#endif + + #if NRNLONGSGID #if NRNMPI -static void sgid_alltoallv(sgid_t* s, int* scnt, int* sdispl, sgid_t* r, int* rcnt, int* rdispl) { +static void sgid_alltoallv(Data& s, Data& r) { if (nrn_sparse_partrans > 0) { - nrnmpi_long_alltoallv_sparse(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_long_alltoallv_sparse(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } else { - nrnmpi_long_alltoallv(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_long_alltoallv(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } } #endif // NRNMPI #else // not NRNLONGSGID #if NRNMPI -static void sgid_alltoallv(sgid_t* s, int* scnt, int* sdispl, sgid_t* r, int* rcnt, int* rdispl) { +static void sgid_alltoallv(Data& s, Data& r) { if (nrn_sparse_partrans > 0) { - nrnmpi_int_alltoallv_sparse(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_int_alltoallv_sparse(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } else { - nrnmpi_int_alltoallv(s, scnt, sdispl, r, rcnt, rdispl); + nrnmpi_int_alltoallv(s.data.data(), + s.cnt.data(), + s.displ.data(), + r.data.data(), + r.cnt.data(), + r.displ.data()); } } #endif // NRNMPI @@ -152,8 +177,12 @@ static std::vector> poutsrc_; // prior t // to proper place in // outsrc_buf_ static int* poutsrc_indices_; // for recalc pointers -static int insrc_buf_size_, *insrccnt_, *insrcdspl_; -static int outsrc_buf_size_, *outsrccnt_, *outsrcdspl_; +static int insrc_buf_size_; +static std::vector insrccnt_; +static std::vector insrcdspl_; +static int outsrc_buf_size_; +static std::vector outsrccnt_; +static std::vector outsrcdspl_; static MapSgid2Int sid2insrc_; // received interprocessor sid data is // associated with which insrc_buf index. Created by nrnmpi_setup_transfer // and used by mk_ttd @@ -522,11 +551,19 @@ static void mpi_transfer() { if (nrnmpi_numprocs > 1) { double wt = nrnmpi_wtime(); if (nrn_sparse_partrans > 0) { - nrnmpi_dbl_alltoallv_sparse( - outsrc_buf_, outsrccnt_, outsrcdspl_, insrc_buf_, insrccnt_, insrcdspl_); + nrnmpi_dbl_alltoallv_sparse(outsrc_buf_, + outsrccnt_.data(), + outsrcdspl_.data(), + insrc_buf_, + insrccnt_.data(), + insrcdspl_.data()); } else { - nrnmpi_dbl_alltoallv( - outsrc_buf_, outsrccnt_, outsrcdspl_, insrc_buf_, insrccnt_, insrcdspl_); + nrnmpi_dbl_alltoallv(outsrc_buf_, + outsrccnt_.data(), + outsrcdspl_.data(), + insrc_buf_, + insrccnt_.data(), + insrcdspl_.data()); } nrnmpi_transfer_wait_ += nrnmpi_wtime() - wt; errno = 0; @@ -577,15 +614,6 @@ static void thread_transfer(NrnThread* _nt) { // " But this was a mistake as many mpi implementations do not allow overlap // of send and receive buffers. -// 22-08-2014 For setup of the All2allv pattern, use the rendezvous rank -// idiom. -#define HAVEWANT_t sgid_t -#define HAVEWANT_alltoallv sgid_alltoallv -#define HAVEWANT2Int MapSgid2Int -#if NRNMPI -#include "have2want.cpp" -#endif - void nrnmpi_setup_transfer() { #if !NRNMPI if (nrnmpi_numprocs > 1) { @@ -611,10 +639,14 @@ void nrnmpi_setup_transfer() { return; } if (nrnmpi_numprocs > 1) { - delete[] std::exchange(insrccnt_, nullptr); - delete[] std::exchange(insrcdspl_, nullptr); - delete[] std::exchange(outsrccnt_, nullptr); - delete[] std::exchange(outsrcdspl_, nullptr); + insrccnt_.clear(); + insrccnt_.shrink_to_fit(); + insrcdspl_.clear(); + insrcdspl_.shrink_to_fit(); + outsrccnt_.clear(); + outsrccnt_.shrink_to_fit(); + outsrcdspl_.clear(); + outsrcdspl_.shrink_to_fit(); // This is an old comment prior to using the want_to_have rendezvous // rank function in want2have.cpp. The old method did not scale // to more sgids than could fit on a single rank, because @@ -649,78 +681,54 @@ void nrnmpi_setup_transfer() { // sids needed by this machine. The 'seen' table values are unused // but the keys are all the (unique) sgid needed by this process. // At the end seen is in fact what we want for sid2insrc_. - int needsrc_cnt = 0; int szalloc = targets_.size(); szalloc = szalloc ? szalloc : 1; // At the moment sid2insrc_ is serving as 'seen' sid2insrc_.clear(); - sid2insrc_.reserve(szalloc); // for single counting - sgid_t* needsrc = new sgid_t[szalloc]; // more than we need + sid2insrc_.reserve(szalloc); // for single counting + std::vector needsrc{}; for (size_t i = 0; i < sgid2targets_.size(); ++i) { sgid_t sid = sgid2targets_[i]; auto search = sid2insrc_.find(sid); if (search == sid2insrc_.end()) { sid2insrc_[sid] = 0; // at the moment, value does not matter - needsrc[needsrc_cnt++] = sid; + needsrc.push_back(sid); } } // 1 continued) Create an array of sources this rank owns. // This already exists as a vector in the SgidList sgids_ but // that is private so go ahead and copy. - sgid_t* ownsrc = new sgid_t[sgids_.size() + 1]; // not 0 length if count is 0 - for (size_t i = 0; i < sgids_.size(); ++i) { - ownsrc[i] = sgids_[i]; - } + std::vector ownsrc = sgids_; // 2) Call the have_to_want function. - sgid_t* send_to_want; - int *send_to_want_cnt, *send_to_want_displ; - sgid_t* recv_from_have; - int *recv_from_have_cnt, *recv_from_have_displ; - - have_to_want(ownsrc, - sgids_.size(), - needsrc, - needsrc_cnt, - send_to_want, - send_to_want_cnt, - send_to_want_displ, - recv_from_have, - recv_from_have_cnt, - recv_from_have_displ, - default_rendezvous); + auto [send_to_want, recv_from_have] = have_to_want(ownsrc, needsrc, sgid_alltoallv); // sanity check. all the sgids we are asked to send, we actually have - int n = send_to_want_displ[nhost]; + int n = send_to_want.displ[nhost]; // sanity check. all the sgids we receive, we actually need. // also set the sid2insrc_ value to the proper recv_from_have index. - n = recv_from_have_displ[nhost]; + n = recv_from_have.displ[nhost]; for (int i = 0; i < n; ++i) { - sgid_t sgid = recv_from_have[i]; + sgid_t sgid = recv_from_have.data[i]; nrn_assert(sid2insrc_.find(sgid) != sid2insrc_.end()); sid2insrc_[sgid] = i; } - // clean up a little - delete[] std::exchange(ownsrc, nullptr); - delete[] std::exchange(needsrc, nullptr); - delete[] std::exchange(recv_from_have, nullptr); - // 3) First return triple creates the proper outsrc_buf_. // Now that we know what machines are interested in our sids... // construct outsrc_buf, outsrc_buf_size, outsrccnt_, outsrcdspl_ // and poutsrc_; - outsrccnt_ = send_to_want_cnt; - outsrcdspl_ = send_to_want_displ; + std::swap(outsrccnt_, send_to_want.cnt); + std::swap(outsrcdspl_, send_to_want.displ); outsrc_buf_size_ = outsrcdspl_[nrnmpi_numprocs]; szalloc = std::max(1, outsrc_buf_size_); outsrc_buf_ = new double[szalloc]; poutsrc_.resize(szalloc); poutsrc_indices_ = new int[szalloc]; for (int i = 0; i < outsrc_buf_size_; ++i) { - sgid_t sid = send_to_want[i]; + sgid_t sid = send_to_want.data[i]; auto search = sgid2srcindex_.find(sid); nrn_assert(search != sgid2srcindex_.end()); Node* nd = visources_[search->second]; @@ -735,12 +743,11 @@ void nrnmpi_setup_transfer() { poutsrc_indices_[i] = search->second; outsrc_buf_[i] = double(sid); // see step 5 } - delete[] send_to_want; // 4) The second triple is creates the insrc_buf_. // From the recv_from_have and sid2insrc_ table, construct the insrc... - insrccnt_ = recv_from_have_cnt; - insrcdspl_ = recv_from_have_displ; + std::swap(insrccnt_, recv_from_have.cnt); + std::swap(insrcdspl_, recv_from_have.displ); insrc_buf_size_ = insrcdspl_[nrnmpi_numprocs]; szalloc = insrc_buf_size_ ? insrc_buf_size_ : 1; insrc_buf_ = new double[szalloc]; From a0f1848c88f0520987679cf8be41675268e239f9 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 9 Jan 2024 10:59:09 +0100 Subject: [PATCH 27/37] Remove __MWERKS__ and FILE_OPEN_RETRY code (#2655) __MWERKS__ is a macro to detect a compiler names CodeWarrior develop by MetroWerks. The last release of this compiler was in 2004. FILE_OPEN_RETRY is unset and unused --- src/coreneuron/mpi/nrnmpiuse.h | 3 --- src/ivoc/ocfile.cpp | 21 --------------------- src/ivoc/pwman.cpp | 8 -------- src/ivoc/rect.h | 15 --------------- src/mswin/mwprefix.h | 4 ---- src/nrniv/kschan.cpp | 4 ---- src/oc/mswinprt.cpp | 2 -- src/oc/nrnmpiuse.h.in | 3 --- 8 files changed, 60 deletions(-) diff --git a/src/coreneuron/mpi/nrnmpiuse.h b/src/coreneuron/mpi/nrnmpiuse.h index 72bcd90092..eec7609094 100644 --- a/src/coreneuron/mpi/nrnmpiuse.h +++ b/src/coreneuron/mpi/nrnmpiuse.h @@ -24,6 +24,3 @@ /* define to the dll path if you want to load automatically */ #undef DLL_DEFAULT_FNAME - -/* Number of times to retry a failed open */ -#undef FILE_OPEN_RETRY diff --git a/src/ivoc/ocfile.cpp b/src/ivoc/ocfile.cpp index ee5b24da7b..61216c9694 100644 --- a/src/ivoc/ocfile.cpp +++ b/src/ivoc/ocfile.cpp @@ -325,14 +325,7 @@ void OcFile::binary_mode() { Use File.seek(0) after opening or use a binary style read/write as first\n\ access to file."); } -#if defined(__MWERKS__) - // printf("can't switch to binary mode. No setmode\n"); - mode_[1] = 'b'; - mode_[2] = '\0'; - file_ = freopen(filename_.c_str(), mode_, file()); -#else setmode(fileno(file()), O_BINARY); -#endif binary_ = true; } } @@ -345,20 +338,6 @@ bool OcFile::open(const char* name, const char* type) { strcpy(mode_, type); #endif file_ = fopen(expand_env_var(name), type); -#if defined(FILE_OPEN_RETRY) && FILE_OPEN_RETRY > 0 - int i; - for (i = 0; !file_ && i < FILE_OPEN_RETRY; ++i) { - // retry occasionally needed on BlueGene - file_ = fopen(expand_env_var(name), type); - } - if (i > 0) { - if (file_) { - printf("%d opened %s after %d retries\n", nrnmpi_myid_world, name, i); - } else { - printf("%d open %s failed after %d retries\n", nrnmpi_myid_world, name, i); - } - } -#endif return is_open(); } diff --git a/src/ivoc/pwman.cpp b/src/ivoc/pwman.cpp index c393eb47a7..898ece83f1 100644 --- a/src/ivoc/pwman.cpp +++ b/src/ivoc/pwman.cpp @@ -3340,13 +3340,6 @@ char* ivoc_get_temp_file() { if (!tdir) { tdir = "/tmp"; } -#if defined(WIN32) && defined(__MWERKS__) - char tname[L_tmpnam + 1]; - tmpnam(tname); - auto const length = strlen(tdir) + 1 + strlen(tname) + 1; - tmpfile = new char[length]; - std::snprintf(tmpfile, length, "%s/%s", tdir, tname); -#else auto const length = strlen(tdir) + 1 + 9 + 1; tmpfile = new char[length]; std::snprintf(tmpfile, length, "%s/nrnXXXXXX", tdir); @@ -3359,7 +3352,6 @@ char* ivoc_get_temp_file() { #else mktemp(tmpfile); #endif -#endif #if defined(WIN32) tmpfile = hoc_back2forward(tmpfile); #endif diff --git a/src/ivoc/rect.h b/src/ivoc/rect.h index 37715f3d0e..74865e1db7 100644 --- a/src/ivoc/rect.h +++ b/src/ivoc/rect.h @@ -36,11 +36,6 @@ class Appear: public Glyph { static const Brush* db_; }; -#if defined(__MWERKS__) -#undef Rect -#define Rect ivoc_Rect -#endif - class Rect: public Appear { public: Rect(Coord left, @@ -63,11 +58,6 @@ class Rect: public Appear { Coord l_, b_, w_, h_; }; -#if defined(__MWERKS__) -#undef Line -#define Line ivoc_Line -#endif - class Line: public Appear { public: Line(Coord dx, Coord dy, const Color* color = NULL, const Brush* brush = NULL); @@ -117,11 +107,6 @@ class Triangle: public Appear { bool filled_; }; -#if defined(__MWERKS__) -#undef Rectangle -#define Rectangle ivoc_Rectangle -#endif - class Rectangle: public Appear { public: Rectangle(float height, diff --git a/src/mswin/mwprefix.h b/src/mswin/mwprefix.h index b28199b663..877b31f9ec 100644 --- a/src/mswin/mwprefix.h +++ b/src/mswin/mwprefix.h @@ -1,9 +1,5 @@ /* auto include file for metrowerks codewarrior for all nrn */ -#if __MWERKS__ >= 7 -#define _MSL_DIRENT_H -#else #include -#endif #pragma once off #ifndef __WIN32__ #define __WIN32__ 1 diff --git a/src/nrniv/kschan.cpp b/src/nrniv/kschan.cpp index f5e2fe0ac4..b3ab59ddc1 100644 --- a/src/nrniv/kschan.cpp +++ b/src/nrniv/kschan.cpp @@ -11,10 +11,6 @@ #include "nrniv_mf.h" #define NSingleIndex 0 -#if defined(__MWERKS__) && !defined(_MSC_VER) -#include -#define strdup _strdup -#endif using KSChanList = std::vector; static KSChanList* channels; diff --git a/src/oc/mswinprt.cpp b/src/oc/mswinprt.cpp index 202a86db82..905aa101b2 100644 --- a/src/oc/mswinprt.cpp +++ b/src/oc/mswinprt.cpp @@ -148,11 +148,9 @@ void hoc_win_exec(void) { void hoc_winio_show(int b) {} -#if !defined(__MWERKS__) int getpid() { return 1; } -#endif void hoc_Plt() { TRY_GUI_REDIRECT_DOUBLE("plt", NULL); diff --git a/src/oc/nrnmpiuse.h.in b/src/oc/nrnmpiuse.h.in index 7f17e5c0b2..4db1c75155 100755 --- a/src/oc/nrnmpiuse.h.in +++ b/src/oc/nrnmpiuse.h.in @@ -16,7 +16,4 @@ /* define if needed */ #undef ALWAYS_CALL_MPI_INIT -/* Number of times to retry a failed open */ -#undef FILE_OPEN_RETRY - #endif From 67f6bedcfc172ea5545af3731c2b636cfd04e922 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 9 Jan 2024 20:32:22 +0100 Subject: [PATCH 28/37] Add python 3.12 for CircleCI for aarch64 wheel (#2656) --- .circleci/config.yml | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index c2d9b65aa7..f17eef0de5 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -49,10 +49,11 @@ jobs: # choose available python versions from pyenv pyenv_py_ver="" case << parameters.NRN_PYTHON_VERSION >> in - 38) pyenv_py_ver="3.8.7" ;; - 39) pyenv_py_ver="3.9.1" ;; - 310) pyenv_py_ver="3.10.1" ;; - 311) pyenv_py_ver="3.11.0" ;; + 38) pyenv_py_ver="3.8.10" ;; + 39) pyenv_py_ver="3.9.13" ;; + 310) pyenv_py_ver="3.10.11" ;; + 311) pyenv_py_ver="3.11.7" ;; + 312) pyenv_py_ver="3.12.0" ;; *) echo "Error: pyenv python version not specified!" && exit 1;; esac @@ -95,7 +96,7 @@ workflows: - /circleci\/.*/ matrix: parameters: - NRN_PYTHON_VERSION: ["311"] + NRN_PYTHON_VERSION: ["312"] NRN_NIGHTLY_UPLOAD: ["false"] nightly: @@ -110,5 +111,5 @@ workflows: - manylinux2014-aarch64: matrix: parameters: - NRN_PYTHON_VERSION: ["38", "39", "310", "311"] + NRN_PYTHON_VERSION: ["38", "39", "310", "311", "312"] NRN_NIGHTLY_UPLOAD: ["true"] From 45cdea7e74b0d695d1597589b47353d728082454 Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Fri, 12 Jan 2024 16:10:37 +0100 Subject: [PATCH 29/37] Fix Apple M1 gitlab CI : ensure pip is setup correctly (#2660) * Fix Apple M1 GitLab CI: ensure pip is setup correctly On MacOS, when we create venv without shims (for sanitizers build) the pip may not be set up. Use `ensurepip` module to make sure the pip is is set up. * Help music for MPI library detection --- .github/workflows/neuron-ci.yml | 5 +++-- .gitlab-ci.yml | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/neuron-ci.yml b/.github/workflows/neuron-ci.yml index 9472e60686..0d0d05cc46 100644 --- a/.github/workflows/neuron-ci.yml +++ b/.github/workflows/neuron-ci.yml @@ -40,7 +40,7 @@ jobs: PY_MIN_VERSION: ${{ matrix.config.python_min_version || '3.8' }} PY_MAX_VERSION: ${{ matrix.config.python_max_version || '3.11' }} MUSIC_INSTALL_DIR: /opt/MUSIC - MUSIC_VERSION: 1.2.0 + MUSIC_VERSION: 1.2.1 strategy: matrix: @@ -173,7 +173,8 @@ jobs: curl -L -o MUSIC.zip https://github.com/INCF/MUSIC/archive/refs/tags/${MUSIC_VERSION}.zip unzip MUSIC.zip && mv MUSIC-* MUSIC && cd MUSIC ./autogen.sh - ./configure --with-python-sys-prefix --prefix=$MUSIC_INSTALL_DIR --disable-anysource + # on some systems MPI library detection fails, provide exact flags/compilers + ./configure --with-python-sys-prefix --prefix=$MUSIC_INSTALL_DIR --disable-anysource MPI_CXXFLAGS="-g -O3" MPI_CFLAGS="-g -O3" MPI_LDFLAGS=" " CC=mpicc CXX=mpicxx make -j install deactivate working-directory: ${{runner.temp}} diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d6baa64fb0..bd944ee463 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -53,6 +53,7 @@ mac_m1_cmake_build: - real_python=$(python3 resolve_shim.py) - echo "python3=$(command -v python3) is really ${real_python}" - PYTHONEXECUTABLE=${real_python} ${real_python} -mvenv venv + - venv/bin/python3 -m ensurepip --upgrade --default-pip - venv/bin/pip install --upgrade pip -r nrn_requirements.txt - git submodule update --init --recursive --force --depth 1 -- external/nmodl - venv/bin/pip install --upgrade -r external/nmodl/requirements.txt From d6b99558d6ae10904d5009abd0950f1195bf9803 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 16 Jan 2024 15:26:25 +0100 Subject: [PATCH 30/37] nonlinz.cpp: use eigen instead of sparse13 (#2623) * diag_ have been removed. Those was pointer for direct access to the diagonal elements, but Eigen can relocate matrix and invalidate pointers. * double* rv_ (real v_) and double* jv_ (imag v_), have been changed to std::vector v_. * char* m_ (the sparse13 matrix) have been changed to a Eigen::SparseMatrix> m_. The solver Eigen::SparseLU lu_ have been added. One of the main difference is that sparse13 is 1-indexed and Eigen: 0-indexed. * The whole complex version of sparse13 have been removed because this was the last use. * transfer_amp, input_amp, transfer_phase, input_phase and ratio_amp have all been simplified by using standard functions of std::complex instead of computing by ourself. * v_index have been removed because this was a vector i -> i + 1. Now that Eigen is 0-indexed it becomes really useless. * In the file partrans.cpp the function pargag_jacob_rhs has been modified to support std::complex, so know we do only one call of the function (before one with real and one with complex part), but inside we loop twice on real and complex because the algorithm is hard to understand to do everything in once. --- cmake/NeuronFileLists.cmake | 2 - src/nrniv/nonlinz.cpp | 308 ++++++++++++++---------------------- src/nrniv/nonlinz.h | 12 +- src/nrniv/partrans.cpp | 87 +++++----- src/sparse13/CMakeLists.txt | 16 +- src/sparse13/cspalloc.cpp | 2 - src/sparse13/cspbuild.cpp | 2 - src/sparse13/cspfactor.cpp | 2 - src/sparse13/cspmatrix.h | 2 - src/sparse13/cspoutput.cpp | 2 - src/sparse13/cspredef.h | 49 ------ src/sparse13/cspsolve.cpp | 2 - src/sparse13/csputils.cpp | 2 - 13 files changed, 183 insertions(+), 305 deletions(-) delete mode 100644 src/sparse13/cspalloc.cpp delete mode 100644 src/sparse13/cspbuild.cpp delete mode 100644 src/sparse13/cspfactor.cpp delete mode 100644 src/sparse13/cspmatrix.h delete mode 100644 src/sparse13/cspoutput.cpp delete mode 100644 src/sparse13/cspredef.h delete mode 100644 src/sparse13/cspsolve.cpp delete mode 100644 src/sparse13/csputils.cpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 8cd9a06553..a3e88a85f9 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -68,8 +68,6 @@ set(HEADER_FILES_TO_INSTALL scopmath/sparse_thread.hpp scopmath/ssimplic.hpp scopmath/ssimplic_thread.hpp - sparse13/cspmatrix.h - sparse13/cspredef.h sparse13/spconfig.h sparse13/spmatrix.h) diff --git a/src/nrniv/nonlinz.cpp b/src/nrniv/nonlinz.cpp index 65bb00dd5f..f0a24eac8d 100644 --- a/src/nrniv/nonlinz.cpp +++ b/src/nrniv/nonlinz.cpp @@ -6,56 +6,58 @@ #include "nrniv_mf.h" #include "nrnoc2iv.h" #include "nrnmpi.h" -#include "cspmatrix.h" #include "membfunc.h" +#include +#include + +using namespace std::complex_literals; + extern void v_setup_vectors(); extern int nrndae_extra_eqn_count(); extern Symlist* hoc_built_in_symlist; extern void (*nrnthread_v_transfer_)(NrnThread*); -extern spREAL* spGetElement(char*, int, int); -extern void pargap_jacobi_rhs(double*, double*); +extern void pargap_jacobi_rhs(std::vector>&, + const std::vector>&); extern void pargap_jacobi_setup(int mode); class NonLinImpRep { public: NonLinImpRep(); - virtual ~NonLinImpRep(); void delta(double); + + // Functions to fill the matrix void didv(); void dids(); void dsdv(); void dsds(); + int gapsolve(); - char* m_; + // Matrix containing the non linear system to solve. + Eigen::SparseMatrix> m_{}; + // The solver of the matrix using the LU decomposition method. + Eigen::SparseLU>> lu_{}; int scnt_; // structure_change int n_v_, n_ext_, n_lin_, n_ode_, neq_v_, neq_; std::vector> pv_, pvdot_; - int* v_index_; - double* rv_; - double* jv_; - double** diag_; - double* deltavec_; // just like cvode.atol*cvode.atolscale for ode's - double delta_; // slightly more efficient and easier for v. + std::vector> v_; + std::vector deltavec_; // just like cvode.atol*cvode.atolscale for ode's + double delta_; // slightly more efficient and easier for v. void current(int, Memb_list*, int); void ode(int, Memb_list*); double omega_; int iloc_; // current injection site of last solve - float* vsymtol_; - int maxiter_; + float* vsymtol_{}; + int maxiter_{500}; }; -NonLinImp::NonLinImp() { - rep_ = NULL; -} NonLinImp::~NonLinImp() { - if (rep_) { - delete rep_; - } + delete rep_; } + double NonLinImp::transfer_amp(int curloc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) { hoc_execerror( @@ -64,9 +66,7 @@ double NonLinImp::transfer_amp(int curloc, int vloc) { if (curloc != rep_->iloc_) { solve(curloc); } - double x = rep_->rv_[vloc]; - double y = rep_->jv_[vloc]; - return sqrt(x * x + y * y); + return std::abs(rep_->v_[vloc]); } double NonLinImp::input_amp(int curloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -78,9 +78,7 @@ double NonLinImp::input_amp(int curloc) { if (curloc < 0) { return 0.0; } - double x = rep_->rv_[curloc]; - double y = rep_->jv_[curloc]; - return sqrt(x * x + y * y); + return std::abs(rep_->v_[curloc]); } double NonLinImp::transfer_phase(int curloc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_ && curloc != rep_->iloc_) { @@ -90,9 +88,7 @@ double NonLinImp::transfer_phase(int curloc, int vloc) { if (curloc != rep_->iloc_) { solve(curloc); } - double x = rep_->rv_[vloc]; - double y = rep_->jv_[vloc]; - return atan2(y, x); + return std::arg(rep_->v_[vloc]); } double NonLinImp::input_phase(int curloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -104,9 +100,7 @@ double NonLinImp::input_phase(int curloc) { if (curloc < 0) { return 0.0; } - double x = rep_->rv_[curloc]; - double y = rep_->jv_[curloc]; - return atan2(y, x); + return std::arg(rep_->v_[curloc]); } double NonLinImp::ratio_amp(int clmploc, int vloc) { if (nrnmpi_numprocs > 1 && nrnthread_v_transfer_) { @@ -118,22 +112,14 @@ double NonLinImp::ratio_amp(int clmploc, int vloc) { if (clmploc != rep_->iloc_) { solve(clmploc); } - double ax, bx, cx, ay, by, cy, bb; - ax = rep_->rv_[vloc]; - ay = rep_->jv_[vloc]; - bx = rep_->rv_[clmploc]; - by = rep_->jv_[clmploc]; - bb = bx * bx + by * by; - cx = (ax * bx + ay * by) / bb; - cy = (ay * bx - ax * by) / bb; - return sqrt(cx * cx + cy * cy); + return std::abs(rep_->v_[vloc] * std::conj(rep_->v_[clmploc]) / std::norm(rep_->v_[clmploc])); } void NonLinImp::compute(double omega, double deltafac, int maxiter) { v_setup_vectors(); nrn_rhs(nrn_ensure_model_data_are_sorted(), nrn_threads[0]); if (rep_ && rep_->scnt_ != structure_change_cnt) { delete rep_; - rep_ = NULL; + rep_ = nullptr; } if (!rep_) { rep_ = new NonLinImpRep(); @@ -151,8 +137,10 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->omega_ = 1000. * omega; rep_->delta(deltafac); + + rep_->m_.setZero(); + // fill matrix - cmplx_spClear(rep_->m_); rep_->didv(); rep_->dsds(); #if 1 // when 0 equivalent to standard method @@ -160,18 +148,32 @@ void NonLinImp::compute(double omega, double deltafac, int maxiter) { rep_->dsdv(); #endif - // cmplx_spPrint(rep_->m_, 0, 1, 1); - // for (int i=0; i < rep_->neq_; ++i) { - // printf("i=%d %g %g\n", i, rep_->diag_[i][0], rep_->diag_[i][1]); - // } - int e = cmplx_spFactor(rep_->m_); - switch (e) { - case spZERO_DIAG: - hoc_execerror("cmplx_spFactor error:", "Zero Diagonal"); - case spNO_MEMORY: - hoc_execerror("cmplx_spFactor error:", "No Memory"); - case spSINGULAR: - hoc_execerror("cmplx_spFactor error:", "Singular"); + // Now that the matrix is filled we can compress it (mandatory for SparseLU) + rep_->m_.makeCompressed(); + + // Factorize the matrix so this is ready to solve + rep_->lu_.compute(rep_->m_); + switch (rep_->lu_.info()) { + case Eigen::Success: + // Everything fine + break; + case Eigen::NumericalIssue: + hoc_execerror( + "Eigen Sparse LU factorization failed with Eigen::NumericalIssue, please check the " + "input matrix:", + rep_->lu_.lastErrorMessage().c_str()); + break; + case Eigen::NoConvergence: + hoc_execerror( + "Eigen Sparse LU factorization reports Eigen::NonConvergence after calling compute():", + rep_->lu_.lastErrorMessage().c_str()); + break; + case Eigen::InvalidInput: + hoc_execerror( + "Eigen Sparse LU factorization failed with Eigen::InvalidInput, the input matrix seems " + "invalid:", + rep_->lu_.lastErrorMessage().c_str()); + break; } rep_->iloc_ = -2; @@ -184,20 +186,18 @@ int NonLinImp::solve(int curloc) { hoc_execerror("Must call Impedance.compute first", 0); } if (rep_->iloc_ != curloc) { - int i; rep_->iloc_ = curloc; - for (i = 0; i < rep_->neq_; ++i) { - rep_->rv_[i] = 0; - rep_->jv_[i] = 0; - } + rep_->v_ = std::vector>(rep_->neq_); if (curloc >= 0) { - rep_->rv_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); + rep_->v_[curloc] = 1.e2 / NODEAREA(_nt->_v_node[curloc]); } if (nrnthread_v_transfer_) { rval = rep_->gapsolve(); } else { - assert(rep_->m_); - cmplx_spSolve(rep_->m_, rep_->rv_ - 1, rep_->rv_ - 1, rep_->jv_ - 1, rep_->jv_ - 1); + auto v = + Eigen::Map, Eigen::Dynamic>>(rep_->v_.data(), + rep_->v_.size()); + v = rep_->lu_.solve(v); } } return rval; @@ -207,13 +207,8 @@ int NonLinImp::solve(int curloc) { // mapping is already done there. NonLinImpRep::NonLinImpRep() { - int err; - int i, j, ieq, cnt; NrnThread* _nt = nrn_threads; - maxiter_ = 500; - m_ = NULL; - vsymtol_ = NULL; Symbol* vsym = hoc_table_lookup("v", hoc_built_in_symlist); if (vsym->extra) { vsymtol_ = &vsym->extra->tolerance; @@ -233,10 +228,10 @@ NonLinImpRep::NonLinImpRep() { n_ode_ = 0; for (NrnThreadMembList* tml = _nt->tml; tml; tml = tml->next) { Memb_list* ml = tml->ml; - i = tml->index; + int i = tml->index; nrn_ode_count_t s = memb_func[i].ode_count; if (s) { - cnt = (*s)(i); + int cnt = (*s)(i); n_ode_ += cnt * ml->nodecount; } } @@ -245,46 +240,21 @@ NonLinImpRep::NonLinImpRep() { if (neq_ == 0) { return; } - m_ = cmplx_spCreate(neq_, 1, &err); - assert(err == spOKAY); + m_.resize(neq_, neq_); pv_.resize(neq_); pvdot_.resize(neq_); - v_index_ = new int[n_v_]; - rv_ = new double[neq_ + 1]; - rv_ += 1; - jv_ = new double[neq_ + 1]; - jv_ += 1; - diag_ = new double*[neq_]; - deltavec_ = new double[neq_]; - - for (i = 0; i < n_v_; ++i) { + v_.resize(neq_); + deltavec_.resize(neq_); + + for (int i = 0; i < n_v_; ++i) { // utilize nd->eqn_index in case of use_sparse13 later Node* nd = _nt->_v_node[i]; pv_[i] = nd->v_handle(); pvdot_[i] = nd->rhs_handle(); - v_index_[i] = i + 1; - } - for (i = 0; i < n_v_; ++i) { - diag_[i] = cmplx_spGetElement(m_, v_index_[i], v_index_[i]); - } - for (i = neq_v_; i < neq_; ++i) { - diag_[i] = cmplx_spGetElement(m_, i + 1, i + 1); } scnt_ = structure_change_cnt; } -NonLinImpRep::~NonLinImpRep() { - if (!m_) { - return; - } - cmplx_spDestroy(m_); - delete[] v_index_; - delete[](rv_ - 1); - delete[](jv_ - 1); - delete[] diag_; - delete[] deltavec_; -} - void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for ode int i, nc, cnt, ieq; NrnThread* nt = nrn_threads; @@ -299,8 +269,12 @@ void NonLinImpRep::delta(double deltafac) { // also defines pv_,pvdot_ map for if (nrn_ode_count_t s = memb_func[i].ode_count; s && (cnt = s(i)) > 0) { nrn_ode_map_t ode_map = memb_func[i].ode_map; for (auto j = 0; j < nc; ++j) { - ode_map( - ml->prop[j], ieq, pv_.data() + ieq, pvdot_.data() + ieq, deltavec_ + ieq, i); + ode_map(ml->prop[j], + ieq, + pv_.data() + ieq, + pvdot_.data() + ieq, + deltavec_.data() + ieq, + i); ieq += cnt; } } @@ -317,19 +291,17 @@ void NonLinImpRep::didv() { for (i = _nt->ncell; i < n_v_; ++i) { nd = _nt->_v_node[i]; ip = _nt->_v_parent[i]->v_node_index; - double* a = cmplx_spGetElement(m_, v_index_[ip], v_index_[i]); - double* b = cmplx_spGetElement(m_, v_index_[i], v_index_[ip]); - *a += NODEA(nd); - *b += NODEB(nd); - *diag_[i] -= NODEB(nd); - *diag_[ip] -= NODEA(nd); + m_.coeffRef(ip, i) += NODEA(nd); + m_.coeffRef(i, ip) += NODEB(nd); + m_.coeffRef(i, i) -= NODEB(nd); + m_.coeffRef(ip, ip) -= NODEA(nd); } // jwC term Memb_list* mlc = _nt->tml->ml; int n = mlc->nodecount; for (i = 0; i < n; ++i) { j = mlc->nodelist[i]->v_node_index; - diag_[v_index_[j] - 1][1] += .001 * mlc->data(i, 0) * omega_; + m_.coeffRef(j, j) += .001 * mlc->data(i, 0) * omega_ * 1i; } // di/dv terms // because there may be several point processes of the same type @@ -370,7 +342,7 @@ void NonLinImpRep::didv() { current(i, ml, j); // conductance // add to matrix - *diag_[v_index_[nd->v_node_index] - 1] -= (x2 - NODERHS(nd)) / delta_; + m_.coeffRef(nd->v_node_index, nd->v_node_index) -= (x2 - NODERHS(nd)) / delta_; } } } @@ -390,8 +362,6 @@ void NonLinImpRep::dids() { nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); if (memb_func[i].current) { - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; // zero rhs @@ -399,22 +369,20 @@ void NonLinImpRep::dids() { // compute rhs current(i, ml, in); // save rhs - x2[in] = NODERHS(nd); + v_[in].imag(NODERHS(nd)); // each state incremented separately and restored for (iis = 0; iis < cnt; ++iis) { is = ieq + in * cnt + iis; // save s - x1[is] = *pv_[is]; + v_[is].real(*pv_[is]); // increment s and zero rhs *pv_[is] += deltavec_[is]; NODERHS(nd) = 0; current(i, ml, in); - *pv_[is] = x1[is]; // restore s - double g = (NODERHS(nd) - x2[in]) / deltavec_[is]; + *pv_[is] = v_[is].real(); // restore s + double g = (NODERHS(nd) - v_[in].imag()) / deltavec_[is]; if (g != 0.) { - double* elm = - cmplx_spGetElement(m_, v_index_[nd->v_node_index], is + 1); - elm[0] = -g; + m_.coeffRef(nd->v_node_index, is) = -g; } } // don't know if this is necessary but make sure last @@ -439,22 +407,20 @@ void NonLinImpRep::dsdv() { nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); if (memb_func[i].current) { - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; // zero rhs, save v for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { *pvdot_[is] = 0.; } - x1[in] = NODEV(nd); + v_[in].real(NODEV(nd)); } // increment v only once in case there are multiple // point processes at the same location for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; auto const v = nd->v(); - if (x1[in] == v) { + if (v_[in].real() == v) { nd->v() = v + delta_; } } @@ -464,10 +430,10 @@ void NonLinImpRep::dsdv() { for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - x2[is] = *pvdot_[is]; + v_[is].imag(*pvdot_[is]); *pvdot_[is] = 0; } - nd->v() = x1[in]; + nd->v() = v_[in].real(); } // compute the rhs(v) ode(i, ml); @@ -475,11 +441,9 @@ void NonLinImpRep::dsdv() { for (in = 0; in < ml->nodecount; ++in) { Node* nd = ml->nodelist[in]; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - double ds = (x2[is] - *pvdot_[is]) / delta_; + double ds = (v_[is].imag() - *pvdot_[is]) / delta_; if (ds != 0.) { - double* elm = - cmplx_spGetElement(m_, is + 1, v_index_[nd->v_node_index]); - elm[0] = -ds; + m_.coeffRef(is, nd->v_node_index) = -ds; } } } @@ -494,7 +458,7 @@ void NonLinImpRep::dsds() { NrnThread* nt = nrn_threads; // jw term for (i = neq_v_; i < neq_; ++i) { - diag_[i][1] += omega_; + m_.coeffRef(i, i) += omega_ * 1i; } ieq = neq_v_; for (NrnThreadMembList* tml = nt->tml; tml; tml = tml->next) { @@ -504,13 +468,11 @@ void NonLinImpRep::dsds() { int nc = ml->nodecount; nrn_ode_count_t s = memb_func[i].ode_count; int cnt = (*s)(i); - double* x1 = rv_; // use as temporary storage - double* x2 = jv_; // zero rhs, save s for (in = 0; in < ml->nodecount; ++in) { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { *pvdot_[is] = 0.; - x1[is] = *pv_[is]; + v_[is].real(*pv_[is]); } } // compute rhs. this is the rhs(s) @@ -518,7 +480,7 @@ void NonLinImpRep::dsds() { // save rhs for (in = 0; in < ml->nodecount; ++in) { for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - x2[is] = *pvdot_[is]; + v_[is].imag(*pvdot_[is]); } } // iterate over the states @@ -538,12 +500,11 @@ void NonLinImpRep::dsds() { Node* nd = ml->nodelist[in]; ks = ieq + in * cnt + kks; for (is = ieq + in * cnt, iis = 0; iis < cnt; ++iis, ++is) { - double ds = (*pvdot_[is] - x2[is]) / deltavec_[is]; + double ds = (*pvdot_[is] - v_[is].imag()) / deltavec_[is]; if (ds != 0.) { - double* elm = cmplx_spGetElement(m_, is + 1, ks + 1); - elm[0] = -ds; + m_.coeffRef(is, ks) = -ds; } - *pv_[ks] = x1[ks]; + *pv_[ks] = v_[ks].real(); } } // perhaps not necessary but ensures the last computation is with @@ -573,9 +534,11 @@ void NonLinImpRep::ode(int im, Memb_list* ml) { // assume there is in fact an o memb_func[im].ode_spec(nrn_ensure_model_data_are_sorted(), nrn_threads, ml, im); } +// This function compute a solution of a converging system by iteration. +// The value returned is the number of iterations to reach a precision of "tol" (1e-9). int NonLinImpRep::gapsolve() { - // On entry, rv_ and jv_ contain the complex b for A*x = b. - // On return rv_ and jv_ contain complex solution, x. + // On entry, v_ contains the complex b for A*x = b. + // On return v_ contains complex solution, x. // m_ is the factored matrix for the trees without gap junctions // Jacobi method (easy for parallel) // A = D + R @@ -592,48 +555,40 @@ int NonLinImpRep::gapsolve() { } #endif - pargap_jacobi_setup(0); + pargap_jacobi_setup(0); // 0 means 'setup' - double *rx, *jx, *rx1, *jx1, *rb, *jb; - if (neq_) { - rx = new double[neq_]; - jx = new double[neq_]; - rx1 = new double[neq_]; - jx1 = new double[neq_]; - rb = new double[neq_]; - jb = new double[neq_]; - } - - // initialize for first iteration - for (int i = 0; i < neq_; ++i) { - rx[i] = jx[i] = 0.0; - rb[i] = rv_[i]; - jb[i] = jv_[i]; - } + std::vector> x_old(neq_); + std::vector> x(neq_); + std::vector> b(v_); // iterate till change in x is small double tol = 1e-9; - double delta; + double delta{}; int success = 0; int iter; for (iter = 1; iter <= maxiter_; ++iter) { if (neq_) { - cmplx_spSolve(m_, rb - 1, rx1 - 1, jb - 1, jx1 - 1); + auto b_ = Eigen::Map, Eigen::Dynamic>>(b.data(), + b.size()); + auto x_ = Eigen::Map, Eigen::Dynamic>>(x.data(), + x.size()); + x_ = lu_.solve(b_); } // if any change in x > tol, then do another iteration. success = 1; delta = 0.0; + // Do the substraction of the previous result (x_old) and current result (x). + // If all differences are < tol stop the loop, otherwise continue to iterate for (int i = 0; i < neq_; ++i) { - double err = fabs(rx1[i] - rx[i]) + fabs(jx1[i] - jx[i]); + auto diff = x[i] - x_old[i]; + double err = std::abs(diff.real()) + std::abs(diff.imag()); if (err > tol) { success = 0; } - if (delta < err) { - delta = err; - } + delta = std::max(err, delta); } #if NRNMPI if (nrnmpi_numprocs > 1) { @@ -641,34 +596,17 @@ int NonLinImpRep::gapsolve() { } #endif if (success) { - for (int i = 0; i < neq_; ++i) { - rv_[i] = rx1[i]; - jv_[i] = jx1[i]; - } + v_ = x; break; } // setup for next iteration - for (int i = 0; i < neq_; ++i) { - rx[i] = rx1[i]; - jx[i] = jx1[i]; - rb[i] = rv_[i]; - jb[i] = jv_[i]; - } - pargap_jacobi_rhs(rb, rx); - pargap_jacobi_rhs(jb, jx); + x_old = x; + b = v_; + pargap_jacobi_rhs(b, x_old); } - pargap_jacobi_setup(1); // tear down - - if (neq_) { - delete[] rx; - delete[] jx; - delete[] rx1; - delete[] jx1; - delete[] rb; - delete[] jb; - } + pargap_jacobi_setup(1); // 1 means 'tear down' if (!success) { char buf[256]; @@ -678,7 +616,7 @@ int NonLinImpRep::gapsolve() { maxiter_, delta, tol); - execerror(buf, 0); + hoc_execerror(buf, nullptr); } return iter; } diff --git a/src/nrniv/nonlinz.h b/src/nrniv/nonlinz.h index 115a191b27..b477bf783a 100644 --- a/src/nrniv/nonlinz.h +++ b/src/nrniv/nonlinz.h @@ -3,20 +3,24 @@ class NonLinImpRep; +// A solver for non linear equation of complex numbers. +// Matrix should be squared. class NonLinImp { public: - NonLinImp(); - virtual ~NonLinImp(); + ~NonLinImp(); + // Prepare the matrix before solving it. void compute(double omega, double deltafac, int maxiter); - double transfer_amp(int curloc, int vloc); // v_node[arg] is the node + + double transfer_amp(int curloc, int vloc); double transfer_phase(int curloc, int vloc); double input_amp(int curloc); double input_phase(int curloc); double ratio_amp(int clmploc, int vloc); + int solve(int curloc); private: - NonLinImpRep* rep_; + NonLinImpRep* rep_{}; }; #endif diff --git a/src/nrniv/partrans.cpp b/src/nrniv/partrans.cpp index 1ebfdf335d..a6779bad9e 100644 --- a/src/nrniv/partrans.cpp +++ b/src/nrniv/partrans.cpp @@ -13,6 +13,7 @@ #include #include +#include #include // Replaces NrnHash for MapSgid2Int #include #include @@ -880,47 +881,59 @@ void pargap_jacobi_setup(int mode) { } } -void pargap_jacobi_rhs(double* b, double* x) { - // helper for complex impedance with parallel gap junctions - // b = b - R*x R are the off diagonal gap elements of the jacobian. - // we presume 1 thread. First nrn_thread[0].end equations are in node order. - if (!nrnthread_v_transfer_) { - return; - } +void pargap_jacobi_rhs(std::vector>& b, + const std::vector>& x) { + // First loop for real, second for imag + for (int real_imag = 0; real_imag < 2; ++real_imag) { + // helper for complex impedance with parallel gap junctions + // b = b - R*x R are the off diagonal gap elements of the jacobian. + // we presume 1 thread. First nrn_thread[0].end equations are in node order. + if (!nrnthread_v_transfer_) { + return; + } - NrnThread* _nt = nrn_threads; + NrnThread* _nt = nrn_threads; - // transfer gap node voltages to gap vpre - for (size_t i = 0; i < visources_.size(); ++i) { - Node* nd = visources_[i]; - nd->v() = x[nd->v_node_index]; - } - mpi_transfer(); - thread_transfer(_nt); + // transfer gap node voltages to gap vpre + for (size_t i = 0; i < visources_.size(); ++i) { + Node* nd = visources_[i]; + if (real_imag == 0) { + nd->v() = x[nd->v_node_index].real(); + } else { + nd->v() = x[nd->v_node_index].imag(); + } + } + mpi_transfer(); + thread_transfer(_nt); - // set gap node voltages to 0 so we can use nrn_cur to set rhs - for (size_t i = 0; i < visources_.size(); ++i) { - Node* nd = visources_[i]; - nd->v() = 0.0; - } - auto const sorted_token = nrn_ensure_model_data_are_sorted(); - auto* const vec_rhs = _nt->node_rhs_storage(); - // Initialize rhs to 0. - for (int i = 0; i < _nt->end; ++i) { - vec_rhs[i] = 0.0; - } - for (int k = 0; k < imped_current_type_count_; ++k) { - int type = imped_current_type_[k]; - Memb_list* ml = imped_current_ml_[k]; - memb_func[type].current(sorted_token, _nt, ml, type); - } + // set gap node voltages to 0 so we can use nrn_cur to set rhs + for (size_t i = 0; i < visources_.size(); ++i) { + Node* nd = visources_[i]; + nd->v() = 0.0; + } + auto const sorted_token = nrn_ensure_model_data_are_sorted(); + auto* const vec_rhs = _nt->node_rhs_storage(); + // Initialize rhs to 0. + for (int i = 0; i < _nt->end; ++i) { + vec_rhs[i] = 0.0; + } + for (int k = 0; k < imped_current_type_count_; ++k) { + int type = imped_current_type_[k]; + Memb_list* ml = imped_current_ml_[k]; + memb_func[type].current(sorted_token, _nt, ml, type); + } - // possibly many gap junctions in same node (and possible even different - // types) but rhs is the accumulation of all those instances at each node - // so ... The only thing that can go wrong is if there are intances of - // gap junctions that are not being used (not in the target list). - for (int i = 0; i < _nt->end; ++i) { - b[i] += vec_rhs[i]; + // possibly many gap junctions in same node (and possible even different + // types) but rhs is the accumulation of all those instances at each node + // so ... The only thing that can go wrong is if there are intances of + // gap junctions that are not being used (not in the target list). + for (int i = 0; i < _nt->end; ++i) { + if (real_imag == 0) { + b[i] += vec_rhs[i]; + } else { + b[i] += std::complex(0, vec_rhs[i]); + } + } } } diff --git a/src/sparse13/CMakeLists.txt b/src/sparse13/CMakeLists.txt index 69f9e1a955..b2adc31dc2 100644 --- a/src/sparse13/CMakeLists.txt +++ b/src/sparse13/CMakeLists.txt @@ -1,15 +1,3 @@ -add_library( - sparse13 STATIC - spalloc.cpp - spbuild.cpp - spfactor.cpp - spoutput.cpp - spsolve.cpp - sputils.cpp - cspalloc.cpp - cspbuild.cpp - cspfactor.cpp - cspoutput.cpp - cspsolve.cpp - csputils.cpp) +add_library(sparse13 STATIC spalloc.cpp spbuild.cpp spfactor.cpp spoutput.cpp spsolve.cpp + sputils.cpp) set_property(TARGET sparse13 PROPERTY POSITION_INDEPENDENT_CODE ON) diff --git a/src/sparse13/cspalloc.cpp b/src/sparse13/cspalloc.cpp deleted file mode 100644 index 29b7ee44c8..0000000000 --- a/src/sparse13/cspalloc.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spalloc.cpp" diff --git a/src/sparse13/cspbuild.cpp b/src/sparse13/cspbuild.cpp deleted file mode 100644 index b002c1e5c2..0000000000 --- a/src/sparse13/cspbuild.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spbuild.cpp" diff --git a/src/sparse13/cspfactor.cpp b/src/sparse13/cspfactor.cpp deleted file mode 100644 index 2a9d91aaaa..0000000000 --- a/src/sparse13/cspfactor.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spfactor.cpp" diff --git a/src/sparse13/cspmatrix.h b/src/sparse13/cspmatrix.h deleted file mode 100644 index 7200d07b55..0000000000 --- a/src/sparse13/cspmatrix.h +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spmatrix.h" diff --git a/src/sparse13/cspoutput.cpp b/src/sparse13/cspoutput.cpp deleted file mode 100644 index 59325c189a..0000000000 --- a/src/sparse13/cspoutput.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spoutput.cpp" diff --git a/src/sparse13/cspredef.h b/src/sparse13/cspredef.h deleted file mode 100644 index 2df7687583..0000000000 --- a/src/sparse13/cspredef.h +++ /dev/null @@ -1,49 +0,0 @@ -/* mostly generated from -cat temp | /usr/bin/tr -cs "[:alpha:]" "[\n*]" | sort | uniq | grep '^sp' > sp_redef.h -where temp is the last part of spmatrix.h -*/ -#define spClear cmplx_spClear -#define spCondition cmplx_spCondition -#define spCreate cmplx_spCreate -#define spDeleteRowAndCol cmplx_spDeleteRowAndCol -#define spDestroy cmplx_spDestroy -#define spDeterminant cmplx_spDeterminant -#define spElementCount cmplx_spElementCount -#define spError cmplx_spError -#define spFactor cmplx_spFactor -#define spFileMatrix cmplx_spFileMatrix -#define spFileStats cmplx_spFileStats -#define spFileVector cmplx_spFileVector -#define spFillinCount cmplx_spFillinCount -#define spGetAdmittance cmplx_spGetAdmittance -#define spGetElement cmplx_spGetElement -#define spGetInitInfo cmplx_spGetInitInfo -#define spGetOnes cmplx_spGetOnes -#define spGetQuad cmplx_spGetQuad -#define spGetSize cmplx_spGetSize -#define spInitialize cmplx_spInitialize -#define spInstallInitInfo cmplx_spInstallInitInfo -#define spLargestElement cmplx_spLargestElement -#define spMNA_Preorder cmplx_spMNA_Preorder -#define spMultTransposed cmplx_spMultTransposed -#define spMultiply cmplx_spMultiply -#define spNorm cmplx_spNorm -#define spOrderAndFactor cmplx_spOrderAndFactor -#define spPartition cmplx_spPartition -#define spPrint cmplx_spPrint -#define spPseudoCondition cmplx_spPseudoCondition -#define spRoundoff cmplx_spRoundoff -#define spScale cmplx_spScale -#define spSetComplex cmplx_spSetComplex -#define spSetReal cmplx_spSetReal -#define spSolve cmplx_spSolve -#define spSolveTransposed cmplx_spSolveTransposed -#define spStripFills cmplx_spStripFills -#define spWhereSingular cmplx_spWhereSingular -#define spcGetFillin cmplx_spcGetFillin -#define spcGetElement cmplx_spcGetElement -#define spcLinkRows cmplx_spcLinkRows -#define spcFindElementInCol cmplx_spcFindElementInCol -#define spcCreateElement cmplx_spcCreateElement -#define spcRowExchange cmplx_spcRowExchange -#define spcColExchange cmplx_spcColExchange diff --git a/src/sparse13/cspsolve.cpp b/src/sparse13/cspsolve.cpp deleted file mode 100644 index f1aeba7d67..0000000000 --- a/src/sparse13/cspsolve.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "spsolve.cpp" diff --git a/src/sparse13/csputils.cpp b/src/sparse13/csputils.cpp deleted file mode 100644 index 0f0363b8ae..0000000000 --- a/src/sparse13/csputils.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#define cmplx_spPrefix -#include "sputils.cpp" From 3b0ffc9ff8d32d384f628c75dd630f396a9e3664 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 16 Jan 2024 17:32:27 +0100 Subject: [PATCH 31/37] Remove useless geometry3d (#2666) --- cmake/NeuronFileLists.cmake | 1 - src/nrniv/geometry3d.cpp | 824 ------------------------------------ 2 files changed, 825 deletions(-) delete mode 100644 src/nrniv/geometry3d.cpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index a3e88a85f9..30f36f5052 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -220,7 +220,6 @@ set(NRNIV_FILE_LIST cxprop.cpp datapath.cpp finithnd.cpp - geometry3d.cpp glinerec.cpp hocmech.cpp impedanc.cpp diff --git a/src/nrniv/geometry3d.cpp b/src/nrniv/geometry3d.cpp deleted file mode 100644 index e0ffe7839a..0000000000 --- a/src/nrniv/geometry3d.cpp +++ /dev/null @@ -1,824 +0,0 @@ -/* - This file contains code adapted from p3d.py in - http://code.google.com/p/pythonisosurfaces/source/checkout - which was released under the new BSD license. - accessed 31 July 2012 -*/ - -#include -#include -#include -#include -#include -//#include - -int geometry3d_find_triangles(double value0, - double value1, - double value2, - double value3, - double value4, - double value5, - double value6, - double value7, - double x0, - double x1, - double y0, - double y1, - double z0, - double z1, - double* out, - int offset); - -double geometry3d_llgramarea(double* p0, double* p1, double* p2); -double geometry3d_sum_area_of_triangles(double* tri_vec, int len); - -const int edgeTable[] = { - 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, - 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, - 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, - 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0xaa, - 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, - 0x569, 0x663, 0x76a, 0x66, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, - 0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, - 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c, - 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, - 0x2cf, 0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9, - 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, - 0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, - 0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, - 0x3f5, 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, - 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3, - 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, - 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, - 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, - 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, - 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}; - -/* CTNG:tritable */ -const int triTable[256][16] = {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, - {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, - {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, - {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, - {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, - {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, - {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, - {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, - {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, - {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, - {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, - {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, - {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, - {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, - {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, - {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, - {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, - {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, - {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, - {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, - {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, - {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, - {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, - {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, - {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, - {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, - {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, - {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, - {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, - {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, - {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, - {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, - {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, - {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, - {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, - {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, - {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, - {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, - {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, - {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, - {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, - {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, - {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, - {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, - {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, - {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, - {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, - {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, - {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, - {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, - {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, - {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, - {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, - {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, - {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, - {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, - {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, - {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, - {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, - {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, - {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, - {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, - {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, - {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, - {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, - {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, - {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, - {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, - {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, - {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, - {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, - {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, - {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, - {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, - {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, - {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, - {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, - {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, - {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, - {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, - {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, - {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, - {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, - {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, - {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, - {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, - {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, - {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, - {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, - {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, - {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, - {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, - {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, - {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, - {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, - {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, - {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, - {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, - {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, - {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, - {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, - {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, - {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, - {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, - {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, - {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, - {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, - {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, - {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, - {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, - {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, - {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, - {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, - {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, - {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, - {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, - {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, - {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, - {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, - {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, - {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, - {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, - {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, - {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, - {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, - {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, - {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, - {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, - {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, - {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, - {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, - {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, - {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, - {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, - {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, - {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, - {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, - {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, - {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, - {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, - {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, - {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, - {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, - {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, - {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, - {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, - {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, - {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, - {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, - {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, - {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, - {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, - {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, - {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, - {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, - {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, - {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, - {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, - {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, - {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, - {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, - {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, - {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, - {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, - {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, - {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, - {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, - {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, - {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, - {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, - {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, - {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, - {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, - {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, - {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}; - - -/* CTNG:interpedge */ -void geometry3d_vi(double* p1, double* p2, double v1, double v2, double* out) { - if (fabs(v1) < 0.000000000001) { - out[0] = p1[0]; - out[1] = p1[1]; - out[2] = p1[2]; - return; - } - if (fabs(v2) < 0.000000000001) { - out[0] = p2[0]; - out[1] = p2[1]; - out[2] = p2[2]; - return; - } - double delta_v = v1 - v2; - if (fabs(delta_v) < 0.0000000001) { - out[0] = p1[0]; - out[1] = p1[1]; - out[2] = p1[2]; - return; - } - double mu = v1 / delta_v; - if (std::isnan(mu)) { - printf("geometry3d_vi error. delta_v = %g, v1 = %g, v2 = %g\n", delta_v, v1, v2); - } - out[0] = p1[0] + mu * (p2[0] - p1[0]); - out[1] = p1[1] + mu * (p2[1] - p1[1]); - out[2] = p1[2] + mu * (p2[2] - p1[2]); -} - -int geometry3d_find_triangles(double value0, - double value1, - double value2, - double value3, - double value4, - double value5, - double value6, - double value7, - double x0, - double x1, - double y0, - double y1, - double z0, - double z1, - double* out, - int offset) { - out = out + offset; - double position[8][3] = {{x0, y0, z0}, - {x1, y0, z0}, - {x1, y1, z0}, - {x0, y1, z0}, - {x0, y0, z1}, - {x1, y0, z1}, - {x1, y1, z1}, - {x0, y1, z1}}; - - /* CTNG:domarch */ - int cubeIndex = 0; - if (value0 < 0) - cubeIndex |= 1; - if (value1 < 0) - cubeIndex |= 2; - if (value2 < 0) - cubeIndex |= 4; - if (value3 < 0) - cubeIndex |= 8; - if (value4 < 0) - cubeIndex |= 16; - if (value5 < 0) - cubeIndex |= 32; - if (value6 < 0) - cubeIndex |= 64; - if (value7 < 0) - cubeIndex |= 128; - - int et = edgeTable[cubeIndex]; - - if (et == 0) - return 0; - - double vertexList[12][3]; - if (et & 1) - geometry3d_vi(position[0], position[1], value0, value1, vertexList[0]); - if (et & 2) - geometry3d_vi(position[1], position[2], value1, value2, vertexList[1]); - if (et & 4) - geometry3d_vi(position[2], position[3], value2, value3, vertexList[2]); - if (et & 8) - geometry3d_vi(position[3], position[0], value3, value0, vertexList[3]); - if (et & 16) - geometry3d_vi(position[4], position[5], value4, value5, vertexList[4]); - if (et & 32) - geometry3d_vi(position[5], position[6], value5, value6, vertexList[5]); - if (et & 64) - geometry3d_vi(position[6], position[7], value6, value7, vertexList[6]); - if (et & 128) - geometry3d_vi(position[7], position[4], value7, value4, vertexList[7]); - if (et & 256) - geometry3d_vi(position[0], position[4], value0, value4, vertexList[8]); - if (et & 512) - geometry3d_vi(position[1], position[5], value1, value5, vertexList[9]); - if (et & 1024) - geometry3d_vi(position[2], position[6], value2, value6, vertexList[10]); - if (et & 2048) - geometry3d_vi(position[3], position[7], value3, value7, vertexList[11]); - - int const* const tt = triTable[cubeIndex]; - int i, j, k, count; - for (i = 0, count = 0; i < 16; i += 3, count++) { - if (tt[i] == -1) - break; - for (k = 0; k < 3; k++) { - for (j = 0; j < 3; j++) - out[j] = vertexList[tt[i + k]][j]; - out += 3; - } - } - return count; -} - - -double geometry3d_llgramarea(double* p0, double* p1, double* p2) { - /* setup the vectors */ - double a[] = {p0[0] - p1[0], p0[1] - p1[1], p0[2] - p1[2]}; - double b[] = {p0[0] - p2[0], p0[1] - p2[1], p0[2] - p2[2]}; - - /* take the cross-product */ - double cpx = a[1] * b[2] - a[2] * b[1]; - double cpy = a[2] * b[0] - a[0] * b[2]; - double cpz = a[0] * b[1] - a[1] * b[0]; - return std::sqrt(cpx * cpx + cpy * cpy + cpz * cpz); -} - - -double geometry3d_sum_area_of_triangles(double* tri_vec, int len) { - double sum = 0; - for (int i = 0; i < len; i += 9) { - sum += geometry3d_llgramarea(tri_vec + i, tri_vec + i + 3, tri_vec + i + 6); - } - return sum / 2.; -} - - -/***************************************************************************** - * this contains cone, and cylinder code translated and modified - * from Barbier and Galin 2004's example code - * see /u/ramcd/spatial/experiments/one_time_tests/2012-06-28/cone.cpp - * on 7 june 2012, the original code was available online at - * http://jgt.akpeters.com/papers/BarbierGalin04/Cone-Sphere.zip - ****************************************************************************/ - - -class geometry3d_Cylinder { - public: - geometry3d_Cylinder(double x0, double y0, double z0, double x1, double y1, double z1, double r); - //~geometry3d_Cylinder(); - double signed_distance(double px, double py, double pz); - - private: - // double x0, y0, z0, x1, y1, z1, r; - double r, rr, axisx, axisy, axisz, cx, cy, cz, h; -}; - -geometry3d_Cylinder::geometry3d_Cylinder(double x0, - double y0, - double z0, - double x1, - double y1, - double z1, - double r) - : r(r) - , cx((x0 + x1) / 2.) - , cy((y0 + y1) / 2.) - , cz((z0 + z1) / 2.) - , rr(r * r) { - axisx = x1 - x0; - axisy = y1 - y0; - axisz = z1 - z0; - double axislength = std::sqrt(axisx * axisx + axisy * axisy + axisz * axisz); - axisx /= axislength; - axisy /= axislength; - axisz /= axislength; - h = axislength / 2.; -} - -double geometry3d_Cylinder::signed_distance(double px, double py, double pz) { - double const nx{px - cx}; - double const ny{py - cy}; - double const nz{pz - cz}; - double y{std::abs(axisx * nx + axisy * ny + axisz * nz)}; - double yy{y * y}; - double const xx{nx * nx + ny * ny + nz * nz - yy}; - if (y < h) { - return std::max(-std::abs(h - y), std::sqrt(xx) - r); - } else { - y -= h; - yy = y * y; - if (xx < rr) { - return std::abs(y); - } else { - double const x{std::sqrt(xx) - r}; - return std::sqrt(yy + x * x); - } - } -} - -void* geometry3d_new_Cylinder(double x0, - double y0, - double z0, - double x1, - double y1, - double z1, - double r) { - return new geometry3d_Cylinder(x0, y0, z0, x1, y1, z1, r); -} -void geometry3d_delete_Cylinder(void* ptr) { - delete (geometry3d_Cylinder*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Cylinder_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Cylinder*) ptr)->signed_distance(px, py, pz); -} - - -class geometry3d_Cone { - public: - geometry3d_Cone(double x0, - double y0, - double z0, - double r0, - double x1, - double y1, - double z1, - double r1); - double signed_distance(double px, double py, double pz); - - private: - double axisx, axisy, axisz, h, rra, rrb, conelength; - double side1, side2, x0, y0, z0, r0, axislength; -}; - -geometry3d_Cone::geometry3d_Cone(double x0, - double y0, - double z0, - double r0, - double x1, - double y1, - double z1, - double r1) - : rra(r0 * r0) - , rrb(r1 * r1) - , x0(x0) - , y0(y0) - , z0(z0) - , r0(r0) { - // TODO: these are preconditions; the python assures them, but could/should - // take care of that here - assert(r1 <= r0); - assert(r1 >= 0); - axisx = x1 - x0; - axisy = y1 - y0; - axisz = z1 - z0; - axislength = std::sqrt(axisx * axisx + axisy * axisy + axisz * axisz); - axisx /= axislength; - axisy /= axislength; - axisz /= axislength; - h = axislength / 2.; - rra = r0 * r0; - rrb = r1 * r1; - conelength = std::sqrt((r1 - r0) * (r1 - r0) + axislength * axislength); - side1 = (r1 - r0) / conelength; - side2 = axislength / conelength; -} - -double geometry3d_Cone::signed_distance(double px, double py, double pz) { - double nx, ny, nz, y, yy, xx, x, rx, ry; - nx = px - x0; - ny = py - y0; - nz = pz - z0; - y = axisx * nx + axisy * ny + axisz * nz; - yy = y * y; - xx = nx * nx + ny * ny + nz * nz - yy; - // in principle, xx >= 0, but roundoff errors may cause trouble - if (xx < 0) - xx = 0; - if (y < 0) { - if (xx < rra) - return -y; - x = std::sqrt(xx) - r0; - return std::sqrt(x * x + yy); - } else if (xx < rrb) { - return y - axislength; - } else { - x = std::sqrt(xx) - r0; - if (y < 0) { - if (x < 0) - return y; - return std::sqrt(x * x + yy); - } else { - ry = x * side1 + y * side2; - if (ry < 0) - return std::sqrt(x * x + yy); - rx = x * side2 - y * side1; - if (ry > conelength) { - ry -= conelength; - return std::sqrt(rx * rx + ry * ry); - } else { - return rx; - } - } - } -} - - -void* geometry3d_new_Cone(double x0, - double y0, - double z0, - double r0, - double x1, - double y1, - double z1, - double r1) { - return new geometry3d_Cone(x0, y0, z0, r0, x1, y1, z1, r1); -} -void geometry3d_delete_Cone(void* ptr) { - delete (geometry3d_Cone*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Cone_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Cone*) ptr)->signed_distance(px, py, pz); -} - - -class geometry3d_Sphere { - public: - geometry3d_Sphere(double x, double y, double z, double r); - double signed_distance(double px, double py, double pz); - - private: - double x, y, z, r; -}; - -geometry3d_Sphere::geometry3d_Sphere(double x, double y, double z, double r) - : x(x) - , y(y) - , z(z) - , r(r) {} - -double geometry3d_Sphere::signed_distance(double px, double py, double pz) { - return std::sqrt(std::pow(x - px, 2) + std::pow(y - py, 2) + std::pow(z - pz, 2)) - r; -} - -void* geometry3d_new_Sphere(double x, double y, double z, double r) { - return new geometry3d_Sphere(x, y, z, r); -} -void geometry3d_delete_Sphere(void* ptr) { - delete (geometry3d_Sphere*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Sphere_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Sphere*) ptr)->signed_distance(px, py, pz); -} - -class geometry3d_Plane { - public: - geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz); - double signed_distance(double px, double py, double pz); - - private: - double nx, ny, nz; - double d, mul; -}; - -geometry3d_Plane::geometry3d_Plane(double x, double y, double z, double nx, double ny, double nz) - : nx(nx) - , ny(ny) - , nz(nz) - , d(-(nx * x + ny * y + nz * z)) - , mul(1. / std::sqrt(nx * nx + ny * ny + nz * nz)) {} - -double geometry3d_Plane::signed_distance(double px, double py, double pz) { - return (nx * px + ny * py + nz * pz + d) * mul; -} - -void* geometry3d_new_Plane(double x, double y, double z, double nx, double ny, double nz) { - return new geometry3d_Plane(x, y, z, nx, ny, nz); -} -void geometry3d_delete_Plane(void* ptr) { - delete (geometry3d_Plane*) ptr; -} -// TODO: add validation for ptr -double geometry3d_Plane_signed_distance(void* ptr, double px, double py, double pz) { - return ((geometry3d_Plane*) ptr)->signed_distance(px, py, pz); -} - -/* - void print_numbers(PyObject *p) { - for (Py_ssize_t i = 0; i< PyList_Size(p); i++) { - PyObject* obj = PyList_GetItem(p, i); - printf("%g ", PyFloat_AsDouble(obj)); - } - printf("\n"); - } - - // TODO: it would be nice to remove the python dependence, because that - // limits us to mostly single threaded due to the global interpreter - // lock - int geometry3d_process_cell(int i, int j, int k, PyObject* objects, double* xs, double* ys, -double* zs, double* tridata, int start) { double x, y, z, x1, y1, z1, xx, yy, zz; x = xs[i]; y = -ys[j]; z = zs[k]; x1 = xs[i + 1]; y1 = ys[j + 1]; z1 = zs[k + 1]; double value[8], current_value; - PyObject* result; - PyObject* obj; - printf("inside process_cell\n"); - - // march around the cube - for (int m = 0; m < 8; m++) { - // NOTE: only describing changes from previous case - switch(m) { - case 0: xx = x; yy = y; zz = z; break; - case 1: xx = x1; break; - case 2: yy = y1; break; - case 3: xx = x; break; - case 4: yy = y; zz = z1; break; - case 5: xx = x1; break; - case 6: yy = y1; break; - case 7: xx = x; break; - } - printf("phase 0, len(objects) = %ld\n", PyList_Size(objects)); - obj = PyList_GetItem(objects, 0); - printf("phase 0a, obj = %x\n", (void*) obj); - result = PyEval_CallMethod(obj, "distance", "ddd", xx, yy, zz); - //Py_DECREF(obj); - printf("phase 1\n"); - current_value = PyFloat_AsDouble(result); - //Py_DECREF(result); - for (Py_ssize_t n = 1; n < PyList_Size(objects); n++) { - printf("phase 2, n = %ld\n", n); - obj = PyList_GetItem(objects, n); - result = PyEval_CallMethod(obj, "distance", "ddd", xx, yy, zz); - Py_DECREF(obj); - current_value = min(current_value, PyFloat_AsDouble(result)); - //Py_DECREF(result); - } - value[m] = current_value; - } - printf("finishing up; start = %d\n", start); - return start + 9 * geometry3d_find_triangles(value[0], value[1], value[2], value[3], -value[4], value[5], value[6], value[7], x, x1, y, y1, z, z1, tridata, start); - } - - - int geometry3d_test_call_function(PyObject* obj) { - printf("inside\n"); - if (obj == NULL) printf("obj is NULL\n"); - Py_INCREF(obj); - PyEval_CallObject(obj, obj); - return 0; - } - - int geometry3d_test_call_function3(int (*func) (void)) { - return func(); - } - - int geometry3d_test_call_method(PyObject* list, PyObject* obj) { - PyEval_CallMethod(list, "insert", "O", obj); - return 0; - } - - // this works, but is slower than the python version - int geometry3d_contains_surface(int i, int j, int k, double (*objdist)(double, double, double), -double* xs, double* ys, double* zs, double dx, double r_inner, double r_outer) { bool has_neg = -false, has_pos = false; double xbar = xs[i] + dx / 2; double ybar = ys[j] + dx / 2; double zbar = -zs[k] + dx / 2; - - double d = fabs(objdist(xbar, ybar, zbar)); - //if (i == 586 && j == 2169 && k == 83) {printf("at magic point, d=%g\n", d);} - //printf("sphere test: d = %g, r_inner = %g, r_outer = %g\n", d, r_inner, r_outer); - if (d <= r_inner) return 1; - if (d >= r_outer) return 0; -// // spheres alone are indeterminant; check corners -// for (int di = 0; di < 2; di++) { -// for (int dj = 0; dj < 2; dj++) { -// for (int dk = 0; dk < 2; dk++) { -// d = objdist(xs[i + di], ys[j + dj], zs[k + dk]); -// //printf("d = %g\n", d); -// if (d <= 0) has_neg = true; -// if (d >= 0) has_pos = true; -// if (has_neg && has_pos) return 1; -// } -// } -// } - - // spheres alone are indeterminant; check corners - for (double* x = xs + i; x < xs + i + 2; x++) { - for (double* y = ys + j; y < ys + j + 2; y++) { - for (double* z = zs + k; z < zs + k + 2; z++) { - d = objdist(*x, *y, *z); - if (d <= 0) has_neg = true; - if (d >= 0) has_pos = true; - if (has_neg && has_pos) return 1; - } - } - } - - return 0; - } - -} -*/ From 4a821d0adda1ee1a9842d3ed2d2a3c89ad6637a7 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Wed, 17 Jan 2024 22:50:20 +0100 Subject: [PATCH 32/37] Windows ci: add support for python 3.12 (#2671) --- ci/win_build_cmake.sh | 2 +- ci/win_download_deps.cmd | 1 + ci/win_install_deps.cmd | 3 +++ nrn_requirements.txt | 3 --- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/ci/win_build_cmake.sh b/ci/win_build_cmake.sh index 3bf55bc842..c3c84ad37e 100755 --- a/ci/win_build_cmake.sh +++ b/ci/win_build_cmake.sh @@ -32,7 +32,7 @@ cd $BUILD_SOURCESDIRECTORY/build -DNRN_BINARY_DIST_BUILD=ON \ -DPYTHON_EXECUTABLE=/c/Python38/python.exe \ -DNRN_ENABLE_PYTHON_DYNAMIC=ON \ - -DNRN_PYTHON_DYNAMIC='c:/Python38/python.exe;c:/Python39/python.exe;c:/Python310/python.exe;c:/Python311/python.exe' \ + -DNRN_PYTHON_DYNAMIC='c:/Python38/python.exe;c:/Python39/python.exe;c:/Python310/python.exe;c:/Python311/python.exe;c:/Python312/python.exe' \ -DCMAKE_INSTALL_PREFIX='/c/nrn-install' \ -DMPI_CXX_LIB_NAMES:STRING=msmpi \ -DMPI_C_LIB_NAMES:STRING=msmpi \ diff --git a/ci/win_download_deps.cmd b/ci/win_download_deps.cmd index fc58ed4254..a006ffc0d2 100644 --- a/ci/win_download_deps.cmd +++ b/ci/win_download_deps.cmd @@ -7,6 +7,7 @@ pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.8.exe htt pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.9.exe https://www.python.org/ftp/python/3.9.0/python-3.9.0-amd64.exe || goto :error pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.10.exe https://www.python.org/ftp/python/3.10.0/python-3.10.0-amd64.exe || goto :error pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.11.exe https://www.python.org/ftp/python/3.11.1/python-3.11.1-amd64.exe || goto :error +pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile python-3.12.exe https://www.python.org/ftp/python/3.12.1/python-3.12.1-amd64.exe || goto :error :: mpi pwsh -command Invoke-WebRequest -MaximumRetryCount 4 -OutFile msmpisetup.exe https://download.microsoft.com/download/a/5/2/a5207ca5-1203-491a-8fb8-906fd68ae623/msmpisetup.exe || goto :error diff --git a/ci/win_install_deps.cmd b/ci/win_install_deps.cmd index 189287b966..2284b911c0 100644 --- a/ci/win_install_deps.cmd +++ b/ci/win_install_deps.cmd @@ -7,6 +7,7 @@ python-3.8.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustFo python-3.9.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python39 || goto :error python-3.10.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python310 || goto :error python-3.11.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python311 || goto :error +python-3.12.exe /passive Include_pip=1 Include_test=0 PrependPath=1 DefaultJustForMeTargetDir=C:\Python312 || goto :error :: fix msvcc version for all python3 pwsh -command "(Get-Content C:\Python38\Lib\distutils\cygwinccompiler.py) -replace 'elif msc_ver == ''1600'':', 'elif msc_ver == ''1916'':' | Out-File C:\Python38\Lib\distutils\cygwinccompiler.py" @@ -26,6 +27,8 @@ C:\Python38\python.exe -m pip install numpy==1.17.5 "cython < 3" || goto :error C:\Python39\python.exe -m pip install numpy==1.19.3 "cython < 3" || goto :error C:\Python310\python.exe -m pip install numpy==1.21.3 "cython < 3" || goto :error C:\Python311\python.exe -m pip install numpy==1.23.5 "cython < 3" || goto :error +C:\Python312\python.exe -m pip install numpy==1.26.3 "cython < 3" || goto :error +C:\Python312\python.exe -m pip install setuptools || goto :error :: install nsis nsis-3.05-setup.exe /S || goto :error diff --git a/nrn_requirements.txt b/nrn_requirements.txt index 06c1f1920a..77ccb32ddc 100644 --- a/nrn_requirements.txt +++ b/nrn_requirements.txt @@ -13,6 +13,3 @@ pytest-cov mpi4py numpy find_libpython - -# For python 3.12 -setuptools From fbacf4f850543aef7695ea78ffc6a9c8c23c9bfa Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Thu, 18 Jan 2024 10:19:06 +0100 Subject: [PATCH 33/37] CircleCI: use 'default' as machine (#2673) --- .circleci/config.yml | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index f17eef0de5..6ec6f56ba6 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,7 +1,7 @@ version: 2.1 orbs: - python: circleci/python@0.3.2 + python: circleci/python@2.1.1 jobs: manylinux2014-aarch64: @@ -13,7 +13,7 @@ jobs: type: string machine: - image: ubuntu-2004:202101-01 + image: default resource_class: arm.medium @@ -49,26 +49,18 @@ jobs: # choose available python versions from pyenv pyenv_py_ver="" case << parameters.NRN_PYTHON_VERSION >> in - 38) pyenv_py_ver="3.8.10" ;; - 39) pyenv_py_ver="3.9.13" ;; - 310) pyenv_py_ver="3.10.11" ;; - 311) pyenv_py_ver="3.11.7" ;; - 312) pyenv_py_ver="3.12.0" ;; + 38) pyenv_py_ver="3.8" ;; + 39) pyenv_py_ver="3.9" ;; + 310) pyenv_py_ver="3.10" ;; + 311) pyenv_py_ver="3.11" ;; + 312) pyenv_py_ver="3.12" ;; *) echo "Error: pyenv python version not specified!" && exit 1;; esac - # install python dependencies: .10 is not available pyenv - if [ "<< parameters.NRN_PYTHON_VERSION >>" == "310" ]; then - sudo apt install software-properties-common -y - sudo add-apt-repository ppa:deadsnakes/ppa -y - sudo apt install python3.10 libpython3.10 python3.10-venv - export PYTHON_EXE=$(which python3.10) - else - cd /opt/circleci/.pyenv/plugins/python-build/../.. && git pull && cd - - env PYTHON_CONFIGURE_OPTS="--enable-shared" pyenv install $pyenv_py_ver --force - pyenv global $pyenv_py_ver - export PYTHON_EXE=$(which python) - fi + cd /opt/circleci/.pyenv/plugins/python-build/../.. && git fetch --all && git checkout -B master origin/master && cd - + env PYTHON_CONFIGURE_OPTS="--enable-shared" pyenv install $pyenv_py_ver --force + pyenv global $pyenv_py_ver + export PYTHON_EXE=$(which python) # test wheel packaging/python/test_wheels.sh $PYTHON_EXE $(ls -t wheelhouse/*.whl) From c164100ef00d67c23933c1fce8700239aca07c69 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Fri, 19 Jan 2024 10:24:07 +0100 Subject: [PATCH 34/37] Move memory related functions to memory.[ch]pp (#2677) --- cmake/NeuronFileLists.cmake | 2 + src/oc/hocdec.h | 4 -- src/oc/memory.cpp | 112 ++++++++++++++++++++++++++++++++++++ src/oc/memory.hpp | 18 ++++++ src/oc/oc_ansi.h | 9 +-- src/oc/symbol.cpp | 101 -------------------------------- 6 files changed, 134 insertions(+), 112 deletions(-) create mode 100644 src/oc/memory.cpp create mode 100644 src/oc/memory.hpp diff --git a/cmake/NeuronFileLists.cmake b/cmake/NeuronFileLists.cmake index 30f36f5052..252f3e9ff6 100644 --- a/cmake/NeuronFileLists.cmake +++ b/cmake/NeuronFileLists.cmake @@ -39,6 +39,7 @@ set(HEADER_FILES_TO_INSTALL oc/hocparse.h oc/mcran4.h oc/mech_api.h + oc/memory.hpp oc/nrnapi.h oc/nrnassrt.h oc/nrnisaac.h @@ -101,6 +102,7 @@ set(OC_FILE_LIST hoc_oop.cpp list.cpp math.cpp + memory.cpp mswinprt.cpp nonlin.cpp ocerf.cpp diff --git a/src/oc/hocdec.h b/src/oc/hocdec.h index 113defaf53..944d4e9d72 100644 --- a/src/oc/hocdec.h +++ b/src/oc/hocdec.h @@ -224,10 +224,6 @@ struct HocParmUnits { /* units for symbol values */ #include "oc_ansi.h" -void* emalloc(size_t n); -void* ecalloc(size_t n, size_t size); -void* erealloc(void* ptr, size_t n); - extern Inst *hoc_progp, *hoc_progbase, *hoc_prog, *hoc_prog_parse_recover; extern Inst* hoc_pc; diff --git a/src/oc/memory.cpp b/src/oc/memory.cpp new file mode 100644 index 0000000000..44a6b93e51 --- /dev/null +++ b/src/oc/memory.cpp @@ -0,0 +1,112 @@ +#include "memory.hpp" + +#include +#include + +// For hoc_warning and hoc_execerror +#include "oc_ansi.h" + +#if HAVE_POSIX_MEMALIGN +#define HAVE_MEMALIGN 1 +#endif +#if defined(DARWIN) /* posix_memalign seems not to work on Darwin 10.6.2 */ +#undef HAVE_MEMALIGN +#endif +#if HAVE_MEMALIGN +#undef _XOPEN_SOURCE /* avoid warnings about redefining this */ +#define _XOPEN_SOURCE 600 +#endif + +static bool emalloc_error = false; + +void* hoc_Emalloc(std::size_t n) { /* check return from malloc */ + void* p = std::malloc(n); + if (p == nullptr) { + emalloc_error = true; + } + return p; +} + +void* hoc_Ecalloc(std::size_t n, std::size_t size) { /* check return from calloc */ + if (n == 0) { + return nullptr; + } + void* p = std::calloc(n, size); + if (p == nullptr) { + emalloc_error = true; + } + return p; +} + +void* hoc_Erealloc(void* ptr, std::size_t size) { /* check return from realloc */ + if (!ptr) { + return hoc_Emalloc(size); + } + void* p = std::realloc(ptr, size); + if (p == nullptr) { + std::free(ptr); + emalloc_error = true; + } + return p; +} + +void hoc_malchk(void) { + if (emalloc_error) { + emalloc_error = false; + hoc_execerror("out of memory", nullptr); + } +} + +void* emalloc(std::size_t n) { + void* p = hoc_Emalloc(n); + if (emalloc_error) { + hoc_malchk(); + } + return p; +} + +void* ecalloc(std::size_t n, std::size_t size) { + void* p = hoc_Ecalloc(n, size); + if (emalloc_error) { + hoc_malchk(); + } + return p; +} + +void* erealloc(void* ptr, std::size_t size) { + void* p = hoc_Erealloc(ptr, size); + if (emalloc_error) { + hoc_malchk(); + } + return p; +} + +void* nrn_cacheline_alloc(void** memptr, std::size_t size) { +#if HAVE_MEMALIGN + static bool memalign_is_working = true; + if (memalign_is_working) { + if (posix_memalign(memptr, 64, size) != 0) { + hoc_warning("posix_memalign not working, falling back to using malloc\n", nullptr); + memalign_is_working = false; + *memptr = hoc_Emalloc(size); + hoc_malchk(); + } + } else +#endif + { + *memptr = hoc_Emalloc(size); + hoc_malchk(); + } + return *memptr; +} + +void* nrn_cacheline_calloc(void** memptr, std::size_t nmemb, std::size_t size) { +#if HAVE_MEMALIGN + nrn_cacheline_alloc(memptr, nmemb * size); + std::memset(*memptr, 0, nmemb * size); +#else + *memptr = hoc_Ecalloc(nmemb, size); + hoc_malchk(); +#endif + return *memptr; +} diff --git a/src/oc/memory.hpp b/src/oc/memory.hpp new file mode 100644 index 0000000000..38960e0eba --- /dev/null +++ b/src/oc/memory.hpp @@ -0,0 +1,18 @@ +#pragma once + +// Some functions here are prepend with 'hoc_' but they are unrelated to hoc itself. +#include + +/* check return from malloc */ +void* hoc_Emalloc(std::size_t n); +void* hoc_Ecalloc(std::size_t n, std::size_t size); +void* hoc_Erealloc(void* ptr, std::size_t size); + +void hoc_malchk(void); + +void* emalloc(std::size_t n); +void* ecalloc(std::size_t n, std::size_t size); +void* erealloc(void* ptr, std::size_t size); + +void* nrn_cacheline_alloc(void** memptr, std::size_t size); +void* nrn_cacheline_calloc(void** memptr, std::size_t nmemb, std::size_t size); diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index 399a86e408..34cc37062f 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -8,6 +8,8 @@ #include #include #include + +#include "memory.hpp" /** * \dir * \brief HOC Interpreter @@ -49,9 +51,6 @@ void ivoc_help(const char*); Symbol* hoc_lookup(const char*); -void* hoc_Ecalloc(std::size_t nmemb, std::size_t size); -void* hoc_Emalloc(size_t size); -void hoc_malchk(); [[noreturn]] void hoc_execerror(const char*, const char*); [[noreturn]] void hoc_execerr_ext(const char* fmt, ...); char* hoc_object_name(Object*); @@ -320,10 +319,6 @@ void hoc_obj_set(int i, Object*); void nrn_hoc_lock(); void nrn_hoc_unlock(); -void* hoc_Erealloc(void* ptr, std::size_t size); - -void* nrn_cacheline_alloc(void** memptr, std::size_t size); -void* nrn_cacheline_calloc(void** memptr, std::size_t nmemb, std::size_t size); [[noreturn]] void nrn_exit(int); void hoc_free_list(Symlist**); int hoc_errno_check(); diff --git a/src/oc/symbol.cpp b/src/oc/symbol.cpp index 7257a91aa3..605a7ba2a6 100644 --- a/src/oc/symbol.cpp +++ b/src/oc/symbol.cpp @@ -2,17 +2,6 @@ /* /local/src/master/nrn/src/oc/symbol.cpp,v 1.9 1999/02/25 18:01:58 hines Exp */ /* version 7.2.1 2-jan-89 */ -#if HAVE_POSIX_MEMALIGN -#define HAVE_MEMALIGN 1 -#endif -#if defined(DARWIN) /* posix_memalign seems not to work on Darwin 10.6.2 */ -#undef HAVE_MEMALIGN -#endif -#if HAVE_MEMALIGN -#undef _XOPEN_SOURCE /* avoid warnings about redefining this */ -#define _XOPEN_SOURCE 600 -#endif - #include "hoc.h" #include "hocdec.h" #include "hoclist.h" @@ -173,96 +162,6 @@ void hoc_link_symbol(Symbol* sp, Symlist* list) { sp->next = nullptr; } -static int emalloc_error = 0; - -void hoc_malchk(void) { - if (emalloc_error) { - emalloc_error = 0; - execerror("out of memory", nullptr); - } -} - -void* hoc_Emalloc(size_t n) { /* check return from malloc */ - void* p = malloc(n); - if (p == nullptr) - emalloc_error = 1; - return p; -} - -void* emalloc(size_t n) { - void* p = hoc_Emalloc(n); - if (emalloc_error) { - hoc_malchk(); - } - return p; -} - -void* hoc_Ecalloc(size_t n, size_t size) { /* check return from calloc */ - if (n == 0) { - return nullptr; - } - void* p = calloc(n, size); - if (p == nullptr) - emalloc_error = 1; - return p; -} - -void* ecalloc(size_t n, size_t size) { - void* p = hoc_Ecalloc(n, size); - if (emalloc_error) { - hoc_malchk(); - } - return p; -} - -void* nrn_cacheline_alloc(void** memptr, size_t size) { -#if HAVE_MEMALIGN - static int memalign_is_working = 1; - if (memalign_is_working) { - if (posix_memalign(memptr, 64, size) != 0) { - fprintf(stderr, "posix_memalign not working, falling back to using malloc\n"); - memalign_is_working = 0; - *memptr = hoc_Emalloc(size); - hoc_malchk(); - } - } else -#endif - *memptr = hoc_Emalloc(size); - hoc_malchk(); - return *memptr; -} - -void* nrn_cacheline_calloc(void** memptr, size_t nmemb, size_t size) { -#if HAVE_MEMALIGN - nrn_cacheline_alloc(memptr, nmemb * size); - memset(*memptr, 0, nmemb * size); -#else - *memptr = hoc_Ecalloc(nmemb, size); - hoc_malchk(); -#endif - return *memptr; -} - -void* hoc_Erealloc(void* ptr, size_t size) { /* check return from realloc */ - if (!ptr) { - return hoc_Emalloc(size); - } - void* p = realloc(ptr, size); - if (p == nullptr) { - free(ptr); - emalloc_error = 1; - } - return p; -} - -void* erealloc(void* ptr, size_t size) { - void* p = hoc_Erealloc(ptr, size); - if (emalloc_error) { - hoc_malchk(); - } - return p; -} - void hoc_free_symspace(Symbol* s1) { /* frees symbol space. Marks it UNDEF */ if (s1 && s1->cpublic != 2) { switch (s1->type) { From d59d45865bd66ea162d01100567697f4f69a8f38 Mon Sep 17 00:00:00 2001 From: JCGoran Date: Fri, 19 Jan 2024 11:49:48 +0100 Subject: [PATCH 35/37] Fix RX3D test (#2676) --- share/lib/python/neuron/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/share/lib/python/neuron/__init__.py b/share/lib/python/neuron/__init__.py index df12031798..c48fdd56ef 100644 --- a/share/lib/python/neuron/__init__.py +++ b/share/lib/python/neuron/__init__.py @@ -1217,7 +1217,7 @@ def _get_color(variable, val, cmap, lo, hi, val_range): elif val > hi: col = color_to_hex(cmap(255)) else: - val = color_to_hex(128) + col = color_to_hex(cmap(128)) else: col = color_to_hex( cmap(int(255 * (min(max(val, lo), hi) - lo) / (val_range))) From d7750508d713877bd55e25186fa3b56df615ae16 Mon Sep 17 00:00:00 2001 From: Pramod Kumbhar Date: Mon, 22 Jan 2024 08:31:27 +0100 Subject: [PATCH 36/37] Avoid issue with GPU builds with the latest NVHPC (>23.1) (#2680) * @matz-e pointed out that the issue #2563 would be solved by initializing the key with {0, 0}. * As value initialization is supposed to work, we believe this is a regression in NVHPC. --- src/coreneuron/utils/randoms/nrnran123.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreneuron/utils/randoms/nrnran123.cpp b/src/coreneuron/utils/randoms/nrnran123.cpp index a618312a34..2f1b12cb3d 100644 --- a/src/coreneuron/utils/randoms/nrnran123.cpp +++ b/src/coreneuron/utils/randoms/nrnran123.cpp @@ -88,7 +88,7 @@ namespace random123_global { #else #define g_k_qualifiers #endif -g_k_qualifiers philox4x32_key_t g_k{}; +g_k_qualifiers philox4x32_key_t g_k{{0, 0}}; // Cannot refer to g_k directly from a nrn_pragma_acc(routine seq) method like // coreneuron_random123_philox4x32_helper, and cannot have this inlined there at From 9db73b4a50336546a8e8f2f766207074619c47dc Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Tue, 23 Jan 2024 20:08:59 +0100 Subject: [PATCH 37/37] Change badge for windows installer (#2684) --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 96d5b467f2..f6cba70ead 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![Build Status](https://dev.azure.com/neuronsimulator/nrn/_apis/build/status/neuronsimulator.nrn?branchName=master)](https://dev.azure.com/neuronsimulator/nrn/_build/latest?definitionId=1&branchName=master) [![Actions Status](https://github.com/neuronsimulator/nrn/workflows/Windows%20Installer/badge.svg)](https://github.com/neuronsimulator/nrn/actions) [![Actions Status](https://github.com/neuronsimulator/nrn/workflows/NEURON%20CI/badge.svg)](https://github.com/neuronsimulator/nrn/actions) [![codecov](https://codecov.io/gh/neuronsimulator/nrn/branch/master/graph/badge.svg?token=T7PIDw6LrC)](https://codecov.io/gh/neuronsimulator/nrn) [![Documentation Status](https://readthedocs.org/projects/nrn/badge/?version=latest)](http://nrn.readthedocs.io/?badge=latest) +[![Build Status](https://dev.azure.com/neuronsimulator/nrn/_apis/build/status/neuronsimulator.nrn?branchName=master)](https://dev.azure.com/neuronsimulator/nrn/_build/latest?definitionId=1&branchName=master) [![Actions Status](https://github.com/neuronsimulator/nrn/actions/workflows/windows.yml/badge.svg?branch=master)](https://github.com/neuronsimulator/nrn/actions) [![Actions Status](https://github.com/neuronsimulator/nrn/workflows/NEURON%20CI/badge.svg)](https://github.com/neuronsimulator/nrn/actions) [![codecov](https://codecov.io/gh/neuronsimulator/nrn/branch/master/graph/badge.svg?token=T7PIDw6LrC)](https://codecov.io/gh/neuronsimulator/nrn) [![Documentation Status](https://readthedocs.org/projects/nrn/badge/?version=latest)](http://nrn.readthedocs.io/?badge=latest) # NEURON NEURON is a simulator for models of neurons and networks of neuron. See [http://neuron.yale.edu](http://neuron.yale.edu) for installers, source code, documentation, tutorials, announcements of