Examples with ReservedSymbol cbit.vcell.model.Model.ReservedSymbol used on opensource projects

Example 11 with ReservedSymbol

use of cbit.vcell.model.Model.ReservedSymbol 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 12 with ReservedSymbol

use of cbit.vcell.model.Model.ReservedSymbol in project vcell by virtualcell.

the class SimulationRepresentation method getParameters.

private static ParameterRepresentation[] getParameters(BioModel bioModel, SimulationRep simulationRep) {
    SimulationContext simContext = null;
    for (SimulationContext sc : bioModel.getSimulationContexts()) {
        if (sc.getMathDescription().getKey().equals(simulationRep.getMathKey())) {
            simContext = sc;
            break;
        }
    }
    if (simContext == null) {
        return null;
    }
    // initialize to old mathDescription in case error generating math
    MathDescription mathDesc = simContext.getMathDescription();
    MathMapping mathMapping = simContext.createNewMathMapping();
    MathSymbolMapping mathSymbolMapping = null;
    try {
        mathDesc = mathMapping.getMathDescription();
        mathSymbolMapping = mathMapping.getMathSymbolMapping();
    } catch (Exception e1) {
        System.err.println(e1.getMessage());
    }
    ArrayList<ParameterRepresentation> parameterReps = new ArrayList<ParameterRepresentation>();
    Enumeration<Constant> enumMath = mathDesc.getConstants();
    while (enumMath.hasMoreElements()) {
        Constant constant = enumMath.nextElement();
        if (constant.getExpression().isNumeric()) {
            SymbolTableEntry biologicalSymbolTableEntry = null;
            if (mathSymbolMapping != null) {
                SymbolTableEntry[] stes = mathSymbolMapping.getBiologicalSymbol(constant);
                if (stes != null && stes.length >= 1) {
                    biologicalSymbolTableEntry = stes[0];
                }
            }
            if (biologicalSymbolTableEntry instanceof ReservedSymbol) {
                continue;
            }
            try {
                parameterReps.add(new ParameterRepresentation(constant.getName(), constant.getExpression().evaluateConstant(), biologicalSymbolTableEntry));
            } catch (ExpressionException e) {
                // can't happen, because constant expression is numeric
                e.printStackTrace();
            }
        }
    }
    return parameterReps.toArray(new ParameterRepresentation[0]);
}

Example 13 with ReservedSymbol

use of cbit.vcell.model.Model.ReservedSymbol in project vcell by virtualcell.

the class Microscopic_IRRKinetics method updateGeneratedExpressions.

/**
 * Insert the method's description here.
 * Creation date: (10/19/2003 12:05:14 AM)
 * @exception cbit.vcell.parser.ExpressionException The exception description.
 */
protected void updateGeneratedExpressions() throws cbit.vcell.parser.ExpressionException, PropertyVetoException {
    KineticsParameter rateParm = getKineticsParameterFromRole(ROLE_ReactionRate);
    KineticsParameter currentParm = getKineticsParameterFromRole(ROLE_CurrentDensity);
    KineticsParameter kOnParam = getKineticsParameterFromRole(ROLE_KOn);
    KineticsParameter bindingRadiusParam = getKineticsParameterFromRole(ROLE_Binding_Radius);
    KineticsParameter diff_react1Param = getKineticsParameterFromRole(ROLE_Diffusion_Reactant1);
    KineticsParameter diff_react2Param = getKineticsParameterFromRole(ROLE_Diffusion_Reactant2);
    KineticsParameter conc_react1Param = getKineticsParameterFromRole(ROLE_Concentration_Reactant1);
    KineticsParameter conc_react2Param = getKineticsParameterFromRole(ROLE_Concentration_Reactant2);
    if (currentParm == null && rateParm == null) {
        return;
    }
    // rate prameter expr.
    ReactionParticipant[] rp_Array = getReactionStep().getReactionParticipants();
    Expression kOn_exp = getSymbolExpression(kOnParam);
    Expression newRateExp = null;
    int reactantCount = 0;
    for (int i = 0; i < rp_Array.length; i++) {
        Expression term = null;
        Expression speciesContext = getSymbolExpression(rp_Array[i].getSpeciesContext());
        int stoichiometry = rp_Array[i].getStoichiometry();
        if (rp_Array[i] instanceof Reactant) {
            reactantCount++;
            if (stoichiometry < 1) {
                throw new ExpressionException("reactant must have stoichiometry of at least 1");
            } else if (stoichiometry == 1) {
                term = speciesContext;
            } else {
                term = Expression.power(speciesContext, new Expression(stoichiometry));
            }
            kOn_exp = Expression.mult(kOn_exp, term);
        }
    }
    if (reactantCount > 0) {
        newRateExp = kOn_exp;
    } else {
        newRateExp = new Expression(0.0);
    }
    rateParm.setExpression(newRateExp);
    // current Parameter. set to 0??
    currentParm.setExpression(new Expression(0.0));
    // Kon parameter
    ReservedSymbol pi_ReservedSymbol = getReactionStep().getModel().getPI_CONSTANT();
    Expression Pa = Expression.max(getSymbolExpression(conc_react1Param), getSymbolExpression(conc_react2Param));
    // sqrt(Pa*PI)
    Expression sqrt_Pa_PI = Expression.sqrt(Expression.mult(Pa, getSymbolExpression(pi_ReservedSymbol)));
    // 1/sqrt(Pa*PI)
    Expression b = Expression.div(new Expression(1), sqrt_Pa_PI);
    Expression ln_b = Expression.log(b);
    Expression ln_Radius = Expression.log(getSymbolExpression(bindingRadiusParam));
    Expression sumD = Expression.add(getSymbolExpression(diff_react1Param), getSymbolExpression(diff_react2Param));
    // 2*PI*D
    Expression exp2_PI_D = Expression.mult(new Expression(2.0), getSymbolExpression(pi_ReservedSymbol), sumD);
    // Lnb-LnR
    Expression expLnb_LnR = Expression.add(ln_b, Expression.negate(ln_Radius));
    // Kon = 2*PI*D/Ln(b/R)
    Expression kOnExp = Expression.div(exp2_PI_D, expLnb_LnR);
    if (kOnParam != null && kOn_exp != null) {
        kOnParam.setExpression(kOnExp);
    }
    // SECONDARY CURRENT DENSITY
    // update from reaction rate
    updateInwardCurrentDensityFromReactionRate();
}