Examples with ErrorTolerance cbit.vcell.solver.ErrorTolerance used on opensource projects

Example 11 with ErrorTolerance

use of cbit.vcell.solver.ErrorTolerance in project vcell by virtualcell.

the class SbmlVcmlConverter method main.

/**
 * @param args -import or -export
 */
public static void main(String[] args) {
    Logging.init();
    ResourceUtil.setNativeLibraryDirectory();
    if (args.length < 2 || args.length > 3) {
        System.out.println("Usage:\n\t -export path_of_input_file\n\tOR\n\t -import path_of_input_file [-simulate]");
        System.exit(1);
    }
    if (args.length > 2 && args[0].equals("-export")) {
        System.out.println("Export cannot have arguments other than input file.");
        System.out.println("Usage:\n\t -export path_of_input_file\n\tOR\n\t-import path_of_input_file [-simulate]");
        System.exit(1);
    }
    try {
        String pathName = args[1];
        // Read xml file (Sbml or Vcml) into stringbuffer, pass the string into VCell's importer/exporter
        String xmlString = getXMLString(pathName);
        if (args[0].equals("-import")) {
            if (args.length > 3 || (args.length == 3 && !args[2].equals("-simulate"))) {
                System.out.println(args[2] + " : Unknown arguments for import; please check and re-run import command");
                System.out.println("Usage:\n\t -import path_of_input_file [-simulate]");
                System.exit(1);
            }
            boolean bSimulate = false;
            if (args.length == 3 && args[2].equals("-simulate")) {
                bSimulate = true;
            }
            // Create a default VCLogger - SBMLImporter needs it
            cbit.util.xml.VCLogger logger = new cbit.util.xml.VCLogger() {

                @Override
                public void sendMessage(Priority p, ErrorType et, String message) {
                    System.err.println("LOGGER: msgLevel=" + p + ", msgType=" + et + ", " + message);
                    if (p == VCLogger.Priority.HighPriority) {
                        throw new RuntimeException("Import failed : " + message);
                    }
                }

                public void sendAllMessages() {
                }

                public boolean hasMessages() {
                    return false;
                }
            };
            // invoke SBMLImporter, which returns a Biomodel, convert that to VCML using XMLHelper
            try {
                // import SBML model into VCML, store VCML string in 'fileName.vcml'
                // 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(pathName, logger, false);
                BioModel bioModel = sbmlImporter.getBioModel();
                // Hashtable<String, SBMLImporter.SBVCConcentrationUnits> speciesUnitsHash = sbmlImporter.getSpeciesUnitsHash();
                // double timeFactor = sbmlImporter.getSBMLTimeUnitsFactor();
                // convert biomodel to vcml and save to file.
                String vcmlString = XmlHelper.bioModelToXML(bioModel);
                String fileExtensionStr;
                if (pathName.endsWith(".xml")) {
                    fileExtensionStr = ".xml";
                } else if (pathName.endsWith(".sbml")) {
                    fileExtensionStr = ".sbml";
                } else {
                    throw new RuntimeException("Unknown file extension for SBML file name; Exiting SbmlConverter.");
                }
                String vcmlFileName = pathName.replace(fileExtensionStr, ".vcml");
                File vcmlFile = new File(vcmlFileName);
                XmlUtil.writeXMLStringToFile(vcmlString, vcmlFile.getAbsolutePath(), true);
                // If user doesn't choose to simulate, you are done.
                if (!bSimulate) {
                    return;
                }
                // Generate math for lone simContext
                SimulationContext simContext = (SimulationContext) bioModel.getSimulationContext(0);
                MathDescription mathDesc = simContext.createNewMathMapping().getMathDescription();
                simContext.setMathDescription(mathDesc);
                // Create basic simulation, with IDA solver (set in solve method) and other defaults, and end time 'Te'
                org.vcell.util.document.SimulationVersion simVersion = new org.vcell.util.document.SimulationVersion(new KeyValue("12345"), "name", new org.vcell.util.document.User("user", new KeyValue("123")), new org.vcell.util.document.GroupAccessNone(), // versionBranchPointRef
                null, // branchID
                new java.math.BigDecimal(1.0), new java.util.Date(), org.vcell.util.document.VersionFlag.Archived, "", null);
                Simulation sim1 = new Simulation(simVersion, simContext.getMathDescription());
                simContext.addSimulation(sim1);
                sim1.setName("sim1");
                // using a default end time of 10.0 secs, and a forcing default output time step of 0.01.
                // If user wants anything different, user can modify the .idaInput file and re-run simulation separately.
                // double newEndTime = endTime * timeFactor;
                double newEndTime = endTime;
                sim1.getSolverTaskDescription().setTimeBounds(new cbit.vcell.solver.TimeBounds(0, newEndTime));
                TimeStep timeStep_1 = new TimeStep();
                sim1.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep_1.getMinimumTimeStep(), timeStep_1.getDefaultTimeStep(), newEndTime / 10000));
                sim1.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((newEndTime - 0) / numTimeSteps));
                sim1.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
                // save the new vcml with generated math and new sim in file named : fileName_IDA_simulation.vcml
                String newVcmlString = XmlHelper.bioModelToXML(bioModel);
                String newVcmlFileName = vcmlFile.getPath().replace(".vcml", "_IDA_simulation.vcml");
                File newVcmlFile = new File(newVcmlFileName);
                XmlUtil.writeXMLStringToFile(newVcmlString, newVcmlFile.getAbsolutePath(), true);
                // Solve simulation, which also generates and saves the .idainput file and .csv
                // SimSpec to get vars to solve and output in .csv
                SimSpec simSpec = SimSpec.fromSBML(xmlString);
                solveSimulation(new SimulationJob(sim1, 0, null), vcmlFile.getPath(), simSpec);
            } catch (Exception e) {
                e.printStackTrace(System.err);
                // and an annotation saying why it failed.
                try {
                    String fileExtensionStr;
                    if (pathName.endsWith(".xml")) {
                        fileExtensionStr = ".xml";
                    } else if (pathName.endsWith(".sbml")) {
                        fileExtensionStr = ".sbml";
                    } else {
                        throw new RuntimeException("Unknown file extension for SBML file name; Exiting SbmlConverter.");
                    }
                    String dummyVcmlFileName = pathName.replace(fileExtensionStr, ".vcml");
                    File dummyVcmlFile = new File(dummyVcmlFileName);
                    BioModel dummy_biomodel = new BioModel(null);
                    String dummyName = dummyVcmlFile.getName().substring(0, dummyVcmlFile.getName().indexOf(".vcml"));
                    dummy_biomodel.setName(dummyName);
                    dummy_biomodel.setDescription("SBML Model could not be automatically converted to VCML : " + e.getMessage());
                    String vcmlString = XmlHelper.bioModelToXML(dummy_biomodel);
                    XmlUtil.writeXMLStringToFile(vcmlString, dummyVcmlFile.getAbsolutePath(), true);
                } catch (Exception e1) {
                    e.printStackTrace(System.err);
                }
            }
        } else if (args[0].equals("-export")) {
            try {
                BioModel bioModel = XmlHelper.XMLToBioModel(new XMLSource(xmlString));
                for (int i = 0; i < bioModel.getSimulationContexts().length; i++) {
                    SimulationContext simContext = bioModel.getSimulationContext(i);
                    if (simContext.getGeometry().getDimension() == 0 && !simContext.isStoch()) {
                        if (!simContext.getGeometryContext().isAllSizeSpecifiedPositive()) {
                            Structure structure = simContext.getModel().getStructure(0);
                            double structureSize = 1.0;
                            StructureMapping structMapping = simContext.getGeometryContext().getStructureMapping(structure);
                            StructureSizeSolver.updateAbsoluteStructureSizes(simContext, structure, structureSize, structMapping.getSizeParameter().getUnitDefinition());
                        }
                        // Export the application itself, with default overrides
                        String sbmlString = XmlHelper.exportSBML(bioModel, 2, 4, 0, false, simContext, null);
                        String filePath = pathName.substring(0, pathName.lastIndexOf("\\") + 1);
                        String sbmlFileName = TokenMangler.mangleToSName(bioModel.getName() + "_" + i);
                        File sbmlFile = new File(filePath + sbmlFileName + ".xml");
                        XmlUtil.writeXMLStringToFile(sbmlString, sbmlFile.getAbsolutePath(), true);
                        // Now export each simulation in the application that has overrides
                        Simulation[] simulations = bioModel.getSimulations(simContext);
                        for (int j = 0; j < simulations.length; j++) {
                            if (simulations[j].getMathOverrides().hasOverrides()) {
                                // Check for parameter scan and create simJob to pass into exporter
                                for (int k = 0; k < simulations[j].getScanCount(); k++) {
                                    SimulationJob simJob = new SimulationJob(simulations[j], k, null);
                                    sbmlString = XmlHelper.exportSBML(bioModel, 2, 4, 0, false, simContext, simJob);
                                    String fileName = TokenMangler.mangleToSName(sbmlFileName + "_" + j + "_" + k);
                                    sbmlFile = new File(filePath + fileName + ".xml");
                                    XmlUtil.writeXMLStringToFile(sbmlString, sbmlFile.getAbsolutePath(), true);
                                }
                            }
                        }
                    }
                }
                System.out.println("Successfully translated model : " + pathName);
            } catch (Exception e) {
                e.printStackTrace(System.err);
            }
        }
        System.exit(0);
    } catch (IOException e) {
        e.printStackTrace(System.err);
        System.exit(1);
    }
}

