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

Example 11 with TimeStep

use of cbit.vcell.solver.TimeStep 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 12 with TimeStep

use of cbit.vcell.solver.TimeStep 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 13 with TimeStep

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

the class ComsolModelBuilder method getVCCModel.

public static VCCModel getVCCModel(SimulationJob vcellSimJob) throws ExpressionException {
    MathDescription vcellMathDesc = vcellSimJob.getSimulation().getMathDescription();
    Geometry vcellGeometry = vcellMathDesc.getGeometry();
    GeometrySpec vcellGeometrySpec = vcellGeometry.getGeometrySpec();
    int vcellDim = vcellGeometrySpec.getDimension();
    VCCModel model = new VCCModel("Model", vcellDim);
    model.modelpath = "D:\\Developer\\eclipse\\workspace_refactor\\comsol_java\\src";
    model.comments = "Untitled\n\n";
    VCCModelNode comp1 = new VCCModelNode("comp1");
    model.modelnodes.add(comp1);
    // if (vcellDim != 2){
    // throw new RuntimeException("expecting 2D simulation");
    // }
    // 
    // assume initial geometry is circle centered at 0.5, 0.5 of radius 0.3
    // 
    // String comsolOutsideDomainName = "dif1";
    // String comsolInsideDomainName = "c1";
    VCCGeomSequence geom1 = new VCCGeomSequence("geom1", vcellDim);
    model.geometrysequences.add(geom1);
    VCCMeshSequence mesh1 = new VCCMeshSequence("mesh1", geom1);
    model.meshes.add(mesh1);
    VCCStudy std1 = new VCCStudy("std1");
    model.study = std1;
    TimeBounds timeBounds = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeBounds();
    TimeStep timeStep = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeStep();
    String beginTime = Double.toString(timeBounds.getStartingTime());
    String endTime = Double.toString(timeBounds.getEndingTime());
    String step = Double.toString(timeStep.getDefaultTimeStep());
    VCCStudyFeature time = new VCCTransientStudyFeature("time", beginTime, step, endTime);
    std1.features.add(time);
    if (vcellGeometrySpec.getImage() != null) {
        throw new RuntimeException("image-based geometries not yet supported by VCell's COMSOL model builder");
    }
    if (vcellGeometrySpec.getNumSubVolumes() == 0) {
        throw new RuntimeException("no subvolumes defined in geometry");
    }
    if (vcellGeometrySpec.getNumAnalyticOrCSGSubVolumes() != vcellGeometrySpec.getNumSubVolumes()) {
        throw new RuntimeException("only analytic and CSG subvolumes currently supported by VCell's COMSOL model builder");
    }
    // 
    // add geometry for all subvolumes
    // 
    HashMap<String, VCCGeomFeature> subvolumeNameFeatureMap = new HashMap<String, VCCGeomFeature>();
    SubVolume[] subVolumes = vcellGeometrySpec.getSubVolumes();
    for (int i = 0; i < subVolumes.length; i++) {
        SubVolume subvolume = subVolumes[i];
        if (subvolume instanceof CSGObject) {
            CSGObject vcellCSGObject = (CSGObject) subvolume;
            CSGNode vcellCSGNode = vcellCSGObject.getRoot();
            ArrayList<VCCGeomFeature> geomFeatureList = new ArrayList<VCCGeomFeature>();
            VCCGeomFeature feature = csgVisitor(vcellCSGNode, geomFeatureList, subvolume.getName());
            geom1.geomfeatures.addAll(geomFeatureList);
            if (i == 0) {
                // first subvolume (on top in ordinals) doesn't need any differencing
                subvolumeNameFeatureMap.put(subvolume.getName(), feature);
            } else {
                // have to subtract union of prior subvolumes
                ArrayList<VCCGeomFeature> priorFeatures = new ArrayList<VCCGeomFeature>();
                for (int j = 0; j < i; j++) {
                    CSGObject priorCSGObject = (CSGObject) subVolumes[j];
                    CSGNode priorCSGNode = priorCSGObject.getRoot();
                    geomFeatureList.clear();
                    VCCGeomFeature priorFeature = csgVisitor(priorCSGNode, geomFeatureList, subvolume.getName());
                    priorFeatures.add(priorFeature);
                    geom1.geomfeatures.addAll(geomFeatureList);
                }
                VCCDifference diff = new VCCDifference("diff" + subvolume.getName(), Keep.off);
                diff.input.add(feature);
                diff.input2.addAll(priorFeatures);
                geom1.geomfeatures.add(diff);
                subvolumeNameFeatureMap.put(subvolume.getName(), diff);
            }
        } else {
            throw new RuntimeException("only CSG subvolumes currently supported by VCell's COMSOL model builder");
        }
    }
    // 
    // add geometry for all surfaceClasses
    // 
    HashMap<String, VCCGeomFeature> surfaceclassNameFeatureMap = new HashMap<String, VCCGeomFeature>();
    SurfaceClass[] surfaceClasses = vcellGeometry.getGeometrySurfaceDescription().getSurfaceClasses();
    for (int i = 0; i < surfaceClasses.length; i++) {
        SurfaceClass surfaceClass = surfaceClasses[i];
        Set<SubVolume> adjacentSubvolumes = surfaceClass.getAdjacentSubvolumes();
        if (adjacentSubvolumes.size() != 2) {
            throw new RuntimeException("expecting two adjacent subvolumes for surface " + surfaceClass.getName() + " in COMSOL model builder");
        }
        // find adjacent Geometry Features (for subvolumes)
        Iterator<SubVolume> svIter = adjacentSubvolumes.iterator();
        SubVolume subvolume0 = svIter.next();
        SubVolume subvolume1 = svIter.next();
        ArrayList<VCCGeomFeature> adjacentFeatures = new ArrayList<VCCGeomFeature>();
        adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume0.getName()));
        adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume1.getName()));
        String name = "inter_" + subvolume0.getName() + "_" + subvolume1.getName();
        // surfaces are dimension N-1
        int entitydim = vcellDim - 1;
        VCCIntersectionSelection intersect_subvolumes = new VCCIntersectionSelection(name, entitydim);
        intersect_subvolumes.input.addAll(adjacentFeatures);
        geom1.geomfeatures.add(intersect_subvolumes);
        surfaceclassNameFeatureMap.put(surfaceClass.getName(), intersect_subvolumes);
    }
    SimulationSymbolTable symbolTable = new SimulationSymbolTable(vcellSimJob.getSimulation(), vcellSimJob.getJobIndex());
    // 
    for (SubDomain subDomain : Collections.list(vcellMathDesc.getSubDomains())) {
        for (Equation equ : subDomain.getEquationCollection()) {
            if (equ instanceof PdeEquation || equ instanceof OdeEquation) {
                VCCGeomFeature geomFeature = null;
                final int dim;
                if (subDomain instanceof CompartmentSubDomain) {
                    geomFeature = subvolumeNameFeatureMap.get(subDomain.getName());
                    dim = vcellDim;
                } else if (subDomain instanceof MembraneSubDomain) {
                    geomFeature = surfaceclassNameFeatureMap.get(subDomain.getName());
                    dim = vcellDim - 1;
                } else {
                    throw new RuntimeException("subdomains of type '" + subDomain.getClass().getSimpleName() + "' not yet supported in COMSOL model builder");
                }
                if (geomFeature == null) {
                    throw new RuntimeException("cannot find COMSOL geometry feature named " + subDomain.getName() + " in COMSOL model builder");
                }
                VCCConvectionDiffusionEquation cdeq = new VCCConvectionDiffusionEquation("cdeq_" + equ.getVariable().getName(), geom1, geomFeature, dim);
                cdeq.fieldName = equ.getVariable().getName();
                cdeq.initial = MathUtilities.substituteModelParameters(equ.getInitialExpression(), symbolTable).flatten().infix();
                cdeq.sourceTerm_f = MathUtilities.substituteModelParameters(equ.getRateExpression(), symbolTable).flatten().infix();
                if (equ instanceof PdeEquation) {
                    PdeEquation pde = (PdeEquation) equ;
                    cdeq.diffTerm_c = MathUtilities.substituteModelParameters(pde.getDiffusionExpression(), symbolTable).flatten().infix();
                    if (subDomain instanceof CompartmentSubDomain) {
                        CompartmentSubDomain compartmentSubdomain = (CompartmentSubDomain) subDomain;
                        ArrayList<String> be = new ArrayList<String>();
                        if (pde.getVelocityX() != null) {
                            be.add(MathUtilities.substituteModelParameters(pde.getVelocityX(), symbolTable).flatten().infix());
                        } else {
                            be.add("0");
                        }
                        if (vcellDim >= 2) {
                            if (pde.getVelocityY() != null) {
                                be.add(MathUtilities.substituteModelParameters(pde.getVelocityY(), symbolTable).flatten().infix());
                            } else {
                                be.add("0");
                            }
                        }
                        if (vcellDim == 3) {
                            if (pde.getVelocityY() != null) {
                                be.add(MathUtilities.substituteModelParameters(pde.getVelocityZ(), symbolTable).flatten().infix());
                            } else {
                                be.add("0");
                            }
                        }
                        cdeq.advection_be = be.toArray(new String[vcellDim]);
                        // 
                        // look for membrane boundary conditions for this variable
                        // 
                        MembraneSubDomain[] membraneSubdomains = vcellMathDesc.getMembraneSubDomains(compartmentSubdomain);
                        for (MembraneSubDomain membraneSubdomain : membraneSubdomains) {
                            JumpCondition jumpCondition = membraneSubdomain.getJumpCondition((VolVariable) pde.getVariable());
                            if (jumpCondition != null) {
                                Expression fluxExpr = null;
                                if (membraneSubdomain.getInsideCompartment() == compartmentSubdomain) {
                                    fluxExpr = jumpCondition.getInFluxExpression();
                                } else if (membraneSubdomain.getOutsideCompartment() == compartmentSubdomain) {
                                    fluxExpr = jumpCondition.getOutFluxExpression();
                                }
                                String name = equ.getVariable().getName() + "_flux_" + membraneSubdomain.getName();
                                VCCGeomFeature selection = surfaceclassNameFeatureMap.get(membraneSubdomain.getName());
                                VCCFluxBoundary fluxBoundary = new VCCFluxBoundary(name, selection, vcellDim - 1);
                                fluxBoundary.flux_g = MathUtilities.substituteModelParameters(fluxExpr, symbolTable).flatten().infix();
                                cdeq.features.add(fluxBoundary);
                            }
                        }
                    }
                }
                model.physics.add(cdeq);
            }
        }
    }
    // 
    return model;
}

