Search in sources :

Example 26 with ODESolverResultSet

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);
}
Also used : BNGException(cbit.vcell.server.bionetgen.BNGException) HelpFormatter(org.apache.commons.cli.HelpFormatter) TreeSet(java.util.TreeSet) SBMLSolver(org.vcell.sbml.SBMLSolver) VariableComparisonSummary(cbit.vcell.solver.test.VariableComparisonSummary) BioModelsWebServicesServiceLocator(uk.ac.ebi.www.biomodels_main.services.BioModelsWebServices.BioModelsWebServicesServiceLocator) Collection(java.util.Collection) Option(org.apache.commons.cli.Option) File(java.io.File) Options(org.apache.commons.cli.Options) SortedSet(java.util.SortedSet) Set(java.util.Set) ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet) TreeSet(java.util.TreeSet) HashSet(java.util.HashSet) Category(org.vcell.sbml.vcell.SBMLImportException.Category) OperatingSystemInfo(cbit.vcell.resource.OperatingSystemInfo) FileWriter(java.io.FileWriter) Element(org.jdom.Element) NullWriter(org.apache.commons.io.output.NullWriter) OptionGroup(org.apache.commons.cli.OptionGroup) ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet) SimulationComparisonSummary(cbit.vcell.solver.test.SimulationComparisonSummary) CommandLineParser(org.apache.commons.cli.CommandLineParser) DefaultParser(org.apache.commons.cli.DefaultParser) PrintWriter(java.io.PrintWriter) PrintStream(java.io.PrintStream) SBMLImportException(org.vcell.sbml.vcell.SBMLImportException) SimSpec(org.vcell.sbml.SimSpec) SbmlException(org.vcell.sbml.SbmlException) FileNotFoundException(java.io.FileNotFoundException) BNGException(cbit.vcell.server.bionetgen.BNGException) SBMLImportException(org.vcell.sbml.vcell.SBMLImportException) ParseException(org.apache.commons.cli.ParseException) SolverException(cbit.vcell.solver.SolverException) IOException(java.io.IOException) BioModelsWebServices(uk.ac.ebi.www.biomodels_main.services.BioModelsWebServices.BioModelsWebServices) CommandLine(org.apache.commons.cli.CommandLine) ParseException(org.apache.commons.cli.ParseException)

Example 27 with ODESolverResultSet

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;
}
Also used : StringTokenizer(java.util.StringTokenizer) ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet) ODESolverResultSetColumnDescription(cbit.vcell.math.ODESolverResultSetColumnDescription)

Example 28 with ODESolverResultSet

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);
    }
}
Also used : KeyValue(org.vcell.util.document.KeyValue) SimulationTask(cbit.vcell.messaging.server.SimulationTask) Variable(cbit.vcell.math.Variable) MathDescription(cbit.vcell.math.MathDescription) ReservedSymbol(cbit.vcell.model.Model.ReservedSymbol) CVodeFileWriter(cbit.vcell.solver.ode.CVodeFileWriter) TimeBounds(cbit.vcell.solver.TimeBounds) TimeStep(cbit.vcell.solver.TimeStep) SimulationVersion(org.vcell.util.document.SimulationVersion) ErrorTolerance(cbit.vcell.solver.ErrorTolerance) ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet) Executable(org.vcell.util.exe.Executable) SimulationJob(cbit.vcell.solver.SimulationJob) PrintWriter(java.io.PrintWriter) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) VCLogger(cbit.util.xml.VCLogger) IOException(java.io.IOException) SimulationContext(cbit.vcell.mapping.SimulationContext) ExecutableException(org.vcell.util.exe.ExecutableException) XMLStreamException(javax.xml.stream.XMLStreamException) XmlParseException(cbit.vcell.xml.XmlParseException) SolverException(cbit.vcell.solver.SolverException) SbmlException(org.vcell.sbml.SbmlException) IOException(java.io.IOException) SBMLImportException(org.vcell.sbml.vcell.SBMLImportException) Simulation(cbit.vcell.solver.Simulation) BioModel(cbit.vcell.biomodel.BioModel) MathMapping(cbit.vcell.mapping.MathMapping) BioModel(cbit.vcell.biomodel.BioModel) Model(cbit.vcell.model.Model) SolverException(cbit.vcell.solver.SolverException) File(java.io.File) XMLSource(cbit.vcell.xml.XMLSource)

