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

Example 6 with ReservedSymbol

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

the class Macroscopic_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 bindingRadiusParam = getKineticsParameterFromRole(ROLE_Binding_Radius);
    KineticsParameter kOnParam = getKineticsParameterFromRole(ROLE_KOn);
    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));
    // binding radius, computed by Kon = 2*PI*D/Ln(b/R), b = 1/sqrt(Pa*PI)
    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 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);
    // 2*PI*D/kon
    Expression exponentNumExp = Expression.div(exp2_PI_D, getSymbolExpression(kOnParam));
    // exp(-2*PI*D/Kon)
    Expression exponentExp = Expression.exp(Expression.negate(exponentNumExp));
    // b*exp(-2*PI*D/Kon)
    Expression radius = Expression.mult(b, exponentExp);
    if (bindingRadiusParam != null && radius != null) {
        bindingRadiusParam.setExpression(radius);
    }
    // SECONDARY CURRENT DENSITY
    // update from reaction rate
    updateInwardCurrentDensityFromReactionRate();
}

Example 7 with ReservedSymbol

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

the class ModelOptimizationMapping method getRemappedReferenceData.

/**
 * Gets the constraintData property (cbit.vcell.opt.ConstraintData) value.
 * @return The constraintData property value.
 * @see #setConstraintData
 */
private ReferenceData getRemappedReferenceData(MathMapping mathMapping) throws MappingException {
    if (modelOptimizationSpec.getReferenceData() == null) {
        return null;
    }
    // 
    // make sure time is mapped
    // 
    ReferenceData refData = modelOptimizationSpec.getReferenceData();
    ReferenceDataMappingSpec[] refDataMappingSpecs = modelOptimizationSpec.getReferenceDataMappingSpecs();
    RowColumnResultSet rowColResultSet = new RowColumnResultSet();
    Vector<SymbolTableEntry> modelObjectList = new Vector<SymbolTableEntry>();
    Vector<double[]> dataList = new Vector<double[]>();
    // 
    // find bound columns, (time is always mapped to the first column)
    // 
    int mappedColumnCount = 0;
    for (int i = 0; i < refDataMappingSpecs.length; i++) {
        SymbolTableEntry modelObject = refDataMappingSpecs[i].getModelObject();
        if (modelObject != null) {
            int mappedColumnIndex = mappedColumnCount;
            if (modelObject instanceof Model.ReservedSymbol && ((ReservedSymbol) modelObject).isTime()) {
                mappedColumnIndex = 0;
            }
            String origRefDataColumnName = refDataMappingSpecs[i].getReferenceDataColumnName();
            int origRefDataColumnIndex = refData.findColumn(origRefDataColumnName);
            if (origRefDataColumnIndex < 0) {
                throw new RuntimeException("reference data column named '" + origRefDataColumnName + "' not found");
            }
            double[] columnData = refData.getDataByColumn(origRefDataColumnIndex);
            if (modelObjectList.contains(modelObject)) {
                throw new RuntimeException("multiple reference data columns mapped to same model object '" + modelObject.getName() + "'");
            }
            modelObjectList.insertElementAt(modelObject, mappedColumnIndex);
            dataList.insertElementAt(columnData, mappedColumnIndex);
            mappedColumnCount++;
        }
    }
    // 
    if (modelObjectList.size() == 0) {
        throw new RuntimeException("reference data was not associated with model");
    }
    if (modelObjectList.size() == 1) {
        throw new RuntimeException("reference data was not associated with model, must map time and at least one other column");
    }
    boolean bFoundTimeVar = false;
    for (SymbolTableEntry ste : modelObjectList) {
        if (ste instanceof Model.ReservedSymbol && ((ReservedSymbol) ste).isTime()) {
            bFoundTimeVar = true;
            break;
        }
    }
    if (!bFoundTimeVar) {
        throw new RuntimeException("must map time column of reference data to model");
    }
    // 
    for (int i = 0; i < modelObjectList.size(); i++) {
        SymbolTableEntry modelObject = (SymbolTableEntry) modelObjectList.elementAt(i);
        try {
            // Find by name because MathSybolMapping has different 'objects' than refDataMapping 'objects'
            Variable variable = mathMapping.getMathSymbolMapping().findVariableByName(modelObject.getName());
            if (variable != null) {
                String symbol = variable.getName();
                rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
            } else if (modelObject instanceof Model.ReservedSymbol && ((Model.ReservedSymbol) modelObject).isTime()) {
                Model.ReservedSymbol time = (Model.ReservedSymbol) modelObject;
                String symbol = time.getName();
                rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
            }
        } catch (MathException | MatrixException | ExpressionException | ModelException e) {
            e.printStackTrace();
            throw new MappingException(e.getMessage(), e);
        }
    }
    // 
    // populate data columns (time and rest)
    // 
    double[] weights = new double[rowColResultSet.getColumnDescriptionsCount()];
    weights[0] = 1.0;
    int numRows = ((double[]) dataList.elementAt(0)).length;
    int numColumns = modelObjectList.size();
    for (int j = 0; j < numRows; j++) {
        double[] row = new double[numColumns];
        for (int i = 0; i < numColumns; i++) {
            row[i] = ((double[]) dataList.elementAt(i))[j];
            if (i > 0) {
                weights[i] += row[i] * row[i];
            }
        }
        rowColResultSet.addRow(row);
    }
    for (int i = 0; i < numColumns; i++) {
        if (weights[i] == 0) {
            weights[i] = 1;
        } else {
            weights[i] = 1 / weights[i];
        }
    }
    SimpleReferenceData remappedRefData = new SimpleReferenceData(rowColResultSet, weights);
    return remappedRefData;
}

