From d7ebbe792121758eaf4b98f438f65970e892a520 Mon Sep 17 00:00:00 2001 From: Dilawar Singh Date: Wed, 22 Jun 2016 10:08:43 +0530 Subject: [PATCH 1/4] No need to have header name in inverted commas. --- builtins/StreamerBase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/builtins/StreamerBase.cpp b/builtins/StreamerBase.cpp index 807cc49e8b..90904c34b7 100644 --- a/builtins/StreamerBase.cpp +++ b/builtins/StreamerBase.cpp @@ -98,7 +98,7 @@ void StreamerBase::writeToCSVFile( const string& filepath, const string& openmod string headerText = ""; for( vector::const_iterator it = columns.begin(); it != columns.end(); it++ ) - headerText += "\"" + *it + "\"" + delimiter_; + headerText += ( *it + delimiter_ ); headerText += eol; fprintf( fp, "%s", headerText.c_str() ); } From 002edb63982a6d6f0367f75aea0166d9f24ae70a Mon Sep 17 00:00:00 2001 From: Dilawar Singh Date: Wed, 22 Jun 2016 10:54:51 +0530 Subject: [PATCH 2/4] Fix to BhallaLab/moose-core#130. --- builtins/Streamer.cpp | 17 +++++++++++++++-- scheduling/Clock.cpp | 19 ++++++++++--------- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/builtins/Streamer.cpp b/builtins/Streamer.cpp index e6170967c9..31eedecb38 100644 --- a/builtins/Streamer.cpp +++ b/builtins/Streamer.cpp @@ -191,6 +191,7 @@ void Streamer::reinit(const Eref& e, ProcPtr p) tableDt_.push_back( clk->getTickDt( tickNum ) ); } + // Make sure all tables have same dt_ else disable the streamer. vector invalidTables; for (size_t i = 1; i < tableTick_.size(); i++) @@ -238,10 +239,9 @@ void Streamer::process(const Eref& e, ProcPtr p) // Prepare data. zipWithTime( data_, currTime_ ); StreamerBase::writeToOutFile( outfilePath_, format_, "a", data_, columns_ ); + // clean the arrays data_.clear(); - for(size_t i = 0; i < tables_.size(); i++ ) - tables_[i]->clearVec(); } @@ -252,6 +252,15 @@ void Streamer::process(const Eref& e, ProcPtr p) */ void Streamer::addTable( Id table ) { + Clock* clk = reinterpret_cast( Id(1).eref().data() ); + int tickNum = table.element()->getTick(); + double tick = clk->getTickDt( tickNum ); + tableDt_.push_back( tick ); + + // Set tick of stramer to 100 times of table tick or 10 seconds whichever is + // higher. + clk->setTickDt( 19, max( 10.0, 100 * tick ) ); + // If this table is not already in the vector, add it. for( size_t i = 0; i < tableIds_.size(); i++) if( table.path() == tableIds_[i].path() ) @@ -371,4 +380,8 @@ void Streamer::zipWithTime( vector& data, double currTime) for( size_t i = 0; i < tables_.size(); i++) data.push_back( tables_[i]->getVec()[i] ); } + + // clear the data from tables now. + for(size_t i = 0; i < tables_.size(); i++ ) + tables_[i]->clearVec(); } diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp index f532d9c2ea..61444f2938 100644 --- a/scheduling/Clock.cpp +++ b/scheduling/Clock.cpp @@ -413,9 +413,9 @@ const Cinfo* Clock::initCinfo() " Gsolve 16 0.1\n" " Ksolve 16 0.1\n" " Stats 17 0.1\n" - " Table2 18 1\n" - " Streamer 29 10\n" + " Streamer 19 10\n" + " HDF5DataWriter 30 1\n" " HDF5WriterBase 30 1\n" " NSDFWriter 30 1\n" @@ -897,7 +897,7 @@ void Clock::buildDefaultTick() defaultTick_["Stats"] = 17; defaultTick_["Table2"] = 18; - defaultTick_["Streamer"] = 29; + defaultTick_["Streamer"] = 19; defaultTick_["HDF5DataWriter"] = 30; defaultTick_["HDF5WriterBase"] = 30; defaultTick_["NSDFWriter"] = 30; @@ -962,9 +962,9 @@ void Clock::buildDefaultTick() defaultDt_[5] = 50.0e-6; defaultDt_[6] = 50.0e-6; defaultDt_[7] = 50.0e-6; - defaultDt_[8] = 1.0e-4; // For the tables for electrical calculations - defaultDt_[9] = 0.0; // Not assigned - defaultDt_[10] = 0.01; // For diffusion. + defaultDt_[8] = 1.0e-4; // For the tables for electrical calculations + defaultDt_[9] = 0.0; // Not assigned + defaultDt_[10] = 0.01; // For diffusion. defaultDt_[11] = 0.1; defaultDt_[12] = 0.1; defaultDt_[13] = 0.1; @@ -972,9 +972,10 @@ void Clock::buildDefaultTick() defaultDt_[15] = 0.1; defaultDt_[16] = 0.1; defaultDt_[17] = 0.1; - defaultDt_[18] = 1; // For tables for chemical calculations. - // 19-28 are not assigned. - defaultDt_[29] = 10; // For Streamer + defaultDt_[18] = 1; // For tables for chemical calculations. + defaultDt_[19] = 10; // For Streamer + + // 20-29 are not assigned. defaultDt_[30] = 1; // For the HDF writer defaultDt_[31] = 0.01; // For the postmaster. } From 8e877070c8b815c3581fd7194c57680f9d0aef0f Mon Sep 17 00:00:00 2001 From: Dilawar Singh Date: Fri, 24 Jun 2016 16:43:53 +0530 Subject: [PATCH 3/4] Fixed the streamer test. --- tests/python/test_streamer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/python/test_streamer.py b/tests/python/test_streamer.py index 97884ea7e6..3c16b76c97 100644 --- a/tests/python/test_streamer.py +++ b/tests/python/test_streamer.py @@ -113,7 +113,7 @@ def test( ): # print( data.dtype ) time = data['time'] print( time ) - assert data.shape >= (58,4), data.shape + assert data.shape >= (58,), data.shape print( '[INFO] Test 2 passed' ) def main( ): From 2ad812f070520ec964549cbfc085910895cfecc7 Mon Sep 17 00:00:00 2001 From: Dilawar Singh Date: Sat, 25 Jun 2016 17:37:50 +0530 Subject: [PATCH 4/4] Rdesigeur by sarthak, --- builtins/Streamer.cpp | 32 +++-- python/rdesigneur/rdesigneur.py | 227 +++++++++++++++++++++++++++++++- scheduling/Clock.cpp | 1 + shell/Shell.cpp | 10 ++ 4 files changed, 255 insertions(+), 15 deletions(-) diff --git a/builtins/Streamer.cpp b/builtins/Streamer.cpp index 31eedecb38..93d4419411 100644 --- a/builtins/Streamer.cpp +++ b/builtins/Streamer.cpp @@ -175,7 +175,7 @@ void Streamer::cleanUp( void ) */ void Streamer::reinit(const Eref& e, ProcPtr p) { - Clock* clk = reinterpret_cast( Id(1).eref().data() ); + if( tables_.size() == 0 ) { moose::showWarn( "Zero tables in streamer. Disabling Streamer" ); @@ -183,6 +183,27 @@ void Streamer::reinit(const Eref& e, ProcPtr p) return; } + Clock* clk = reinterpret_cast( Id(1).eref().data() ); + for (size_t i = 0; i < tableIds_.size(); i++) + { + int tickNum = tableIds_[i].element()->getTick(); + double tick = clk->getTickDt( tickNum ); + tableDt_.push_back( tick ); + // Make sure that all tables have the same tick. + if( i > 0 ) + { + if( tick != tableDt_[0] ) + { + moose::showWarn( "Table " + tableIds_[i].path() + " has " + " different clock dt. " + " Make sure all tables added to Streamer have the same " + " dt value." + ); + } + } + } + + // Push each table dt_ into vector of dt for( size_t i = 0; i < tables_.size(); i++) { @@ -252,15 +273,6 @@ void Streamer::process(const Eref& e, ProcPtr p) */ void Streamer::addTable( Id table ) { - Clock* clk = reinterpret_cast( Id(1).eref().data() ); - int tickNum = table.element()->getTick(); - double tick = clk->getTickDt( tickNum ); - tableDt_.push_back( tick ); - - // Set tick of stramer to 100 times of table tick or 10 seconds whichever is - // higher. - clk->setTickDt( 19, max( 10.0, 100 * tick ) ); - // If this table is not already in the vector, add it. for( size_t i = 0; i < tableIds_.size(); i++) if( table.path() == tableIds_[i].path() ) diff --git a/python/rdesigneur/rdesigneur.py b/python/rdesigneur/rdesigneur.py index cb46f26b51..803d3144b5 100644 --- a/python/rdesigneur/rdesigneur.py +++ b/python/rdesigneur/rdesigneur.py @@ -1,5 +1,5 @@ ######################################################################### -## rdesigneur0_4.py --- +## rdesigneur0_5.py --- ## This program is part of 'MOOSE', the ## Messaging Object Oriented Simulation Environment. ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS @@ -27,6 +27,10 @@ from rdesigneurProtos import * from moose.neuroml.NeuroML import NeuroML from moose.neuroml.ChannelML import ChannelML +import lxml +from lxml import etree +import h5py as h5 +import csv #EREST_ACT = -70e-3 @@ -76,7 +80,8 @@ def __init__(self, adaptorList= [], stimList = [], plotList = [], - moogList = [] + moogList = [], + params = None ): """ Constructor of the rdesigner. This just sets up internal fields for the model building, it doesn't actually create any objects. @@ -108,14 +113,20 @@ def __init__(self, self.chanDistrib = chanDistrib self.chemDistrib = chemDistrib + self.params = params + self.adaptorList = adaptorList self.stimList = stimList self.plotList = plotList + self.saveList = plotList #ADDED BY Sarthak + self.saveAs = [] self.moogList = moogList self.plotNames = [] + self.saveNames = [] self.moogNames = [] self.cellPortionElist = [] self.spineComptElist = [] + self.tabForXML = [] if not moose.exists( '/library' ): library = moose.Neutral( '/library' ) @@ -172,6 +183,7 @@ def buildModel( self, modelPath = '/model' ): self._buildStims() self._configureClocks() self._printModelStats() + self._savePlots() except BuildError as msg: print("Error: rdesigneur: model build failed:", msg) @@ -485,6 +497,7 @@ def _parseComptField( self, comptList, plotSpec, knownFields ): kf = knownFields[field] # Find the field to decide type. if ( kf[0] == 'CaConcBase' or kf[0] == 'ChanBase' ): objList = self._collapseElistToPathAndClass( comptList, plotSpec[2], kf[0] ) + # print ("objList: ", len(objList), kf[1]) return objList, kf[1] elif (field == 'n' or field == 'conc' ): path = plotSpec[2] @@ -538,8 +551,6 @@ def _buildPlots( self ): pair = i[0] + " " + i[1] dendCompts = self.elecid.compartmentsFromExpression[ pair ] spineCompts = self.elecid.spinesFromExpression[ pair ] - #print( "DENDENDENDNEDN = ", len(dendCompts), pair ) - #print( "SPINESPINESPINE = ", len(spineCompts), pair ) plotObj, plotField = self._parseComptField( dendCompts, i, knownFields ) plotObj2, plotField2 = self._parseComptField( spineCompts, i, knownFields ) assert( plotField == plotField2 ) @@ -613,7 +624,212 @@ def display( self ): pylab.plot( t, j.vector * i[3] ) if len( self.moogList ) > 0: pylab.ion() - pylab.show() + pylab.show(block=True) + self._save() #This calls the _save function which saves only if the filenames have been specified + + ################################################################ + # Here we get the time-series data and write to various formats + ################################################################ + #[TO DO] Add NSDF output function + ''' + The author of the functions -- [_savePlots(), _getTimeSeriesTable(), _writeXML(), _writeCSV(), _saveFormats(), _save()] is + Sarthak Sharma. + Email address: sarthaks442@gmail.com + ''' + + def _savePlots( self ): + + knownFields = { + 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ), + 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ), + 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ), + 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ), + 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)' ), + 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)' ), + 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ), + 'n':('PoolBase', 'getN', 1, '# of molecules'), + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ) + } + + save_graphs = moose.Neutral( self.modelPath + '/save_graphs' ) + dummy = moose.element( '/' ) + k = 0 + + for i in self.saveList: + pair = i[0] + " " + i[1] + dendCompts = self.elecid.compartmentsFromExpression[ pair ] + spineCompts = self.elecid.spinesFromExpression[ pair ] + plotObj, plotField = self._parseComptField( dendCompts, i, knownFields ) + plotObj2, plotField2 = self._parseComptField( spineCompts, i, knownFields ) + assert( plotField == plotField2 ) + plotObj3 = plotObj + plotObj2 + numPlots = sum( i != dummy for i in plotObj3 ) + if numPlots > 0: + save_tabname = save_graphs.path + '/save_plot' + str(k) + scale = knownFields[i[3]][2] + units = knownFields[i[3]][3] + self.saveNames.append( ( save_tabname, i[4], k, scale, units ) ) + k += 1 + if i[3] == 'n' or i[3] == 'conc': + save_tabs = moose.Table2( save_tabname, numPlots ) + else: + save_tabs = moose.Table( save_tabname, numPlots ) + save_vtabs = moose.vec( save_tabs ) + q = 0 + for p in [ x for x in plotObj3 if x != dummy ]: + moose.connect( save_vtabs[q], 'requestOut', p, plotField ) + q += 1 + + def _getTimeSeriesTable( self ): + + ''' + This function gets the list with all the details of the simulation + required for plotting. + This function adds flexibility in terms of the details + we wish to store. + ''' + + knownFields = { + 'Vm':('CompartmentBase', 'getVm', 1000, 'Memb. Potential (mV)' ), + 'Im':('CompartmentBase', 'getIm', 1e9, 'Memb. current (nA)' ), + 'inject':('CompartmentBase', 'getInject', 1e9, 'inject current (nA)' ), + 'Gbar':('ChanBase', 'getGbar', 1e9, 'chan max conductance (nS)' ), + 'Gk':('ChanBase', 'getGk', 1e9, 'chan conductance (nS)' ), + 'Ik':('ChanBase', 'getIk', 1e9, 'chan current (nA)' ), + 'Ca':('CaConcBase', 'getCa', 1e3, 'Ca conc (uM)' ), + 'n':('PoolBase', 'getN', 1, '# of molecules'), + 'conc':('PoolBase', 'getConc', 1000, 'Concentration (uM)' ) + } + + ''' + This takes data from plotList + saveList is exactly like plotList but with a few additional arguments: + ->It will have a resolution option, i.e., the number of decimal figures to which the value should be rounded + ->There is a list of "saveAs" formats + With saveList, the user will able to set what all details he wishes to be saved. + ''' + + for i,ind in enumerate(self.saveNames): + pair = self.saveList[i][0] + " " + self.saveList[i][1] + dendCompts = self.elecid.compartmentsFromExpression[ pair ] + spineCompts = self.elecid.spinesFromExpression[ pair ] + # Here we get the object details from plotList + savePlotObj, plotField = self._parseComptField( dendCompts, self.saveList[i], knownFields ) + savePlotObj2, plotField2 = self._parseComptField( spineCompts, self.saveList[i], knownFields ) + savePlotObj3 = savePlotObj + savePlotObj2 + + rowList = list(ind) + save_vtab = moose.vec( ind[0] ) + t = np.arange( 0, save_vtab[0].vector.size, 1 ) * save_vtab[0].dt + + rowList.append(save_vtab[0].dt) + rowList.append(t) + rowList.append([jvec.vector * ind[3] for jvec in save_vtab]) #get values + rowList.append(self.saveList[i][3]) + rowList.append(filter(lambda obj: obj.path != '/', savePlotObj3)) #this filters out dummy elements + + if (type(self.saveList[i][-1])==int): + rowList.append(self.saveList[i][-1]) + else: + rowList.append(12) + + self.tabForXML.append(rowList) + rowList = [] + + timeSeriesTable = self.tabForXML # the list with all the details of plot + return timeSeriesTable + + def _writeXML( self, filename, timeSeriesData ): #to write to XML file + + plotData = timeSeriesData + print("[CAUTION] The '%s' file might be very large if all the compartments are to be saved." % filename) + root = etree.Element("TimeSeriesPlot") + parameters = etree.SubElement( root, "parameters" ) + if self.params == None: + parameters.text = "None" + else: + assert(isinstance(self.params, dict)), "'params' should be a dictionary." + for pkey, pvalue in self.params.items(): + parameter = etree.SubElement( parameters, str(pkey) ) + parameter.text = str(pvalue) + + #plotData contains all the details of a single plot + title = etree.SubElement( root, "timeSeries" ) + title.set( 'title', str(plotData[1])) + title.set( 'field', str(plotData[8])) + title.set( 'scale', str(plotData[3])) + title.set( 'units', str(plotData[4])) + title.set( 'dt', str(plotData[5])) + p = [] + assert(len(plotData[7]) == len(plotData[9])) + + res = plotData[10] + for ind, jvec in enumerate(plotData[7]): + p.append( etree.SubElement( title, "data")) + p[-1].set( 'path', str(plotData[9][ind].path)) + p[-1].text = ''.join( str(round(value,res)) + ' ' for value in jvec ) + + tree = etree.ElementTree(root) + tree.write(filename) + + def _writeCSV(self, filename, timeSeriesData): + + plotData = timeSeriesData + dataList = [] + header = [] + time = plotData[6] + res = plotData[10] + + for ind, jvec in enumerate(plotData[7]): + header.append(plotData[9][ind].path) + dataList.append([round(value,res) for value in jvec.tolist()]) + dl = [tuple(lst) for lst in dataList] + rows = zip(tuple(time), *dl) + header.insert(0, "time") + + with open(filename, 'wb') as f: + writer = csv.writer(f, quoting=csv.QUOTE_MINIMAL) + writer.writerow(header) + for row in rows: + writer.writerow(row) + + ##########****SAVING*****############### + def _saveFormats(self, timeSeriesData, k, *filenames): + "This takes in the filenames and writes to corresponding format." + if filenames: + for filename in filenames: + for name in filename: + print (name) + if name[-4:] == '.xml': + self._writeXML(name, timeSeriesData) + print(name, " written") + elif name[-4:] == '.csv': + self._writeCSV(name, timeSeriesData) + print(name, " written") + else: + print("not possible") + pass + else: + pass + + + def _save( self ): + timeSeriesTable = self._getTimeSeriesTable() + for i,sList in enumerate(self.saveList): + + if (len(sList) >= 6) and (type(sList[5]) != int): + self.saveAs.extend(filter(lambda fmt: type(fmt)!=int, sList[5:])) + try: + timeSeriesData = timeSeriesTable[i] + except IndexError: + print("The object to be plotted has all dummy elements.") + pass + self._saveFormats(timeSeriesData, i, self.saveAs) + self.saveAs=[] + else: + pass + else: + pass ################################################################ # Here we set up the stims @@ -1053,3 +1269,4 @@ def _buildAdaptor( self, meshName, elecRelPath, elecField, \ for j in range( i[1], i[2] ): moose.connect( i[3], 'requestOut', chemVec[j], chemFieldSrc) msg = moose.connect( i[3], 'output', elObj, elecFieldDest ) + diff --git a/scheduling/Clock.cpp b/scheduling/Clock.cpp index 61444f2938..a90cb7c9db 100644 --- a/scheduling/Clock.cpp +++ b/scheduling/Clock.cpp @@ -591,6 +591,7 @@ void Clock::setTickDt( unsigned int i, double v ) } for ( unsigned int j = 0; j < numTicks; ++j ) numUsed += ( ticks_[j] != 0 ); + if ( numUsed == 0 ) { dt_ = v; diff --git a/shell/Shell.cpp b/shell/Shell.cpp index dc6da32d34..85aca3e444 100644 --- a/shell/Shell.cpp +++ b/shell/Shell.cpp @@ -8,6 +8,8 @@ **********************************************************************/ #include +#include + using namespace std; @@ -396,6 +398,14 @@ void Shell::doStop( ) void Shell::doSetClock( unsigned int tickNum, double dt ) { LookupField< unsigned int, double >::set( ObjId( 1 ), "tickDt", tickNum, dt ); + + // FIXME: + // HACK: If clock 18 is being updated, make sure that clock 19 (streamer is also + // updated with correct dt (10 or 100*dt). This is bit hacky. + if( tickNum == 18 ) + LookupField< unsigned int, double >::set( ObjId( 1 ), "tickDt" + , tickNum + 1, max( 100 * dt, 10.0 ) + ); } void Shell::doUseClock( string path, string field, unsigned int tick )