use of cbit.vcell.solver.ode.ODESolverResultSet in project vcell by virtualcell.
the class BiomodelsDB_TestSuite method main.
public static void main(String[] args) {
try {
Logging.init();
/*sanity check -- we currently only have a copasi_java dll for 32-bit, so let's make sure
* we're running on the right JVM. (Note we can run this on 64 bit machine, just have to install
* a 32 bit JVM...)
*/
OperatingSystemInfo osi = OperatingSystemInfo.getInstance();
if (!osi.isWindows() || osi.is64bit()) {
System.err.println("run on 32 bit JVM");
System.exit(99);
}
// following are set in command line processing
SortedSet<BiomodelsNetEntry> modelIDs = new TreeSet<BiomodelsNetEntry>();
BioModelsWebServices service = null;
File outDir = null;
boolean isDetailed;
{
// scope for command line processing
Options commandOptions = new Options();
Option help = new Option("help", false, "show help");
commandOptions.addOption(help);
Option detailed = new Option("detailed", false, "record detailed information");
commandOptions.addOption(detailed);
Option output = new Option("output", true, "output directory");
output.setRequired(true);
commandOptions.addOption(output);
OptionGroup primary = new OptionGroup();
LOption min = new LOption("min", true, "minimum number of model to import", true);
min.setType(Number.class);
LOption only = new LOption("only", true, "run only this model", false);
only.setType(Number.class);
LOption include = new LOption("include", true, "run on models in specified file", false);
include.setType(String.class);
LOption exclude = new LOption("exclude", true, "exclude models in specified file", true);
exclude.setType(String.class);
primary.addOption(min).addOption(only).addOption(include).addOption(exclude);
for (Object obj : primary.getOptions()) {
// CLI old, non-generic API
commandOptions.addOption((Option) obj);
}
CommandLine cmdLine = null;
try {
CommandLineParser parser = new DefaultParser();
cmdLine = parser.parse(commandOptions, args);
} catch (ParseException e1) {
e1.printStackTrace();
HelpFormatter hf = new HelpFormatter();
hf.printHelp("BiomodelsDB_TestSuite", commandOptions);
System.exit(2);
}
Option[] present = cmdLine.getOptions();
Set<Option> optionSet = new HashSet<Option>(Arrays.asList(present));
if (optionSet.contains(help)) {
HelpFormatter hf = new HelpFormatter();
hf.printHelp("BiomodelsDB_TestSuite", commandOptions);
System.exit(0);
}
// placeholder, avoid messing with nulls
LOption primaryOpt = new LOption("", false, "", true);
@SuppressWarnings("unchecked") Collection<? extends Option> priOpts = primary.getOptions();
priOpts.retainAll(optionSet);
if (!priOpts.isEmpty()) {
assert (priOpts.size() == 1);
primaryOpt = (LOption) priOpts.iterator().next();
}
outDir = new File(cmdLine.getOptionValue(output.getOpt()));
if (!outDir.exists()) {
outDir.mkdirs();
}
isDetailed = optionSet.contains(detailed);
BioModelsWebServicesServiceLocator locator = new BioModelsWebServicesServiceLocator();
service = locator.getBioModelsWebServices();
if (primaryOpt.downloads) {
String[] modelStrings = service.getAllCuratedModelsId();
for (String s : modelStrings) {
modelIDs.add(new BiomodelsNetEntry(s));
}
}
if (primaryOpt == only) {
Long only1 = (Long) cmdLine.getParsedOptionValue(only.getOpt());
modelIDs.add(new BiomodelsNetEntry(only1.intValue()));
} else if (primaryOpt == min) {
Long lv = (Long) cmdLine.getParsedOptionValue(min.getOpt());
int minimumModel = lv.intValue();
for (int m = 0; m < minimumModel; m++) {
modelIDs.remove(new BiomodelsNetEntry(m));
}
} else if (primaryOpt == include) {
FileBaseFilter fbf = new FileBaseFilter(cmdLine.getOptionValue(include.getOpt()));
modelIDs.addAll(fbf.models);
} else if (primaryOpt == exclude) {
FileBaseFilter fbf = new FileBaseFilter(cmdLine.getOptionValue(exclude.getOpt()));
modelIDs.removeAll(fbf.models);
}
}
// end command line processing
WriterFlusher flusher = new WriterFlusher(10);
PrintWriter detailWriter;
PrintWriter bngWriter;
PrintWriter sbmlWriter;
SBMLExceptionSorter sbmlExceptions;
if (isDetailed) {
detailWriter = new PrintWriter(new File(outDir, "compareDetail.txt"));
bngWriter = new PrintWriter(new File(outDir, "bngErrors.txt"));
sbmlWriter = new PrintWriter(new File(outDir, "sbmlErrors.txt"));
sbmlExceptions = new LiveSorter();
flusher.add(detailWriter);
flusher.add(bngWriter);
flusher.add(sbmlWriter);
} else {
detailWriter = new PrintWriter(new NullWriter());
bngWriter = new PrintWriter(new NullWriter());
sbmlWriter = new PrintWriter(new NullWriter());
sbmlExceptions = new NullSorter();
}
PropertyLoader.loadProperties();
/**
* example properties
*
* vcell.COPASI.executable = "C:\\Program Files\\COPASI\\bin\\CopasiSE.exe"
* vcell.mathSBML.directory = "c:\\developer\\eclipse\\workspace\\mathsbml\\"
* vcell.mathematica.kernel.executable = "C:\\Program Files\\Wolfram Research\\Mathematica\\7.0\\MathKernel.exe"
*/
ResourceUtil.setNativeLibraryDirectory();
PrintWriter printWriter = new PrintWriter(new FileWriter(new File(outDir, "summary.log"), true));
flusher.add(printWriter);
listModels(printWriter, modelIDs);
// SBMLSolver copasiSBMLSolver = new CopasiSBMLSolver();
SBMLSolver copasiSBMLSolver = new SBMLSolver() {
@Override
public File solve(String filePrefix, File outDir, String sbmlText, SimSpec testSpec) throws IOException, SolverException, SbmlException {
throw new RuntimeException("COPASI SOLVER REMOVED, REIMPLEMENT");
}
@Override
public String getResultsFileColumnDelimiter() {
throw new RuntimeException("COPASI SOLVER REMOVED, REIMPLEMENT");
}
};
try {
printWriter.println(" | *BIOMODEL ID* | *BioModel name* | *PASS* | *Rel Error (VC/COP)(VC/MSBML)(COP/MSBML)* | *Exception* | ");
removeToxic(modelIDs, printWriter);
for (BiomodelsNetEntry modelID : modelIDs) {
String modelName = service.getModelNameById(modelID.toString());
String modelSBML = service.getModelById(modelID.toString());
Element bioModelInfo = new Element(BioModelsNetPanel.BIOMODELINFO_ELEMENT_NAME);
bioModelInfo.setAttribute(BioModelsNetPanel.ID_ATTRIBUTE_NAME, modelID.toString());
bioModelInfo.setAttribute(BioModelsNetPanel.SUPPORTED_ATTRIBUTE_NAME, "false");
bioModelInfo.setAttribute("vcell_ran", "false");
bioModelInfo.setAttribute("COPASI_ran", "false");
bioModelInfo.setAttribute("mSBML_ran", "false");
bioModelInfo.setAttribute(BioModelsNetPanel.MODELNAME_ATTRIBUTE_NAME, modelName);
boolean bUseUTF8 = true;
File sbmlFile = new File(outDir, modelID + ".sbml");
XmlUtil.writeXMLStringToFile(modelSBML, sbmlFile.getAbsolutePath(), bUseUTF8);
PrintStream saved_sysout = System.out;
PrintStream saved_syserr = System.err;
PrintStream new_sysout = null;
PrintStream new_syserr = null;
try {
String filePrefix = modelID.toString();
String sbmlText = modelSBML;
File logFile = new File(outDir, filePrefix + ".log");
new_sysout = new PrintStream(logFile);
new_syserr = new_sysout;
System.setOut(new_sysout);
System.setErr(new_syserr);
StringBuffer combinedErrorBuffer = new StringBuffer();
SimSpec simSpec = SimSpec.fromSBML(modelSBML);
String[] varsToTest = simSpec.getVarsList();
printWriter.print("ModelId: " + modelID + ": ");
// if a model crashes out the libSBML dll, we will terminate abruptly. Flush
// summary log before that happens
printWriter.flush();
try {
//
// get VCell solution with an embedded "ROUND TRIP" (time and species concentrations)
//
ODESolverResultSet vcellResults_RT = null;
try {
VCellSBMLSolver vcellSBMLSolver_RT = new VCellSBMLSolver();
vcellSBMLSolver_RT.setRoundTrip(false);
// TODO try with round-trip later.
String columnDelimiter = vcellSBMLSolver_RT.getResultsFileColumnDelimiter();
File resultFile = vcellSBMLSolver_RT.solve(filePrefix, outDir, sbmlFile.getAbsolutePath(), simSpec);
vcellResults_RT = readResultFile(resultFile, columnDelimiter);
bioModelInfo.setAttribute("vcell_ran", "true");
} catch (BNGException e) {
bngWriter.println(modelID + " " + e.getMessage());
throw e;
} catch (SBMLImportException e) {
ModelException me = new ModelException(modelID, e);
write(sbmlWriter, me);
sbmlExceptions.add(me);
throw e;
} catch (Exception e) {
printWriter.println("vcell solve(roundtrip=true) failed");
e.printStackTrace(printWriter);
System.out.println("vcell solve(roundtrip=true) failed");
e.printStackTrace(System.out);
combinedErrorBuffer.append(" *VCELL_RT* _" + e.getMessage() + "_ ");
Element vcellSolverReport = new Element("SolverReport");
vcellSolverReport.setAttribute("solverName", "vcell");
vcellSolverReport.setAttribute("errorMessage", e.getMessage());
bioModelInfo.addContent(vcellSolverReport);
bioModelInfo.setAttribute("vcell_ran", "false");
}
//
// get COPASI solution (time and species concentrations)
//
ODESolverResultSet copasiResults = null;
try {
String columnDelimiter = copasiSBMLSolver.getResultsFileColumnDelimiter();
File resultFile = copasiSBMLSolver.solve(filePrefix, outDir, sbmlText, simSpec);
copasiResults = readResultFile(resultFile, columnDelimiter);
bioModelInfo.setAttribute("COPASI_ran", "true");
} catch (Exception e) {
printWriter.println("Copasi solve() failed");
e.printStackTrace(printWriter);
System.out.println("Copasi solve() failed");
e.printStackTrace(System.out);
combinedErrorBuffer.append(" *COPASI* _" + e.getMessage() + "_ ");
Element copasiSolverReport = new Element("SolverReport");
copasiSolverReport.setAttribute("solverName", "COPASI");
copasiSolverReport.setAttribute("errorMessage", e.getMessage());
bioModelInfo.addContent(copasiSolverReport);
bioModelInfo.setAttribute("COPASI_ran", "false");
}
//
// get mSBML solution (time and species concentrations)
//
/*
ODESolverResultSet mSBMLResults = null;
if (idInt != 246 && idInt != 250 && idInt != 285 && idInt != 301) {
try {
MathSBMLSolver mSBMLSolver = new MathSBMLSolver();
String columnDelimiter = mSBMLSolver.getResultsFileColumnDelimiter();
org.sbml.libsbml.SBMLDocument sbmlDocument = new SBMLReader().readSBML(sbmlFile.getAbsolutePath());
long level = sbmlDocument.getLevel();
long version = sbmlDocument.getVersion();
String mathsbmlFilePrefix = filePrefix;
if (level!=2 || (level==2 && version>3)){
// sbmlDocument.setConsistencyChecksForConversion(libsbmlConstants.LIBSBML_CAT_MODELING_PRACTICE, false);
long numErrors = sbmlDocument.checkL2v3Compatibility();
if (numErrors==0){
boolean bConversionWorked = sbmlDocument.setLevelAndVersion(2, 3, false);
SBMLDocument doc = new SBMLDocument(sbmlDocument);
doc.setLevelAndVersion(2,3,false);
if (bConversionWorked){
mathsbmlFilePrefix = filePrefix+"_L2V3";
long newVersion = doc.getVersion();
SBMLWriter sbmlWriter = new SBMLWriter();
sbmlText = sbmlWriter.writeToString(doc);
try {
XmlUtil.writeXMLStringToFile(sbmlText, mathsbmlFilePrefix+".sbml", true);
} catch (IOException e1) {
e1.printStackTrace(System.out);
}
}else{
throw new RuntimeException("couldn't convert SBML from L"+level+"V"+version+" to L2V3");
}
}else{
throw new RuntimeException("couldn't convert SBML from L"+level+"V"+version+" to L2V3, not compatible: "+sbmlDocument.getError(0));
}
}
File resultFile = mSBMLSolver.solve(mathsbmlFilePrefix, outDir, sbmlText, simSpec);
mSBMLResults = readResultFile(resultFile, columnDelimiter);
bioModelInfo.setAttribute("mSBML_ran","true");
}catch (Exception e){
printWriter.println("mSBML solve() failed");
e.printStackTrace(printWriter);
System.out.println("mSBML solve() failed");
e.printStackTrace(System.out);
combinedErrorBuffer.append(" *mSBML* _"+e.getMessage()+"_ ");
Element mSBMLSolverReport = new Element("SolverReport");
mSBMLSolverReport.setAttribute("solverName","mSBML");
mSBMLSolverReport.setAttribute("errorMessage",e.getMessage());
bioModelInfo.addContent(mSBMLSolverReport);
bioModelInfo.setAttribute("mSBML_ran","false");
}
}
*/
//
// compare results from COPASI and VCELL_RT
//
Boolean bCOPASI_VCELL_matched = null;
if (copasiResults != null && vcellResults_RT != null) {
try {
SimulationComparisonSummary summary = MathTestingUtilities.compareResultSets(copasiResults, vcellResults_RT, varsToTest, TestCaseNew.REGRESSION, 1e-5, 1e-5);
double maxRelError = summary.getMaxRelativeError();
bioModelInfo.setAttribute("COPASI_VCELL_maxRelErr", Double.toString(maxRelError));
if (maxRelError < 1) {
bCOPASI_VCELL_matched = true;
} else {
detailWriter.println(modelID + " " + modelName);
bCOPASI_VCELL_matched = false;
Element solverComparison = new Element("SolverComparison");
solverComparison.setAttribute("solver1", "vcell");
solverComparison.setAttribute("solver2", "COPASI");
VariableComparisonSummary[] failedVCSummaries = summary.getFailingVariableComparisonSummaries(1e-5, 1e-5);
for (VariableComparisonSummary vcSummary : failedVCSummaries) {
Element failedVariableSummary = getVariableSummary(vcSummary);
solverComparison.addContent(failedVariableSummary);
String ss = vcSummary.toShortString();
detailWriter.print('\t');
detailWriter.println(ss);
System.out.println(ss);
}
detailWriter.println();
bioModelInfo.addContent(solverComparison);
}
} catch (Exception e) {
printWriter.println("COMPARE VC/COPASI failed");
e.printStackTrace(printWriter);
System.out.println("COMPARE VC/COPASI failed");
e.printStackTrace(System.out);
combinedErrorBuffer.append(" *COMPARE VC/COPASI* _" + e.getMessage() + "_ ");
Element solverComparison = new Element("SolverComparison");
solverComparison.setAttribute("solver1", "vcell");
solverComparison.setAttribute("solver2", "COPASI");
solverComparison.setAttribute("error", e.getMessage());
bioModelInfo.addContent(solverComparison);
}
}
/*
Boolean bmSBML_VCELL_matched = null;
if (mSBMLResults!=null && vcellResults_RT!=null){
try {
SimulationComparisonSummary summary = MathTestingUtilities.compareUnEqualResultSets(mSBMLResults, vcellResults_RT, varsToTest, 1e-5, 1e-5, 1);
double maxRelError = summary.getMaxRelativeError();
bioModelInfo.setAttribute("mSBML_VCELL_maxRelErr", Double.toString(maxRelError));
if (maxRelError<1){
bmSBML_VCELL_matched = true;
}else{
bmSBML_VCELL_matched = false;
Element solverComparison = new Element("SolverComparison");
solverComparison.setAttribute("solver1","vcell");
solverComparison.setAttribute("solver2","mSBML");
VariableComparisonSummary[] failedVCSummaries = summary.getFailingVariableComparisonSummaries(1e-5, 1e-5);
for (VariableComparisonSummary vcSummary : failedVCSummaries){
Element failedVariableSummary = getVariableSummary(vcSummary);
solverComparison.addContent(failedVariableSummary);
System.out.println(vcSummary.toShortString());
}
bioModelInfo.addContent(solverComparison);
}
} catch (Exception e) {
printWriter.println("COMPARE VCRT/mSBML failed");
e.printStackTrace(printWriter);
System.out.println("COMPARE VCRT/mSBML failed");
e.printStackTrace(System.out);
combinedErrorBuffer.append(" *COMPARE VCRT/mSBML* _"+e.getMessage()+"_ ");
Element solverComparison = new Element("SolverComparison");
solverComparison.setAttribute("solver1","vcell");
solverComparison.setAttribute("solver2","mSBML");
solverComparison.setAttribute("error",e.getMessage());
bioModelInfo.addContent(solverComparison);
}
}
//
// compare results from COPASI and mSBML
//
Boolean bCOPASI_mSBML_matched = null;
if (copasiResults!=null && mSBMLResults!=null){
try {
SimulationComparisonSummary summary = MathTestingUtilities.compareUnEqualResultSets(copasiResults, mSBMLResults, varsToTest, 1e-5, 1e-5, 1);
double maxRelError = summary.getMaxRelativeError();
bioModelInfo.setAttribute("COPASI_mSBML_maxRelErr", Double.toString(maxRelError));
if (maxRelError<1){
bCOPASI_mSBML_matched = true;
}else{
bCOPASI_mSBML_matched = false;
Element solverComparison = new Element("SolverComparison");
solverComparison.setAttribute("solver1","COPASI");
solverComparison.setAttribute("solver2","mSBML");
solverComparison.setAttribute("result","failed");
VariableComparisonSummary[] vcSummaries = summary.getVariableComparisonSummaries();
for (VariableComparisonSummary vcSummary : vcSummaries){
Element failedVariableSummary = getVariableSummary(vcSummary);
solverComparison.addContent(failedVariableSummary);
System.out.println(vcSummary.toShortString());
}
bioModelInfo.addContent(solverComparison);
}
} catch (Exception e) {
printWriter.println("COMPARE COPASI/mSBML failed");
e.printStackTrace(printWriter);
System.out.println("COMPARE COPASI/mSBML failed");
e.printStackTrace(System.out);
combinedErrorBuffer.append(" *COMPARE COPASI/mSBML* _"+e.getMessage()+"_ ");
Element solverComparison = new Element("SolverComparison");
solverComparison.setAttribute("solver1","COPASI");
solverComparison.setAttribute("solver2","mSBML");
solverComparison.setAttribute("error",e.getMessage());
}
}
*/
if ((bCOPASI_VCELL_matched != null && bCOPASI_VCELL_matched.booleanValue())) {
// || /*(bmSBML_VCELL_matched!=null && bmSBML_VCELL_matched.booleanValue()) */ )
bioModelInfo.setAttribute(BioModelsNetPanel.SUPPORTED_ATTRIBUTE_NAME, "true");
printWriter.println("PASS");
} else {
bioModelInfo.setAttribute(BioModelsNetPanel.SUPPORTED_ATTRIBUTE_NAME, "false");
printWriter.println("FAIL");
}
} catch (Exception e) {
e.printStackTrace(printWriter);
combinedErrorBuffer.append(" *UNKNOWN* _" + e.getMessage() + "_ ");
bioModelInfo.setAttribute(BioModelsNetPanel.SUPPORTED_ATTRIBUTE_NAME, "false");
bioModelInfo.setAttribute("exception", e.getMessage());
}
printWriter.flush();
// write for each model just in case files get corrupted (it happened).
write(bioModelInfo, new File(outDir, modelID + "_report.xml"));
} finally {
if (new_sysout != null) {
new_sysout.close();
new_sysout = null;
}
if (new_syserr != null) {
new_syserr.close();
new_syserr = null;
}
System.setOut(saved_sysout);
System.setOut(saved_syserr);
}
}
// this writes out the SBML import exceptions grouped by type
if (!sbmlExceptions.isEmpty()) {
Map<Category, Collection<ModelException>> map = sbmlExceptions.getMap();
try (PrintWriter pw = new PrintWriter(new File(outDir, "sbmlSorted.txt"))) {
// SBMLImportException.Category
for (Category c : Category.values()) {
Collection<ModelException> meCollection = map.get(c);
for (ModelException me : meCollection) {
write(pw, me);
}
}
}
}
} finally {
printWriter.close();
detailWriter.close();
bngWriter.close();
}
} catch (Throwable e) {
e.printStackTrace(System.out);
e.printStackTrace(System.err);
}
System.exit(0);
}
use of cbit.vcell.solver.ode.ODESolverResultSet in project vcell by virtualcell.
the class BiomodelsDB_TestSuite method readResultFile.
public static ODESolverResultSet readResultFile(File resultFile, String delimiter) throws IOException {
String reportText = XmlUtil.getXMLString(resultFile.getAbsolutePath());
// System.out.println(reportText);
ODESolverResultSet odeResultSet = new ODESolverResultSet();
StringTokenizer lineTokenizer = new StringTokenizer(reportText, "\r\n", false);
int lineCount = 0;
while (lineTokenizer.hasMoreElements()) {
lineCount++;
String line = lineTokenizer.nextToken();
StringTokenizer columnTokenizer = new StringTokenizer(line, delimiter, false);
if (lineCount == 1) {
int wordCount = 0;
while (columnTokenizer.hasMoreTokens()) {
String label = columnTokenizer.nextToken();
if (wordCount == 0) {
label = "t";
}
odeResultSet.addDataColumn(new ODESolverResultSetColumnDescription(label));
wordCount++;
}
} else {
double[] values = new double[odeResultSet.getDataColumnCount()];
for (int i = 0; i < values.length; i++) {
values[i] = Double.parseDouble(columnTokenizer.nextToken());
}
odeResultSet.addRow(values);
}
}
return odeResultSet;
}
use of cbit.vcell.solver.ode.ODESolverResultSet in project vcell by virtualcell.
the class VCellSBMLSolver method solve.
public File solve(String filePrefix, File outDir, String sbmlFileName, SimSpec testSpec) throws IOException, SolverException, SbmlException {
try {
cbit.util.xml.VCLogger sbmlImportLogger = new LocalLogger();
//
// Instantiate an SBMLImporter to get the speciesUnitsHash - to compute the conversion factor from VC->SB species units.
// and import SBML (sbml->bioModel)
BioModel bioModel = importSBML(sbmlFileName, sbmlImportLogger, false);
// Hashtable<String, SBMLImporter.SBVCConcentrationUnits> speciesUnitsHash = sbmlImporter.getSpeciesUnitsHash();
// double timeFactor = sbmlImporter.getSBMLTimeUnitsFactor();
String vcml_1 = XmlHelper.bioModelToXML(bioModel);
SBMLUtils.writeStringToFile(vcml_1, new File(outDir, filePrefix + ".vcml").getAbsolutePath(), true);
if (bRoundTrip) {
// Round trip the bioModel (bioModel->sbml->bioModel).
// save imported "bioModel" as VCML
// String vcml_1 = XmlHelper.bioModelToXML(bioModel);
// SBMLUtils.writeStringToFile(vcml_1, new File(outDir,filePrefix+".vcml").getAbsolutePath());
// export bioModel as sbml and save
// String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, bioModel.getSimulationContexts(0).getName());
// SimulationJob simJob = new SimulationJob(bioModel.getSimulations(bioModel.getSimulationContexts(0))[0], null, 0);
String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, 0, false, bioModel.getSimulationContext(0), null);
SBMLUtils.writeStringToFile(vcml_sbml, new File(outDir, filePrefix + ".vcml.sbml").getAbsolutePath(), true);
// re-import bioModel from exported sbml
XMLSource vcml_sbml_Src = new XMLSource(vcml_sbml);
BioModel newBioModel = (BioModel) XmlHelper.importSBML(sbmlImportLogger, vcml_sbml_Src, false);
String vcml_sbml_vcml = XmlHelper.bioModelToXML(newBioModel);
SBMLUtils.writeStringToFile(vcml_sbml_vcml, new File(outDir, filePrefix + ".vcml.sbml.vcml").getAbsolutePath(), true);
// have rest of code use the round-tripped biomodel
bioModel = newBioModel;
}
//
// select only Application, generate math, and create a single Simulation.
//
SimulationContext simContext = bioModel.getSimulationContext(0);
MathMapping mathMapping = simContext.createNewMathMapping();
MathDescription mathDesc = mathMapping.getMathDescription();
String vcml = mathDesc.getVCML();
try (PrintWriter pw = new PrintWriter("vcmlTrace.txt")) {
pw.println(vcml);
}
simContext.setMathDescription(mathDesc);
SimulationVersion simVersion = new SimulationVersion(new KeyValue("100"), "unnamed", null, null, null, null, null, null, null, null);
Simulation sim = new Simulation(simVersion, mathDesc);
sim.setName("unnamed");
// if time factor from SBML is not 1 (i.e., it is not in secs but in minutes or hours), convert endTime to min/hr as : endTime*timeFactor
// double endTime = testSpec.getEndTime()*timeFactor;
double endTime = testSpec.getEndTime();
sim.getSolverTaskDescription().setTimeBounds(new TimeBounds(0, endTime));
TimeStep timeStep = new TimeStep();
sim.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep.getMinimumTimeStep(), timeStep.getDefaultTimeStep(), endTime / 10000));
sim.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((endTime - 0) / testSpec.getNumTimeSteps()));
sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
// sim.getSolverTaskDescription().setErrorTolerance(new cbit.vcell.solver.ErrorTolerance(1e-10, 1e-12));
// Generate .idaInput string
/* IDAFileWriter idaFileWriter = new IDAFileWriter(sim);
File idaInputFile = new File(filePathName.replace(".vcml", ".idaInput"));
PrintWriter idaPW = new java.io.PrintWriter(idaInputFile);
idaFileWriter.writeInputFile(idaPW);
idaPW.close();
// use the idastandalone solver
File idaOutputFile = new File(filePathName.replace(".vcml", ".ida"));
Executable executable = new Executable("IDAStandalone " + idaInputFile + " " + idaOutputFile);
executable.start();
*/
// Generate .cvodeInput string
File cvodeFile = new File(outDir, filePrefix + SimDataConstants.CVODEINPUT_DATA_EXTENSION);
PrintWriter cvodePW = new java.io.PrintWriter(cvodeFile);
SimulationJob simJob = new SimulationJob(sim, 0, null);
SimulationTask simTask = new SimulationTask(simJob, 0);
CVodeFileWriter cvodeFileWriter = new CVodeFileWriter(cvodePW, simTask);
cvodeFileWriter.write();
cvodePW.close();
// use the cvodeStandalone solver
File cvodeOutputFile = new File(outDir, filePrefix + SimDataConstants.IDA_DATA_EXTENSION);
String executableName = null;
try {
executableName = SolverUtilities.getExes(SolverDescription.CVODE)[0].getAbsolutePath();
} catch (IOException e) {
throw new RuntimeException("failed to get executable for solver " + SolverDescription.CVODE.getDisplayLabel() + ": " + e.getMessage(), e);
}
Executable executable = new Executable(new String[] { executableName, cvodeFile.getAbsolutePath(), cvodeOutputFile.getAbsolutePath() });
executable.start();
// get the result
ODESolverResultSet odeSolverResultSet = getODESolverResultSet(simJob, cvodeOutputFile.getPath());
//
// print header
//
File outputFile = new File(outDir, filePrefix + ".vcell.csv");
java.io.PrintStream outputStream = new java.io.PrintStream(new java.io.BufferedOutputStream(new java.io.FileOutputStream(outputFile)));
outputStream.print("time");
for (int i = 0; i < testSpec.getVarsList().length; i++) {
outputStream.print("," + testSpec.getVarsList()[i]);
}
outputStream.println();
//
// extract data for time and species
//
double[][] data = new double[testSpec.getVarsList().length + 1][];
int column = odeSolverResultSet.findColumn("t");
data[0] = odeSolverResultSet.extractColumn(column);
int origDataLength = data[0].length;
for (int i = 0; i < testSpec.getVarsList().length; i++) {
column = odeSolverResultSet.findColumn(testSpec.getVarsList()[i]);
if (column == -1) {
Variable var = simJob.getSimulationSymbolTable().getVariable(testSpec.getVarsList()[i]);
data[i + 1] = new double[data[0].length];
if (var instanceof cbit.vcell.math.Constant) {
double value = ((cbit.vcell.math.Constant) var).getExpression().evaluateConstant();
for (int j = 0; j < data[i + 1].length; j++) {
data[i + 1][j] = value;
}
} else {
throw new RuntimeException("Did not find " + testSpec.getVarsList()[i] + " in simulation");
}
} else {
data[i + 1] = odeSolverResultSet.extractColumn(column);
}
}
//
// for each time, print row
//
int index = 0;
double[] sampleTimes = new double[testSpec.getNumTimeSteps() + 1];
for (int i = 0; i <= testSpec.getNumTimeSteps(); i++) {
sampleTimes[i] = endTime * i / testSpec.getNumTimeSteps();
}
Model vcModel = bioModel.getModel();
ReservedSymbol kMole = vcModel.getKMOLE();
for (int i = 0; i < sampleTimes.length; i++) {
//
while (true) {
//
if (index == odeSolverResultSet.getRowCount() - 1) {
if (data[0][index] == sampleTimes[i]) {
break;
} else {
throw new RuntimeException("sampleTime does not match at last time point");
}
}
//
if (data[0][index + 1] > sampleTimes[i]) {
break;
}
//
// sampleTime must be later in our data list.
//
index++;
}
// if data[0][index] == sampleTime no need to interpolate
if (data[0][index] == sampleTimes[i]) {
// if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
// if (timeFactor != 1.0) {
// outputStream.print(data[0][index]/timeFactor);
// } else {
outputStream.print(data[0][index]);
// }
for (int j = 0; j < testSpec.getVarsList().length; j++) {
// SBMLImporter.SBVCConcentrationUnits spConcUnits = speciesUnitsHash.get(testSpec.getVarsList()[j]);
// if (spConcUnits != null) {
// VCUnitDefinition sbunits = spConcUnits.getSBConcentrationUnits();
// VCUnitDefinition vcunits = spConcUnits.getVCConcentrationUnits();
// SBMLUnitParameter unitFactor = SBMLUtils.getConcUnitFactor("spConcParam", vcunits, sbunits, kMole);
// outputStream.print("," + data[j + 1][index] * unitFactor.getExpression().evaluateConstant()); //earlier, hack unitfactor = 0.000001
// earlier, hack unitfactor = 0.000001
outputStream.print("," + data[j + 1][index]);
// }
}
// System.out.println("No interpolation needed!");
outputStream.println();
} else {
// if data[0][index] < sampleTime, must interpolate
double fraction = (sampleTimes[i] - data[0][index]) / (data[0][index + 1] - data[0][index]);
// if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
// if (timeFactor != 1.0) {
// outputStream.print(sampleTimes[i]/timeFactor);
// } else {
outputStream.print(sampleTimes[i]);
// }
for (int j = 0; j < testSpec.getVarsList().length; j++) {
double interpolatedValue = 0.0;
double[] speciesVals = null;
double[] times = null;
// Currently using 2nd order interpolation
if (index == 0) {
// can only do 1st order interpolation
times = new double[] { data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1] };
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
} else if (index >= 1 && index <= origDataLength - 3) {
double val_1 = Math.abs(sampleTimes[i] - data[0][index - 1]);
double val_2 = Math.abs(sampleTimes[i] - data[0][index + 2]);
if (val_1 < val_2) {
times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
} else {
times = new double[] { data[0][index], data[0][index + 1], data[0][index + 2] };
speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1], data[j + 1][index + 2] };
}
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
} else {
times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
}
// // Currently using 1st order interpolation
// times = new double[] { data[0][index], data[0][index+1] };
// speciesVals = new double[] { data[j+1][index], data[j+1][index+1] };
// interpolatedValue = taylorInterpolation(sampleTimes[i], times, speciesVals);
// interpolatedValue = interpolatedValue * unitFactor.getExpression().evaluateConstant(); //earlier, hack unitfactor = 0.000001
// System.out.println("Sample time: " + sampleTimes[i] + ", between time[" + index + "]=" + data[0][index]+" and time["+(index+1)+"]="+(data[0][index+1])+", interpolated = "+interpolatedValue);
outputStream.print("," + interpolatedValue);
}
outputStream.println();
}
}
outputStream.close();
return outputFile;
} catch (RuntimeException e) {
e.printStackTrace(System.out);
// rethrow without losing context
throw e;
} catch (Exception e) {
e.printStackTrace(System.out);
throw new SolverException(e.getMessage(), e);
}
}
use of cbit.vcell.solver.ode.ODESolverResultSet in project vcell by virtualcell.
the class VCellSBMLSolver method solveVCell.
public File solveVCell(String filePrefix, File outDir, String sbmlFileName, SimSpec testSpec) throws IOException, SolverException, SbmlException {
try {
cbit.util.xml.VCLogger logger = new LocalLogger();
//
// Instantiate an SBMLImporter to get the speciesUnitsHash - to compute the conversion factor from VC->SB species units.
// and import SBML (sbml->bioModel)
org.vcell.sbml.vcell.SBMLImporter sbmlImporter = new org.vcell.sbml.vcell.SBMLImporter(sbmlFileName, logger, false);
BioModel bioModel = sbmlImporter.getBioModel();
if (bRoundTrip) {
// Round trip the bioModel (bioModel->sbml->bioModel).
// export bioModel as sbml and save
String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, 0, false, bioModel.getSimulationContext(0), null);
// re-import bioModel from exported sbml
XMLSource vcml_sbml_Src = new XMLSource(vcml_sbml);
BioModel newBioModel = (BioModel) XmlHelper.importSBML(logger, vcml_sbml_Src, false);
// have rest of code use the round-tripped biomodel
bioModel = newBioModel;
}
//
// select only Application, generate math, and create a single Simulation.
//
SimulationContext simContext = bioModel.getSimulationContext(0);
MathMapping mathMapping = simContext.createNewMathMapping();
MathDescription mathDesc = mathMapping.getMathDescription();
simContext.setMathDescription(mathDesc);
SimulationVersion simVersion = new SimulationVersion(new KeyValue("100"), "unnamed", null, null, null, null, null, null, null, null);
Simulation sim = new Simulation(simVersion, mathDesc);
sim.setName("unnamed");
// if time factor from SBML is not 1 (i.e., it is not in secs but in minutes or hours), convert endTime to min/hr as : endTime*timeFactor
// double endTime = testSpec.getEndTime()*timeFactor;
double endTime = testSpec.getEndTime();
sim.getSolverTaskDescription().setTimeBounds(new TimeBounds(0, endTime));
TimeStep timeStep = new TimeStep();
sim.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep.getMinimumTimeStep(), timeStep.getDefaultTimeStep(), endTime / 10000));
sim.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((endTime - 0) / testSpec.getNumTimeSteps()));
sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(testSpec.getAbsTolerance(), testSpec.getRelTolerance()));
// sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
// Generate .idaInput string
File idaInputFile = new File(outDir, filePrefix + SimDataConstants.IDAINPUT_DATA_EXTENSION);
PrintWriter idaPW = new java.io.PrintWriter(idaInputFile);
SimulationJob simJob = new SimulationJob(sim, 0, null);
SimulationTask simTask = new SimulationTask(simJob, 0);
IDAFileWriter idaFileWriter = new IDAFileWriter(idaPW, simTask);
idaFileWriter.write();
idaPW.close();
// use the idastandalone solver
File idaOutputFile = new File(outDir, filePrefix + SimDataConstants.IDA_DATA_EXTENSION);
// String sundialsSolverExecutable = "C:\\Developer\\Eclipse\\workspace\\VCell 4.8\\SundialsSolverStandalone_NoMessaging.exe";
String executableName = null;
try {
executableName = SolverUtilities.getExes(SolverDescription.IDA)[0].getAbsolutePath();
} catch (IOException e) {
throw new RuntimeException("failed to get executable for solver " + SolverDescription.IDA.getDisplayLabel() + ": " + e.getMessage(), e);
}
Executable executable = new Executable(new String[] { executableName, idaInputFile.getAbsolutePath(), idaOutputFile.getAbsolutePath() });
executable.start();
/* // Generate .cvodeInput string
File cvodeFile = new File(outDir,filePrefix+SimDataConstants.CVODEINPUT_DATA_EXTENSION);
PrintWriter cvodePW = new java.io.PrintWriter(cvodeFile);
SimulationJob simJob = new SimulationJob(sim, 0, null);
CVodeFileWriter cvodeFileWriter = new CVodeFileWriter(cvodePW, simJob);
cvodeFileWriter.write();
cvodePW.close();
// use the cvodeStandalone solver
File cvodeOutputFile = new File(outDir,filePrefix+SimDataConstants.IDA_DATA_EXTENSION);
String sundialsSolverExecutable = PropertyLoader.getRequiredProperty(PropertyLoader.sundialsSolverExecutableProperty);
Executable executable = new Executable(new String[]{sundialsSolverExecutable, cvodeFile.getAbsolutePath(), cvodeOutputFile.getAbsolutePath()});
executable.start();
*/
// get the result
ODESolverResultSet odeSolverResultSet = getODESolverResultSet(simJob, idaOutputFile.getPath());
// remove CVOde input and output files ??
idaInputFile.delete();
idaOutputFile.delete();
//
// print header
//
File outputFile = new File(outDir, "results" + filePrefix + ".csv");
java.io.PrintStream outputStream = new java.io.PrintStream(new java.io.BufferedOutputStream(new java.io.FileOutputStream(outputFile)));
outputStream.print("time");
for (int i = 0; i < testSpec.getVarsList().length; i++) {
outputStream.print("," + testSpec.getVarsList()[i]);
}
outputStream.println();
//
// extract data for time and species
//
double[][] data = new double[testSpec.getVarsList().length + 1][];
int column = odeSolverResultSet.findColumn("t");
data[0] = odeSolverResultSet.extractColumn(column);
int origDataLength = data[0].length;
for (int i = 0; i < testSpec.getVarsList().length; i++) {
column = odeSolverResultSet.findColumn(testSpec.getVarsList()[i]);
if (column == -1) {
Variable var = simJob.getSimulationSymbolTable().getVariable(testSpec.getVarsList()[i]);
data[i + 1] = new double[data[0].length];
if (var instanceof cbit.vcell.math.Constant) {
double value = ((cbit.vcell.math.Constant) var).getExpression().evaluateConstant();
for (int j = 0; j < data[i + 1].length; j++) {
data[i + 1][j] = value;
}
} else {
throw new RuntimeException("Did not find " + testSpec.getVarsList()[i] + " in simulation");
}
} else {
data[i + 1] = odeSolverResultSet.extractColumn(column);
}
}
//
// for each time, print row
//
int index = 0;
double[] sampleTimes = new double[testSpec.getNumTimeSteps() + 1];
for (int i = 0; i <= testSpec.getNumTimeSteps(); i++) {
sampleTimes[i] = endTime * i / testSpec.getNumTimeSteps();
}
Model vcModel = bioModel.getModel();
ReservedSymbol kMole = vcModel.getKMOLE();
for (int i = 0; i < sampleTimes.length; i++) {
//
while (true) {
//
if (index == odeSolverResultSet.getRowCount() - 1) {
if (data[0][index] == sampleTimes[i]) {
break;
} else {
throw new RuntimeException("sampleTime does not match at last time point");
}
}
//
if (data[0][index + 1] > sampleTimes[i]) {
break;
}
//
// sampleTime must be later in our data list.
//
index++;
}
// if data[0][index] == sampleTime no need to interpolate
if (data[0][index] == sampleTimes[i]) {
// if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
// if (timeFactor != 1.0) {
// outputStream.print(data[0][index]/timeFactor);
// } else {
outputStream.print(data[0][index]);
// }
for (int j = 0; j < testSpec.getVarsList().length; j++) {
// SBMLImporter.SBVCConcentrationUnits spConcUnits = speciesUnitsHash.get(testSpec.getVarsList()[j]);
// if (spConcUnits != null) {
// VCUnitDefinition sbunits = spConcUnits.getSBConcentrationUnits();
// VCUnitDefinition vcunits = spConcUnits.getVCConcentrationUnits();
// SBMLUnitParameter unitFactor = SBMLUtils.getConcUnitFactor("spConcParam", vcunits, sbunits, kMole);
// outputStream.print("," + data[j + 1][index] * unitFactor.getExpression().evaluateConstant()); //earlier, hack unitfactor = 0.000001
// earlier, hack unitfactor = 0.000001
outputStream.print("," + data[j + 1][index]);
// }
}
// System.out.println("No interpolation needed!");
outputStream.println();
} else {
// if data[0][index] < sampleTime, must interpolate
double fraction = (sampleTimes[i] - data[0][index]) / (data[0][index + 1] - data[0][index]);
// if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
// if (timeFactor != 1.0) {
// outputStream.print(sampleTimes[i]/timeFactor);
// } else {
outputStream.print(sampleTimes[i]);
// }
for (int j = 0; j < testSpec.getVarsList().length; j++) {
double interpolatedValue = 0.0;
double[] speciesVals = null;
double[] times = null;
// Currently using 2nd order interpolation
if (index == 0) {
// can only do 1st order interpolation
times = new double[] { data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1] };
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
} else if (index >= 1 && index <= origDataLength - 3) {
double val_1 = Math.abs(sampleTimes[i] - data[0][index - 1]);
double val_2 = Math.abs(sampleTimes[i] - data[0][index + 2]);
if (val_1 < val_2) {
times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
} else {
times = new double[] { data[0][index], data[0][index + 1], data[0][index + 2] };
speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1], data[j + 1][index + 2] };
}
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
} else {
times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
}
// // Currently using 1st order interpolation
// times = new double[] { data[0][index], data[0][index+1] };
// speciesVals = new double[] { data[j+1][index], data[j+1][index+1] };
// interpolatedValue = taylorInterpolation(sampleTimes[i], times, speciesVals);
// interpolatedValue = interpolatedValue * unitFactor.getExpression().evaluateConstant(); //earlier, hack unitfactor = 0.000001
// System.out.println("Sample time: " + sampleTimes[i] + ", between time[" + index + "]=" + data[0][index]+" and time["+(index+1)+"]="+(data[0][index+1])+", interpolated = "+interpolatedValue);
outputStream.print("," + interpolatedValue);
}
outputStream.println();
}
}
outputStream.close();
return outputFile;
} catch (Exception e) {
e.printStackTrace(System.out);
// File outputFile = new File(outDir,"results" + filePrefix + ".csv");
throw new SolverException(e.getMessage());
}
}
use of cbit.vcell.solver.ode.ODESolverResultSet in project vcell by virtualcell.
the class VCellSBMLSolver method getODESolverResultSet.
public static ODESolverResultSet getODESolverResultSet(SimulationJob argSimJob, String idaFileName) {
// read .ida file
ODESolverResultSet odeSolverResultSet = new ODESolverResultSet();
FileInputStream inputStream = null;
try {
inputStream = new FileInputStream(idaFileName);
InputStreamReader inputStreamReader = new InputStreamReader(inputStream);
BufferedReader bufferedReader = new BufferedReader(inputStreamReader);
// Read header
String line = bufferedReader.readLine();
if (line == null) {
// throw exception
}
while (line.indexOf(':') > 0) {
String name = line.substring(0, line.indexOf(':'));
odeSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(name));
line = line.substring(line.indexOf(':') + 1);
}
// Read data
while ((line = bufferedReader.readLine()) != null) {
line = line + "\t";
double[] values = new double[odeSolverResultSet.getDataColumnCount()];
boolean bCompleteRow = true;
for (int i = 0; i < odeSolverResultSet.getDataColumnCount(); i++) {
if (line.indexOf('\t') == -1) {
bCompleteRow = false;
break;
} else {
String value = line.substring(0, line.indexOf('\t')).trim();
values[i] = Double.valueOf(value).doubleValue();
line = line.substring(line.indexOf('\t') + 1);
}
}
if (bCompleteRow) {
odeSolverResultSet.addRow(values);
} else {
break;
}
}
//
} catch (Exception e) {
e.printStackTrace(System.out);
} finally {
try {
if (inputStream != null) {
inputStream.close();
}
} catch (Exception ex) {
ex.printStackTrace(System.out);
}
}
// add appropriate Function columns to result set
cbit.vcell.math.Function[] functions = argSimJob.getSimulationSymbolTable().getFunctions();
for (int i = 0; i < functions.length; i++) {
if (SimulationSymbolTable.isFunctionSaved(functions[i])) {
Expression exp1 = new Expression(functions[i].getExpression());
try {
exp1 = argSimJob.getSimulationSymbolTable().substituteFunctions(exp1);
} catch (cbit.vcell.math.MathException e) {
e.printStackTrace(System.out);
throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
} catch (cbit.vcell.parser.ExpressionException e) {
e.printStackTrace(System.out);
throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
}
try {
FunctionColumnDescription cd = new FunctionColumnDescription(exp1.flatten(), functions[i].getName(), null, functions[i].getName(), false);
odeSolverResultSet.addFunctionColumn(cd);
} catch (cbit.vcell.parser.ExpressionException e) {
e.printStackTrace(System.out);
}
}
}
return odeSolverResultSet;
}
Aggregations