Example 8 with ReservedSymbol

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

the class XmlReader method getKinetics.

/**
 * This method returns a Kinetics object from a XML Element based on the value of the kinetics type attribute.
 * Creation date: (3/19/2001 4:42:04 PM)
 * @return cbit.vcell.model.Kinetics
 * @param param org.jdom.Element
 */
private Kinetics getKinetics(Element param, ReactionStep reaction, Model model) throws XmlParseException {
    VariableHash varHash = new VariableHash();
    addResevedSymbols(varHash, model);
    String type = param.getAttributeValue(XMLTags.KineticsTypeAttrTag);
    Kinetics newKinetics = null;
    try {
        if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralKinetics)) {
            // create a general kinetics
            newKinetics = new GeneralKinetics(reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralCurrentKinetics)) {
            // Create GeneralCurrentKinetics
            newKinetics = new GeneralCurrentKinetics(reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeMassAction) && reaction instanceof SimpleReaction) {
            // create a Mass Action kinetics
            newKinetics = new MassActionKinetics((SimpleReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeNernst) && reaction instanceof FluxReaction) {
            // create NernstKinetics
            newKinetics = new NernstKinetics((FluxReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGHK) && reaction instanceof FluxReaction) {
            // create GHKKinetics
            newKinetics = new GHKKinetics((FluxReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeHMM_Irr) && reaction instanceof SimpleReaction) {
            // create HMM_IrrKinetics
            newKinetics = new HMM_IRRKinetics((SimpleReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeHMM_Rev) && reaction instanceof SimpleReaction) {
            // create HMM_RevKinetics
            newKinetics = new HMM_REVKinetics((SimpleReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralTotal_oldname)) {
            // create GeneralTotalKinetics
            newKinetics = new GeneralLumpedKinetics(reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralLumped)) {
            // create GeneralLumpedKinetics
            newKinetics = new GeneralLumpedKinetics(reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralCurrentLumped)) {
            // create GeneralCurrentLumpedKinetics
            newKinetics = new GeneralCurrentLumpedKinetics(reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralPermeability) && reaction instanceof FluxReaction) {
            // create GeneralPermeabilityKinetics
            newKinetics = new GeneralPermeabilityKinetics((FluxReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeMacroscopic_Irr) && reaction instanceof SimpleReaction) {
            // create Macroscopic_IRRKinetics
            newKinetics = new Macroscopic_IRRKinetics((SimpleReaction) reaction);
        } else if (type.equalsIgnoreCase(XMLTags.KineticsTypeMicroscopic_Irr) && reaction instanceof SimpleReaction) {
            // create Microscopic_IRRKinetics
            newKinetics = new Microscopic_IRRKinetics((SimpleReaction) reaction);
        } else {
            throw new XmlParseException("Unknown kinetics type: " + type);
        }
    } catch (ExpressionException e) {
        e.printStackTrace();
        throw new XmlParseException("Error creating the kinetics for reaction: " + reaction.getName(), e);
    }
    try {
        // transaction begin flag ... yeah, this is a hack
        newKinetics.reading(true);
        // Read all of the parameters
        List<Element> list = param.getChildren(XMLTags.ParameterTag, vcNamespace);
        // add constants that may be used in kinetics.
        // VariableHash varHash = getVariablesHash();
        ArrayList<String> reserved = new ArrayList<String>();
        ReservedSymbol[] reservedSymbols = reaction.getModel().getReservedSymbols();
        for (ReservedSymbol rs : reservedSymbols) {
            reserved.add(rs.getName());
        }
        try {
            if (reaction.getStructure() instanceof Membrane) {
                Membrane membrane = (Membrane) reaction.getStructure();
                varHash.addVariable(new Constant(membrane.getMembraneVoltage().getName(), new Expression(0.0)));
                reserved.add(membrane.getMembraneVoltage().getName());
            }
            // 
            // add Reactants, Products, and Catalysts (ReactionParticipants)
            // 
            ReactionParticipant[] rp = reaction.getReactionParticipants();
            for (int i = 0; i < rp.length; i++) {
                varHash.addVariable(new Constant(rp[i].getName(), new Expression(0.0)));
            }
        } catch (MathException e) {
            e.printStackTrace(System.out);
            throw new XmlParseException("error reordering parameters according to dependencies: ", e);
        }
        // 
        for (Element xmlParam : list) {
            String paramName = unMangle(xmlParam.getAttributeValue(XMLTags.NameAttrTag));
            String role = xmlParam.getAttributeValue(XMLTags.ParamRoleAttrTag);
            String paramExpStr = xmlParam.getText();
            Expression paramExp = unMangleExpression(paramExpStr);
            try {
                if (varHash.getVariable(paramName) == null) {
                    varHash.addVariable(new Function(paramName, paramExp, null));
                } else {
                    if (reserved.contains(paramName)) {
                        varHash.removeVariable(paramName);
                        varHash.addVariable(new Function(paramName, paramExp, null));
                    }
                }
            } catch (MathException e) {
                e.printStackTrace(System.out);
                throw new XmlParseException("error reordering parameters according to dependencies: ", e);
            }
            Kinetics.KineticsParameter tempParam = null;
            if (!role.equals(XMLTags.ParamRoleUserDefinedTag)) {
                tempParam = newKinetics.getKineticsParameterFromRole(Kinetics.getParamRoleFromDefaultDesc(role));
            } else {
                continue;
            }
            // hack for bringing in General Total kinetics without breaking.
            if (tempParam == null && newKinetics instanceof GeneralLumpedKinetics) {
                if (role.equals(Kinetics.GTK_AssumedCompartmentSize_oldname) || role.equals(Kinetics.GTK_ReactionRate_oldname) || role.equals(Kinetics.GTK_CurrentDensity_oldname)) {
                    continue;
                } else if (role.equals(VCMODL.TotalRate_oldname)) {
                    tempParam = newKinetics.getKineticsParameterFromRole(Kinetics.ROLE_LumpedReactionRate);
                }
            }
            // hack from bringing in chargeValence parameters without breaking
            if (tempParam == null && Kinetics.getParamRoleFromDefaultDesc(role) == Kinetics.ROLE_ChargeValence) {
                tempParam = newKinetics.getChargeValenceParameter();
            }
            if (tempParam == null) {
                throw new XmlParseException("parameter with role '" + role + "' not found in kinetics type '" + type + "'");
            }
            // 
            if (!tempParam.getName().equals(paramName)) {
                Kinetics.KineticsParameter multNameParam = newKinetics.getKineticsParameter(paramName);
                int n = 0;
                while (multNameParam != null) {
                    String tempName = paramName + "_" + n++;
                    newKinetics.renameParameter(paramName, tempName);
                    multNameParam = newKinetics.getKineticsParameter(tempName);
                }
                newKinetics.renameParameter(tempParam.getName(), paramName);
            }
        }
        // 
        // create unresolved parameters for all unresolved symbols
        // 
        String unresolvedSymbol = varHash.getFirstUnresolvedSymbol();
        while (unresolvedSymbol != null) {
            try {
                // will turn into an UnresolvedParameter.
                varHash.addVariable(new Function(unresolvedSymbol, new Expression(0.0), null));
            } catch (MathException e) {
                e.printStackTrace(System.out);
                throw new XmlParseException(e);
            }
            newKinetics.addUnresolvedParameter(unresolvedSymbol);
            unresolvedSymbol = varHash.getFirstUnresolvedSymbol();
        }
        Variable[] sortedVariables = varHash.getTopologicallyReorderedVariables();
        ModelUnitSystem modelUnitSystem = reaction.getModel().getUnitSystem();
        for (int i = sortedVariables.length - 1; i >= 0; i--) {
            if (sortedVariables[i] instanceof Function) {
                Function paramFunction = (Function) sortedVariables[i];
                Element xmlParam = null;
                for (int j = 0; j < list.size(); j++) {
                    Element tempParam = (Element) list.get(j);
                    if (paramFunction.getName().equals(unMangle(tempParam.getAttributeValue(XMLTags.NameAttrTag)))) {
                        xmlParam = tempParam;
                        break;
                    }
                }
                if (xmlParam == null) {
                    // must have been an unresolved parameter
                    continue;
                }
                String symbol = xmlParam.getAttributeValue(XMLTags.VCUnitDefinitionAttrTag);
                VCUnitDefinition unit = null;
                if (symbol != null) {
                    unit = modelUnitSystem.getInstance(symbol);
                }
                Kinetics.KineticsParameter tempParam = newKinetics.getKineticsParameter(paramFunction.getName());
                if (tempParam == null) {
                    newKinetics.addUserDefinedKineticsParameter(paramFunction.getName(), paramFunction.getExpression(), unit);
                } else {
                    newKinetics.setParameterValue(tempParam, paramFunction.getExpression());
                    tempParam.setUnitDefinition(unit);
                }
            }
        }
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace(System.out);
        throw new XmlParseException("Exception while setting parameters for Reaction : " + reaction.getName(), e);
    } catch (ExpressionException e) {
        e.printStackTrace(System.out);
        throw new XmlParseException("Exception while settings parameters for Reaction : " + reaction.getName(), e);
    } finally {
        newKinetics.reading(false);
    }
    return newKinetics;
}

Example 9 with ReservedSymbol

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

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

the class TriggerTemplatePanel method init.

public void init(SimulationContext simulationContext, AutoCompleteSymbolFilter autoCompleteSymbolFilter, BioEvent existingBioEvent) {
    if (simulationContext != null) {
        // mathOpComboBox.removeAllItems();
        // mathOpComboBox.addItem(BioEvent.TriggerComparison.greaterThan);
        // mathOpComboBox.addItem(BioEvent.TriggerComparison.lessThan);
        // mathOpComboBox.addItem(BioEvent.TriggerComparison.greaterThanOrEqual);
        // mathOpComboBox.addItem(BioEvent.TriggerComparison.lessThanOrEqual);
        varAboveComboBox.removeAllItems();
        varBelowComboBox.removeAllItems();
        SymbolTableEntry[] triggerExprSymbols = ModelOptimizationSpec.calculateTimeDependentModelObjects(simulationContext);
        Arrays.sort(triggerExprSymbols, new Comparator<SymbolTableEntry>() {

            @Override
            public int compare(SymbolTableEntry o1, SymbolTableEntry o2) {
                return o1.getName().compareToIgnoreCase(o2.getName());
            }
        });
        for (SymbolTableEntry ste : triggerExprSymbols) {
            if (!(ste instanceof ReservedSymbol)) {
                varAboveComboBox.addItem(ste);
                varBelowComboBox.addItem(ste);
            }
        }
        // setVerifier(simulationContext);
        Map<String, SymbolTableEntry> entryMap = new HashMap<String, SymbolTableEntry>();
        simulationContext.getEntries(entryMap);
        textFieldAboveVarVal.setAutoCompleteSymbolFilter(autoCompleteSymbolFilter);
        textFieldAboveVarVal.setAutoCompletionWords(entryMap.keySet());
        textFieldBelowVarVal.setAutoCompleteSymbolFilter(autoCompleteSymbolFilter);
        textFieldBelowVarVal.setAutoCompletionWords(entryMap.keySet());
        textFieldSingleTime.setAutoCompleteSymbolFilter(autoCompleteSymbolFilter);
        textFieldSingleTime.setAutoCompletionWords(entryMap.keySet());
        textFieldGeneral.setAutoCompleteSymbolFilter(autoCompleteSymbolFilter);
        textFieldGeneral.setAutoCompletionWords(entryMap.keySet());
        if (existingBioEvent != null) {
            // Edit Trigger mode
            switch(existingBioEvent.getTriggerType()) {
                case GeneralTrigger:
                    {
                        textFieldGeneral.setText(existingBioEvent.getParameter(BioEventParameterType.GeneralTriggerFunction).getExpression().infix());
                        rdbtnGeneral.setSelected(true);
                        break;
                    }
                case LogRangeTimes:
                    {
                        LocalParameter rangeMinParam = existingBioEvent.getParameter(BioEventParameterType.RangeMinTime);
                        if (rangeMinParam != null && rangeMinParam.getExpression() != null) {
                            textFieldLogMin.setText(rangeMinParam.getExpression().infix());
                        } else {
                            textFieldLogMin.setText("");
                        }
                        LocalParameter rangeMaxParam = existingBioEvent.getParameter(BioEventParameterType.RangeMaxTime);
                        if (rangeMaxParam != null && rangeMaxParam.getExpression() != null) {
                            textFieldLogMax.setText(rangeMaxParam.getExpression().infix());
                        } else {
                            textFieldLogMax.setText("");
                        }
                        LocalParameter rangeNumParam = existingBioEvent.getParameter(BioEventParameterType.RangeNumTimes);
                        if (rangeNumParam != null && rangeNumParam.getExpression() != null) {
                            textFieldLogNum.setText(rangeNumParam.getExpression().infix());
                        } else {
                            textFieldLogNum.setText("");
                        }
                        rdbtnLogTimeRange.setSelected(true);
                        break;
                    }
                case LinearRangeTimes:
                    {
                        LocalParameter rangeMinParam = existingBioEvent.getParameter(BioEventParameterType.RangeMinTime);
                        if (rangeMinParam != null && rangeMinParam.getExpression() != null) {
                            textFieldLinearMin.setText(rangeMinParam.getExpression().infix());
                        } else {
                            textFieldLinearMin.setText("");
                        }
                        LocalParameter rangeMaxParam = existingBioEvent.getParameter(BioEventParameterType.RangeMaxTime);
                        if (rangeMaxParam != null && rangeMaxParam.getExpression() != null) {
                            textFieldLinearMax.setText(rangeMaxParam.getExpression().infix());
                        } else {
                            textFieldLinearMax.setText("");
                        }
                        LocalParameter rangeNumParam = existingBioEvent.getParameter(BioEventParameterType.RangeNumTimes);
                        if (rangeNumParam != null && rangeNumParam.getExpression() != null) {
                            textFieldLinearNum.setText(rangeNumParam.getExpression().infix());
                        } else {
                            textFieldLinearNum.setText("");
                        }
                        rdbtnLinearTimeRange.setSelected(true);
                        break;
                    }
                case ListOfTimes:
                    {
                        StringBuffer sb = new StringBuffer();
                        for (LocalParameter p : existingBioEvent.getEventParameters()) {
                            if (p.getRole() == BioEventParameterType.TimeListItem) {
                                if (sb.length() > 0) {
                                    sb.append(",");
                                }
                                sb.append(p.getExpression().infix());
                            }
                        }
                        textFieldListTimes.setText(sb.toString());
                        rdbtnListTimes.setSelected(true);
                        break;
                    }
                case ObservableAboveThreshold:
                    {
                        LocalParameter observableParameter = existingBioEvent.getParameter(BioEventParameterType.Observable);
                        SymbolTableEntry ste = null;
                        if (observableParameter != null && observableParameter.getExpression() != null && observableParameter.getExpression().isIdentifier()) {
                            ste = observableParameter.getExpression().getSymbolBinding(observableParameter.getExpression().getSymbols()[0]);
                            varAboveComboBox.setSelectedItem(ste);
                        }
                        textFieldAboveVarVal.setText(existingBioEvent.getParameter(BioEventParameterType.Threshold).getExpression().infix());
                        rdbtnVarAboveThreshold.setSelected(true);
                        break;
                    }
                case ObservableBelowThreshold:
                    {
                        LocalParameter observableParameter = existingBioEvent.getParameter(BioEventParameterType.Observable);
                        SymbolTableEntry ste = null;
                        if (observableParameter != null && observableParameter.getExpression() != null && observableParameter.getExpression().isIdentifier()) {
                            ste = observableParameter.getExpression().getSymbolBinding(observableParameter.getExpression().getSymbols()[0]);
                            varBelowComboBox.setSelectedItem(ste);
                        }
                        textFieldBelowVarVal.setText(existingBioEvent.getParameter(BioEventParameterType.Threshold).getExpression().infix());
                        rdbtnVarBelowThreshold.setSelected(true);
                        break;
                    }
                case SingleTriggerTime:
                    {
                        LocalParameter parameter = existingBioEvent.getParameter(BioEventParameterType.SingleTriggerTime);
                        textFieldSingleTime.setText(parameter.getExpression().infix());
                        rdbtnSingleTime.setSelected(true);
                        break;
                    }
            }
        }
    }
}