diff --git a/brain/detail/compartmentReport.h b/brain/detail/compartmentReport.h index e639d84..6b0b2fd 100644 --- a/brain/detail/compartmentReport.h +++ b/brain/detail/compartmentReport.h @@ -72,34 +72,34 @@ struct CompartmentReportView void CompartmentReportView::_initIndices() { size_t indicesCount = 0; - size_t index = 0; - auto gidCount = report->getGIDs().size(); + const auto& gids = report->getGIDs(); - while (gidCount--) - indicesCount += report->getOffsets()[index++].size(); + const auto gidCount = gids.size(); + std::vector indexPositions(gidCount); + for (size_t i = 0; i < gids.size(); ++i) + { + indexPositions[i] = indicesCount; + indicesCount += report->getOffsets()[i].size(); + } - indices.reserve(indicesCount); - index = 0; + indices.resize(indicesCount); - for (auto gid : report->getGIDs()) + const std::vector gidList(report->getGIDs().begin(), + report->getGIDs().end()); + for (size_t i = 0; i < gids.size(); ++i) { const brion::uint16_ts& compartments = - report->getCompartmentCounts()[index]; - const brion::uint64_ts& offsets = report->getOffsets()[index]; + report->getCompartmentCounts()[i]; + const brion::uint64_ts& offsets = report->getOffsets()[i]; uint16_t section = 0; + auto pos = indexPositions[i]; for (auto offset : offsets) { - indices.push_back({offset, gid, section, compartments[section]}); + indices[pos + section] = {offset, gidList[i], section, + compartments[section]}; ++section; } - ++index; } } - -struct CompartmentReportFrame -{ - double timeStamp = -1.0; - brion::floats data; -}; } } // namespaces diff --git a/brion/plugin/.#-emacsD9s7NS b/brion/plugin/.#-emacsD9s7NS new file mode 100644 index 0000000..e69de29 diff --git a/brion/plugin/compartmentReportBinary.cpp b/brion/plugin/compartmentReportBinary.cpp index 0437bbc..f6c1f14 100644 --- a/brion/plugin/compartmentReportBinary.cpp +++ b/brion/plugin/compartmentReportBinary.cpp @@ -42,6 +42,8 @@ namespace { +const size_t _mappingItemSize = 4; + // Offsets of the header information in the file. enum HeaderPositions { @@ -104,6 +106,7 @@ struct CellInfo { int32_t gid; int32_t numCompartments; + uint64_t accumCompartments; uint64_t mappingOffset; uint64_t extraMappingOffset; uint64_t dataOffset; @@ -637,6 +640,25 @@ bool CompartmentReportBinary::_parseMapping() size_t offset = _header.headerSize; + // Prefetching the mapping data + uint64_t dataOffset = get(ptr, DATA_INFO + offset); + if (_header.byteswap) + lunchbox::byteswap(dataOffset); + _header.dataBlockOffset = dataOffset; + + std::unique_ptr buffer(new uint8_t[_header.dataBlockOffset]); + std::copy(ptr, ptr + dataOffset, buffer.get()); + ptr = buffer.get(); + + // According to Garik Suess, all compartments of a cell in a frame are next + // to each other, and all compartments of a section are next to each other. + // However, the sections are not necessarily sorted by their index in the + // frame. It could be that for cell with GID 50 while all data is + // contiguous, the sections are out of order so compartments for section 3 6 + // 22 8 could be next to each other, while the compartments inside these + // sections are in order. + + uint64_t totalCompartments = 0; CellInfos cells(_header.numCells); for (int32_t i = 0; i < _header.numCells; ++i) { @@ -649,10 +671,19 @@ bool CompartmentReportBinary::_parseMapping() cell.extraMappingOffset = get(ptr, EXTRA_MAPPING_INFO + offset); cell.dataOffset = get(ptr, DATA_INFO + offset); + if (cell.dataOffset < dataOffset) + { + LBERROR << "Bad offset in report mapping" << std::endl; + return false; + } offset += SIZE_CELL_INFO_LENGTH; if (_header.byteswap) lunchbox::byteswap(cell); + + cell.accumCompartments = totalCompartments; + totalCompartments += cell.numCompartments; + _originalGIDs.insert(cell.gid); } _header.dataBlockOffset = cells[0].dataOffset; @@ -666,76 +697,70 @@ bool CompartmentReportBinary::_parseMapping() perCellOffsets.resize(cells.size()); _perCellCounts.resize(cells.size()); - // According to Garik Suess, all compartments of a cell in a frame are next - // to each other, and all compartments of a section are next to each other. - // However, the sections are not necessarily sorted by their index in the - // frame. It could be that for cell with GID 50 while all data is - // contiguous, the sections are out of order so compartments for section 3 6 - // 22 8 could be next to each other, while the compartments inside these - // sections are in order. - size_t idx = 0; - size_t totalCompartments = 0; - for (const CellInfo& info : cells) +// With 8 threads there's some speed up, but very small. We stick with 6 +// to avoid using more than 1 socket. +#pragma omp parallel for num_threads(6) schedule(dynamic) + for (size_t idx = 0; idx < cells.size(); ++idx) { + const CellInfo& info = cells[idx]; uint16_t current = LB_UNDEFINED_UINT16; uint16_t previous = LB_UNDEFINED_UINT16; uint16_t count = 0; // < sectionID, < frameIndex, numCompartments > > - typedef std::map> - SectionMapping; - SectionMapping sectionsMapping; + std::vector>> + sectionsMapping; + sectionsMapping.reserve(info.numCompartments); _perCellCounts[idx] = info.numCompartments; - perCellOffsets[idx] = totalCompartments; - totalCompartments += info.numCompartments; + perCellOffsets[idx] = info.accumCompartments; for (int32_t j = 0; j < info.numCompartments; ++j) { previous = current; const size_t pos(info.mappingOffset + - j * sizeof(float) * _header.mappingSize); + j * _mappingItemSize * _header.mappingSize); float value = get(ptr, pos); if (_header.byteswap) lunchbox::byteswap(value); + assert(value < 65536); current = value; // in case this is the start of a new section if (current != previous) { const uint64_t frameIndex = - j + ((info.dataOffset - _header.dataBlockOffset) / - sizeof(float)); + j + ((info.dataOffset - dataOffset) / _mappingItemSize); - sectionsMapping.insert( + sectionsMapping.push_back( std::make_pair(current, std::make_pair(frameIndex, 0))); if (previous != LB_UNDEFINED_UINT16) - sectionsMapping[previous].second = count; + sectionsMapping[sectionsMapping.size() - 2].second.second = + count; count = 0; } ++count; } - sectionsMapping[current].second = count; + sectionsMapping.back().second.second = count; + std::sort(sectionsMapping.begin(), sectionsMapping.end()); // now convert the maps into the desired mapping format - uint64_ts& sectionOffsets = offsetMapping[idx]; - uint16_ts& sectionCompartmentCounts = perSectionCounts[idx]; - ++idx; + auto& sectionOffsets = offsetMapping[idx]; + auto& sectionCompartmentCounts = perSectionCounts[idx]; // get maximum section id const uint16_t maxID = sectionsMapping.rbegin()->first; sectionOffsets.resize(maxID + 1, LB_UNDEFINED_UINT64); sectionCompartmentCounts.resize(maxID + 1, 0); - for (auto k : sectionsMapping) + for (auto sectionInfo : sectionsMapping) { - sectionOffsets[k.first] = k.second.first; - sectionCompartmentCounts[k.first] = k.second.second; + const auto id = sectionInfo.first; + sectionOffsets[id] = sectionInfo.second.first; + sectionCompartmentCounts[id] = sectionInfo.second.second; } - - _originalGIDs.insert(info.gid); } return true; diff --git a/tests/compartmentReport.cpp b/tests/compartmentReport.cpp index 180fdba..1977b65 100644 --- a/tests/compartmentReport.cpp +++ b/tests/compartmentReport.cpp @@ -370,6 +370,19 @@ void testReadSoma(const char* relativePath) BOOST_CHECK_EQUAL(report.getTimestep(), 0.1); BOOST_CHECK_EQUAL(report.getFrameSize(), 1); + const auto& offsets = report.getOffsets(); + for (size_t i = 0; i != offsets.size(); ++i) + { + BOOST_CHECK_EQUAL(offsets[i].size(), 1); + BOOST_CHECK_EQUAL(offsets[i][0], i); + } + const auto& counts = report.getCompartmentCounts(); + for (auto&& c : counts) + { + BOOST_CHECK_EQUAL(c.size(), 1); + BOOST_CHECK_EQUAL(c[0], 1); + } + brion::floatsPtr frame = report.loadFrame(report.getStartTime()).get(); BOOST_CHECK(frame); BOOST_CHECK_EQUAL((*frame)[0], -65); @@ -401,6 +414,14 @@ void testReadAllCompartments(const char* relativePath) BOOST_CHECK_EQUAL(report.getTimestep(), 0.1); BOOST_CHECK_EQUAL(report.getFrameSize(), 20360); + const auto& offsets = report.getOffsets(); + BOOST_CHECK_EQUAL(offsets.size(), report.getGIDs().size()); + BOOST_CHECK_EQUAL(offsets[0][0], 0); + BOOST_CHECK_EQUAL(offsets[0][1], 1); + BOOST_CHECK_EQUAL(offsets[1][0], 629); + BOOST_CHECK_EQUAL(offsets[1][1], 630); + BOOST_CHECK_EQUAL(offsets.back().back(), 20359); + brion::floatsPtr frame = report.loadFrame(.8).get(); BOOST_CHECK(frame); BOOST_CHECK_CLOSE((*frame)[0], -65.2919388, .000001f);