diff --git a/basecode/Eref.cpp b/basecode/Eref.cpp index 5c046458eb..4763902ef9 100644 --- a/basecode/Eref.cpp +++ b/basecode/Eref.cpp @@ -16,6 +16,12 @@ Eref::Eref() ; } +Eref::Eref( const Eref& other ) + : e_( other.e_ ), i_( other.i_ ), f_( other.f_ ) +{ + ; +} + Eref::Eref( Element* e, unsigned int index, unsigned int field ) : e_( e ), i_( index ), f_( field ) { diff --git a/basecode/Eref.h b/basecode/Eref.h index bb71ce251c..d03da1cb23 100644 --- a/basecode/Eref.h +++ b/basecode/Eref.h @@ -16,6 +16,7 @@ class Eref friend ostream& operator <<( ostream& s, const Eref& e ); Eref(); + Eref( const Eref& other ); Eref( Element* e, unsigned int index, unsigned int field = 0 ); /** diff --git a/basecode/ObjId.h b/basecode/ObjId.h index c487f86c8d..1d50b42e4a 100644 --- a/basecode/ObjId.h +++ b/basecode/ObjId.h @@ -33,6 +33,12 @@ class ObjId ; } + ObjId( const ObjId& obj ) + : id( obj.id ), dataIndex( obj.dataIndex ), fieldIndex( obj.fieldIndex ) + { + ; + } + /** * Creates a ObjId using specified Id and DataIndex */ diff --git a/biophysics/CompartmentBase.cpp b/biophysics/CompartmentBase.cpp index 2085589cc6..8c87ceb009 100644 --- a/biophysics/CompartmentBase.cpp +++ b/biophysics/CompartmentBase.cpp @@ -267,6 +267,11 @@ const Cinfo* CompartmentBase::initCinfo() &CompartmentBase::setZ, &CompartmentBase::getZ ); + static ValueFinfo< CompartmentBase, vector< double > > coords( "coords", + "Vector with all coords: [x0 y0 z0 x y z dia]", + &CompartmentBase::setCoords, + &CompartmentBase::getCoords + ); ////////////////////////////////////////////////////////////////// // DestFinfo definitions @@ -322,6 +327,7 @@ const Cinfo* CompartmentBase::initCinfo() &x, // Value &y, // Value &z, // Value + &coords, // Value &injectMsg, // DestFinfo &randInject, // DestFinfo &injectMsg, // DestFinfo @@ -571,6 +577,28 @@ double CompartmentBase::getZ() const return z_; } +void CompartmentBase::setCoords( vector< double > value ) +{ + if (value.size() < 7 ) { + cout << "Warning: CompartmentBase:setCoords. Vector size = " << value.size() << ", require 7 = [ x0 y0 z0 x y z dia ] \n"; + return; + } + x0_ = value[0]; + y0_ = value[1]; + z0_ = value[2]; + x_ = value[3]; + y_ = value[4]; + z_ = value[5]; + diameter_ = value[6]; + updateLength(); +} + +vector< double > CompartmentBase::getCoords() const +{ + vector< double > ret = { x0_, y0_, z0_, x_, y_, z_, diameter_ }; + return ret; +} + ////////////////////////////////////////////////////////////////// // CompartmentBase::Dest function definitions. ////////////////////////////////////////////////////////////////// diff --git a/biophysics/CompartmentBase.h b/biophysics/CompartmentBase.h index 2be0fa05a7..298989cc24 100644 --- a/biophysics/CompartmentBase.h +++ b/biophysics/CompartmentBase.h @@ -57,6 +57,8 @@ class CompartmentBase double getY() const; void setZ( double value ); double getZ() const; + void setCoords( vector< double > value ); + vector< double > getCoords() const; // Dest function definitions. /** diff --git a/builtins/CMakeLists.txt b/builtins/CMakeLists.txt index aefa3964a1..5bde664429 100644 --- a/builtins/CMakeLists.txt +++ b/builtins/CMakeLists.txt @@ -74,6 +74,7 @@ if(WITH_NSDF AND HDF5_FOUND) list(APPEND SRCS HDF5WriterBase.cpp NSDFWriter.cpp + NSDFWriter2.cpp HDF5DataWriter.cpp SpikeStats.cpp testBuiltins.cpp diff --git a/builtins/InputVariable.cpp b/builtins/InputVariable.cpp index 3e4c9a8eee..10b0daadef 100644 --- a/builtins/InputVariable.cpp +++ b/builtins/InputVariable.cpp @@ -50,6 +50,7 @@ #include "hdf5.h" #include "NSDFWriter.h" +#include "NSDFWriter2.h" #include "InputVariable.h" @@ -81,7 +82,7 @@ const Cinfo * InputVariable::initCinfo() static const Cinfo *InputVariableCinfo = InputVariable::initCinfo(); -InputVariable::InputVariable(): owner_(0) +InputVariable::InputVariable(): owner_(0), owner2_(0) { ; } @@ -96,12 +97,18 @@ void InputVariable::setOwner( NSDFWriter * owner) owner_ = owner; } +void InputVariable::setOwner( NSDFWriter2 * owner) +{ + owner2_ = owner; +} + void InputVariable::epSetValue( const Eref& eref, double value) { - if (owner_) - { + if (owner_) { owner_->setInput(eref.fieldIndex(), value); - } + } else if ( owner2_ ) { + owner2_->setInput(eref.fieldIndex(), value); + } } #endif diff --git a/builtins/InputVariable.h b/builtins/InputVariable.h index 552973e275..5b4e416e89 100644 --- a/builtins/InputVariable.h +++ b/builtins/InputVariable.h @@ -52,6 +52,7 @@ #include "Variable.h" class NSDFWriter; +class NSDFWriter2; /** This class is for collecting data and updating a handler object @@ -65,9 +66,11 @@ class InputVariable: public Variable ~InputVariable(); void epSetValue(const Eref &eref, double value); void setOwner(NSDFWriter * owner); + void setOwner(NSDFWriter2 * owner); static const Cinfo * initCinfo(); protected: NSDFWriter * owner_; + NSDFWriter2 * owner2_; }; #endif diff --git a/builtins/NSDFWriter.cpp b/builtins/NSDFWriter.cpp index 413b2b411e..271ba367d8 100644 --- a/builtins/NSDFWriter.cpp +++ b/builtins/NSDFWriter.cpp @@ -6,9 +6,9 @@ // Maintainer: // Created: Thu Jun 18 23:16:11 2015 (-0400) // Version: -// Last-Updated: Sun Dec 20 23:20:19 2015 (-0500) -// By: subha -// Update #: 49 +// Last-Updated: Sat Jan 29 2022 +// By: bhalla +// Update #: 50 // URL: // Keywords: // Compatibility: @@ -22,7 +22,7 @@ // // Change log: -// +// Jan 2022: Many changes added // // // @@ -49,9 +49,9 @@ #include "hdf5.h" #include "hdf5_hl.h" +#include #include #include - #include "../basecode/header.h" #include "../utility/utility.h" #include "../utility/strutil.h" @@ -67,8 +67,11 @@ extern template herr_t writeScalarAttr(hid_t file_id, string path, double value) const char* const EVENTPATH = "/data/event"; const char* const UNIFORMPATH = "/data/uniform"; +const char* const STATICPATH = "/data/static"; const char* const MODELTREEPATH = "/model/modeltree"; +const char* const MODELFILEPATH = "/model/modelfile"; const char* const MAPUNIFORMSRC = "/map/uniform"; +const char* const MAPSTATICSRC = "/map/static"; const char* const MAPEVENTSRC = "/map/event"; string iso_time(time_t * t) @@ -99,10 +102,16 @@ const Cinfo * NSDFWriter::initCinfo() static ValueFinfo modelRoot( "modelRoot", - "The moose element tree root to be saved under /model/modeltree", + "The moose element tree root to be saved under /model/modeltree. If blank, nothing is saved. Default: root object, '/'", &NSDFWriter::setModelRoot, &NSDFWriter::getModelRoot); + static ValueFinfo modelFileNames( + "modelFileNames", + "Comma separated list of model files to save into the NSDF file.", + &NSDFWriter::setModelFiles, + &NSDFWriter::getModelFiles); + static DestFinfo process( "process", "Handle process calls. Collects data in buffer and if number of steps" @@ -125,7 +134,8 @@ const Cinfo * NSDFWriter::initCinfo() processShared, sizeof( processShared ) / sizeof( Finfo* )); static Finfo * finfos[] = { - &eventInputFinfo, + &eventInputFinfo, // FieldElementFinfo + &modelFileNames, // ValueFinfo &proc, }; @@ -267,12 +277,28 @@ void NSDFWriter::openUniformData(const Eref &eref) create the DS for uniform data. */ void NSDFWriter::createUniformMap() +{ + innerCreateMaps( MAPUNIFORMSRC ); +} + +/** + create the DS for static data. + */ +void NSDFWriter::createStaticMap() +{ + innerCreateMaps( MAPSTATICSRC ); +} + + +/** + Generic call for create the DS for static/uniform data + */ +void NSDFWriter::innerCreateMaps( const char* const mapSrcStr ) { // Create the container for all the DS - // TODO: make a common function like `mkdir -p` to avoid repeating this htri_t exists; herr_t status; - hid_t uniformMapContainer = require_group(filehandle_, MAPUNIFORMSRC); + hid_t uniformMapContainer = require_group(filehandle_, mapSrcStr ); // Create the DS themselves for (map< string, vector < unsigned int > >::iterator ii = classFieldToSrcIndex_.begin(); ii != classFieldToSrcIndex_.end(); ++ii){ @@ -280,6 +306,8 @@ void NSDFWriter::createUniformMap() moose::tokenize(ii->first, "/", pathTokens); string className = pathTokens[0]; string fieldName = pathTokens[1]; + if (mapSrcStr == MAPSTATICSRC ) //Hack. for now only static field is coords + fieldName = "coords"; hid_t container = require_group(uniformMapContainer, className); char ** sources = (char **)calloc(ii->second.size(), sizeof(char*)); for (unsigned int jj = 0; jj < ii->second.size(); ++jj){ @@ -291,9 +319,6 @@ void NSDFWriter::createUniformMap() status = H5Tset_size(memtype, H5T_VARIABLE); assert(status >= 0); status = H5Dwrite(ds, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, sources); -#ifndef NDEBUG - cout << "Write dataset: status=" << status << endl; -#endif assert(status >= 0); for (unsigned int jj = 0; jj < ii->second.size(); ++jj){ free(sources[jj]); @@ -534,9 +559,12 @@ void NSDFWriter::reinit(const Eref& eref, const ProcPtr proc) writeScalarAttr< double >(it->second, "dt", proc->dt); } openEventData(eref); + writeModelFiles(); writeModelTree(); createUniformMap(); + createStaticMap(); createEventMap(); + writeStaticCoords(); steps_ = 0; } @@ -627,9 +655,97 @@ string NSDFWriter::getModelRoot() const return modelRoot_; } +void NSDFWriter::setModelFiles(string value) +{ + modelFileNames_.clear(); + moose::tokenize( value, ", ", modelFileNames_); +} + +string NSDFWriter::getModelFiles() const +{ + string ret = ""; + string spacer = ""; + for( auto s = modelFileNames_.begin(); s!= modelFileNames_.end(); ++s) { + ret += spacer + *s; + spacer = ","; + } + return ret; +} + +void NSDFWriter::writeStaticCoords() +{ + hid_t staticObjContainer = require_group(filehandle_, STATICPATH ); + for (map< string, vector < unsigned int > >::iterator ii = classFieldToSrcIndex_.begin(); ii != classFieldToSrcIndex_.end(); ++ii){ + vector < string > pathTokens; + moose::tokenize(ii->first, "/", pathTokens); + string className = pathTokens[0]; + string fieldName = "coords"; // pathTokens[1] is not relevant. + hid_t container = require_group(staticObjContainer, className); + double * buffer = (double*)calloc(ii->second.size() * 7, sizeof(double)); + // Ugly class checking stuff here: Both have a coord field + if ( className.find( "Pool" ) != string::npos || + className.find( "Compartment" ) != string::npos ) { + for (unsigned int jj = 0; jj < ii->second.size(); ++jj) { + vector< double > coords = Field< vector< double > >::get( src_[ii->second[jj]], fieldName.c_str() ); + if ( coords.size() == 11 ) { // For SpineMesh + for ( unsigned int kk = 0; kk < 6; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + } + buffer[jj * 7 + 6] = coords[9]; // head Dia + } else if ( coords.size() == 4 ) { // for EndoMesh + for ( unsigned int kk = 0; kk < 3; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + buffer[jj * 7 + kk+3] = coords[kk]; + } + buffer[jj * 7 + 6] = coords[3]; + } else if ( coords.size() >= 7 ) { // For NeuroMesh + for ( unsigned int kk = 0; kk < 7; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + } + } + } + } else { // Want to check for things like Ca in an elec compt... + for (unsigned int jj = 0; jj < ii->second.size(); ++jj) { + for ( unsigned int kk = 0; kk < 7; ++kk) { + buffer[jj * 7 + kk] = 0.0; + } + } + } + hsize_t dims[2]; + dims[0] = ii->second.size(); + dims[1] = 7; + hid_t memspace = H5Screate_simple(2, dims, NULL); + hid_t dataspace = H5Screate_simple(2, dims, NULL); + hid_t dataset = H5Dcreate2(container, fieldName.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + hid_t filespace = H5Dget_space(dataset); + herr_t status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace, filespace, H5P_DEFAULT, buffer); + if ( status < 0 ) { + cout << "Error: Failed to write coords as static entry\n"; + } + } +} + +void NSDFWriter::writeModelFiles() +{ + for ( const string& fName : modelFileNames_ ) { + // string fPath = MODELFILEPATH + string("/") + fName; + string fPath = MODELFILEPATH; + std::ifstream f( fName ); + auto ss = ostringstream{}; + if ( f.is_open() ) { + ss << f.rdbuf(); + hid_t fGroup = require_group(filehandle_, fPath); + writeScalarAttr(fGroup, fName, ss.str()); + } else { + cout << "Warning: NSDFWriter::writeModelFiles Could not open file '" << fName << "'/n"; + } + } +} void NSDFWriter::writeModelTree() { + if (modelRoot_ == "") + return; vector< string > tokens; ObjId mRoot(modelRoot_); string rootPath = MODELTREEPATH + string("/") + mRoot.element()->getName(); diff --git a/builtins/NSDFWriter.h b/builtins/NSDFWriter.h index 8998ed951a..a3654ff612 100644 --- a/builtins/NSDFWriter.h +++ b/builtins/NSDFWriter.h @@ -80,6 +80,8 @@ class NSDFWriter: public HDF5DataWriter NSDFWriter(); ~NSDFWriter(); virtual void flush(); + void setModelFiles(string value); + string getModelFiles() const; // set the environment specs like title, author, tstart etc. void setEnvironment(string key, string value); // the model tree rooted here is to be copied to NSDF file @@ -95,8 +97,12 @@ class NSDFWriter: public HDF5DataWriter void closeEventData(); virtual void close(); void createUniformMap(); + void createStaticMap(); void createEventMap(); void writeModelTree(); + void writeModelFiles(); + void writeStaticCoords(); + void innerCreateMaps( const char* const mapSrcStr ); // Sort the incoming data lines according to source object/field. void process(const Eref &e, ProcPtr p); void reinit(const Eref &e, ProcPtr p); @@ -116,6 +122,7 @@ class NSDFWriter: public HDF5DataWriter vector < InputVariable > eventInputs_; vector < string > eventSrcFields_; vector < string > eventSrc_; + vector < string > modelFileNames_; map < string, hid_t > eventSrcDataset_; hid_t eventGroup_; // handle for event container hid_t uniformGroup_; // handle for uniform container diff --git a/builtins/NSDFWriter2.cpp b/builtins/NSDFWriter2.cpp new file mode 100644 index 0000000000..fdde494b50 --- /dev/null +++ b/builtins/NSDFWriter2.cpp @@ -0,0 +1,958 @@ +// NSDFWriter2.cpp --- +// +// Filename: NSDFWriter2.cpp +// Description: +// Author: subha +// Maintainer: +// Created: Thu Jun 18 23:16:11 2015 (-0400) +// Version: +// Last-Updated: Sat Jan 29 2022 +// By: bhalla +// Update #: 50 +// URL: +// Keywords: +// Compatibility: +// + +// Commentary: +// + +// Change log: +// Jan 2022: Many changes added +// +// +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; see the file COPYING. If not, write to +// the Free Software Foundation, Inc., 51 Franklin Street, Fifth +// Floor, Boston, MA 02110-1301, USA. +// +// + +// Code: +#ifdef USE_HDF5 + +#include "hdf5.h" +#include "hdf5_hl.h" + +#include +#include +#include +#include +#include "../basecode/header.h" +#include "../utility/utility.h" +#include "../utility/strutil.h" +#include "../shell/Wildcard.h" +#include "../shell/Shell.h" + +#include "HDF5WriterBase.h" +#include "HDF5DataWriter.h" + +#include "NSDFWriter.h" +#include "NSDFWriter2.h" +#include "InputVariable.h" + +extern template herr_t writeScalarAttr(hid_t file_id, string path, string value); +extern template herr_t writeScalarAttr(hid_t file_id, string path, double value); + +const char* const EVENTPATH = "/data/event"; +const char* const UNIFORMPATH = "/data/uniform"; +const char* const STATICPATH = "/data/static"; +const char* const MODELTREEPATH = "/model/modeltree"; +const char* const MODELFILEPATH = "/model/modelfile"; +const char* const MAPUNIFORMSRC = "/map/uniform"; +const char* const MAPSTATICSRC = "/map/static"; +const char* const MAPEVENTSRC = "/map/event"; + +string iso_time(time_t * t); +/* +string iso_time(time_t * t) +{ + struct tm * timeinfo; + if (t == NULL){ + time_t current; + std::time(¤t); + timeinfo = std::gmtime(¤t); + } else { + timeinfo = std::gmtime(t); + } + assert(timeinfo != NULL); + char buf[32]; + strftime(buf, 32, "%FT%T", timeinfo); + return string(buf); +} +*/ + +const Cinfo * NSDFWriter2::initCinfo() +{ + static FieldElementFinfo< NSDFWriter2, InputVariable > eventInputFinfo( + "eventInput", + "Sets up field elements for event inputs", + InputVariable::initCinfo(), + &NSDFWriter2::getEventInput, + &NSDFWriter2::setNumEventInputs, + &NSDFWriter2::getNumEventInputs); + + static ValueFinfo modelRoot( + "modelRoot", + "The moose element tree root to be saved under /model/modeltree. If blank, nothing is saved. Default: root object, '/'", + &NSDFWriter2::setModelRoot, + &NSDFWriter2::getModelRoot); + + static ValueFinfo modelFileNames( + "modelFileNames", + "Comma separated list of model files to save into the NSDF file.", + &NSDFWriter2::setModelFiles, + &NSDFWriter2::getModelFiles); + + static ValueFinfo > blocks( + "blocks", + "Vector of strings to specify data blocks. Format: path.field" + "The path is a wildcard path. It will be split into a single path" + "to a container such as a Neuron or a Mesh, and below this a " + "wildcard path to the actual objects", + &NSDFWriter2::setBlocks, + &NSDFWriter2::getBlocks); + + static DestFinfo process( + "process", + "Handle process calls. Collects data in buffer and if number of steps" + " since last write exceeds flushLimit, writes to file.", + new ProcOpFunc( &NSDFWriter2::process)); + + static DestFinfo reinit( + "reinit", + "Reinitialize the object. If the current file handle is valid, it tries" + " to close that and open the file specified in current filename field.", + new ProcOpFunc( &NSDFWriter2::reinit )); + + static Finfo * processShared[] = { + &process, &reinit + }; + + static SharedFinfo proc( + "proc", + "Shared message to receive process and reinit", + processShared, sizeof( processShared ) / sizeof( Finfo* )); + + static Finfo * finfos[] = { + &eventInputFinfo, // FieldElementFinfo + &modelRoot, // ValueFinfo + &modelFileNames, // ValueFinfo + &blocks, // ValueFinfo + &proc, + }; + + static string doc[] = { + "Name", "NSDFWriter2", + "Author", "Upi Bhalla", + "Description", "NSDF file writer for saving data." + }; + + static Dinfo< NSDFWriter2 > dinfo; + static Cinfo cinfo( + "NSDFWriter2", + HDF5DataWriter::initCinfo(), + finfos, + sizeof(finfos)/sizeof(Finfo*), + &dinfo, + doc, sizeof( doc ) / sizeof( string )); + + return &cinfo; +} + +static const Cinfo * nsdfWriterCinfo = NSDFWriter2::initCinfo(); + +NSDFWriter2::NSDFWriter2(): eventGroup_(-1), uniformGroup_(-1), dataGroup_(-1), modelGroup_(-1), mapGroup_(-1), modelRoot_("/") +{ + ; +} + +NSDFWriter2::~NSDFWriter2() +{ + close(); +} + +void NSDFWriter2::close() +{ + if (filehandle_ < 0){ + return; + } + flush(); + closeUniformData(); + if (uniformGroup_ >= 0){ + H5Gclose(uniformGroup_); + uniformGroup_ = -1; + } + closeEventData(); + if (eventGroup_ >= 0){ + H5Gclose(eventGroup_); + eventGroup_ = -1; + } + if (dataGroup_ >= 0){ + H5Gclose(dataGroup_); + dataGroup_ = -1; + } + HDF5DataWriter::close(); + for ( auto bb = blocks_.begin(); bb != blocks_.end(); ++bb ) { + bb->hasContainer = false; + } +} + +void NSDFWriter2::closeUniformData() +{ + for ( vector< Block >::iterator ii = blocks_.begin(); ii != blocks_.end(); ++ii ) { + if ( ii->dataset >= 0 ) { + H5Dclose( ii->dataset ); + } + /** + ii->hasMsg = false; + ii->hasContainer = false; + ii->objVec.clear(); + ii->objPathList.clear(); + for ( auto jj = ii->data.begin(); jj != ii->data.end(); jj++ ) { + jj.clear(); + } + ii->data.clear(); + */ + } + vars_.clear(); + data_.clear(); + src_.clear(); + func_.clear(); + datasets_.clear(); + +} + +void NSDFWriter2::sortMsgs(const Eref& eref) +{ + const Finfo* tmp = eref.element()->cinfo()->findFinfo("requestOut"); + const SrcFinfo1< vector < double > *>* requestOut = static_cast * > * >(tmp); + vector< ObjId > tgts = eref.element()->getMsgTargets( eref.dataIndex(), requestOut ); + // Make a map from ObjId to index of obj in objVec. + map< ObjId, unsigned int > mapMsgs; + for (unsigned int tgtIdx = 0; tgtIdx < tgts.size(); ++tgtIdx) + mapMsgs[ tgts[ tgtIdx ] ] = tgtIdx; + + mapMsgIdx_.resize( tgts.size() ); + unsigned int consolidatedBlockMsgIdx = 0; + for (unsigned int blockIdx = 0; blockIdx < blocks_.size(); ++blockIdx) { + vector< ObjId >& objVec = blocks_[blockIdx].objVec; + for ( auto obj = objVec.begin(); obj != objVec.end(); ++obj ) { + mapMsgIdx_[consolidatedBlockMsgIdx] = mapMsgs[*obj]; + consolidatedBlockMsgIdx++; + } + } + assert( tgts.size() == consolidatedBlockMsgIdx ); + // make a vector where tgtMsgIdx = mapMsgIdx_[consolidatedBlockMsgIdx] +} + +void NSDFWriter2::buildUniformSources(const Eref& eref) +{ + Shell* shell = reinterpret_cast(Id().eref().data()); + for ( auto bb = blocks_.begin(); bb != blocks_.end(); ++bb ) { + if ( bb->hasMsg ) + continue; + const vector< ObjId >& objVec = bb->objVec; + for( vector< ObjId >::const_iterator obj = objVec.begin(); obj != objVec.end(); ++obj ) { + + ObjId ret = shell->doAddMsg( "single", eref.objId(), "requestOut", *obj, bb->getField ); + if (ret == ObjId() ) { + cout << "Error: NSDFWriter2::buildUniformSources: Failed to build msg from '" << eref.id().path() << "' to '" << bb->containerPath << "/" << bb->relPath << "." << bb->field << endl; + return; + } + } + + bb->hasMsg = true; + } + sortMsgs( eref ); +} + +/** + Handle the datasets for the requested fields (connected to + requestOut). This is is similar to what HDF5DataWriter does. + */ +void NSDFWriter2::openUniformData(const Eref &eref) +{ + buildUniformSources(eref); + htri_t exists; + herr_t status; + if (uniformGroup_ < 0){ + uniformGroup_ = require_group(filehandle_, UNIFORMPATH); + } + for ( auto bb = blocks_.begin(); bb != blocks_.end(); ++bb ) { + if ( bb->hasContainer ) + continue; + // From the documentation: + // https://support.hdfgroup.org/HDF5/doc1.6/UG/09_Groups.html + // "Component link names may be any string of ASCII characters not containing a slash or a dot (/ and ., which are reserved as noted above)." + // So I need to replace path with a string with the slashes + bb->container = require_group(uniformGroup_, bb->nsdfContainerPath); + bb->relPathContainer = require_group(bb->container,bb->nsdfRelPath); + hid_t dataset = createDataset2D(bb->relPathContainer, bb->field.c_str(), bb->data.size()); + bb->dataset = dataset; + writeScalarAttr(dataset, "field", bb->field); + H5Gclose(bb->container); + H5Gclose(bb->relPathContainer); + bb->hasContainer = true; + } +} + +/** + create the DS for uniform data. + */ +void NSDFWriter2::createUniformMap() +{ + innerCreateMaps( MAPUNIFORMSRC ); +} + +/** + create the DS for static data. + */ +void NSDFWriter2::createStaticMap() +{ + innerCreateMaps( MAPSTATICSRC ); +} + + +/** + Generic call for create the DS for static/uniform data + */ +void NSDFWriter2::innerCreateMaps( const char* const mapSrcStr ) +{ + // Create the container for all the DS + htri_t exists; + herr_t status; + hid_t uniformMapContainer = require_group(filehandle_, mapSrcStr ); + // Create the DS themselves + for( auto bit = blocks_.begin(); bit != blocks_.end(); bit++ ) { + //for (map< string, vector < unsigned int > >::iterator ii = classFieldToSrcIndex_.begin(); ii != classFieldToSrcIndex_.end(); ++ii){ + /* + vector < string > pathTokens; + moose::tokenize(ii->first, "/", pathTokens); + string className = pathTokens[0]; + string fieldName = pathTokens[1]; + */ + string className = bit->nsdfContainerPath + "/" + bit->nsdfRelPath; + string fieldName = bit->field; + if (mapSrcStr == MAPSTATICSRC ) //Hack. for now only static field is coords + fieldName = "coords"; + hid_t container = require_group(uniformMapContainer, className); + char ** sources = (char **)calloc(bit->objVec.size(), sizeof(char*)); + for (unsigned int jj = 0; jj < bit->objVec.size(); ++jj){ + sources[jj] = (char*)calloc(bit->objVec[jj].path().length()+1, sizeof(char)); + strcpy(sources[jj],bit->objVec[jj].path().c_str()); + } + hid_t ds = createStringDataset(container, fieldName, (hsize_t)bit->objVec.size(), (hsize_t)bit->objVec.size()); + hid_t memtype = H5Tcopy(H5T_C_S1); + status = H5Tset_size(memtype, H5T_VARIABLE); + assert(status >= 0); + status = H5Dwrite(ds, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, sources); + assert(status >= 0); + for (unsigned int jj = 0; jj < bit->objVec.size(); ++jj){ + free(sources[jj]); + } + free(sources); + /* + status = H5DSset_scale(ds, "source"); + status = H5DSattach_scale(classFieldToUniform_[ii->first], ds, 0); + status = H5DSset_label(classFieldToUniform_[ii->first], 0, "source"); + */ + status = H5Dclose(ds); + status = H5Tclose(memtype); + } +} + +void NSDFWriter2::closeEventData() +{ + for (unsigned int ii = 0; ii < eventDatasets_.size(); ++ii){ + if (eventDatasets_[ii] >= 0){ + H5Dclose(eventDatasets_[ii]); + } + } + events_.clear(); + eventInputs_.clear(); + eventDatasets_.clear(); + eventSrc_.clear(); + eventSrcFields_.clear(); +} + +/** + Populates the vector of event data buffers (vectors), vector of + event source objects, vector of event source fields and the vector + of event datasets by querying the messages on InputVariables. + */ +void NSDFWriter2::openEventData(const Eref &eref) +{ + if (filehandle_ <= 0){ + return; + } + for (unsigned int ii = 0; ii < eventInputs_.size(); ++ii){ + stringstream path; + path << eref.objId().path() << "/" << "eventInput[" << ii << "]"; + ObjId inputObj = ObjId(path.str()); + Element * el = inputObj.element(); + const DestFinfo * dest = static_cast(el->cinfo()->findFinfo("input")); + vector < ObjId > src; + vector < string > srcFields; + el->getMsgSourceAndSender(dest->getFid(), src, srcFields); + if (src.size() > 1){ + cerr << "NSDFWriter2::openEventData - only one source can be connected to an eventInput" < >::iterator ii = classFieldToEventSrc_.begin(); + ii != classFieldToEventSrc_.end(); + ++ii){ + vector < string > pathTokens; + moose::tokenize(ii->first, "/", pathTokens); + string className = pathTokens[0]; + string fieldName = pathTokens[1]; + hid_t classGroup = require_group(eventMapContainer, className); + hid_t strtype = H5Tcopy(H5T_C_S1); + status = H5Tset_size(strtype, H5T_VARIABLE); + // create file space + hid_t ftype = H5Tcreate(H5T_COMPOUND, sizeof(hvl_t) +sizeof(hobj_ref_t)); + status = H5Tinsert(ftype, "source", 0, strtype); + status = H5Tinsert(ftype, "data", sizeof(hvl_t), H5T_STD_REF_OBJ); + hsize_t dims[1] = {ii->second.size()}; + hid_t space = H5Screate_simple(1, dims, NULL); + // The dataset for mapping is named after the field + hid_t ds = H5Dcreate2(classGroup, fieldName.c_str(), ftype, space, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + status = H5Sclose(space); + map_type * buf = (map_type*)calloc(ii->second.size(), sizeof(map_type)); + // Populate the buffer entries with source uid and data + // reference + for (unsigned int jj = 0; jj < ii->second.size(); ++jj){ + buf->source = ii->second[jj].c_str(); + char * dsname = (char*)calloc(256, sizeof(char)); + hsize_t size = H5Iget_name(classFieldToEvent_[ii->first][jj], dsname, 255); + if (size > 255){ + free(dsname); + dsname = (char*)calloc(size, sizeof(char)); + size = H5Iget_name(classFieldToEvent_[ii->first][jj], dsname, 255); + } + status = H5Rcreate(&(buf->data), filehandle_, dsname, H5R_OBJECT, -1); + free(dsname); + assert(status >= 0); + } + // create memory space + hid_t memtype = H5Tcreate(H5T_COMPOUND, sizeof(map_type)); + status = H5Tinsert(memtype, "source", + HOFFSET(map_type, source), strtype); + status = H5Tinsert(memtype, "data", + HOFFSET(map_type, data), H5T_STD_REF_OBJ); + status = H5Dwrite(ds, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf); + free(buf); + status = H5Tclose(strtype); + status = H5Tclose(ftype); + status = H5Tclose(memtype); + status = H5Dclose(ds); + } +} + +/** + Create or retrieve a dataset for an event input. The dataset path + will be /data/event/{class}/{srcFinfo}/{id}_{dataIndex}_{fieldIndex}. + + path : {source_object_id}.{source_field_name} + + TODO: check the returned hid_t and show appropriate error messages. +*/ +hid_t NSDFWriter2::getEventDataset(string srcPath, string srcField) +{ + string eventSrcPath = srcPath + string("/") + srcField; + map< string, hid_t >::iterator it = eventSrcDataset_.find(eventSrcPath); + if (it != eventSrcDataset_.end()){ + return it->second; + } + ObjId source(srcPath); + herr_t status; + htri_t exists = -1; + string className = Field::get(source, "className"); + string path = EVENTPATH + string("/") + className + string("/") + srcField; + hid_t container = require_group(filehandle_, path); + stringstream dsetname; + dsetname << source.id.value() <<"_" << source.dataIndex << "_" << source.fieldIndex; + hid_t dataset = createDoubleDataset(container, dsetname.str().c_str()); + classFieldToEvent_[className + "/" + srcField].push_back(dataset); + classFieldToEventSrc_[className + "/" + srcField].push_back(srcPath); + status = writeScalarAttr(dataset, "source", srcPath); + assert(status >= 0); + status = writeScalarAttr(dataset, "field", srcField); + assert(status >= 0); + eventSrcDataset_[eventSrcPath] = dataset; + return dataset; +} + +void NSDFWriter2::flush() +{ + // We need to update the tend on each write since we do not know + // when the simulation is getting over and when it is just paused. + writeScalarAttr(filehandle_, "tend", iso_time(NULL)); + + // append all uniform data + for ( vector< Block >::iterator bit = blocks_.begin(); (steps_ > 0) && (bit != blocks_.end()); bit++ ) { + assert( steps_ == bit->data[0].size() ); + double* buffer = (double*)calloc(bit->data.size() * steps_, sizeof(double)); + for (unsigned int ii = 0; ii < bit->data.size(); ++ii){ + for (unsigned int jj = 0; jj < steps_; ++jj){ + buffer[ii * steps_ + jj] = bit->data[ii][jj]; + } + bit->data[ii].clear(); + } + hid_t filespace = H5Dget_space(bit->dataset); + if (filespace < 0){ + cout << "Error: NSDFWriter2::flush(): Failed to open filespace\n"; + break; + } + hsize_t dims[2]; + hsize_t maxdims[2]; + // retrieve current datset dimensions + herr_t status = H5Sget_simple_extent_dims(filespace, dims, maxdims); + hsize_t newdims[] = {dims[0], dims[1] + steps_}; // new column count + status = H5Dset_extent(bit->dataset, newdims); // extend dataset to new column count + if ( status < 0 ) { + cout << "Error: NSDFWriter2::flush(): Fail to extend dataset\n"; + break; + } + H5Sclose(filespace); + filespace = H5Dget_space(bit->dataset); // get the updated filespace + hsize_t start[2] = {0, dims[1]}; + dims[1] = steps_; // change dims for memspace & hyperslab + hid_t memspace = H5Screate_simple(2, dims, NULL); + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, start, NULL, dims, NULL); + status = H5Dwrite(bit->dataset, H5T_NATIVE_DOUBLE, memspace, filespace, H5P_DEFAULT, buffer); + if ( status < 0 ) { + cout << "Error: NSDFWriter2::flush(): Failed to write data\n"; + break; + } + H5Sclose(memspace); + H5Sclose(filespace); + free(buffer); + } + steps_ = 0; + + // append all event data + for (unsigned int ii = 0; ii < eventSrc_.size(); ++ii){ + appendToDataset(getEventDataset(eventSrc_[ii], eventSrcFields_[ii]), + events_[ii]); + events_[ii].clear(); + } + // flush HDF5 nodes. + HDF5DataWriter::flush(); +} + +void NSDFWriter2::reinit(const Eref& eref, const ProcPtr proc) +{ + // write environment + // write model + // write map + if (filehandle_ >0){ + close(); + } + // TODO: what to do when reinit is called? Close the existing file + // and open a new one in append mode? Or keep adding to the + // current file? + if (filename_.empty()){ + filename_ = "moose_data.nsdf.h5"; + } + openFile(); + writeScalarAttr(filehandle_, "created", iso_time(0)); + writeScalarAttr(filehandle_, "tstart", iso_time(0)); + writeScalarAttr(filehandle_, "nsdf_version", "1.0"); + openUniformData(eref); + for (vector< Block >::iterator bi = blocks_.begin(); bi != blocks_.end(); ++bi) { + writeScalarAttr< double >(bi->dataset, "tstart", 0.0); + // dt is same for all requested data - that of the NSDFWriter2 + writeScalarAttr< double >(bi->dataset, "dt", proc->dt); + } + openEventData(eref); + writeModelFiles(); + writeModelTree(); + createUniformMap(); + createStaticMap(); + createEventMap(); + writeStaticCoords(); + steps_ = 0; +} + +void NSDFWriter2::process(const Eref& eref, ProcPtr proc) +{ + if (filehandle_ < 0){ + return; + } + vector < double > uniformData; + const Finfo* tmp = eref.element()->cinfo()->findFinfo("requestOut"); + const SrcFinfo1< vector < double > *>* requestOut = static_cast * > * >(tmp); + requestOut->send(eref, &uniformData); + assert( uniformData.size() == mapMsgIdx_.size() ); + // Note that uniformData is ordered by msg tgt order. We want to store + // data in block_->objVec order. + unsigned int ii = 0; + for (unsigned int blockIdx = 0; blockIdx < blocks_.size(); ++blockIdx) { + vector< vector< double > >& bjd = blocks_[blockIdx].data; + for ( auto jj = bjd.begin(); jj != bjd.end(); ++jj ) { + jj->push_back( uniformData[ mapMsgIdx_[ii] ] ); + ii++; + } + } + ++steps_; + if (steps_ < flushLimit_){ + return; + } + NSDFWriter2::flush(); + } + +NSDFWriter2& NSDFWriter2::operator=( const NSDFWriter2& other) +{ + eventInputs_ = other.eventInputs_; + for ( vector< InputVariable >::iterator i = eventInputs_.begin(); i != eventInputs_.end(); ++i ) { + i->setOwner( this ); + } + for (unsigned int ii = 0; ii < getNumEventInputs(); ++ii){ + events_[ii].clear(); + } + return *this; +} + +void NSDFWriter2::setNumEventInputs(unsigned int num) +{ + unsigned int prevSize = eventInputs_.size(); + eventInputs_.resize(num); + for (unsigned int ii = prevSize; ii < num; ++ii){ + eventInputs_[ii].setOwner(this); + } +} + +unsigned int NSDFWriter2::getNumEventInputs() const +{ + return eventInputs_.size(); +} + +void NSDFWriter2::setEnvironment(string key, string value) +{ + env_[key] = value; +} + + +void NSDFWriter2::setInput(unsigned int index, double value) +{ + events_[index].push_back(value); +} + +InputVariable* NSDFWriter2::getEventInput(unsigned int index) +{ + static InputVariable dummy; + if (index < eventInputs_.size()){ + return &eventInputs_[index]; + } + cout << "Warning: NSDFWriter2::getEventInput: index: " << index << + " is out of range: " << eventInputs_.size() << endl; + return &dummy; +} + +void NSDFWriter2::setModelRoot(string value) +{ + modelRoot_ = value; +} + +string NSDFWriter2::getModelRoot() const +{ + return modelRoot_; +} + +void NSDFWriter2::setModelFiles(string value) +{ + modelFileNames_.clear(); + moose::tokenize( value, ", ", modelFileNames_); +} + +string NSDFWriter2::getModelFiles() const +{ + string ret = ""; + string spacer = ""; + for( auto s = modelFileNames_.begin(); s!= modelFileNames_.end(); ++s) { + ret += spacer + *s; + spacer = ","; + } + return ret; +} + +bool parseBlockString( const string& val, Block& block ) +{ + string s = val; + auto ff = s.find_last_of("."); + if (ff == string::npos ) + return false; + block.hasMsg = false; + block.hasContainer = false; + block.dataset = -1; + block.field = s.substr( ff + 1 ); + string temp = block.field; + temp[0] = toupper( temp[0] ); + block.getField = "get" + temp; + s = s.substr( 0, ff ); + vector< string > svec; + moose::tokenize(s, "/", svec); + string path = ""; + unsigned int containerIdx = 0; + block.containerPath = ""; + block.nsdfContainerPath = ""; + string pct = ""; + for ( unsigned int ii = 0; ii < svec.size(); ii++ ) { + path += "/" + svec[ii]; + Id id( path ); + if ( id != Id() ) { + if ( id.element()->cinfo()->isA( "Neuron" ) || + id.element()->cinfo()->isA( "ChemCompt" ) ) + { + containerIdx = ii; + block.containerPath = path; + block.nsdfContainerPath += pct + svec[ii]; + } + pct = "%"; + } else { + cout << "Error: NSDFWriter2: parseBlockString: No object found on '" << path << "'. Ignoring block.\n"; + } + } + if( block.containerPath == "" ) + return false; + block.relPath = ""; + block.nsdfRelPath = ""; + string slash = ""; + pct = ""; + for ( auto jj = containerIdx+1; jj < svec.size(); ++jj) { + block.relPath += slash + svec[jj]; + slash = "/"; + block.nsdfRelPath += pct + svec[jj]; + pct = "%"; + } + if( block.relPath == "" ) + return false; + string objWildcardPath = block.containerPath + '/' + block.relPath; + block.objVec.clear(); + simpleWildcardFind( objWildcardPath, block.objVec ); + if ( block.objVec.size() == 0 ) { + cout << "Error: NSDFWriter2:parseBlockString: No objects found on path '" << objWildcardPath << "'\n"; + return false; + } + block.className = block.objVec[0].element()->cinfo()->name(); + for ( auto obj : block.objVec ) { + // Nasty workaround for different ways of handling CaConcs in Hsolve + if (obj.element()->cinfo()->name().find("CaConc") != string::npos && + block.className.find( "CaConc" ) != string::npos ) + continue; + if (obj.element()->cinfo()->name() != block.className ) { + cout << "Error: NSDFWriter2:parseBlockString: different classes found on path '" << objWildcardPath << "': '" << block.className << "' vs. '" << obj.element()->cinfo()->name() << "'\n"; + return false; + } + } + + block.data.resize( block.objVec.size() ); + for( auto i = block.data.begin(); i != block.data.end(); i++ ) + i->clear(); + return true; +} + +void NSDFWriter2::setBlocks(vector< string > value) +{ + if ( value.size() == 0 ) + blocks_.clear(); + blocks_.resize( value.size() ); + for ( unsigned int i = 0; i < value.size(); ++i ) { + if ( !parseBlockString( value[i], blocks_[i] ) ) { + cout << "Error: NSDFWriter2::setBlocks: block[" << i << "] = '" + << value[i] << "' failed\n"; + return; + } + } + blockStrVec_ = value; +} + +vector< string > NSDFWriter2::getBlocks() const +{ + return blockStrVec_; +} + +//////////////////////////////////////////////////////////////////////// + +ObjId findParentElecCompt( ObjId obj ) +{ + for (ObjId pa = Field< ObjId >::get( obj, "parent" ); pa != ObjId(); + pa = Field< ObjId >::get( pa, "parent" ) ) { + if ( pa.element()->cinfo()->isA( "CompartmentBase" ) ) + return pa; + } + return ObjId(); +} + +void NSDFWriter2::writeStaticCoords() +{ + hid_t staticObjContainer = require_group(filehandle_, STATICPATH ); + for( auto bit = blocks_.begin(); bit != blocks_.end(); bit++ ) { + string coordContainer = bit->nsdfContainerPath + "/" + bit->nsdfRelPath; + string fieldName = "coords"; // pathTokens[1] is not relevant. + hid_t container = require_group(staticObjContainer, coordContainer); + double * buffer = + (double*)calloc(bit->data.size() * 7, sizeof(double)); + if ( bit->className.find( "Pool" ) != string::npos || + bit->className.find( "Compartment" ) != string::npos ) { + for (unsigned int jj = 0; jj < bit->data.size(); ++jj) { + ObjId obj = bit->objVec[jj]; + vector< double > coords = Field< vector< double > >::get( obj, fieldName ); + if ( coords.size() == 11 ) { // For SpineMesh + for ( unsigned int kk = 0; kk < 6; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + } + buffer[jj * 7 + 6] = coords[9]; // head Dia + } else if ( coords.size() == 4 ) { // for EndoMesh + for ( unsigned int kk = 0; kk < 3; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + buffer[jj * 7 + kk+3] = coords[kk]; + } + buffer[jj * 7 + 6] = coords[3]; + } else if ( coords.size() >= 7 ) { // For NeuroMesh + for ( unsigned int kk = 0; kk < 7; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + } + } + } + } else { // Check for things like Ca or chans in an elec compt + for (unsigned int jj = 0; jj < bit->objVec.size(); ++jj) { + ObjId pa = findParentElecCompt( bit->objVec[jj] ); + vector< double > coords( 7, 0.0 ); + if (pa != ObjId()) { + coords = Field< vector< double > >::get(pa, fieldName); + } + for ( unsigned int kk = 0; kk < 7; ++kk) { + buffer[jj * 7 + kk] = coords[kk]; + } + } + } + hsize_t dims[2]; + dims[0] = bit->data.size(); + dims[1] = 7; + hid_t memspace = H5Screate_simple(2, dims, NULL); + hid_t dataspace = H5Screate_simple(2, dims, NULL); + hid_t dataset = H5Dcreate2(container, fieldName.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + hid_t filespace = H5Dget_space(dataset); + herr_t status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace, filespace, H5P_DEFAULT, buffer); + if ( status < 0 ) { + cout << "Error: Failed to write coords as static entry\n"; + } + } +} + +void NSDFWriter2::writeModelFiles() +{ + // These can be large, exceed 64K limit of attributes. So write as + // datasets, not attributes. + for ( const string& fName : modelFileNames_ ) { + // string fPath = MODELFILEPATH + string("/") + fName; + string fPath = MODELFILEPATH; + std::ifstream f( fName ); + auto ss = ostringstream{}; + if ( f.is_open() ) { + ss << f.rdbuf(); + string fstr = ss.str(); + char* filebuf = (char*)calloc( ss.str().length()+1, sizeof(char) ); + char** sources = (char**) calloc( 1, sizeof(char* ) ); + sources[0] = filebuf; + strcpy( filebuf, fstr.c_str() ); + hid_t fGroup = require_group(filehandle_, fPath); + hid_t ds = createStringDataset(fGroup, fName, (hsize_t)1, (hsize_t)1 ); + hid_t memtype = H5Tcopy(H5T_C_S1); + int status = H5Tset_size(memtype, H5T_VARIABLE ); + assert(status >= 0); + // status = H5Tclose(memtype); + status = H5Dwrite(ds, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, sources ); + assert(status >= 0); + free( filebuf ); + free( sources ); + } else { + cout << "Warning: NSDFWriter2::writeModelFiles Could not open file '" << fName << "'/n"; + } + } +} + +void NSDFWriter2::writeModelTree() +{ + if (modelRoot_ == "") + return; + vector< string > tokens; + ObjId mRoot(modelRoot_); + string rootPath = MODELTREEPATH + string("/") + mRoot.element()->getName(); + hid_t rootGroup = require_group(filehandle_, rootPath); + hid_t tmp; + htri_t exists; + herr_t status; + deque nodeQueue; + deque h5nodeQueue; + nodeQueue.push_back(mRoot); + h5nodeQueue.push_back(rootGroup); + // TODO: need to clarify what happens with array elements. We can + // have one node per vec and set a count field for the number of + // elements + while (nodeQueue.size() > 0){ + ObjId node = nodeQueue.front(); + nodeQueue.pop_front(); + hid_t prev = h5nodeQueue.front();; + h5nodeQueue.pop_front(); + vector < Id > children; + Neutral::children(node.eref(), children); + for ( unsigned int ii = 0; ii < children.size(); ++ii){ + string name = children[ii].element()->getName(); + // skip the system elements + if (children[ii].path() == "/Msgs" + || children[ii].path() == "/clock" + || children[ii].path() == "/classes" + || children[ii].path() == "/postmaster"){ + continue; + } + exists = H5Lexists(prev, name.c_str(), H5P_DEFAULT); + if (exists > 0){ + tmp = H5Gopen2(prev, name.c_str(), H5P_DEFAULT); + } else { + tmp = H5Gcreate2(prev, name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + } + writeScalarAttr< string >(tmp, "uid", children[ii].path()); + nodeQueue.push_back(children[ii]); + h5nodeQueue.push_back(tmp); + } + status = H5Gclose(prev); + } +} +#endif // USE_HDF5 + +// +// NSDFWriter2.cpp ends here diff --git a/builtins/NSDFWriter2.h b/builtins/NSDFWriter2.h new file mode 100644 index 0000000000..d5b3d0e8b6 --- /dev/null +++ b/builtins/NSDFWriter2.h @@ -0,0 +1,214 @@ +/* NSDFWriter2.h --- + * + * Filename: NSDFWriter2.h + * Description: + * Author: bhalla + * Maintainer: + * Created: Sun 13 Feb 2022 + * Version: + * Last-Updated: Sun 13 Feb 2022 (-0500) + * By: bhalla + * Update #: 2 + * URL: + * Keywords: + * Compatibility: + * + */ + +/* Commentary: + * + * + * + */ + +/* Change log: + * + * + */ + +/* This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation; either version 3, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, Fifth + * Floor, Boston, MA 02110-1301, USA. + */ + +/* Code: */ + +#ifdef USE_HDF5 +#ifndef _NSDFWRITER2_H +#define _NSDFWRITER2_H + +#include "HDF5DataWriter.h" + +class InputVariable; + +/** + compound data type for storing source->data mapping for event data +typedef struct { + const char * source; + hobj_ref_t data; +} map_type; +*/ + +typedef struct { + string containerPath; + // For elec models this is the path of the Neuron object. + // For chem models this is the path of the Mesh object. + // This is used to provide a single set of coords. + string relPath; // Relative to container + // for elec models it is wildcard path of all objects wrt container + // for chem models it is rel path to pool[]. Braces for indices. + string nsdfContainerPath; // NSDF makes a group of the container. + // Sits on /data/uniform. + // Converts / to % + string nsdfRelPath; // nsdf makes subgroups identifed by nsdfRelPath, + // with respect to thei container. + // Sits on /data/uniform/container + // It is the relpath with the / replaced by % + vector< string > objPathList; // relative paths of all objects here. + // For elec compts, it is the compt names + // For chem obj it is just the indices. + string field; // Regular MOOSE value field. + string getField; // name of call to get the field. + string className; // All obj in a block should be of same class. + vector< vector< double > > data; // data_[objIdx][timeStep] + vector< ObjId > objVec; + hid_t container; // reference to container + hid_t relPathContainer; // reference to objects on nsdfRelPath + // hid_t filespace; + hid_t dataset; + bool hasMsg; + bool hasContainer; +} Block; + +/** + NSDFWriter dumps data in NSDF file format. + + - As of June 2015, MOOSE uses fixed time steps for updating field + values. So we support only Uniform data. + + - This class has one SrcFinfo to request data. Multiple get{Field} + DestFinfos can be connected to this and the NSDFWriter class will put + together the ones with same name in one 2D dataset. + + - One DestFinfo where SrcFinfos sending out event times can be + connected. These will go under Event data. + + */ +class NSDFWriter2: public HDF5DataWriter +{ + public: + NSDFWriter2(); + ~NSDFWriter2(); + virtual void flush(); + void setModelFiles(string value); + string getModelFiles() const; + // set the environment specs like title, author, tstart etc. + void setEnvironment(string key, string value); + // the model tree rooted here is to be copied to NSDF file + void setModelRoot(string root); + string getModelRoot() const; + void setBlocks(vector< string > blocks); + vector< string > getBlocks() const; + InputVariable *getEventInput(unsigned int index); + void setNumEventInputs(unsigned int num); + unsigned int getNumEventInputs() const; + void setInput(unsigned int index, double value); + void openUniformData(const Eref &eref); + void closeUniformData(); + void openEventData(const Eref &eref); + void closeEventData(); + virtual void close(); + void createUniformMap(); + void createStaticMap(); + void createEventMap(); + void writeModelTree(); + void writeModelFiles(); + void writeStaticCoords(); + void innerCreateMaps( const char* const mapSrcStr ); + // Sort the incoming data lines according to source object/field. + void process(const Eref &e, ProcPtr p); + void reinit(const Eref &e, ProcPtr p); + NSDFWriter2& operator=(const NSDFWriter2& other); + + static const Cinfo *initCinfo(); + + protected: + hid_t getEventDataset(string srcPath, string srcField); + // void sortOutUniformSources(const Eref& eref); + void buildUniformSources(const Eref& eref); + void sortMsgs(const Eref& eref); + /* hid_t getUniformDataset(string srcPath, string srcField); */ + map env_; // environment attributes + vector < hid_t > eventDatasets_; + // event times data_ and datasets_ inherited from HDF5DataWriter + // are still attached to requestOut message + vector < vector < double > > events_; + vector < InputVariable > eventInputs_; + vector < string > eventSrcFields_; + vector < string > eventSrc_; + vector < string > modelFileNames_; + map < string, hid_t > eventSrcDataset_; + hid_t eventGroup_; // handle for event container + hid_t uniformGroup_; // handle for uniform container + hid_t dataGroup_; // handle for data container. + hid_t modelGroup_; // handle for model container + hid_t mapGroup_; // handle for map container + vector< string > blockStrVec_; + vector< Block > blocks_; + vector< unsigned int > mapMsgIdx_; // Look up tgt idx from consolidated block idx. + + map< string, vector< hid_t > > classFieldToEvent_; + map< string, vector< string > > classFieldToEventSrc_; + map< string, hid_t > classFieldToUniform_; + // maps a path.srcFinfo to the for + // storing uniform data in dataset. + map< string, pair< hid_t, unsigned int > > uniformSrcDatasetIndex_; + /** + Storing uniform data: + + The data coming in return for requestOut gets stored in + dataBuffer in the same sequence as the connections, + irrespective of class and field. + + The messages are to be made in the correct order by the build + command, which takes specified basepath, relpath wildcard, and field. + A separate block is built for each such command. + In each block we have + The full path of the source object i + + For each source object find the class and field and create map + < class.field, vector > that maps the class.field + to a vector of indices in dataBuffer. the index vector contents + should be in the same order as the sequence of objects + + */ + // The last number of rows in each dataset + // (/data/uniform/className/fieldName -> size) + + map< string, vector< unsigned int > > classFieldToSrcIndex_; + // the index of the recorded field in its dataset (objectPath/fieldName -> index) + map< string, unsigned int > objectFieldToIndex_; + // objectPath -> field - just for reuse later - avoid repeating + // extraction of field name from func_. + vector < pair< string, string > > objectField_; + map< string, vector < string > > classFieldToObjectField_; + vector < string > vars_; + string modelRoot_; + +}; +#endif // _NSDFWRITER2_H +#endif // USE_HDF5 + + +/* NSDFWriter2.h ends here */ diff --git a/kinetics/ReadKkit.cpp b/kinetics/ReadKkit.cpp index bb62e5f763..f67cbc35ca 100644 --- a/kinetics/ReadKkit.cpp +++ b/kinetics/ReadKkit.cpp @@ -225,6 +225,9 @@ void setMethod( Shell* s, Id mgr, double simdt, double plotdt, { vector< ObjId > ret; simpleWildcardFind( mgr.path() + "/#[ISA=ChemCompt]", ret ); + if ( ret.size() == 0 ) { // We may have converted /kinetics to a group + simpleWildcardFind( mgr.path() + "/#/#[ISA=ChemCompt]", ret ); + } assert( ret.size() > 0 ); Id compt( mgr.path() + "/kinetics" ); @@ -271,6 +274,49 @@ void setMethod( Shell* s, Id mgr, double simdt, double plotdt, s->doSetClock( 17, plotdt ); // Stats objects s->doSetClock( 18, plotdt ); // Table2 objects. } + + +bool kineticsHasReactions( Id mgr ) { + static const string choices[] = { "Pool", "PoolBase", "Reac", "Enz", "MMEnz", "Neutral" }; + Id kinId = Neutral::child( mgr.eref(), "kinetics" ); + assert( kinId != Id() ); + vector< Id > kids; + Neutral::children( kinId.eref(), kids ); + for( vector< Id >::iterator k = kids.begin(); k != kids.end(); ++k ) { + string childClass = Field< string >::get( *k, "className" ); + auto const last = std::end( choices ); + auto const pos = std::find( std::begin(choices), std::end(choices), childClass ); + // cout << k->path() << " " << childClass << endl; + if ( pos!= last) + return true; + } + return false; +} + +Id ReadKkit::convertKineticsToGroup( Id mgr ) { + Id kinId = Neutral::child( mgr.eref(), "kinetics" ); + kinId.element()->setName( "kinetics_conv_to_group"); + Id newKin = shell_->doCreate( "Neutral", mgr, "kinetics", 1 ); + // Move all child objects onto newKin. + vector< Id > kids; + Neutral::children( kinId.eref(), kids ); + for ( auto k: kids ) + if ( k.element()->getName() == "mesh" ) { + shell_->doDelete( k ); + } else { + shell_->doMove( k, newKin ); + } + auto oldCompt = compartments_; + compartments_.clear(); + for ( auto ii = oldCompt.begin(); ii != oldCompt.end(); ++ii ) { + if ( *ii != kinId ) { + compartments_.push_back( *ii ); + } + } + shell_->doDelete( kinId ); + return newKin; +} + /** * The readcell function implements the old GENESIS cellreader * functionality. Although it is really a parser operation, I @@ -311,6 +357,8 @@ Id ReadKkit::read( assignReacCompartments(); assignEnzCompartments(); assignMMenzCompartments(); + } else if ( !kineticsHasReactions( mgr ) ) { + Id newKinetics = convertKineticsToGroup( mgr ); } convertParametersToConcUnits(); @@ -327,7 +375,6 @@ Id ReadKkit::read( s->doReinit(); return mgr; } - double ReadKkit::childPoolVol( Id gid ) const { vector< ObjId > pools; @@ -370,6 +417,7 @@ int ReadKkit::findCompartmentsFromAnnotation() // SetGet1< double >::set( comptId, "setVolumeNotRates", vol ); // Now tell the mesh to resize. Field< double >::set( comptId, "volume", vol ); + shell_->doDelete( gid ); } } return ret; diff --git a/kinetics/ReadKkit.h b/kinetics/ReadKkit.h index 1aaeb652f3..9f7a4e9b1c 100644 --- a/kinetics/ReadKkit.h +++ b/kinetics/ReadKkit.h @@ -97,6 +97,9 @@ class ReadKkit */ Id findSumTotSrc( const string& src ); + /// If kinetics has no reacs, we convert it into a group + Id convertKineticsToGroup( Id mgr ); + ////////////////////////////////////////////////////////////////// // Special ops in the model definition ////////////////////////////////////////////////////////////////// diff --git a/ksolve/Gsolve.cpp b/ksolve/Gsolve.cpp index 489e06cfae..b58d10755c 100644 --- a/ksolve/Gsolve.cpp +++ b/ksolve/Gsolve.cpp @@ -938,7 +938,7 @@ void Gsolve::setDsolve( Id dsolve ) ////////////////////////////////////////////////////////////// -// Zombie Pool Access functions +// Pool Access functions ////////////////////////////////////////////////////////////// void Gsolve::setN( const Eref& e, double v ) @@ -946,7 +946,7 @@ void Gsolve::setN( const Eref& e, double v ) unsigned int vox = getVoxelIndex( e ); if ( vox != OFFNODE ) { - if ( e.element()->cinfo()->isA( "ZombieBufPool" ) ) + if ( e.element()->cinfo()->isA( "BufPool" ) ) { // Do NOT round it here, it is folded into rate term. pools_[vox].setN( getPoolIndex( e ), v ); diff --git a/python/rdesigneur/moogul.py b/python/rdesigneur/moogul.py index a5bcbb60f0..50d084d42e 100644 --- a/python/rdesigneur/moogul.py +++ b/python/rdesigneur/moogul.py @@ -1,9 +1,5 @@ -# Moogul.py: MOOSE Graphics Using Lines -# This is a fallback graphics interface for displaying neurons using -# regular matplotlib routines. -# Put in because the GL versions like moogli need all sorts of difficult -# libraries and dependencies. -# Copyright (C) Upinder S. Bhalla NCBS 2018 +# Moogul.py: MOOSE Graphics 3D. +# Copyright (C) Upinder S. Bhalla NCBS 2022 # This program is licensed under the GNU Public License version 3. # @@ -20,6 +16,8 @@ bgvector = vp.vector(0.7, 0.8, 0.9) # RGB bgDict = {'default': bgvector, 'black': vp.color.black, 'white': vp.color.white, 'cyan': vp.color.cyan, 'grey': vp.vector( 0.5, 0.5, 0.5 ) } +sleepTimes = [0.0, 0.0005, 0.001, 0.002, 0.003, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1] + def bgLookup( bg ): col = bgDict.get( bg ) if not col: @@ -35,7 +33,7 @@ def __str__( self ): class MooView: ''' The MooView class is a window in which to display one or more - moose cells, using the MooNeuron class.''' + neurons, using the MooNeuron and MooReacSystemclass.''' viewIdx = 0 origScene = None rgb = [] @@ -58,6 +56,7 @@ def __init__( self, swx = 10, swy = 10, hideAxis = True, title = "view", colorma self.colorbar = None self.valMin = 0.0 self.valMmax = 1.0 + self.plotFlag_ = True @staticmethod def replayLoop(): @@ -97,7 +96,8 @@ def toggleReplay( self ): self.replayButton.background = vp.color.white def setSleepTime( self ): - self.sleep = self.sleepSlider.value + idx = int( round( self.sleepSlider.value ) ) + self.sleep = sleepTimes[idx] self.sleepLabel.text = " Frame dt = {:1.3f} sec".format( self.sleep ) def updateAxis( self ): @@ -119,11 +119,12 @@ def updateAxis( self ): self.zAx.axis = vp.vector( z.dot( right ), z.dot( up ), 0.0 ) self.axisLength.text = "{:.2f} um".format( dx * 1e6*self.scene.range * self.colorbar.width / self.scene.width ) - def makeColorbar( self, doOrnaments = True, colorscale = 'jet', bg = 'default' ): - title = None - if doOrnaments: - title = MooView.consolidatedTitle + "\n" + def innerColorbar( self, title, bg ): barWidth = SCALE_SCENE * 1.5 + if ( bgLookup(bg).mag < 1 ): + barTextColor = vp.color.white + else: + barTextColor = vp.color.black self.colorbar = vp.canvas( title = title, width = barWidth, height = self.swy * SCALE_SCENE, background = bgLookup(bg), align = 'left', range = 1, autoscale = False ) #self.colorbar = vp.canvas( title = title, width = barWidth, height = self.swy * SCALE_SCENE, background = vp.color.cyan, align = 'left', range = 1, autoscale = False ) self.colorbar.userzoom = False @@ -135,17 +136,23 @@ def makeColorbar( self, doOrnaments = True, colorscale = 'jet', bg = 'default' ) for idx, rgb in enumerate( self.rgb ): cbox = vp.box( canvas = self.colorbar, pos = vp.vector( 0, height * (idx - 26), 0), width = width, height = height, color = rgb ) barName = self.title.replace( ' ', '\n' ) - self.barName = vp.label( canvas = self.colorbar, align = 'left', pixel_pos = True, pos = vp.vector( 2, (self.swy - 0.32) * SCALE_SCENE, 0), text = barName, height = 15, color = vp.color.black, box = False, opacity = 0 ) - self.barMin = vp.label( canvas = self.colorbar, align = 'center', pixel_pos = True, pos = vp.vector( barWidth/2, self.swy * SCALE_SCENE * 0.22, 0), text = "{:.3f}".format(self.valMin), height = 12, color = vp.color.black, box = False, opacity = 0 ) - self.barMax = vp.label( canvas = self.colorbar, align = 'center', pixel_pos = True, pos = vp.vector( barWidth/2, (self.swy - 1.2) * SCALE_SCENE, 0), text = "{:.3f}".format(self.valMax), height = 12, color = vp.color.black, box = False, opacity = 0 ) + self.barName = vp.label( canvas = self.colorbar, align = 'left', pixel_pos = True, pos = vp.vector( 2, (self.swy - 0.32) * SCALE_SCENE, 0), text = barName, height = 15, color = barTextColor, box = False, opacity = 0 ) + self.barMin = vp.label( canvas = self.colorbar, align = 'center', pixel_pos = True, pos = vp.vector( barWidth/2, self.swy * SCALE_SCENE * 0.22, 0), text = "{:.3f}".format(self.valMin), height = 12, color = barTextColor, box = False, opacity = 0 ) + self.barMax = vp.label( canvas = self.colorbar, align = 'center', pixel_pos = True, pos = vp.vector( barWidth/2, (self.swy - 1.2) * SCALE_SCENE, 0), text = "{:.3f}".format(self.valMax), height = 12, color = barTextColor, box = False, opacity = 0 ) self.xAx = vp.cylinder( canvas = self.colorbar, pos = axOrigin, axis = vp.vector( 0.8, 0, 0 ), radius = 0.04, color = vp.color.red ) self.yAx = vp.cylinder( canvas = self.colorbar, pos = axOrigin, axis = vp.vector( 0, 0.8, 0 ), radius = 0.04, color = vp.color.green ) self.zAx = vp.cylinder( canvas = self.colorbar, pos = axOrigin, axis = vp.vector( 0, 0, 0 ), radius = 0.04, color = vp.color.blue ) - self.axisLength = vp.label( pos = axOrigin + vp.vector(0, 1, 0), text = "1.00 um", color = vp.color.black, box = False ) + self.axisLength = vp.label( pos = axOrigin + vp.vector(0, 1, 0), text = "1.00 um", color = barTextColor, box = False ) + + def makeColorbar( self, doOrnaments = True, colorscale = 'jet', bg = 'default' ): + title = None + if doOrnaments: + title = MooView.consolidatedTitle + "\n" + self.innerColorbar( title, bg ) if doOrnaments: self.timeLabel = vp.wtext( text = "Time = 0.000 sec", pos = self.colorbar.title_anchor ) self.sleepLabel = vp.wtext( text = " Frame dt = 0.005 sec", pos = self.colorbar.title_anchor ) - self.sleepSlider = vp.slider( pos = self.colorbar.title_anchor, length = 200, bind = self.setSleepTime, min = 0.0, max = 0.2, value = self.sleep ) + self.sleepSlider = vp.slider( pos = self.colorbar.title_anchor, length = 200, bind = self.setSleepTime, min = 0, max = len( sleepTimes ) -1, value = min( len( sleepTimes ), 2 ) ) self.replayButton = vp.button( text = "Start Replay", pos = self.colorbar.title_anchor, bind=self.toggleReplay, disabled = True ) self.colorbar.append_to_title("\n") @@ -153,26 +160,33 @@ def pickObj( self ): obj = self.scene.mouse.pick if obj == None: return - elm = self.innerPickObj( obj ) - if elm: - print( elm.path, elm.dataIndex ) + elmPath = self.innerPickObj( obj ) + if elmPath: + self.handlePick( elmPath ) return elif self.viewIdx == 0: for view in MooView.viewList[1:]: if view.colorbar == None: - elm = view.innerPickObj( obj ) - if elm: - print( elm.path, elm.dataIndex ) + elmPath = view.innerPickObj( obj ) + if elmPath: + self.handlePick( elmPath ) return print( "Object {} not found on view {}".format( obj, self.title ) ) def innerPickObj( self, obj ): for dr in self.drawables_: - elm = dr.findDisplayObject( obj ) - if elm: - return elm + elmPath = dr.findDisplayObject( obj ) + if elmPath: + return (elmPath[0], elmPath[1], dr) return None + def handlePick( self, elmPath ): + path, field, drawable = elmPath + if self.plotFlag_: + drawable.plotHistory( path, field, self.graph, self.graphPlot1 ) + else: + print( path, field ) + def makeScene( self, mergeDisplays, bg = 'default' ): @@ -221,10 +235,14 @@ def firstDraw( self, mergeDisplays, rotation=0.0, elev=0.0, azim=0.0, center = [ else: self.doAutoscale() self.updateAxis() + if self.viewIdx == (MooView.viewIdx-1): + self.graph = vp.graph( title = "Graph", xtitle = "Time (s)", ytitle = " Units here", width = 700, fast=False, align = "left" ) + self.graphPlot1 = vp.gcurve( color = vp.color.blue, interval=-1) + #self.graphPlot1.data = [[0,0], [1,1],[2,0],[3,4],[4,0], [5,1]] + #self.graphPlot1.plot( [[0,0], [1,1],[2,0],[3,4],[4,0]] ) + - def updateValues( self ): - simTime = moose.element( '/clock' ).currentTime - #self.timeStr.set_text( "Time= {:.3f}".format( time ) ) + def updateValues( self, simTime ): for i in self.drawables_: i.updateValues( simTime ) if self.doRotation and abs( self.rotation ) < 2.0 * 3.14 / 3.0: @@ -242,12 +260,12 @@ def replaySnapshot( self, idx ): self.updateAxis() def doAutoscale( self ): - if len( self.drawables_[0].activeDia ) == 0: + if self.drawables_[0].dataWrapper_.numObj() == 0: print( "Warning: No values to display in Moogli view ", self.title ) return - cmin = self.drawables_[0].coordMin - cmax = self.drawables_[0].coordMax - diamax = max( self.drawables_[0].activeDia ) + cmin = self.drawables_[0].dataWrapper_.coordMin_ + cmax = self.drawables_[0].dataWrapper_.coordMax_ + diamax = max( self.drawables_[0].dataWrapper_.getCoords()[:,6] ) v0 = vp.vector( cmin[0], cmin[1], cmin[2] ) v1 = vp.vector( cmax[0], cmax[1], cmax[2] ) #self.scene.camera.axis = self.scene.forward * vp.mag(v1 - v0) * 4 @@ -354,34 +372,73 @@ def printMoogulHelp( self ): def list2vec( arg ): return vp.vector( arg[0], arg[1], arg[2] ) +class DataWrapper: + ''' Class for interfacing between moogli and the data source. Currently + implemented for MOOSE and for nsdf reader. + ''' + def __init__( self, field ): + self.coordMin_ = np.zeros( 3 ) + self.coordMax_ = np.ones( 3 ) + self.field_ = field + self.objList_ = [] + + def getValues( self ): + return np.zeros( 1 ) + + def numObj( self ): + return len( self.objList_ ) + + def getCoords( self ): + return np.array( [] ) + + def getMinMax( self ): + nmin = np.amin(self.coords_, axis = 0) + self.coordMin_ = np.amin( np.array( [nmin[0:3], nmin[3:6]] ), axis = 0 ) + nmax = np.amax(self.coords_, axis = 0) + self.coordMax_ = np.amax( np.array( [nmax[0:3], nmax[3:6]] ), axis = 0 ) + def objPathFromIndex( self, idx ): + if idx < len( self.objList_ ): + return self.objList_[idx].path + return None + + def advance( self, simTime ): + # Checks that the simTime has crossed upcomingTime + return True # used for multi timestep cases. + + def getHistory( self, path, field ): + # stub function. Derived classes fill it in and return useful values + return [0, 1, 2, 3], [ 1, 4, 9, 16] + class MooDrawable: ''' Base class for drawing things''' def __init__( self, - fieldInfo, field, relativeObj, + dataWrapper, colormap, - lenScale, diaScale, autoscale, + lenScale, + diaScale, + fieldScale, + autoscale, valMin, valMax ): - self.field = field - self.relativeObj = relativeObj + self.dataWrapper_ = dataWrapper self.lenScale = lenScale self.diaScale = diaScale + self.fieldScale = fieldScale self.colormap = colormap self.autoscale = autoscale self.valMin = valMin self.valMax = valMax - self.fieldInfo = fieldInfo - self.fieldScale = fieldInfo[2] self.segments = [] self.snapshot = [] - self.coordMin = np.zeros( 3 ) - self.coordMax = np.zeros( 3 ) #cmap = plt.get_cmap( self.colormap, lut = NUM_CMAP ) #self.rgb = [ list2vec(cmap(i)[0:3]) for i in range( NUM_CMAP ) ] def updateValues( self, simTime ): - ''' Obtains values from the associated cell''' - self.val = np.array([moose.getField(i, self.field) for i in self.activeObjs]) * self.fieldScale + if self.dataWrapper_.advance( simTime ): + self.val = self.dataWrapper_.getValues() * self.fieldScale + else: + return + if self.autoscale: valMin = min( self.val ) valMax = max( self.val ) @@ -391,7 +448,10 @@ def updateValues( self, simTime ): scaleVal = NUM_CMAP * (self.val - valMin) / (valMax - valMin) #indices = scaleVal.ndarray.astype( int ) indices = np.maximum( np.minimum( scaleVal, NUM_CMAP-0.5), 0.0).astype(int) + + # Have to figure how this will work with multiple update rates. self.snapshot.append( [simTime, indices] ) + self.displayValues( indices ) def displayValues( self, indices ): @@ -406,18 +466,15 @@ def replaySnapshot( self, idx ): return self.snapshot[idx][0] # return frame time def updateDiameter( self ): - for s, w in zip( self.segments, self.activeDia ): + dia = self.dataWrapper_.getCoords()[:,6] + for s, w in zip( self.segments, dia ): s.radius = self.diaScale * w / 2.0 def cylinderDraw( self, _scene ): - for idx, coord in enumerate( self.activeCoords ): - v0 = list2vec( coord[0] ) - v1 = list2vec( coord[1] ) - self.coordMin = np.minimum( self.coordMin, coord[0][0:3] ) - self.coordMin = np.minimum( self.coordMin, coord[1][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[0][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[1][0:3] ) - radius = self.diaScale * self.activeDia[idx] / 2.0 + for idx, coord in enumerate( self.dataWrapper_.getCoords() ): + v0 = list2vec( coord[0:3] ) + v1 = list2vec( coord[3:6] ) + radius = self.diaScale * coord[6] / 2.0 opacity = self.opacity[idx] rod = vp.cylinder( canvas = _scene, pos = v0, axis = v1 - v0, radius = radius, opacity = opacity ) self.segments.append( rod ) @@ -425,62 +482,56 @@ def cylinderDraw( self, _scene ): def findDisplayObject( self, obj ): try: idx = self.segments.index( obj ) + return self.dataWrapper_.objPathFromIndex( idx ), self.dataWrapper_.field_ except ValueError: return None - if idx >= len( self.activeObjs ): - return None - return self.activeObjs[idx] + + def plotHistory( self, path, field, graph, plot ): + t, v = self.dataWrapper_.getHistory( path, field ) + if len( t ) == 0: + print( "No data history for '", path, ".", field ) + return + #self.graph = vp.graph( title = path + "." + field, xtitle = "Time (s)", ytitle = field + " Units here", width = 800, fast=False, pos=self.colorbar.caption_anchor ) + graph.title = path + "." + field + dat = [[x,y] for x, y in zip( t, v ) ] + plot.data = dat + #print (dat) + #print( "IN plotHistory, ", len( dat), len( v ) ) + #plot.data = [[x,y] for x, y in zip( t, v ) ] + #plot.data = [[x,sin(x)] for x in range( 0.0, 10.0, 0.1 ) ] + ''' + fig = plt.figure( 1 ) + plt.ion() + plt.title( path + "." + field ) + plt.xlabel( "Time (s)" ) + plt.ylabel( field + " um, units?" ) + plt.plot( t, v ) + plt.show( block = False ) + fig.canvas.draw() + ''' + ##################################################################### class MooNeuron( MooDrawable ): ''' Draws collection of line segments of defined dia and color''' - def __init__( self, - neuronId, - fieldInfo, + def __init__( self, + dataWrapper, field = 'Vm', - relativeObj = '.', colormap = 'jet', - lenScale = 1.0, diaScale = 1.0, autoscale = False, + lenScale = 1.0, diaScale = 1.0, fieldScale = 1.0, + autoscale = False, valMin = -0.1, valMax = 0.05, ): #self.isFieldOnCompt = #field in ( 'Vm', 'Im', 'Rm', 'Cm', 'Ra', 'inject', 'diameter' ) - MooDrawable.__init__( self, fieldInfo, field = field, - relativeObj = relativeObj, + MooDrawable.__init__( self, dataWrapper, colormap = colormap, lenScale = lenScale, - diaScale = diaScale, autoscale = autoscale, + diaScale = diaScale, fieldScale = fieldScale, + autoscale = autoscale, valMin = valMin, valMax = valMax ) - self.neuronId = neuronId - self.updateCoords() - - def updateCoords( self ): - ''' Obtains coords from the associated cell''' - self.compts_ = moose.wildcardFind( self.neuronId.path + "/#[ISA=CompartmentBase]" ) - coords = np.array([[[i.x0,i.y0,i.z0],[i.x,i.y,i.z]] - for i in self.compts_]) - dia = np.array([i.diameter for i in self.compts_]) - if self.relativeObj == '.': - self.activeCoords = coords - self.activeDia = dia - self.activeObjs = self.compts_ - else: - self.activeObjs = [] - self.activeCoords = [] - self.activeDia = [] - for i,j,k in zip( self.compts_, coords, dia ): - if moose.exists( i.path + '/' + self.relativeObj ): - elm = moose.element( i.path + '/' + self.relativeObj ) - self.activeObjs.append( elm ) - self.activeCoords.append( j ) - self.activeDia.append( k ) - - self.activeCoords = np.array( self.activeCoords ) * self.lenScale - self.opacity = np.ones( len( self.activeDia ) ) * 0.5 - super().updateDiameter() - - return + self.opacity = np.ones( dataWrapper.numObj() ) * 0.5 def drawForTheFirstTime( self, _scene ): self.cylinderDraw( _scene ) @@ -489,104 +540,58 @@ def drawForTheFirstTime( self, _scene ): class MooReacSystem( MooDrawable ): ''' Draws collection of line segments of defined dia and color''' def __init__( self, - mooObj, fieldInfo, - field = 'conc', - relativeObj = '.', + dataWrapper, colormap = 'jet', - lenScale = 1e0, diaScale = 1.0, autoscale = False, + lenScale = 1e0, diaScale = 1.0, fieldScale = 1.0, + autoscale = False, valMin = 0.0, valMax = 1.0 ): - MooDrawable.__init__( self, fieldInfo, field = field, - relativeObj = relativeObj, + MooDrawable.__init__( self, dataWrapper, colormap = colormap, lenScale = lenScale, - diaScale = diaScale, autoscale = autoscale, + diaScale = diaScale, fieldScale = fieldScale, + autoscale = autoscale, valMin = valMin, valMax = valMax ) - self.mooObj = mooObj - self.updateCoords() - - def updateCoords( self ): - activeCoords = [] - self.activeDia = [] - for pool in self.mooObj: - coords = pool.coords - meshType = pool.compartment.className - if meshType in ["NeuroMesh", "CylMesh", "PsdMesh"]: - # Unfortunately at present these return radius rather than - # diameter in argument 6. To fix. - # Make a cylinder - activeCoords.append( [coords[0:3], coords[3:6]] ) - self.activeDia.append( coords[6] * 2 ) - elif meshType == "SpineMesh": - # Spine entry has head[3], shaft[3], root[3], dia. - activeCoords.append( [coords[0:3], coords[3:6]] ) - self.activeDia.append( coords[9] ) - elif meshType == "PresynMesh": - # This returns diameter in argument 6. - # first vec is centre of base, second axis pointing - # toward postsyn - # Hack: make each bouton as a cone with length == dia. - activeCoords.append( [coords[0:3], coords[6]*coords[3:6] + coords[0:3]] ) - self.activeDia.append( coords[6] ) - # Returns centre as args 0,1,2, diameter as argument 3. - # Make a hemisphere - elif meshType == "EndoMesh": - # Make a sphere. - activeCoords.append( [ coords[0:3], coords[0:3] ] ) - self.activeDia.append( coords[3] ) - self.activeCoords = np.array( activeCoords ) * self.lenScale - self.activeDia = np.array( self.activeDia ) * self.diaScale - self.opacity = np.ones( len( self.activeDia ) ) - self.activeObjs = self.mooObj - return + self.opacity = np.ones( dataWrapper.numObj() ) + def drawForTheFirstTime( self, _scene ): - if len( self.mooObj ) == 0: + if self.dataWrapper_.numObj() == 0: return - meshType = self.mooObj[0].compartment.className - if meshType in ["NeuroMesh", "CylMesh", "SpineMesh", "PsdMesh"]: + mt = self.dataWrapper_.meshType() + if mt in ["NeuroMesh", "CylMesh", "SpineMesh", "PsdMesh"]: self.cylinderDraw( _scene ) - elif meshType == "SpineMesh": + elif mt == "SpineMesh": self.spineDraw( _scene ) - elif meshType == "PresynMesh": + elif mt == "PresynMesh": self.presynDraw( _scene ) - elif meshType == "EndoMesh": + elif mt == "EndoMesh": self.endoDraw( _scene ) def spineDraw( self, _scene ): # Spine entry has head[3], shaft[3], root[3], dia. - for idx, coord in enumerate( self.activeCoords ): - v0 = list2vec( coord[0] ) - v1 = list2vec( coord[1] ) - self.coordMin = np.minimum( self.coordMin, coord[0][0:3] ) - self.coordMin = np.minimum( self.coordMin, coord[1][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[0][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[1][0:3] ) - radius = self.diaScale * self.activeDia[idx] / 2.0 + for idx, coord in enumerate( self.dataWrapper_.getCoords() ): + v0 = list2vec( coord[0:3] ) + v1 = list2vec( coord[3:6] ) + radius = self.diaScale * coord[6] / 2.0 opacity = self.opacity[idx] rod = vp.cylinder( canvas = _scene, pos = v0, axis = v1 - v0, radius = radius, opacity = opacity ) self.segments.append( rod ) def presynDraw( self, _scene ): - for idx, coord in enumerate( self.activeCoords ): - v0 = list2vec( coord[0] ) - v1 = list2vec( coord[1] ) - self.coordMin = np.minimum( self.coordMin, coord[0][0:3] ) - self.coordMin = np.minimum( self.coordMin, coord[1][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[0][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[1][0:3] ) - radius = self.diaScale * self.activeDia[idx] / 2.0 + for idx, coord in enumerate( self.dataWrapper_.getCoords() ): + v0 = list2vec( coord[0:3] ) + v1 = list2vec( coord[3:6] ) + radius = self.diaScale * coord[6] / 2.0 opacity = self.opacity[idx] cone = vp.cone( canvas = _scene, pos = v0, axis = v0 - v1, radius = radius, opacity = opacity ) self.segments.append( cone ) def endoDraw( self, _scene ): - for idx, coord in enumerate( self.activeCoords ): - v0 = list2vec( coord[0] ) - v1 = list2vec( coord[1] ) - self.coordMin = np.minimum( self.coordMin, coord[0][0:3] ) - self.coordMax = np.maximum( self.coordMax, coord[0][0:3] ) - radius = self.diaScale * self.activeDia[idx] / 2.0 + for idx, coord in enumerate( self.dataWrapper_.getCoords() ): + v0 = list2vec( coord[0:3] ) + v1 = list2vec( coord[3:6] ) + radius = self.diaScale * coord[6] / 2.0 opacity = self.opacity[idx] sphere = vp.sphere( canvas = _scene, pos = (v0 + v1)/2.0, radius = radius, opacity = opacity ) self.segments.append( sphere ) diff --git a/python/rdesigneur/nsdfview.py b/python/rdesigneur/nsdfview.py new file mode 100644 index 0000000000..efc01e68d1 --- /dev/null +++ b/python/rdesigneur/nsdfview.py @@ -0,0 +1,215 @@ +# nmoogli.py: nsdf Moogli interface +# This program displays NSDF files. +# Copyright (C) Upinder S. Bhalla NCBS 2022 +# This program is licensed under the GNU Public License version 3. + +import argparse +import numpy as np +import h5py +import time +#import rdesigneur.moogul as moogul +import moogul +mooViews = [] + +class ObjHandle: + def __init__( self, path ): + self.path = str( path, "utf-8" ) + +class NsdfNeuronDataWrapper( moogul.DataWrapper ): + def __init__( self, nsdf, path, field ): + moogul.DataWrapper.__init__( self, field ) + self.nsdf_ = nsdf + self.path_ = path + spl = path.split( '/', 1 ) + self.objBase_ = "%elec" + if len( spl ) > 1: # Some relative path tree. + if spl[0] == "#": + self.objRel_ = "##[ISA=CompartmentBase]%" + spl[1].replace( '/', '%' ) + else: + self.objRel_ = spl[0] + "%" + spl[1].replace( '/', '%' ) + else: + if spl[0] == "#": + self.objRel_ = "##[ISA=CompartmentBase]" + else: + self.objRel_ = path # No slashes in the string. + self.field_ = field + self.simTime_ = 0.0 + self.idx_ = 0 + self.dt_= nsdf["/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, field)].attrs['dt'] + self.coords_ = np.array( nsdf['/data/static/{}/{}/coords'.format(self.objBase_, self.objRel_) ] ) + objPaths = nsdf["/map/static/{}/{}/coords".format( self.objBase_, self.objRel_)] + self.objList_ = [ ObjHandle( path ) for path in objPaths ] + #print( "COORDS SHAPE for ", neuronName, " = ", self.coords_.shape ) + self.getMinMax() + + def getValues( self ): + npath = "/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, self.field_ ) + ret=np.array( self.nsdf_[npath][:,self.idx_] ) + #ret = np.array( nsdf["/data/uniform/{}/{}".format( neuronName, field)][self.idx_] ) + self.idx_ += 1 + self.simTime_ = self.idx_ * self.dt_ + return ret + + def getCoords( self ): + ''' Obtains 2-D array [comptidx, coord#] from the associated cell. + There can be any number of rows, but only 7 columns (i.e, coords). + These are x0, y0, z0, x, y, z, dia + ''' + return self.coords_ + + def getShape( self ): + npath = "/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, self.field_ ) + return self.nsdf_[npath].shape + + def numObj( self ): + return len( self.coords_ ) + + def getHistory( self, path, field ): + if field != self.field_: + print( "NsdfNeuronDataWrapper:getHistory Error: field name does not match: ", field, self.field_ ) + return + objIdx = [ idx for idx, val in enumerate( self.objList_ ) if val.path == path ] + if len( objIdx ) == 0: + print( "NsdfNeuronDataWrapper:getHistory Error: path not found: ", path ) + return + + npath = "/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, self.field_ ) + val = self.nsdf_[npath][objIdx[0],:] + t = np.arange( 0.0, len( val ), 1.0 ) * self.dt_ + return t, val + +##################################################################### + +class NsdfChemDataWrapper( moogul.DataWrapper ): + def __init__( self, nsdf, objname, field ): + if len( objname.split('/' ) ) != 2: + print( "Error: objname needs '/' between base and rel parts: ", objname ) + assert( 0 ) + moogul.DataWrapper.__init__( self, field ) + self.nsdf_ = nsdf + spl = objname.split( '/', 1 ) + if len( spl ) > 1: # Can't deal with case where we have /chem as base but multiple levels below it. + self.objBase_ = '%' + spl[0] + self.objRel_ = spl[1].replace( '/', '%' ) + if self.objRel_[-1] != ']': # put in wildcard, if none provided + self.objRel_ += '[]' + else: # obj is direct child of /chem. + self.objBase_ = "%chem" + self.objRel_ = objname # doesn't have any slashes either. + if self.objRel_[-1] != ']': # put in wildcard, if none provided + self.objRel_ += '[]' + + self.field_ = field + self.simTime_ = 0.0 + self.idx_ = 0 + self.dt_= nsdf["/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, field)].attrs['dt'] + self.coords_ = np.array( nsdf['/data/static/{}/{}/coords'.format(self.objBase_, self.objRel_) ] ) + if self.coords_.shape[1] == 10: + self.coords_[:,6] = self.coords_[:,9] # Temp hack to get radius correctly + objPaths = nsdf["/map/static/{}/{}/coords".format( self.objBase_, self.objRel_)] + self.objList_ = [ ObjHandle( path ) for path in objPaths ] + + self.getMinMax() + + def getValues( self ): + npath = "/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, self.field_ ) + ret=np.array( self.nsdf_[npath][:,self.idx_] ) + self.idx_ += 1 + self.simTime_ = self.idx_ * self.dt_ + return ret + + def getCoords( self ): + ''' Obtains 2-D array [comptidx, coord#] from the associated cell. + There can be any number of rows, but only 7 columns (i.e, coords). + These are x0, y0, z0, x, y, z, dia + ''' + return self.coords_ + + def meshType( self ): + #return self.meshType_ + return "NeuroMesh" + + def getShape( self ): + npath = "/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, self.field_ ) + return self.nsdf_[npath].shape + + def numObj( self ): + return len( self.coords_ ) + + def getHistory( self, path, field ): + if field != self.field_: + print( "NsdfChemDataWrapper:getHistory Error: field name does not match: ", field, self.field_ ) + return + objIdx = [ idx for idx, val in enumerate( self.objList_ ) if val.path == path ] + if len( objIdx ) == 0: + print( "NsdfChemDataWrapper:getHistory Error: path not found: ", path ) + return + + npath = "/data/uniform/{}/{}/{}".format( self.objBase_, self.objRel_, self.field_ ) + val = self.nsdf_[npath][objIdx[0],:] + t = np.arange( 0, len( val ), 1 ) * self.dt_ + return t, val + + +##################################################################### +def main(): + parser = argparse.ArgumentParser( description = "NSDF Moogli viewer." ) + parser.add_argument( "NSDF_filename", type = str, help = "Required: name of NSDF format File record of simulation" ) + parser.add_argument( "-v", "--viewspec", nargs = '+', default=[], action='append', help="Specification for each view: [objname, field, min, max]. Any number of views may be specified, each indicated by -v or --viewspec." ) + parser.add_argument( "-r", "--rotation", type = float, default=0.0, help="Rotate display around vertical axis by this angle in radians every step. Default 0.0") + parser.add_argument( "-c", "--colormap", type = str, default="plasma", help="Name of matplotlib colormap to use. Default is 'plasma'") + parser.add_argument( "-bg", "--background", type = str, default="default", help="Name of matplotlib color to use for background. Default is a light blue-gray.") + parser.add_argument( '-m', '--merge_displays', action="store_true", help="Display all data in the same view." ) + parser.add_argument( '-l', '--list_datasets', action="store_true", help="List possible datasets available to view." ) + args = parser.parse_args() + + if len( args.viewspec ) == 0: + print( "warning: No viewpsec defined in command line" ) + quit() + + nsdf = h5py.File( args.NSDF_filename, 'r' ) + + viewer = [] + for vs in args.viewspec: + viewer.append( makeMoogli( nsdf, vs ) ) + dt = viewer[0].drawables_[0].dataWrapper_.dt_ + shape = viewer[0].drawables_[0].dataWrapper_.getShape() + numSteps = shape[1] + for v in viewer: + v.firstDraw( args.merge_displays, rotation = args.rotation, colormap = args.colormap, bg = args.background ) + + simTime = 0.0 + for step in range( numSteps ): + for v in viewer: + v.updateValues( simTime ) + simTime += dt + viewer[0].notifySimulationEnd() + + while True: + time.sleep(1) + + +def makeMoogli( nsdf, args ): + objname, field, symin, symax = args + ymin = float( symin ) + ymax = float( symax ) + fieldScale = 1.0 + diaScale = 1.0 + + viewer = moogul.MooView( title = objname + "." + field ) + if field == 'n' or field == 'conc': + dw = NsdfChemDataWrapper( nsdf, objname, field ) + if field == 'conc': + fieldScale = 1.0e3 + reacSystem = moogul.MooReacSystem( dw, + valMin = ymin, valMax = ymax, diaScale = diaScale, fieldScale = fieldScale ) + viewer.addDrawable( reacSystem ) + else: + dw = NsdfNeuronDataWrapper( nsdf, objname, field ) + neuron = moogul.MooNeuron( dw, + valMin = ymin, valMax = ymax, diaScale = diaScale, fieldScale = fieldScale ) + viewer.addDrawable( neuron ) + return viewer + +if __name__ == "__main__": + main() diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py index 72e8c8de19..4026433b6b 100644 --- a/python/rdesigneur/rdesigneur.py +++ b/python/rdesigneur/rdesigneur.py @@ -46,6 +46,21 @@ meshOrder = ['soma', 'dend', 'spine', 'psd', 'psd_dend', 'presyn_dend', 'presyn_spine', 'endo'] +knownFieldsDefault = { + 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)', -80.0, 40.0 ), + 'initVm':('CompartmentBase', 'getInitVm', 1000, 'Init. Memb. Potl (mV)', -80.0, 40.0 ), + 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)', -10.0, 10.0 ), + 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)', -10.0, 10.0 ), + 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)', 0.0, 1.0 ), + 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)', 0.0, 1.0 ), + 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)', -10.0, 10.0 ), + 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)', -10.0, 10.0 ), + 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)', 0.0, 10.0 ), + 'n':('PoolBase', 'getN', 1, '# of molecules', 0.0, 200.0 ), + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)', 0.0, 2.0 ), + 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ) +} + #EREST_ACT = -70e-3 def _profile(func): @@ -102,6 +117,7 @@ def __init__(self, elecPlotDt = 0.1e-3, # Same default as from MOOSE funcDt = 0.1e-3, # Used when turnOffElec is False. # Otherwise system uses chemDt. + statusDt = 0.0, # Dt to print out status. 0 = no print numWaveFrames = 100, # Number of frames to use for waveplots cellProto = [], spineProto = [], @@ -115,6 +131,8 @@ def __init__(self, stimList = [], plotList = [], # elecpath, geom_expr, object, field, title ['wave' [min max]] moogList = [], + outputFileList = [], # List of all file save specifications. + modelFileNameList = [], # List of any files used to build. ode_method = "gsl", # gsl, lsoda, gssa, gillespie isLegacyMethod = False, params = None @@ -145,6 +163,7 @@ def __init__(self, self.elecDt= elecDt self.elecPlotDt= elecPlotDt self.funcDt= funcDt + self.statusDt= statusDt self.chemPlotDt= chemPlotDt self.numWaveFrames = numWaveFrames self.isLegacyMethod = isLegacyMethod @@ -158,6 +177,7 @@ def __init__(self, self.spineDistrib = spineDistrib self.chanDistrib = chanDistrib self.chemDistrib = chemDistrib + self.modelFileNameList = modelFileNameList self.params = params @@ -166,6 +186,7 @@ def __init__(self, self.stimList = [ rstim.convertArg(i) for i in stimList ] self.plotList = [ rplot.convertArg(i) for i in plotList ] self.moogList = [ rmoog.convertArg(i) for i in moogList ] + self.outputFileList = [ rfile.convertArg(i) for i in outputFileList ] except BuildError as msg: print("Error: rdesigneur: " + msg) quit() @@ -175,10 +196,13 @@ def __init__(self, self.wavePlotNames = [] self.saveNames = [] self.moogNames = [] + self.fileDumpNames = [] self.cellPortionElist = [] self.spineComptElist = [] self.tabForXML = [] self._endos = [] + self.nsdfPathList = [] # List of paths of nsdf objects. + self._finishedSaving = False if not moose.exists( '/library' ): library = moose.Neutral( '/library' ) @@ -215,9 +239,11 @@ def buildModel( self, modelPath = '/model' ): self.model = moose.Neutral( modelPath ) self.modelPath = modelPath funcs = [self.installCellFromProtos, self.buildPassiveDistrib - , self.buildChanDistrib, self.buildSpineDistrib, self.buildChemDistrib + , self.buildChanDistrib, self.buildSpineDistrib + , self.buildChemDistrib , self._configureSolvers, self.buildAdaptors, self._buildStims - , self._buildPlots, self._buildMoogli, self._configureHSolve + , self._buildPlots, self._buildMoogli, self._buildFileOutput + , self._configureHSolve , self._configureClocks, self._printModelStats] for i, _func in enumerate(funcs): @@ -237,6 +263,13 @@ def buildModel( self, modelPath = '/model' ): msg += ' %.3f sec' % t print(msg) sys.stdout.flush() + if self.statusDt > min( self.elecDt, self.chemDt, self.diffDt ): + pr = moose.PyRun( modelPath + '/updateStatus' ) + pr.initString = "_status_t0 = time.time()" + pr.runString = ''' +print( "Wall Clock Time = {:8.2f}, simtime = {:8.3f}".format( time.time() - _status_t0, moose.element( '/clock' ).currentTime ), flush=True ) +''' + moose.setClock( pr.tick, self.statusDt ) def installCellFromProtos( self ): if self.stealCellFromLibrary: @@ -961,22 +994,8 @@ def _buildPlots( self ): q += 1 def _buildMoogli( self ): - knownFields = { - 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)', -80.0, 40.0 ), - 'initVm':('CompartmentBase', 'getInitVm', 1000, 'Init. Memb. Potl (mV)', -80.0, 40.0 ), - 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)', -10.0, 10.0 ), - 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)', -10.0, 10.0 ), - 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)', 0.0, 1.0 ), - 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)', 0.0, 1.0 ), - 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)', -10.0, 10.0 ), - 'ICa':('NMDAChan', 'getICa', 1e9, 'Ca current (nA)', -10.0, 10.0 ), - 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)', 0.0, 10.0 ), - 'n':('PoolBase', 'getN', 1, '# of molecules', 0.0, 200.0 ), - 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)', 0.0, 2.0 ), - 'volume':('PoolBase', 'getVolume', 1e18, 'Volume (um^3)' ) - } + knownFields = knownFieldsDefault moogliBase = moose.Neutral( self.modelPath + '/moogli' ) - k = 0 for i in self.moogList: kf = knownFields[i.field] pair = i.elecpath + " " + i.geom_expr @@ -989,6 +1008,57 @@ def _buildMoogli( self ): numMoogli = len( dendObj ) self.moogNames.append( rmoogli.makeMoogli( self, dendObj, i, kf ) ) + def _buildFileOutput( self ): + fileBase = moose.Neutral( self.modelPath + "/file" ) + knownFields = knownFieldsDefault + + nsdfBlocks = {} + nsdf = None + + for idx, fentry in enumerate( self.outputFileList ): + oname = fentry.fname.split( "." )[0] + if fentry.ftype in ["h5", "nsdf"]: + # Should check for duplication. + nsdfPath = fileBase.path + '/' + oname + if fentry.field in ["n", "conc"]: + modelPath = self.modelPath + "/chem" + basePath = modelPath + "/" + fentry.path + if fentry.path[-1] in [']', '#']: # explicit index + pathStr = basePath + "." + fentry.field + else: + pathStr = basePath + "[]." + fentry.field + else: + modelPath = self.modelPath + "/elec" + spl = fentry.path.split('/') + if spl[0] == "#": + if len( spl ) == 1: + fentry.path = "##[ISA=CompartmentBase]" + else: + fentry.path = "##[ISA=CompartmentBase]" + fentry.path[1:] + # Otherwise we use basepath as is. + basePath = modelPath + "/" + fentry.path + pathStr = basePath + "." + fentry.field + if not nsdfPath in nsdfBlocks: + self.nsdfPathList.append( nsdfPath ) + nsdfBlocks[nsdfPath] = [pathStr] + nsdf = moose.NSDFWriter2( nsdfPath ) + nsdf.modelRoot = "" # Blank means don't save tree. + nsdf.filename = fentry.fname + # Insert the model setup files here. + nsdf.mode = 2 + nsdf.flushLimit = fentry.flushSteps # Number of timesteps between flush + nsdf.tick = 20 + len( nsdfBlocks ) + moose.setClock( nsdf.tick, fentry.dt ) + mfns = sys.argv[0] + for ii in self.modelFileNameList: + mfns += "," + ii + nsdf.modelFileNames = mfns + else: + nsdfBlocks[nsdfPath].append( pathStr ) + for nsdfPath in self.nsdfPathList: + nsdf = moose.element( nsdfPath ) + nsdf.blocks = nsdfBlocks[nsdfPath] + ################################################################ # Here we display the plots and moogli @@ -1006,11 +1076,13 @@ def displayMoogli( self, moogliDt, runtime, rotation = math.pi/500.0, fullscreen moose.setClock( pr.tick, moogliDt ) moose.reinit() moose.start( runtime ) + self._save() rmoogli.notifySimulationEnd() if block: self.display( len( self.moogNames ) + 1) def display( self, startIndex = 0, block=True ): + self._save() for i in self.plotNames: plt.figure( i[2] + startIndex ) plt.title( i[1] ) @@ -1037,7 +1109,6 @@ def display( self, startIndex = 0, block=True ): for i in range( 3 ): self.displayWavePlots() plt.show( block=block ) - self._save() def initWavePlots( self, startIndex ): @@ -1151,6 +1222,10 @@ def _writeCSV( self, plotData, time, vtab ): def _save( self ): + self._finishedSaving = True + for nsdfPath in self.nsdfPathList: + nsdf = moose.element( nsdfPath ) + nsdf.close() for i in self.saveNames: tabname = i[0] idx = i[1] @@ -1186,7 +1261,7 @@ def _buildStims( self ): } stims = moose.Neutral( self.modelPath + '/stims' ) k = 0 - # rstim class has {elecpath, geom_expr, relpath, field, expr} + # rstim class has {fname, path, field, dt, flush_steps } for i in self.stimList: pair = i.elecpath + " " + i.geom_expr dendCompts = self.elecid.compartmentsFromExpression[ pair ] @@ -1393,9 +1468,9 @@ def _decorateWithSpines( self ): ################################################################ def _loadElec( self, efile, elecname ): + self.modelFileNameList.append( efile ) if ( efile[ len( efile ) - 2:] == ".p" ): self.elecid = moose.loadModel( efile, '/library/' + elecname) - print(self.elecid) elif ( efile[ len( efile ) - 4:] == ".swc" ): self.elecid = moose.loadModel( efile, '/library/' + elecname) else: @@ -1429,32 +1504,12 @@ def _loadElec( self, efile, elecname ): # with those names and volumes in decreasing order. def validateChem( self ): cpath = self.chemid.path - comptlist = moose.wildcardFind( cpath + '/#[ISA=ChemCompt]' ) + comptlist = moose.wildcardFind( cpath + '/##[ISA=ChemCompt],' ) if len( comptlist ) == 0: raise BuildError( "validateChem: no compartment on: " + cpath ) return True - ''' - if len( comptlist ) == 1: - return; - - # Sort comptlist in decreasing order of volume - sortedComptlist = sorted( comptlist, key=lambda x: -x.volume ) - if ( len( sortedComptlist ) != 3 ): - print(cpath, sortedComptlist) - raise BuildError( "validateChem: Require 3 chem compartments, have: " + str( len( sortedComptlist ) ) ) - ''' - ''' - if not( sortedComptlist[0].name.lower() == 'dend' and \ - sortedComptlist[1].name.lower() == 'spine' and \ - sortedComptlist[2].name.lower() == 'psd' ): - raise BuildError( "validateChem: Invalid compt names: require dend, spine and PSD.\nActual names = " \ - + sortedComptList[0].name + ", " \ - + sortedComptList[1].name + ", " \ - + sortedComptList[2].name ) - ''' - ################################################################# def _isModelFromKkit_SBML( self ): @@ -1691,13 +1746,14 @@ def _oldConfigureSolvers( self ) : ################################################################ def _loadChem( self, fname, chemName ): + self.modelFileNameList.append( fname ) chem = moose.Neutral( '/library/' + chemName ) pre, ext = os.path.splitext( fname ) if ext == '.xml' or ext == '.sbml': modelId = moose.readSBML( fname, chem.path ) else: modelId = moose.loadModel( fname, chem.path, 'ee' ) - comptlist = moose.wildcardFind( chem.path + '/#[ISA=ChemCompt]' ) + comptlist = moose.wildcardFind( chem.path + '/##[ISA=ChemCompt]' ) if len( comptlist ) == 0: print("loadChem: No compartment found in file: ", fname) return @@ -1905,3 +1961,36 @@ def convertArg( arg ): else: raise BuildError( "rstim initialization failed" ) + +class rfile: + def __init__( self, + fname = 'output.h5', path = 'soma', field = 'Vm', dt = 1e-4, flushSteps = 200, start = 0.0, stop = -1.0, ftype = 'nsdf'): + self.fname = fname + self.path = path + if not field in knownFieldsDefault: + print( "Error: Field '{}' not known.".format( field ) ) + assert( 0 ) + self.field = field + self.dt = dt + self.flushSteps = flushSteps + self.start = start + self.stop = stop + self.ftype = self.fname.split(".")[-1] + if not self.ftype in ["txt", "csv", "h5", "nsdf"]: + print( "Error: output file format for ", fname , " not known") + assert( 0 ) + self.fname = self.fname.split("/")[-1] + + def printme( self ): + print( "{0}, {1}, {2}, {3}".format( + self.fname, self.path, self.field, self.dt) ) + + @staticmethod + def convertArg( arg ): + if isinstance( arg, rfile ): + return arg + elif isinstance( arg, list ): + return rfile( *arg ) + else: + raise BuildError( "rfile initialization failed" ) + diff --git a/python/rdesigneur/rmoogli.py b/python/rdesigneur/rmoogli.py index 05653280c9..c312cca79b 100644 --- a/python/rdesigneur/rmoogli.py +++ b/python/rdesigneur/rmoogli.py @@ -7,9 +7,96 @@ # Copyright (C) Upinder S. Bhalla NCBS 2018 # This program is licensed under the GNU Public License version 3. +import numpy as np import rdesigneur.moogul as moogul +import moose mooViews = [] +class MooseNeuronDataWrapper( moogul.DataWrapper ): + def __init__( self, neuronId, relativeObjPath, field ): + moogul.DataWrapper.__init__( self, field ) + self.neuronId_ = neuronId + self.relativeObjPath_ = relativeObjPath + if field == "Ca": + self.dummyObj = moose.CaConc( "/dummyCa" ) + elif field in ["Ik","Gk", "Ek", "Gbar", "modulation"]: + self.dummyObj = moose.SynChan( "/dummyCa" ) + + compts = moose.wildcardFind( neuronId.path + "/#[ISA=CompartmentBase]" ) + self.coords_ = np.array( [ii.coords for ii in compts] ) + self.getMinMax() + if relativeObjPath == ".": + self.objList_ = compts + else: + self.objList_ = [] + for i in compts: + if moose.exists( i.path + '/' + relativeObjPath ): + elm = moose.element(i.path + '/' + relativeObjPath ) + self.objList_.append( elm ) + else: + self.objList_.append( self.dummyObj ) + + def getValues( self ): + return np.array( [moose.getField(i, self.field_) for i in self.objList_] ) + + def getCoords( self ): + ''' Obtains 2-D array [comptidx, coord#] from the associated cell. + There can be any number of rows, but only 7 columns (i.e, coords). + These are x0, y0, z0, x, y, z, dia + ''' + return self.coords_ + + +class MooseChemDataWrapper( moogul.DataWrapper ): + def __init__( self, objList, field ): + moogul.DataWrapper.__init__( self, field ) + self.objList_ = objList + if len( objList ) == 0: + return + #coords = np.array( [ii.coords for ii in objList] ) + self.coords_ = np.array( [ii.coords for ii in objList] ) + self.meshType_ = objList[0].compartment.className + if self.meshType_ in ["NeuroMesh", "CylMesh", "PsdMesh"]: + # Unfortunately at present these return radius rather than + # diameter in argument 6. To fix. + self.coords_[:,6] *= 2.0 + elif self.meshType_ == "SpineMesh": + # Spine entry has head[3], shaft[3], root[3], dia. + # We want to put dia in idx == 6. + self.coords_[:,6] = self.coords_[:,9] + elif self.meshType_ == "PresynMesh": + # This returns diameter in argument 6. + # first vec is centre of base, second axis pointing + # toward postsyn + # Hack: make each bouton as a cone with length == dia. + diaColumn = np.array(self.coords_[:,6]) + diaColumn.shape = ( len( diaColumn ), 1 ) + self.coords_[:,3:6] = self.coords_[:,0:3] + self.coords_[:,3:6] * diaColumn + elif self.meshType_ == "EndoMesh": + # Returns centre as args 0,1,2, diameter as argument 3. + # Make a sphere. The coords need to be set for 7 columns. + temp = np.array( self.coords_ ) + self.coords_ = np.zeros( (temp.shape[0], 7 ) ) + self.coords_[:,0:3] = temp[:,0:3] + self.coords_[:,3:6] = temp[:,0:3] + self.coords_[:,6] = temp[:,3] + #= np.array( self.coords_[0:3].extend( self.coords_[0:3] ).extend(self.coords_[3] ) ) + self.getMinMax() + + def getValues( self ): + return np.array( [moose.getField(i, self.field_) for i in self.objList_] ) + + def getCoords( self ): + ''' Obtains 2-D array [comptidx, coord#] from the associated cell. + There can be any number of rows, but only 7 columns (i.e, coords). + These are x0, y0, z0, x, y, z, dia + ''' + return self.coords_ + + def meshType( self ): + return self.meshType_ + + def makeMoogli( rd, mooObj, args, fieldInfo ): #mooObj is currently poorly handled. Ideally it should simply be a # vector of objects to plot, each with coords. This could be readily @@ -33,23 +120,21 @@ def makeMoogli( rd, mooObj, args, fieldInfo ): viewer = moogul.MooView( title = args.title ) if mooField == 'n' or mooField == 'conc': - #moogul.updateDiffCoords( mooObj ) - reacSystem = moogul.MooReacSystem( mooObj, fieldInfo, - field = mooField, relativeObj = relObjPath, - valMin = ymin, valMax = ymax, diaScale = args.diaScale ) + dw = MooseChemDataWrapper( mooObj, mooField ) + reacSystem = moogul.MooReacSystem( dw, + valMin = ymin, valMax = ymax, diaScale = args.diaScale, fieldScale = fieldInfo[2] ) viewer.addDrawable( reacSystem ) else: - neuron = moogul.MooNeuron( rd.elecid, fieldInfo, - field = mooField, relativeObj = relObjPath, - valMin = ymin, valMax = ymax, diaScale = args.diaScale ) - #print( "min = {}, max = {}".format(ymin, ymax) ) + dw = MooseNeuronDataWrapper( rd.elecid, relObjPath, mooField ) + neuron = moogul.MooNeuron( dw, + valMin = ymin, valMax = ymax, diaScale = args.diaScale, fieldScale = fieldInfo[2] ) viewer.addDrawable( neuron ) - return viewer def updateMoogliViewer(): + simTime = moose.element( '/clock' ).currentTime for i in mooViews: - i.updateValues() + i.updateValues( simTime ) def notifySimulationEnd(): if len( mooViews ) > 0: diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp index e06e995fc7..0bc3a8970b 100644 --- a/scheduling/Clock.cpp +++ b/scheduling/Clock.cpp @@ -421,6 +421,7 @@ const Cinfo* Clock::initCinfo() " HDF5DataWriter 30 1\n" " HDF5WriterBase 30 1\n" " NSDFWriter 30 1\n" + " NSDFWriter2 30 1\n" " PyRun 30 1\n" " PostMaster 31 0.01\n" " \n" @@ -929,6 +930,7 @@ void Clock::buildDefaultTick() defaultTick_["HDF5DataWriter"] = 30; defaultTick_["HDF5WriterBase"] = 30; defaultTick_["NSDFWriter"] = 30; + defaultTick_["NSDFWriter2"] = 30; defaultTick_["PyRun"] = 30; defaultTick_["PostMaster"] = 31;