Skip to content

Commit

Permalink
Add support for H5 format in output (#14)
Browse files Browse the repository at this point in the history 
This commit adds the option to use h5 as output format. Note that
arbitrary precision is not supported by this (optional) output.
  • Loading branch information
@jbreue16
jbreue16 committed Nov 24, 2024
1 parent d24f0a0 commit d586b69
Showing 1 changed file with 205 additions and 11 deletions.
216 changes: 205 additions & 11 deletions src/chromatogram.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,11 +92,182 @@ class ProgressBarUpdater
int _prevPhase;
};

#ifdef ENABLE_HDF5
void writeMetaAndResultToH5(const ProgramOptions& opts, const casema::model::ModelSystem* sys, const std::vector<mpfr::mpreal>& solutionTimes, const casema::util::SlicedVector<mpfr::mpreal>& output, const double time_sim, const int timeOffset)
{
casema::io::IFileWriter* writer = casema::io::createWriter("h5");

bool split_ports = true; // from CADET file format
bool single_as_multi_port = false; // from CADET file format

if(opts.outFile.empty()) // write to input file
{
casema::io::IFileReader* reader = casema::io::createReader("h5");
reader->openFile(opts.inFile, "r");
reader->pushGroup("input");
reader->pushGroup("return");
if (reader->exists("SPLIT_PORTS_DATA"))
split_ports = reader->getBool("SPLIT_PORTS_DATA");
if (reader->exists("SINGLE_AS_MULTI_PORT"))
single_as_multi_port = reader->getBool("SINGLE_AS_MULTI_PORT");

reader->closeFile();
writer->openFile(opts.inFile, "rw");
}
else if (opts.outFile == opts.inFile)
{
casema::io::IFileReader* reader = casema::io::createReader("h5");
reader->openFile(opts.outFile, "r");
reader->pushGroup("input");
reader->pushGroup("return");
if (reader->exists("SPLIT_PORTS_DATA"))
split_ports = reader->getBool("SPLIT_PORTS_DATA");
if (reader->exists("SINGLE_AS_MULTI_PORT"))
single_as_multi_port = reader->getBool("SINGLE_AS_MULTI_PORT");

reader->closeFile();
writer->openFile(opts.outFile, "rw");
}
else // create output file
writer->openFile(opts.outFile, "co");

// Write meta data
writer->deleteGroup("meta");
writer->pushGroup("meta");
writer->writeString("CASEMA_VERSION", casema::getVersion());
writer->writeString("CASEMA_BRANCH", casema::getBranchRefspec());
writer->writeString("CASEMA_COMMIT", casema::getCommitHash());
writer->writeInt("FILE_FORMAT", casema::getFileFormat());
writer->writeDouble("TIME_SIM", time_sim);
{
std::ostringstream oss;
oss << opts.errorTrunc + opts.errorCons;
writer->writeString("ERROR", oss.str());
}
writer->popGroup();

// Write program options
writer->deleteGroup("program_options");
writer->pushGroup("program_options");
writer->writeBool("KAHAN_SUMMATION", opts.kahan);
writer->writeBool("IGNORE_CSTR", opts.ignoreCSTR);
#ifdef ENABLE_BESSEL_TRUNCATION
writer->writeBool("BESSEL_FUNCTION_TRUNCATION", true);
#else
writer->writeBool("BESSEL_FUNCTION_TRUNCATION", false);
#endif
{
std::ostringstream oss;
oss << mpfr::digits2bits(opts.precision);
writer->writeString("PRECISION", oss.str());
oss.str("");
//writer->writeDouble("OUTPUT_PRECISION", std::numeric_limits<double>::digits10 + 1);

writer->writeBool("EXTRAPOLATION", false);
writer->writeInt("MAX_LAPLACE_SUMMANDS", opts.laplaceSummands);
writer->writeInt("MAX_HANKEL_SUMMANDS", opts.hankelSummands);
oss << opts.abscissa;
writer->writeString("ABSCISSA", oss.str());
oss.str("");
oss << opts.errorCons;
writer->writeString("CONSISTENCY_ERROR", oss.str());
oss.str("");
oss << opts.errorTrunc;
writer->writeString("TRUNCATION_ERROR", oss.str());
oss.str("");
oss << opts.errorTrunc + opts.errorCons;
writer->writeString("ERROR", oss.str());
oss.str("");
oss << opts.error;
writer->writeString("REQUESTED_ERROR", oss.str());
oss.str("");
oss << opts.errorWeight;
writer->writeString("ERROR_WEIGHT", oss.str());
}
writer->popGroup();

writer->deleteGroup("output");
writer->pushGroup("output");

writer->pushGroup("solution");

//writer->writeVectorDouble("SOLUTION_TIMES", length, *buffer, stride, blocksize);
{
const int timePoints = solutionTimes.size() - timeOffset;
std::vector<double> solutionBuffer(timePoints);

const int numUnits = sys->numModels();
const int numOutlets = sys->numOutlets();
int numWrittenOutlets = 0;

if (!split_ports) // resize if ports data should be written together
{
int maxNumPorts = 1;
for (int j = 0; j < numUnits; ++j)
maxNumPorts = std::max(maxNumPorts, sys->unitOperation(j)->numOutletPorts());

solutionBuffer.push_back(maxNumPorts * timePoints);
}

for (int i = 0; i < timePoints; i++)
solutionBuffer[i] = static_cast<double>(solutionTimes[i + timeOffset]);

writer->writeVectorDouble("SOLUTION_TIMES", timePoints, &solutionBuffer[0], 1, timePoints);

std::ostringstream oss;

for (int j = 0; j < numUnits; ++j)
{
oss << std::setw(3) << std::setfill('0') << j;
writer->pushGroup("unit_" + oss.str());
oss.str("");

void writeMeta(std::ostream& os, const casema::model::ModelSystem& model, const ProgramOptions& opts)
const int nPorts = std::max(1, sys->unitOperation(j)->numOutletPorts());

if (split_ports) // write vectors for each port
{
for (int k = 0; k < nPorts; ++k)
{
for (int i = 0; i < timePoints; i++)
solutionBuffer[i] = static_cast<double>(output(k + numWrittenOutlets, i));

oss << std::setw(3) << std::setfill('0') << k;
const std::string outName = single_as_multi_port ? "SOLUTION_OUTLET_PORT_" + oss.str() : "SOLUTION_OUTLET";
writer->writeVectorDouble(outName, timePoints, &solutionBuffer[0], 1, timePoints);
oss.str("");
numWrittenOutlets++;
}
}
else
{
const int nPorts = std::max(1, sys->unitOperation(j)->numOutletPorts());

for (int i = 0; i < timePoints; i++)
{
for (int k = 0; k < nPorts; ++k)
solutionBuffer[i * nPorts + k] = static_cast<double>(output(k + numWrittenOutlets, i));
}

writer->writeMatrixDouble("SOLUTION_OUTLET", timePoints, nPorts, &solutionBuffer[0], 1, nPorts * timePoints);
numWrittenOutlets += nPorts;
}
writer->popGroup();
}
}

writer->closeFile();
}
#endif