Example 12 with ErrorTolerance

use of cbit.vcell.solver.ErrorTolerance 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);
    }
}

Example 13 with ErrorTolerance

use of cbit.vcell.solver.ErrorTolerance 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());
    }
}

Example 14 with ErrorTolerance

use of cbit.vcell.solver.ErrorTolerance in project vcell by virtualcell.

the class OdeFileWriter method write.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
public void write(String[] parameterNames) throws Exception {
    createStateVariables();
    createSymbolTable();
    Simulation simulation = simTask.getSimulation();
    if (simulation.getSolverTaskDescription().getUseSymbolicJacobian()) {
        throw new RuntimeException("symbolic jacobian option not yet supported in interpreted Stiff solver");
    }
    writeJMSParamters();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
    ErrorTolerance errorTolerance = solverTaskDescription.getErrorTolerance();
    printWriter.println("SOLVER " + getSolverName());
    printWriter.println("STARTING_TIME " + timeBounds.getStartingTime());
    printWriter.println("ENDING_TIME " + timeBounds.getEndingTime());
    printWriter.println("RELATIVE_TOLERANCE " + errorTolerance.getRelativeErrorTolerance());
    printWriter.println("ABSOLUTE_TOLERANCE " + errorTolerance.getAbsoluteErrorTolerance());
    printWriter.println("MAX_TIME_STEP " + simulation.getSolverTaskDescription().getTimeStep().getMaximumTimeStep());
    OutputTimeSpec ots = simulation.getSolverTaskDescription().getOutputTimeSpec();
    if (ots.isDefault()) {
        printWriter.println("KEEP_EVERY " + ((DefaultOutputTimeSpec) ots).getKeepEvery());
    } else if (ots.isUniform()) {
        printWriter.println("OUTPUT_TIME_STEP " + ((UniformOutputTimeSpec) ots).getOutputTimeStep());
    } else if (ots.isExplicit()) {
        printWriter.println("OUTPUT_TIMES " + ((ExplicitOutputTimeSpec) ots).getNumTimePoints());
        printWriter.println(((ExplicitOutputTimeSpec) ots).toSpaceSeperatedMultiLinesOfString());
    }
    if (parameterNames != null && parameterNames.length != 0) {
        printWriter.println("NUM_PARAMETERS " + parameterNames.length);
        for (int i = 0; i < parameterNames.length; i++) {
            printWriter.println(parameterNames[i]);
        }
    }
    HashMap<Discontinuity, String> discontinuityNameMap = new HashMap<Discontinuity, String>();
    String eventString = null;
    if (simulation.getMathDescription().hasEvents()) {
        eventString = writeEvents(discontinuityNameMap);
    }
    String equationString = writeEquations(discontinuityNameMap);
    if (discontinuityNameMap.size() > 0) {
        printWriter.println("DISCONTINUITIES " + discontinuityNameMap.size());
        for (Discontinuity od : discontinuityNameMap.keySet()) {
            printWriter.println(discontinuityNameMap.get(od) + " " + od.getDiscontinuityExp().flatten().infix() + "; " + od.getRootFindingExp().flatten().infix() + ";");
        }
    }
    if (eventString != null) {
        printWriter.print(eventString);
    }
    printWriter.println("NUM_EQUATIONS " + getStateVariableCount());
    printWriter.println(equationString);
}