diff --git a/brain/detail/circuit.h b/brain/detail/circuit.h index e9b94cd..f54c25e 100644 --- a/brain/detail/circuit.h +++ b/brain/detail/circuit.h @@ -149,6 +149,7 @@ class Circuit::Impl , _synapseSource( config.getSynapseSource( )) , _targetSources( config.getTargetSources( )) , _cache( keyv::Map::createCache( )) + , _synapsePositionColumns( 0 ) { } @@ -266,11 +267,17 @@ class Circuit::Impl const size_t i = afferent ? 0 : 1; lunchbox::ScopedWrite mutex( _synapsePositions[i] ); - if( !(*_synapsePositions[i] )) - _synapsePositions[i]->reset( + auto& positions = *_synapsePositions[i]; + if( !positions ) + positions.reset( new brion::Synapse( _synapseSource.getPath() + (afferent ? afferentPositionsFilename : efferentPositionsFilename ))); - return **_synapsePositions[i]; + + if( _synapsePositionColumns == 0 ) + _synapsePositionColumns = positions->getNumAttributes(); + assert(_synapsePositionColumns == positions->getNumAttributes()); + + return *positions; } void saveMorphologiesToCache( const std::string& uri, @@ -354,19 +361,21 @@ class Circuit::Impl std::vector< Future > futures; futures.reserve( keys.size( )); - _cache->takeValues( keys, [&futures] ( const std::string& key, - char* data, const size_t size ) + _findSynapsePositionsColumns(); + + _cache->takeValues( keys, [this, &futures]( const std::string& key, + char* data, const size_t size ) { - futures.push_back( std::async( [key, data, size] + futures.push_back( std::async( [this, key, data, size] { // there is no constructor in multi_array which just accepts the // size in bytes (although there's a getter for it used in // saveSynapsePositionsToCache()), so we reconstruct the row and // column count here. - const size_t numColumns = brion::SYNAPSE_POSITION_ALL; + const size_t numColumns = _synapsePositionColumns; const size_t numRows = size / sizeof(float) / numColumns; - brion::SynapseMatrix values( boost::extents[numRows] - [numColumns]); + brion::SynapseMatrix values( + boost::extents[numRows][numColumns]); ::memmove( values.data(), data, size ); std::free( data ); return std::make_pair( key, values ); @@ -382,6 +391,17 @@ class Circuit::Impl return loaded; } + void _findSynapsePositionsColumns() const + { + lunchbox::ScopedWrite mutex( _synapsePositions[0] ); + if( !(*_synapsePositions[0] )) + _synapsePositions[0]->reset( + new brion::Synapse( _synapseSource.getPath() + + afferentPositionsFilename )); + _synapsePositionColumns = (*_synapsePositions[0])->getNumAttributes(); + + } + const brion::URI _circuitSource; const brion::URI _morphologySource; const brion::URI _synapseSource; @@ -396,6 +416,7 @@ class Circuit::Impl mutable LockPtr< brion::Synapse > _synapseAttributes[2]; mutable LockPtr< brion::Synapse > _synapseExtra; mutable LockPtr< brion::Synapse > _synapsePositions[2]; + mutable size_t _synapsePositionColumns; }; class MVD2 : public Circuit::Impl diff --git a/brain/python/synapses.cpp b/brain/python/synapses.cpp index 6bdbb33..5a9d821 100644 --- a/brain/python/synapses.cpp +++ b/brain/python/synapses.cpp @@ -54,8 +54,9 @@ SynapseWrapper Synapses_get( const SynapsesWrapper& synapses, long int index ) } #define GET_SYNAPSES_ARRAY_PROPERTY( type, name ) \ -bp::object Synapses_##name( const SynapsesWrapper& synapses ) \ - { return toNumpy( synapses.name(), synapses.size(), synapses._impl ); } +bp::object Synapses_##name( const SynapsesWrapper& synapses ) \ + { if( !synapses.name() ) return bp::object(); \ + return toNumpy( synapses.name(), synapses.size(), synapses._impl ); } GET_SYNAPSES_ARRAY_PROPERTY( size_t, indices ) GET_SYNAPSES_ARRAY_PROPERTY( uin32_t, preGIDs ) diff --git a/brain/synapse.cpp b/brain/synapse.cpp index e993845..f07ad6e 100644 --- a/brain/synapse.cpp +++ b/brain/synapse.cpp @@ -64,6 +64,9 @@ float Synapse::getPresynapticDistance() const Vector3f Synapse::getPresynapticSurfacePosition() const { + if (!_synapses.preSurfaceXPositions()) + throw std::runtime_error( "Surface synapse positions not available" ); + return Vector3f( _synapses.preSurfaceXPositions()[_index], _synapses.preSurfaceYPositions()[_index], _synapses.preSurfaceZPositions()[_index] ); @@ -98,6 +101,9 @@ float Synapse::getPostsynapticDistance() const Vector3f Synapse::getPostsynapticSurfacePosition() const { + if (!_synapses.postSurfaceXPositions()) + throw std::runtime_error( "Surface synapse positions not available" ); + return Vector3f( _synapses.postSurfaceXPositions()[_index], _synapses.postSurfaceYPositions()[_index], _synapses.postSurfaceZPositions()[_index] ); diff --git a/brain/synapses.cpp b/brain/synapses.cpp index 16fb64d..c859119 100644 --- a/brain/synapses.cpp +++ b/brain/synapses.cpp @@ -181,7 +181,7 @@ struct Synapses::Impl : public Synapses::BaseImpl void _loadPositions( const GIDSet& gids, const GIDSet& filterGIDs ) const { - if( _preSurfacePositionX ) + if( _preCenterPositionX ) return; Strings hashes; @@ -197,7 +197,7 @@ struct Synapses::Impl : public Synapses::BaseImpl hashes.push_back( gidHash ); } - CachedSynapses loaded = + const CachedSynapses loaded = _circuit._impl->loadSynapsePositionsFromCache( hashes ); const bool haveSize = _size > 0; @@ -205,23 +205,23 @@ struct Synapses::Impl : public Synapses::BaseImpl // delay the opening of the synapse file as much as possible, even // though the code looks ugly... As the circuit impl keeps the file // opened, we can safely just get a loose pointer here. - const brion::Synapse* synapsePositions = nullptr; + const brion::Synapse* positions = nullptr; if( !haveSize ) { auto hash = hashes.begin(); for( const auto gid : gids ) { - auto it = loaded.find( *hash ); + const auto it = loaded.find( *hash ); ++hash; if( it != loaded.end() ) _size += it->second.shape()[0]; else { - if( !synapsePositions ) - synapsePositions = + if( !positions ) + positions = &_circuit._impl->getSynapsePositions( _afferent ); - _size += synapsePositions->getNumSynapses( GIDSet{ gid } ); + _size += positions->getNumSynapses( GIDSet{ gid } ); } } } @@ -230,6 +230,7 @@ struct Synapses::Impl : public Synapses::BaseImpl size_t i = 0; auto hash = hashes.begin(); + bool haveSurfacePositions = false; for( const auto gid : gids ) { auto it = loaded.find( *hash ); @@ -237,10 +238,14 @@ struct Synapses::Impl : public Synapses::BaseImpl const auto readFromFile = [&] { - if( !synapsePositions ) - synapsePositions = - &_circuit._impl->getSynapsePositions( _afferent ); - return synapsePositions->read( gid, brion::SYNAPSE_POSITION ); + if( !positions ) + positions = + &_circuit._impl->getSynapsePositions( _afferent ); + if( positions->getNumAttributes() == + brion::SYNAPSE_POSITION_ALL ) + return positions->read( gid, brion::SYNAPSE_POSITION ); + else + return positions->read( gid, brion::SYNAPSE_OLD_POSITION ); }; const brion::SynapseMatrix pos = cached ? it->second @@ -250,7 +255,7 @@ struct Synapses::Impl : public Synapses::BaseImpl _circuit._impl->saveSynapsePositionsToCache( gid, *hash, pos ); ++hash; - for( size_t j = 0; j < pos.shape()[0]; ++j ) + for( size_t j = 0; j < pos.size(); ++j ) { const uint32_t preGid = pos[j][0]; FILTER( preGid ); @@ -261,22 +266,45 @@ struct Synapses::Impl : public Synapses::BaseImpl _postGID.get()[i] = gid; } - _preSurfacePositionX.get()[i] = pos[j][1]; - _preSurfacePositionY.get()[i] = pos[j][2]; - _preSurfacePositionZ.get()[i] = pos[j][3]; - _postSurfacePositionX.get()[i] = pos[j][4]; - _postSurfacePositionY.get()[i] = pos[j][5]; - _postSurfacePositionZ.get()[i] = pos[j][6]; - _preCenterPositionX.get()[i] = pos[j][7]; - _preCenterPositionY.get()[i] = pos[j][8]; - _preCenterPositionZ.get()[i] = pos[j][9]; - _postCenterPositionX.get()[i] = pos[j][10]; - _postCenterPositionY.get()[i] = pos[j][11]; - _postCenterPositionZ.get()[i] = pos[j][12]; + if( pos.shape()[1] == brion::SYNAPSE_POSITION_ALL ) + { + haveSurfacePositions = true; + _preSurfacePositionX.get()[i] = pos[j][1]; + _preSurfacePositionY.get()[i] = pos[j][2]; + _preSurfacePositionZ.get()[i] = pos[j][3]; + _postSurfacePositionX.get()[i] = pos[j][4]; + _postSurfacePositionY.get()[i] = pos[j][5]; + _postSurfacePositionZ.get()[i] = pos[j][6]; + _preCenterPositionX.get()[i] = pos[j][7]; + _preCenterPositionY.get()[i] = pos[j][8]; + _preCenterPositionZ.get()[i] = pos[j][9]; + _postCenterPositionX.get()[i] = pos[j][10]; + _postCenterPositionY.get()[i] = pos[j][11]; + _postCenterPositionZ.get()[i] = pos[j][12]; + } + else + { + _preCenterPositionX.get()[i] = pos[j][1]; + _preCenterPositionY.get()[i] = pos[j][2]; + _preCenterPositionZ.get()[i] = pos[j][3]; + _postCenterPositionX.get()[i] = pos[j][4]; + _postCenterPositionY.get()[i] = pos[j][5]; + _postCenterPositionZ.get()[i] = pos[j][6]; + } ++i; } } + if( !haveSurfacePositions) + { + _preSurfacePositionX.reset(); + _preSurfacePositionY.reset(); + _preSurfacePositionZ.reset(); + _postSurfacePositionX.reset(); + _postSurfacePositionY.reset(); + _postSurfacePositionZ.reset(); + } + if( !haveSize ) { if( !_afferent ) @@ -356,7 +384,7 @@ struct Synapses::Impl : public Synapses::BaseImpl bool _hasPositions() const { lunchbox::ScopedRead mutex( _lock ); - return _preSurfacePositionX.get() != nullptr; + return _preCenterPositionX.get() != nullptr; } const Circuit& _circuit; diff --git a/brain/synapses.h b/brain/synapses.h index 9174a56..d5e69ee 100644 --- a/brain/synapses.h +++ b/brain/synapses.h @@ -101,19 +101,19 @@ class Synapses /** * @return the presynaptic touch position x-coordinates on the surfaces of - * the segments. + * the segments. May be null in old circuits. */ BRAIN_API const float* preSurfaceXPositions() const; /** * @return the presynaptic touch position y-coordinates on the surfaces of - * the segments. + * the segments. May be null in old circuits. */ BRAIN_API const float* preSurfaceYPositions() const; /** * @return the presynaptic touch position z-coordinates on the surfaces of - * the segments. + * the segments. May be null in old circuits. */ BRAIN_API const float* preSurfaceZPositions() const; @@ -152,19 +152,19 @@ class Synapses /** * @return the postsynaptic touch position x-coordinates on the surfaces of - * the segments. + * the segments. May be null in old circuits. */ BRAIN_API const float* postSurfaceXPositions() const; /** * @return the postsynaptic touch position x-coordinates on the surfaces of - * the segments. + * the segments. May be null in old circuits. */ BRAIN_API const float* postSurfaceYPositions() const; /** * @return the postsynaptic touch position x-coordinates on the surfaces of - * the segments. + * the segments. May be null in old circuits. */ BRAIN_API const float* postSurfaceZPositions() const; diff --git a/brion/enums.h b/brion/enums.h index f29e735..ddf4603 100644 --- a/brion/enums.h +++ b/brion/enums.h @@ -244,6 +244,29 @@ enum SynapsePositions SYNAPSE_POSTSYNAPTIC_POSITION }; +enum SynapseOldPositions +{ + SYNAPSE_OLD_PRESYNAPTIC_CENTER_X = 1 << 1, + SYNAPSE_OLD_PRESYNAPTIC_CENTER_Y = 1 << 2, + SYNAPSE_OLD_PRESYNAPTIC_CENTER_Z = 1 << 3, + SYNAPSE_OLD_POSTSYNAPTIC_CENTER_X = 1 << 4, + SYNAPSE_OLD_POSTSYNAPTIC_CENTER_Y = 1 << 5, + SYNAPSE_OLD_POSTSYNAPTIC_CENTER_Z = 1 << 6, + SYNAPSE_OLD_POSITION_ALL = 7, + + SYNAPSE_OLD_PRESYNAPTIC_POSITION = SYNAPSE_OLD_PRESYNAPTIC_CENTER_X | + SYNAPSE_OLD_PRESYNAPTIC_CENTER_Y | + SYNAPSE_OLD_PRESYNAPTIC_CENTER_Z, + + SYNAPSE_OLD_POSTSYNAPTIC_POSITION = SYNAPSE_OLD_POSTSYNAPTIC_CENTER_X | + SYNAPSE_OLD_POSTSYNAPTIC_CENTER_Y | + SYNAPSE_OLD_POSTSYNAPTIC_CENTER_Z, + + SYNAPSE_OLD_POSITION = SYNAPSE_CONNECTED_NEURON | + SYNAPSE_OLD_PRESYNAPTIC_POSITION | + SYNAPSE_OLD_POSTSYNAPTIC_POSITION +}; + /** * Specify the access mode of data. * @version 1.4 diff --git a/brion/synapse.cpp b/brion/synapse.cpp index fdcefa1..13fdde1 100644 --- a/brion/synapse.cpp +++ b/brion/synapse.cpp @@ -44,12 +44,48 @@ namespace brion namespace detail { +namespace +{ struct Dataset { H5::DataSet dataset; H5::DataSpace dataspace; hsize_t dims[2]; }; + +bool _openDataset( const H5::H5File& file, const std::string& name, + Dataset& dataset ) +{ + try + { + SilenceHDF5 silence; + dataset.dataset = file.openDataSet( name ); + } + catch( const H5::Exception& ) + { + LBVERB << "Could not find synapse dataset for " << name << ": " + << std::endl; + return false; + } + + dataset.dataspace = dataset.dataset.getSpace(); + if( dataset.dataspace.getSimpleExtentNdims() != 2 ) + { + LBERROR << "Synapse dataset is not 2 dimensional" << std::endl; + return false; + } + + if( dataset.dataspace.getSimpleExtentDims( dataset.dims ) < 0 ) + { + LBERROR << "Synapse dataset dimensions could not be retrieved" + << std::endl; + return false; + } + + return true; +} + +}; namespace fs = boost::filesystem; using boost::lexical_cast; @@ -74,13 +110,14 @@ class SynapseFile : public boost::noncopyable Dataset dataset; const std::string& datasetName = _file.getObjnameByIdx( 0 ); - if( !_openDataset( datasetName, dataset )) + if( !detail::_openDataset( _file, datasetName, dataset )) LBTHROW( std::runtime_error( "Cannot open dataset in synapse file " + source )); _numAttributes = dataset.dims[1]; if( _numAttributes != SYNAPSE_ALL && _numAttributes != SYNAPSE_POSITION_ALL && + _numAttributes != SYNAPSE_OLD_POSITION_ALL && _numAttributes != 1 /* nrn_extra */) { LBTHROW( std::runtime_error( source + " not a valid synapse file")); @@ -127,6 +164,11 @@ class SynapseFile : public boost::noncopyable return values; } + size_t getNumAttributes() const + { + return _numAttributes; + } + size_t getNumSynapses( const GIDSet& gids ) const { lunchbox::ScopedWrite mutex( detail::_hdf5Lock ); @@ -146,38 +188,7 @@ class SynapseFile : public boost::noncopyable { std::stringstream name; name << "a" << gid; - return _openDataset( name.str(), dataset ); - } - - bool _openDataset( const std::string& name, Dataset& dataset ) const - { - try - { - SilenceHDF5 silence; - dataset.dataset = _file.openDataSet( name ); - } - catch( const H5::Exception& ) - { - LBVERB << "Could not find synapse dataset for " << name << ": " - << std::endl; - return false; - } - - dataset.dataspace = dataset.dataset.getSpace(); - if( dataset.dataspace.getSimpleExtentNdims() != 2 ) - { - LBERROR << "Synapse dataset is not 2 dimensional" << std::endl; - return false; - } - - if( dataset.dataspace.getSimpleExtentDims( dataset.dims ) < 0 ) - { - LBERROR << "Synapse dataset dimensions could not be retrieved" - << std::endl; - return false; - } - - return true; + return detail::_openDataset( _file, name.str(), dataset ); } SynapseMatrix read( const uint32_t gid, const uint32_t attributes ) const @@ -188,6 +199,8 @@ class SynapseFile : public boost::noncopyable return read< SYNAPSE_ALL >( gid, attributes ); case SYNAPSE_POSITION_ALL: return read< SYNAPSE_POSITION_ALL >( gid, attributes ); + case SYNAPSE_OLD_POSITION_ALL: + return read< SYNAPSE_OLD_POSITION_ALL >( gid, attributes ); case 1: // nrn_extra return read< 1 >( gid, 1 ); @@ -225,7 +238,7 @@ class Synapse : public boost::noncopyable const std::string filename = sourcePath.filename().generic_string(); _createIndex( dir, filename ); - _fillFilemap( dir, filename ); + _findCandidateFiles( dir, filename ); } } @@ -256,13 +269,39 @@ class Synapse : public boost::noncopyable return numSynapses; } + size_t getNumAttributes() const + { + if( _file ) + return _file->getNumAttributes(); + + assert( !_fileNames.empty( )); + + const std::string& filename = _fileNames.front(); + H5::H5File file; + try + { + SilenceHDF5 silence; + file.openFile( filename, H5F_ACC_RDONLY ); + } + catch( const H5::Exception& exc ) + { + LBTHROW( std::runtime_error( + "Could not open synapse file " + filename + ": " + + exc.getDetailMsg( ))); + } + Dataset dataset; + if( !_openDataset( file, file.getObjnameByIdx( 0 ), dataset )) + LBTHROW( std::runtime_error( "Cannot open dataset in synapse file " + + filename )); + return dataset.dims[1]; + } + private: typedef boost::unordered_map< uint32_t, std::string > GidFileMap; - typedef boost::unordered_map< std::string, H5::H5File > UnmappedFiles; mutable SynapseFile* _file; mutable uint32_t _gid; // current or 0 for all - mutable UnmappedFiles _unmappedFiles; + Strings _fileNames; mutable GidFileMap _fileMap; void _createIndex( const fs::path& dir, const std::string& filename ) @@ -274,8 +313,8 @@ class Synapse : public boost::noncopyable const std::ifstream mergeFile( merge_nrn.generic_string().c_str( )); if( !mergeFile.is_open( )) { - LBWARN << "No merge file found in " << dir - << " to build lookup index; loading data will be slow" + LBWARN << "No merged file found in " << dir + << " to build lookup index; loading data will be very slow" << std::endl; return; } @@ -303,7 +342,7 @@ class Synapse : public boost::noncopyable } } - void _fillFilemap( const fs::path& dir, const std::string& filename ) + void _findCandidateFiles( const fs::path& dir, const std::string& filename ) { // try to open in individual files const boost::regex filter( filename + "\\.[0-9]+$" ); @@ -313,17 +352,16 @@ class Synapse : public boost::noncopyable const fs::path candidate = i->path().filename(); boost::smatch match; - if( !boost::filesystem::is_regular_file( i->status( )) || - !boost::regex_match( candidate.string(), match, filter )) + // Testing first if the regex matches. We don't want to throw on + // is_regular_file for directory entries that are not valid anyway. + if( boost::regex_match( candidate.string(), match, filter ) && + boost::filesystem::is_regular_file( i->status( ))) { - continue; + _fileNames.push_back( i->path().string()); } - - _unmappedFiles.insert( std::make_pair( i->path().string(), - H5::H5File( ))); } - if( _unmappedFiles.empty( )) + if( _fileNames.empty( )) { LBTHROW( std::runtime_error( "Could not find synapse files " + dir.string() + "/" + filename )); @@ -357,20 +395,17 @@ class Synapse : public boost::noncopyable lunchbox::ScopedWrite mutex( detail::_hdf5Lock ); SilenceHDF5 silence; - BOOST_FOREACH( UnmappedFiles::value_type& entry, _unmappedFiles ) + // this trial-and-error is the 'fastest' path found + BOOST_FOREACH( const std::string& candidate, _fileNames ) { - const std::string& candidate = entry.first; - H5::H5File& file = entry.second; - // keeping the files open 'saves' some time - if( file.getId() <= 0 ) - file.openFile( candidate, H5F_ACC_RDONLY ); - try { + H5::H5File file; + file.openFile( candidate, H5F_ACC_RDONLY ); + std::stringstream name; name << "a" << gid; - // this trial-and-error is the 'fastest' path found file.openDataSet( name.str( )); _fileMap[ gid ] = candidate; return candidate; @@ -404,4 +439,9 @@ size_t Synapse::getNumSynapses( const GIDSet& gids ) const return _impl->getNumSynapses( gids ); } +size_t Synapse::getNumAttributes() const +{ + return _impl->getNumAttributes(); +} + } diff --git a/brion/synapse.h b/brion/synapse.h index 79d1522..e668dcc 100644 --- a/brion/synapse.h +++ b/brion/synapse.h @@ -87,6 +87,11 @@ class Synapse : public boost::noncopyable * @version 1.0 */ BRION_API size_t getNumSynapses( const GIDSet& gids ) const; + + /** Return the number of columns of the synapse file. + * @version 1.9 + */ + BRION_API size_t getNumAttributes() const; //@} private: diff --git a/doc/Changelog.md b/doc/Changelog.md index 25ac381..fa6d03c 100644 --- a/doc/Changelog.md +++ b/doc/Changelog.md @@ -3,6 +3,8 @@ Changelog {#Changelog} # git master +* [113](https://github.com/BlueBrain/Brion/pull/113): + Support for old circuits containing only synapse center positions. * [110](https://github.com/BlueBrain/Brion/pull/110): Break PersistentMap out into keyv::Map * [107](https://github.com/BlueBrain/Brion/pull/107):