void writeMeta(std::ostream& os, const casema::model::ModelSystem& model, const ProgramOptions& opts, const double simTime)
{
os << "# Version: " << casema::getVersion() << "\n";
os << "# Commit: " << casema::getCommitHash() << "\n";
os << "# CASEMA Version: " << casema::getVersion() << "\n";
os << "# CASEMA Branch: " << casema::getBranchRefspec() << "\n";
os << "# CASEMA Commit: " << casema::getCommitHash() << "\n";
//os << "# CASEMA File Format: " << casema::getFileFormat() << "\n";
if(simTime < 0.0)
os << "# Compute time: " << simTime << "\n";

if (opts.kahan)
os << "# Kahan summation: On\n";
else
Expand Down Expand Up @@ -433,23 +604,46 @@ int main(int argc, char** argv)
std::cout << "[DONE]" << std::endl;
}

writeMeta(std::cout, *sys, opts);
writeMeta(std::cout, *sys, opts, -1.0);

std::cout << "Starting Laplace inversion" << std::endl;

const auto start = std::chrono::high_resolution_clock::now();

const casema::util::SlicedVector<mpfr::mpreal> output = invert(*sys, maxTime, st.solutionTimes, timeOffset, opts);

const double time_sim = std::chrono::duration_cast<std::chrono::duration<double>>(std::chrono::high_resolution_clock::now() - start).count();

std::cout << "Writing output" << std::endl;
if (opts.outFile.empty())

// Check if input file is h5 file and if so, write to that. Otherwise, print to console
bool inIsH5 = false;
const std::size_t dotPos = opts.inFile.find_last_of('.');

#ifdef ENABLE_HDF5
if (dotPos != std::string::npos && dotPos != 0)
{
writeResult(std::cout, st.solutionTimes, output, opts.outPrecision, timeOffset, *sys);
if ((opts.inFile.substr(dotPos + 1) == "h5"))
{
inIsH5 = true;

writeMetaAndResultToH5(opts, sys, st.solutionTimes, output, time_sim, timeOffset);
}
}
else
if (!inIsH5)
#endif
{
std::ofstream fs(opts.outFile, std::ofstream::out | std::ofstream::trunc);

writeMeta(fs, *sys, opts);
writeResult(fs, st.solutionTimes, output, opts.outPrecision, timeOffset, *sys);
if (opts.outFile.empty())
{
writeResult(std::cout, st.solutionTimes, output, opts.outPrecision, timeOffset, *sys);
}
else
{
std::ofstream fs(opts.outFile, std::ofstream::out | std::ofstream::trunc);

writeMeta(fs, *sys, opts, time_sim);
writeResult(fs, st.solutionTimes, output, opts.outPrecision, timeOffset, *sys);
}
}

::mpfr_free_cache();
Expand Down

0 comments on commit d586b69

Please sign in to comment.