Example 29 with ODESolverResultSet

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());
    }
}
Also used : KeyValue(org.vcell.util.document.KeyValue) SimulationTask(cbit.vcell.messaging.server.SimulationTask) Variable(cbit.vcell.math.Variable) SBMLImporter(org.vcell.sbml.vcell.SBMLImporter) MathDescription(cbit.vcell.math.MathDescription) ReservedSymbol(cbit.vcell.model.Model.ReservedSymbol) TimeBounds(cbit.vcell.solver.TimeBounds) TimeStep(cbit.vcell.solver.TimeStep) SimulationVersion(org.vcell.util.document.SimulationVersion) ErrorTolerance(cbit.vcell.solver.ErrorTolerance) ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet) Executable(org.vcell.util.exe.Executable) SimulationJob(cbit.vcell.solver.SimulationJob) PrintWriter(java.io.PrintWriter) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) VCLogger(cbit.util.xml.VCLogger) SBMLImporter(org.vcell.sbml.vcell.SBMLImporter) IOException(java.io.IOException) SimulationContext(cbit.vcell.mapping.SimulationContext) ExecutableException(org.vcell.util.exe.ExecutableException) XMLStreamException(javax.xml.stream.XMLStreamException) XmlParseException(cbit.vcell.xml.XmlParseException) SolverException(cbit.vcell.solver.SolverException) SbmlException(org.vcell.sbml.SbmlException) IOException(java.io.IOException) SBMLImportException(org.vcell.sbml.vcell.SBMLImportException) IDAFileWriter(cbit.vcell.solver.ode.IDAFileWriter) Simulation(cbit.vcell.solver.Simulation) BioModel(cbit.vcell.biomodel.BioModel) MathMapping(cbit.vcell.mapping.MathMapping) BioModel(cbit.vcell.biomodel.BioModel) Model(cbit.vcell.model.Model) SolverException(cbit.vcell.solver.SolverException) XMLSource(cbit.vcell.xml.XMLSource) File(java.io.File)

Example 30 with ODESolverResultSet

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;
}
Also used : InputStreamReader(java.io.InputStreamReader) FileInputStream(java.io.FileInputStream) ExecutableException(org.vcell.util.exe.ExecutableException) XMLStreamException(javax.xml.stream.XMLStreamException) XmlParseException(cbit.vcell.xml.XmlParseException) SolverException(cbit.vcell.solver.SolverException) SbmlException(org.vcell.sbml.SbmlException) IOException(java.io.IOException) SBMLImportException(org.vcell.sbml.vcell.SBMLImportException) Expression(cbit.vcell.parser.Expression) BufferedReader(java.io.BufferedReader) ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet) ODESolverResultSetColumnDescription(cbit.vcell.math.ODESolverResultSetColumnDescription) FunctionColumnDescription(cbit.vcell.math.FunctionColumnDescription)

Aggregations

ODESolverResultSet (cbit.vcell.solver.ode.ODESolverResultSet)37 ODESolverResultSetColumnDescription (cbit.vcell.math.ODESolverResultSetColumnDescription)18 FunctionColumnDescription (cbit.vcell.math.FunctionColumnDescription)11 Expression (cbit.vcell.parser.Expression)10 DataSource (cbit.vcell.modelopt.DataSource)9 ReferenceData (cbit.vcell.opt.ReferenceData)8 ExpressionException (cbit.vcell.parser.ExpressionException)7 Simulation (cbit.vcell.solver.Simulation)7 SimulationSymbolTable (cbit.vcell.solver.SimulationSymbolTable)6 SolverException (cbit.vcell.solver.SolverException)6 IOException (java.io.IOException)6 Variable (cbit.vcell.math.Variable)4 SimulationTask (cbit.vcell.messaging.server.SimulationTask)4 SimpleReferenceData (cbit.vcell.opt.SimpleReferenceData)4 SbmlException (org.vcell.sbml.SbmlException)4 SBMLImportException (org.vcell.sbml.vcell.SBMLImportException)4 BioModel (cbit.vcell.biomodel.BioModel)3 Function (cbit.vcell.math.Function)3 MathDescription (cbit.vcell.math.MathDescription)3 MathException (cbit.vcell.math.MathException)3