Example 14 with TimeStep

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

the class ComsolSolver method createLogFileContents.

public static String createLogFileContents(SolverTaskDescription solverTaskDescription) {
    StringBuffer buffer = new StringBuffer();
    buffer.append(SolverDataType.COMSOL.name() + "\n");
    double startingTime = solverTaskDescription.getTimeBounds().getStartingTime();
    double endingTime = solverTaskDescription.getTimeBounds().getEndingTime();
    TimeStep timeStep = solverTaskDescription.getTimeStep();
    int iteration = 0;
    double time = startingTime;
    while (time < (endingTime + timeStep.getDefaultTimeStep() / 2.0)) {
        time = startingTime + iteration * timeStep.getDefaultTimeStep();
        iteration++;
        buffer.append(time + "\n");
    }
    return buffer.toString();
}

Example 15 with TimeStep

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

the class SmoldynFileWriter method writeRuntimeCommands.

// uncomment for debug
/*private void writeGraphicsLegend() throws MathException{
	try {
		java.awt.image.BufferedImage cmapImage = new java.awt.image.BufferedImage(200, particleVariableList.size()*30,java.awt.image.BufferedImage.TYPE_INT_RGB);
		Graphics g = cmapImage.getGraphics();
		for (int i = 0; i < particleVariableList.size(); i ++) {
			Color c = colors[i];
			System.out.println("color for legend: " + "red--"+ c.getRed() + "  green--" + c.getGreen() + "  blue--" + c.getBlue());
			String variableName = getVariableName(particleVariableList.get(i),null);
			g.setColor(c);
			g.drawString(variableName, 5, 30*i + 20);
			g.fillRect(105, 30*i + 10, 20, 10);
		}
		g.dispose();
		File tmpFile = File.createTempFile("legend", ".jpg");

		FileOutputStream fios = null;
		try {
			printWriter.println("# legend file: " + tmpFile.getAbsolutePath());
			fios = new FileOutputStream(tmpFile);
			ImageIO.write(cmapImage,"jpg",fios);
		}  finally {
			if(fios != null) {fios.close();}
		}
	} catch (Exception e) {
		e.printStackTrace();
		throw new MathException(e.getMessage());
	}
}*/
private void writeRuntimeCommands() throws SolverException, DivideByZeroException, DataAccessException, IOException, MathException, ExpressionException {
    printWriter.println("# " + SmoldynVCellMapper.SmoldynKeyword.killmolincmpt + " runtime command to kill molecules misplaced during initial condtions");
    for (ParticleVariable pv : particleVariableList) {
        CompartmentSubDomain varDomain = mathDesc.getCompartmentSubDomain(pv.getDomain().getName());
        if (varDomain == null) {
            continue;
        }
        boolean bkillMol = false;
        ArrayList<ParticleInitialCondition> iniConditionList = varDomain.getParticleProperties(pv).getParticleInitialConditions();
        for (ParticleInitialCondition iniCon : iniConditionList) {
            if (iniCon instanceof ParticleInitialConditionConcentration) {
                try {
                    subsituteFlattenToConstant(((ParticleInitialConditionConcentration) iniCon).getDistribution());
                } catch (// can not be evaluated to a constant
                Exception e) {
                    bkillMol = true;
                    break;
                }
            }
        }
        if (bkillMol) {
            Enumeration<SubDomain> subDomainEnumeration = mathDesc.getSubDomains();
            while (subDomainEnumeration.hasMoreElements()) {
                SubDomain subDomain = subDomainEnumeration.nextElement();
                if (subDomain instanceof CompartmentSubDomain && varDomain != subDomain) {
                    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.B + " " + SmoldynVCellMapper.SmoldynKeyword.killmolincmpt + " " + pv.getName() + "(" + SmoldynVCellMapper.SmoldynKeyword.all + ") " + subDomain.getName());
                }
            }
        }
    }
    printWriter.println();
    // write command to kill molecules on membrane for adsortption to nothing
    printWriter.println("# kill membrane molecues that are absorbed (to nothing)");
    for (String killMolCmd : killMolCommands) {
        printWriter.println(killMolCmd);
    }
    printWriter.println();
    printWriter.println("# runtime command");
    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.E + " " + VCellSmoldynKeyword.vcellPrintProgress);
    if (outputFile != null && cartesianMesh != null) {
        OutputTimeSpec ots = simulation.getSolverTaskDescription().getOutputTimeSpec();
        if (ots.isUniform()) {
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.output_files + " " + outputFile.getName());
            ISize sampleSize = simulation.getMeshSpecification().getSamplingSize();
            TimeStep timeStep = simulation.getSolverTaskDescription().getTimeStep();
            int n = (int) Math.round(((UniformOutputTimeSpec) ots).getOutputTimeStep() / timeStep.getDefaultTimeStep());
            if (simulation.getSolverTaskDescription().getSmoldynSimulationOptions().isSaveParticleLocations()) {
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + SmoldynVCellMapper.SmoldynKeyword.incrementfile + " " + outputFile.getName());
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + SmoldynVCellMapper.SmoldynKeyword.listmols + " " + outputFile.getName());
            }
            // DON'T CHANGE THE ORDER HERE.
            // DataProcess must be before vcellWriteOutput
            writeDataProcessor();
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + VCellSmoldynKeyword.vcellWriteOutput + " begin");
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + VCellSmoldynKeyword.vcellWriteOutput + " " + VCellSmoldynKeyword.dimension + " " + dimension);
            printWriter.print(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + VCellSmoldynKeyword.vcellWriteOutput + " " + VCellSmoldynKeyword.sampleSize + " " + sampleSize.getX());
            if (dimension > 1) {
                printWriter.print(" " + sampleSize.getY());
                if (dimension > 2) {
                    printWriter.print(" " + sampleSize.getZ());
                }
            }
            printWriter.println();
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + VCellSmoldynKeyword.vcellWriteOutput + " " + VCellSmoldynKeyword.numMembraneElements + " " + cartesianMesh.getNumMembraneElements());
            for (ParticleVariable pv : particleVariableList) {
                String type = pv instanceof MembraneParticleVariable ? VCellSmoldynKeyword.membrane.name() : VCellSmoldynKeyword.volume.name();
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + VCellSmoldynKeyword.vcellWriteOutput + " " + VCellSmoldynKeyword.variable + " " + pv.getName() + " " + type + " " + pv.getDomain().getName());
            }
            printWriter.println(SmoldynVCellMapper.SmoldynKeyword.cmd + " " + SmoldynVCellMapper.SmoldynKeyword.N + " " + n + " " + VCellSmoldynKeyword.vcellWriteOutput + " end");
        } else {
            throw new SolverException(SolverDescription.Smoldyn.getDisplayLabel() + " only supports uniform output.");
        }
    }
    printWriter.println();
}