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

Example 21 with SimulationSymbolTable

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

the class FiniteVolumeFileWriter method writeVariables.

/**
 *# Variables : type name unit time_dependent_flag advection_flag solve_whole_mesh_flag solve_regions
 *VARIABLE_BEGIN
 *VOLUME_ODE rB uM
 *VOLUME_PDE rf uM false false
 *VOLUME_PDE r uM false false
 *VOLUME_ODE rfB uM
 *VARIABLE_END
 * @throws MathException
 * @throws ExpressionException
 * @throws IOException
 */
private void writeVariables() throws MathException, ExpressionException, IOException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    printWriter.println("# Variables : type name domain time_dependent_flag advection_flag grad_flag solve_whole_mesh_flag solve_regions");
    printWriter.println(FVInputFileKeyword.VARIABLE_BEGIN);
    MathDescription mathDesc = simSymbolTable.getSimulation().getMathDescription();
    Variable[] vars = simSymbolTable.getVariables();
    ArrayList<RandomVariable> rvList = new ArrayList<RandomVariable>();
    for (int i = 0; i < vars.length; i++) {
        String varName = vars[i].getName();
        String domainName = vars[i].getDomain() == null ? null : vars[i].getDomain().getName();
        if (vars[i] instanceof VolumeRandomVariable || vars[i] instanceof MembraneRandomVariable) {
            rvList.add((RandomVariable) vars[i]);
        } else if (vars[i] instanceof VolVariable) {
            if (bChomboSolver && domainName == null) {
                throw new MathException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " requires that every variable is defined in a single domain");
            }
            VolVariable volVar = (VolVariable) vars[i];
            if (mathDesc.isPDE(volVar)) {
                boolean hasTimeVaryingDiffusionOrAdvection = simSymbolTable.hasTimeVaryingDiffusionOrAdvection(volVar);
                final boolean hasVelocity = mathDesc.hasVelocity(volVar);
                final boolean hasGradient = mathDesc.hasGradient(volVar);
                if (mathDesc.isPdeSteady(volVar)) {
                    printWriter.print("VOLUME_PDE_STEADY ");
                } else {
                    printWriter.print("VOLUME_PDE ");
                }
                printWriter.print(varName + " " + domainName + " " + hasTimeVaryingDiffusionOrAdvection + " " + hasVelocity + " " + hasGradient);
            } else {
                printWriter.print("VOLUME_ODE " + varName + " " + domainName);
            }
            if (domainName == null) {
                Vector<SubDomain> listOfSubDomains = new Vector<SubDomain>();
                int totalNumCompartments = 0;
                Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
                while (subDomainEnum.hasMoreElements()) {
                    SubDomain subDomain = subDomainEnum.nextElement();
                    if (subDomain instanceof CompartmentSubDomain) {
                        CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
                        totalNumCompartments++;
                        Equation varEquation = subDomain.getEquation(vars[i]);
                        if (varEquation != null) {
                            if (!(varEquation instanceof PdeEquation) || !((PdeEquation) varEquation).isDummy(simSymbolTable, compartmentSubDomain)) {
                                listOfSubDomains.add(compartmentSubDomain);
                            }
                        }
                    }
                }
                if ((totalNumCompartments == listOfSubDomains.size()) || (listOfSubDomains.size() == 0 && simTask.getSimulation().getSolverTaskDescription().getSolverDescription().equals(SolverDescription.SundialsPDE))) {
                    printWriter.print(" true");
                } else {
                    printWriter.print(" false");
                    for (int j = 0; j < listOfSubDomains.size(); j++) {
                        CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) listOfSubDomains.elementAt(j);
                        printWriter.print(" " + compartmentSubDomain.getName());
                    }
                }
                printWriter.println();
            } else {
                printWriter.println(" false " + domainName);
            }
        } else if (vars[i] instanceof VolumeParticleVariable) {
            printWriter.println(FVInputFileKeyword.VOLUME_PARTICLE + " " + varName + " " + domainName);
        } else if (vars[i] instanceof MembraneParticleVariable) {
            printWriter.println(FVInputFileKeyword.MEMBRANE_PARTICLE + " " + varName + " " + domainName);
        } else if (vars[i] instanceof VolumeRegionVariable) {
            printWriter.println("VOLUME_REGION " + varName + " " + domainName);
        } else if (vars[i] instanceof MemVariable) {
            if (bChomboSolver && domainName == null) {
                throw new MathException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " requires that every variable is defined in a single domain");
            }
            MemVariable memVar = (MemVariable) vars[i];
            if (mathDesc.isPDE(memVar)) {
                printWriter.println("MEMBRANE_PDE " + varName + " " + domainName + " " + simSymbolTable.hasTimeVaryingDiffusionOrAdvection(memVar));
            } else {
                printWriter.println("MEMBRANE_ODE " + varName + " " + domainName);
            }
        } else if (vars[i] instanceof MembraneRegionVariable) {
            printWriter.println("MEMBRANE_REGION " + varName + " " + domainName);
        } else if (vars[i] instanceof FilamentVariable) {
            throw new RuntimeException("Filament application not supported yet");
        }
    }
    int numRandomVariables = rvList.size();
    if (numRandomVariables > 0) {
        ISize samplingSize = simTask.getSimulation().getMeshSpecification().getSamplingSize();
        String[] varNameArr = new String[numRandomVariables];
        VariableType[] varTypeArr = new VariableType[numRandomVariables];
        double[][] dataArr = new double[numRandomVariables][];
        for (int i = 0; i < numRandomVariables; i++) {
            RandomVariable rv = rvList.get(i);
            varNameArr[i] = rv.getName();
            int numRandomNumbers = 0;
            if (rv instanceof VolumeRandomVariable) {
                printWriter.print("VOLUME_RANDOM");
                varTypeArr[i] = VariableType.VOLUME;
                numRandomNumbers = samplingSize.getXYZ();
            } else if (rv instanceof MembraneRandomVariable) {
                printWriter.print("MEMBRANE_RANDOM");
                varTypeArr[i] = VariableType.MEMBRANE;
                numRandomNumbers = resampledGeometry.getGeometrySurfaceDescription().getSurfaceCollection().getTotalPolygonCount();
            } else {
                throw new RuntimeException("Unknown RandomVariable type");
            }
            printWriter.println(" " + varNameArr[i]);
            dataArr[i] = generateRandomNumbers(rv, numRandomNumbers);
        }
        File rvFile = new File(workingDirectory, simTask.getSimulationJobID() + RANDOM_VARIABLE_FILE_EXTENSION);
        DataSet.writeNew(rvFile, varNameArr, varTypeArr, samplingSize, dataArr);
    }
    printWriter.println(FVInputFileKeyword.VARIABLE_END);
    printWriter.println();
}

Example 22 with SimulationSymbolTable

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

the class OptXmlWriter method getObjectiveFunctionXML.

public static Element getObjectiveFunctionXML(OptimizationSpec optimizationSpec) {
    Element objectiveFunctionElement = new Element(OptXmlTags.ObjectiveFunction_Tag);
    ObjectiveFunction objectiveFunction = optimizationSpec.getObjectiveFunction();
    if (objectiveFunction instanceof ExplicitObjectiveFunction) {
        ExplicitObjectiveFunction explicitObjectiveFunction = (ExplicitObjectiveFunction) objectiveFunction;
        objectiveFunctionElement.setAttribute(OptXmlTags.ObjectiveFunctionType_Attr, OptXmlTags.ObjectiveFunctionType_Attr_Explicit);
        Element expressionElement = new Element(OptXmlTags.Expression_Tag);
        expressionElement.addContent(explicitObjectiveFunction.getExpression().infix());
        objectiveFunctionElement.addContent(expressionElement);
    } else if (objectiveFunction instanceof ExplicitFitObjectiveFunction) {
        ExplicitFitObjectiveFunction explicitFitObjectiveFunction = (ExplicitFitObjectiveFunction) objectiveFunction;
        objectiveFunctionElement.setAttribute(OptXmlTags.ObjectiveFunctionType_Attr, OptXmlTags.ObjectiveFunctionType_Attr_ExplicitFit);
        // add expression list
        Element expressionListElement = new Element(OptXmlTags.ExpressionList_Tag);
        objectiveFunctionElement.addContent(expressionListElement);
        ExplicitFitObjectiveFunction.ExpressionDataPair[] expDataPairs = explicitFitObjectiveFunction.getExpressionDataPairs();
        for (int i = 0; i < expDataPairs.length; i++) {
            Element expressionElement = new Element(OptXmlTags.Expression_Tag);
            expressionElement.setAttribute(OptXmlTags.ExpRefDataColumnIndex_Attr, expDataPairs[i].getReferenceDataIndex() + "");
            expressionElement.addContent(expDataPairs[i].getFitFunction().infix());
            expressionListElement.addContent(expressionElement);
        }
        if (explicitFitObjectiveFunction.getReferenceData() instanceof SimpleReferenceData) {
            Element dataElement = getDataXML(explicitFitObjectiveFunction.getReferenceData());
            objectiveFunctionElement.addContent(dataElement);
        } else {
            throw new OptimizationException("Optimization XML Writer exports SimpleReferenceData only.");
        }
    } else if (objectiveFunction instanceof PdeObjectiveFunction) {
        PdeObjectiveFunction pdeObjectiveFunction = (PdeObjectiveFunction) objectiveFunction;
        objectiveFunctionElement.setAttribute(OptXmlTags.ObjectiveFunctionType_Attr, OptXmlTags.ObjectiveFunctionType_Attr_PDE);
        SpatialReferenceData refData = pdeObjectiveFunction.getReferenceData();
        Element dataElement = getDataXML(refData);
        objectiveFunctionElement.addContent(dataElement);
        Element modelElement = getModelXML(pdeObjectiveFunction, optimizationSpec.getParameterNames());
        objectiveFunctionElement.addContent(modelElement);
        Simulation tempSim = new Simulation(pdeObjectiveFunction.getMathDescription());
        SimulationSymbolTable simSymbolTable = new SimulationSymbolTable(tempSim, 0);
        // 
        for (int i = 1; i < refData.getNumVariables(); i++) {
            Element modelMappingElement = new Element(OptXmlTags.ModelMapping_Tag);
            modelMappingElement.setAttribute(OptXmlTags.ModelMappingDataColumn_Attr, refData.getVariableNames()[i]);
            modelMappingElement.setAttribute(OptXmlTags.ModelMappingWeight_Attr, String.valueOf(refData.getVariableWeights()[i]));
            Expression mappingExpression = null;
            try {
                Variable var = pdeObjectiveFunction.getMathDescription().getVariable(refData.getVariableNames()[i]);
                if (var instanceof Function) {
                    Expression exp = new Expression(var.getExpression());
                    exp.bindExpression(simSymbolTable);
                    exp = MathUtilities.substituteFunctions(exp, simSymbolTable);
                    mappingExpression = exp.flatten();
                } else {
                    mappingExpression = new Expression(var.getName());
                }
            } catch (ExpressionException e) {
                e.printStackTrace();
                throw new OptimizationException(e.getMessage());
            }
            modelMappingElement.addContent(mappingExpression.infix());
            objectiveFunctionElement.addContent(modelMappingElement);
        }
    }
    return objectiveFunctionElement;
}

Example 23 with SimulationSymbolTable

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

the class GibsonSolver method getStochSolverResultSet.

/**
 * Get data columns and function columns to be ready for plot.
 * Creation date: (8/15/2006 11:36:43 AM)
 * @return cbit.vcell.solver.stoch.StochSolverResultSet
 */
public ODESolverResultSet getStochSolverResultSet() {
    // read .stoch file, this funciton here equals to getODESolverRestultSet()+getStateVariableResultSet()  in ODE.
    ODESolverResultSet stSolverResultSet = new ODESolverResultSet();
    FileInputStream inputStream = null;
    try {
        inputStream = new FileInputStream(getBaseName() + IDA_DATA_EXTENSION);
        InputStreamReader inputStreamReader = new InputStreamReader(inputStream);
        BufferedReader bufferedReader = new BufferedReader(inputStreamReader);
        // Read header
        String line = bufferedReader.readLine();
        if (line == null) {
            // throw exception
            System.out.println("There is no data in output file!");
            return null;
        }
        while (line.indexOf(':') > 0) {
            String name = line.substring(0, line.indexOf(':'));
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(name));
            line = line.substring(line.indexOf(':') + 1);
        }
        // Read data
        while ((line = bufferedReader.readLine()) != null) {
            line = line + "\t";
            double[] values = new double[stSolverResultSet.getDataColumnCount()];
            boolean bCompleteRow = true;
            for (int i = 0; i < stSolverResultSet.getDataColumnCount(); i++) {
                if (line.indexOf('\t') == -1) {
                    bCompleteRow = false;
                    break;
                } else {
                    String value = line.substring(0, line.indexOf('\t')).trim();
                    values[i] = Double.valueOf(value).doubleValue();
                    line = line.substring(line.indexOf('\t') + 1);
                }
            }
            if (bCompleteRow) {
                stSolverResultSet.addRow(values);
            } else {
                break;
            }
        }
    // 
    } catch (Exception e) {
        e.printStackTrace(System.out);
        return null;
    } finally {
        try {
            if (inputStream != null) {
                inputStream.close();
            }
        } catch (Exception ex) {
            lg.error(ex.getMessage(), ex);
        }
    }
    /*
	Add appropriate Function columns to result set if the stochastic simulation is to display the trajectory.
	No function columns for the results of multiple stochastic trials
	*/
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (simSymbolTable.getSimulation().getSolverTaskDescription().getStochOpt().getNumOfTrials() == 1) {
        Function[] functions = simSymbolTable.getFunctions();
        for (int i = 0; i < functions.length; i++) {
            if (SimulationSymbolTable.isFunctionSaved(functions[i])) {
                Expression exp1 = new Expression(functions[i].getExpression());
                try {
                    exp1 = simSymbolTable.substituteFunctions(exp1);
                } catch (MathException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
                }
                try {
                    FunctionColumnDescription cd = new FunctionColumnDescription(exp1.flatten(), functions[i].getName(), null, functions[i].getName(), false);
                    stSolverResultSet.addFunctionColumn(cd);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                }
            }
        }
    }
    return stSolverResultSet;
}

Example 24 with SimulationSymbolTable

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

the class NetCDFWriter method writeHybridInputFile.

/**
 * Write the model to a NetCDF file which serves as an input for stoch hybrid simulator.
 * To write to a NetCDF file is a bit complicated. First, we have to create a NetCDF-3
 * file. And then feed in the data.
 * Creation date: (5/22/2007 5:36:03 PM)
 */
public void writeHybridInputFile(String[] parameterNames) throws Exception, cbit.vcell.parser.ExpressionException, IOException, MathException, InvalidRangeException {
    Simulation simulation = simTask.getSimulation();
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (initialize()) {
        // we need to get model and control information first
        NetcdfFileWriteable ncfile = NetcdfFileWriteable.createNew(filename, false);
        // Model info. will be extracted from subDomain of mathDescription
        java.util.Enumeration<SubDomain> e = simulation.getMathDescription().getSubDomains();
        // remember we are dealing with compartmental model here. only 1 subdomain.
        SubDomain subDomain = e.nextElement();
        JumpProcess[] reactions = (JumpProcess[]) subDomain.getJumpProcesses().toArray(new JumpProcess[subDomain.getJumpProcesses().size()]);
        // get species variable names
        Variable[] variables = simSymbolTable.getVariables();
        String[] speciesNames = new String[variables.length];
        for (int i = 0; i < variables.length; i++) speciesNames[i] = variables[i].getName();
        // the probabilities for reactions
        Expression[] probs = new Expression[reactions.length];
        for (int i = 0; i < reactions.length; i++) {
            probs[i] = simSymbolTable.substituteFunctions(reactions[i].getProbabilityRate());
            probs[i] = probs[i].flatten();
        }
        VarIniCondition[] varInis = (VarIniCondition[]) subDomain.getVarIniConditions().toArray(new VarIniCondition[subDomain.getVarIniConditions().size()]);
        // the non-constant stoch variables
        Vector<Variable> vars = new Vector<Variable>();
        for (int i = 0; i < varInis.length; i++) {
            if (varInis[i].getVar() instanceof StochVolVariable) {
                vars.addElement(varInis[i].getVar());
            }
        }
        // get reaction rate law types and rate constants
        ReactionRateLaw[] reactionRateLaws = getReactionRateLaws(probs);
        SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
        TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
        UniformOutputTimeSpec timeSpec = (UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec();
        UniformOutputTimeSpec outputTimeSpec = ((UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec());
        NonspatialStochSimOptions stochOpt = solverTaskDescription.getStochOpt();
        // create an empty NetCDF-3 file
        // define dimensions
        /* these sizes must match the buffers allocated in corresponding Fortran code -- see globalvariables.f90
			in numerics Hy3S/src directory */
        Dimension numTrial = ncfile.addDimension("NumTrials", (int) stochOpt.getNumOfTrials());
        Dimension numSpecies = ncfile.addDimension("NumSpecies", vars.size());
        Dimension numReactions = ncfile.addDimension("NumReactions", subDomain.getJumpProcesses().size());
        int outPoints = ((int) ((timeBounds.getEndingTime() - timeBounds.getStartingTime()) / outputTimeSpec.getOutputTimeStep())) + 1;
        Dimension numTimePoints = ncfile.addDimension("NumTimePoints", outPoints);
        Dimension numModels = ncfile.addDimension("NumModels", 1);
        Dimension numMaxDepList = ncfile.addDimension("NumMaxDepList", 6);
        Dimension numMaxStoichList = ncfile.addDimension("NumMaxStoichList", 25);
        Dimension stringLen = ncfile.addDimension("StringLen", 72);
        // define variables
        // jms info
        ArrayList<Dimension> dims = new ArrayList<Dimension>();
        dims.add(stringLen);
        if (bMessaging) {
            ncfile.addVariable("JMS_BROKER", DataType.CHAR, dims);
            ncfile.addVariable("JMS_USER", DataType.CHAR, dims);
            ncfile.addVariable("JMS_PASSWORD", DataType.CHAR, dims);
            ncfile.addVariable("JMS_QUEUE", DataType.CHAR, dims);
            ncfile.addVariable("JMS_TOPIC", DataType.CHAR, dims);
            ncfile.addVariable("VCELL_USER", DataType.CHAR, dims);
            ncfile.addVariable("SIMULATION_KEY", DataType.INT, new ArrayList<Dimension>());
            ncfile.addVariable("JOB_INDEX", DataType.INT, new ArrayList<Dimension>());
        }
        // scalars
        ncfile.addVariable("TStart", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("TEnd", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("SaveTime", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("Volume", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("CellGrowthTime", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("CellGrowthTimeSD", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("ExpType", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("LastTrial", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("LastModel", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("MaxNumModels", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("NumModels", DataType.INT, new ArrayList<Dimension>());
        // variables with at least 1 dimension
        ArrayList<Dimension> dimspecies = new ArrayList<Dimension>();
        dimspecies.add(numSpecies);
        ArrayList<Dimension> dimreactions = new ArrayList<Dimension>();
        dimreactions.add(numReactions);
        ncfile.addVariable("SpeciesSplitOnDivision", DataType.INT, dimspecies);
        ncfile.addVariable("SaveSpeciesData", DataType.INT, dimspecies);
        ncfile.addVariable("Reaction_Rate_Laws", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_DListLen", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_StoichListLen", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_OptionalData", DataType.INT, dimreactions);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxStoichList);
        ncfile.addVariable("Reaction_StoichCoeff", DataType.INT, dims);
        ncfile.addVariable("Reaction_StoichSpecies", DataType.INT, dims);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxDepList);
        ncfile.addVariable("Reaction_DepList", DataType.INT, dims);
        dims.clear();
        dims.add(numReactions);
        dims.add(stringLen);
        ncfile.addVariable("Reaction_names", DataType.CHAR, dims);
        dims.clear();
        dims.add(numSpecies);
        dims.add(stringLen);
        ncfile.addVariable("Species_names", DataType.CHAR, dims);
        ncfile.addVariable("SpeciesIC", DataType.INT, dimspecies);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxDepList);
        ncfile.addVariable("Reaction_Rate_Constants", DataType.DOUBLE, dims);
        // create the file
        try {
            ncfile.create();
        } catch (IOException ioe) {
            ioe.printStackTrace(System.err);
            throw new IOException("Error creating hybrid file " + filename + ": " + ioe.getMessage());
        }
        // write data to the NetCDF file
        try {
            // write jms info
            if (bMessaging) {
                ArrayChar.D1 jmsString = new ArrayChar.D1(stringLen.getLength());
                String jmshost = PropertyLoader.getRequiredProperty(PropertyLoader.jmsHostExternal);
                // 
                // Used for new REST HTTP messaging api (USE THIS WHEN Hyrbid Solvers are compiled).
                // 
                // String jmsrestport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsRestPortExternal);
                // String jmsurl = jmshost+":"+jmsrestport;
                // 
                // connect to messaging using legacy AMQP protocol instead of new REST api.  Needed for legacy pre-compiled solvers.
                // 
                String jmsport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsPortExternal);
                String jmsurl = "failover:(tcp://" + jmshost + ":" + jmsport + ")";
                jmsString.setString(jmsurl);
                ncfile.write("JMS_BROKER", jmsString);
                jmsString.setString(PropertyLoader.getRequiredProperty(PropertyLoader.jmsUser));
                ncfile.write("JMS_USER", jmsString);
                String jmsPassword = PropertyLoader.getSecretValue(PropertyLoader.jmsPasswordValue, PropertyLoader.jmsPasswordFile);
                jmsString.setString(jmsPassword);
                ncfile.write("JMS_PASSWORD", jmsString);
                jmsString.setString(VCellQueue.WorkerEventQueue.getName());
                ncfile.write("JMS_QUEUE", jmsString);
                jmsString.setString(VCellTopic.ServiceControlTopic.getName());
                ncfile.write("JMS_TOPIC", jmsString);
                jmsString.setString(simulation.getVersion().getOwner().getName());
                ncfile.write("VCELL_USER", jmsString);
                ArrayInt.D0 scalarJMS = new ArrayInt.D0();
                scalarJMS.set(Integer.parseInt(simulation.getVersion().getVersionKey() + ""));
                ncfile.write("SIMULATION_KEY", scalarJMS);
                scalarJMS.set(simTask.getSimulationJob().getJobIndex());
                ncfile.write("JOB_INDEX", scalarJMS);
            }
            ArrayDouble.D0 scalarDouble = new ArrayDouble.D0();
            // TStart, TEnd, SaveTime
            if ((timeBounds.getEndingTime() > timeBounds.getStartingTime()) && (outputTimeSpec.getOutputTimeStep() > 0)) {
                scalarDouble.set(timeBounds.getStartingTime());
                ncfile.write("TStart", scalarDouble);
                scalarDouble.set(timeBounds.getEndingTime());
                ncfile.write("TEnd", scalarDouble);
                scalarDouble.set(outputTimeSpec.getOutputTimeStep());
                ncfile.write("SaveTime", scalarDouble);
            } else {
                System.err.println("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
                throw new RuntimeException("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
            }
            // Volume
            // we set volume to 1. This model file cannot support multi-compartmental sizes.
            // When writting the rate constants, we must take the volume into account according to the reaction type.
            scalarDouble.set(1);
            ncfile.write("Volume", scalarDouble);
            // CellGrowthTime, CellGrowthTimeSD,
            scalarDouble.set(0);
            ncfile.write("CellGrowthTime", scalarDouble);
            ncfile.write("CellGrowthTimeSD", scalarDouble);
            // ExpType, Last Trial, Last Model, MaxNumModels, NumModels
            ArrayInt.D0 scalarInt = new ArrayInt.D0();
            scalarInt.set(0);
            ncfile.write("LastTrial", scalarInt);
            ncfile.write("LastModel", scalarInt);
            scalarInt.set(1);
            ncfile.write("ExpType", scalarInt);
            ncfile.write("MaxNumModels", scalarInt);
            ncfile.write("NumModels", scalarInt);
            // SpeciesSplitOnDivision
            ArrayInt A1 = new ArrayInt.D1(numSpecies.getLength());
            Index idx = A1.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                A1.setInt(idx.set(i), 0);
            }
            ncfile.write("SpeciesSplitOnDivision", new int[1], A1);
            // SaveSpeciesData
            ArrayInt A2 = new ArrayInt.D1(numSpecies.getLength());
            idx = A2.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                A2.setInt(idx.set(i), 1);
            }
            ncfile.write("SaveSpeciesData", new int[1], A2);
            // Reaction_Rate_Laws
            ArrayInt A3 = new ArrayInt.D1(numReactions.getLength());
            idx = A3.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A3.setInt(idx.set(i), reactionRateLaws[i].getLawType());
            }
            ncfile.write("Reaction_Rate_Laws", new int[1], A3);
            // Reaction_DListLen
            ArrayInt A4 = new ArrayInt.D1(numReactions.getLength());
            idx = A4.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_0)
                    A4.setInt(idx.set(i), 0);
                else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_1) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_1substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_1substrate))
                    A4.setInt(idx.set(i), 1);
                else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_2substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_2substrate))
                    A4.setInt(idx.set(i), 2);
                else if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_3substrate)
                    A4.setInt(idx.set(i), 3);
            }
            ncfile.write("Reaction_DListLen", new int[1], A4);
            // Reaction_StoichListLen
            ArrayInt A5 = new ArrayInt.D1(numReactions.getLength());
            idx = A5.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A5.setInt(idx.set(i), reactions[i].getActions().size());
            }
            ncfile.write("Reaction_StoichListLen", new int[1], A5);
            // Reaction_OptionalData
            ArrayInt A6 = new ArrayInt.D1(numReactions.getLength());
            idx = A6.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A6.setInt(idx.set(i), 0);
            }
            ncfile.write("Reaction_OptionalData", new int[1], A6);
            // Reaction_StoichCoeff
            ArrayInt A7 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
            idx = A7.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                Action[] actions = (Action[]) reactions[i].getActions().toArray(new Action[reactions[i].getActions().size()]);
                for (int j = 0; j < actions.length; j++) {
                    try {
                        actions[j].getOperand().evaluateConstant();
                        int coeff = (int) Math.round(actions[j].getOperand().evaluateConstant());
                        A7.setInt(idx.set(i, j), coeff);
                    } catch (ExpressionException ex) {
                        ex.printStackTrace(System.err);
                        throw new ExpressionException(ex.getMessage());
                    }
                }
            }
            ncfile.write("Reaction_StoichCoeff", new int[2], A7);
            // Reaction_StoichSpecies
            ArrayInt A8 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
            idx = A8.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ArrayList<Action> actions = reactions[i].getActions();
                for (int j = 0; j < actions.size(); j++) {
                    A8.setInt(idx.set(i, j), getVariableIndex(((Action) actions.get(j)).getVar().getName(), vars));
                }
            }
            ncfile.write("Reaction_StoichSpecies", new int[2], A8);
            // Reaction_DepList
            ArrayInt A9 = new ArrayInt.D2(numReactions.getLength(), numMaxDepList.getLength());
            idx = A9.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ReactionRateLaw rl = reactionRateLaws[i];
                Hashtable<String, Integer> tem = varInProbOrderHash[i];
                Enumeration<String> varnames = tem.keys();
                if (rl.getLawType() == ReactionRateLaw.order_0) {
                // don't do anything here.
                } else if ((rl.getLawType() == ReactionRateLaw.order_1) || (rl.getLawType() == ReactionRateLaw.order_2_1substrate) || (rl.getLawType() == ReactionRateLaw.order_3_1substrate) || (rl.getLawType() == ReactionRateLaw.order_2_2substrate) || (rl.getLawType() == ReactionRateLaw.order_3_3substrate)) {
                    int j = 0;
                    while (varnames.hasMoreElements()) {
                        String name = varnames.nextElement();
                        A9.setInt(idx.set(i, j), getVariableIndex(name, vars));
                        j++;
                    }
                } else if (rl.getLawType() == ReactionRateLaw.order_3_2substrate) {
                    int order = 0;
                    String highOrderName = "";
                    String lowOrderName = "";
                    // we must make sure to put the higher order species first.
                    while (varnames.hasMoreElements()) {
                        lowOrderName = varnames.nextElement();
                        if (tem.get(lowOrderName) > order) {
                            String s = highOrderName;
                            highOrderName = lowOrderName;
                            lowOrderName = s;
                            order = tem.get(highOrderName);
                        }
                    }
                    A9.setInt(idx.set(i, 0), getVariableIndex(highOrderName, vars));
                    A9.setInt(idx.set(i, 1), getVariableIndex(lowOrderName, vars));
                }
            }
            ncfile.write("Reaction_DepList", new int[2], A9);
            // Reaction_names
            ArrayChar A10 = new ArrayChar.D2(numReactions.getLength(), stringLen.getLength());
            for (int i = 0; i < numReactions.getLength(); i++) {
                String name = reactions[i].getName();
                int diff = stringLen.getLength() - name.length();
                if (diff >= 0) {
                    for (int j = 0; j < diff; j++) {
                        name = name + " ";
                    }
                    A10.setString(i, name);
                } else
                    throw new RuntimeException("Name of Reaction:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
            }
            ncfile.write("Reaction_names", A10);
            // Species_names
            ArrayChar A11 = new ArrayChar.D2(numSpecies.getLength(), stringLen.getLength());
            for (int i = 0; i < numSpecies.getLength(); i++) {
                String name = vars.elementAt(i).getName();
                int diff = stringLen.getLength() - name.length();
                if (diff >= 0) {
                    for (int j = 0; j < diff; j++) {
                        name = name + " ";
                    }
                    A11.setString(i, name);
                } else
                    throw new RuntimeException("Name of Species:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
            }
            ncfile.write("Species_names", A11);
            // Species Initial Condition (in number of molecules).
            // Species iniCondition are sampled from a poisson distribution(which has a mean of the current iniExp value)
            RandomDataGenerator dist = new RandomDataGenerator();
            if (stochOpt.isUseCustomSeed()) {
                Integer randomSeed = stochOpt.getCustomSeed();
                if (randomSeed != null) {
                    dist.reSeed(randomSeed);
                }
            }
            ArrayLong A12 = new ArrayLong.D1(numSpecies.getLength());
            idx = A12.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                try {
                    VarIniCondition varIniCondition = subDomain.getVarIniCondition(vars.elementAt(i));
                    Expression varIniExp = varIniCondition.getIniVal();
                    varIniExp.bindExpression(simSymbolTable);
                    varIniExp = simSymbolTable.substituteFunctions(varIniExp).flatten();
                    double expectedCount = varIniExp.evaluateConstant();
                    long varCount = 0;
                    if (varIniCondition instanceof VarIniCount) {
                        varCount = (long) expectedCount;
                    } else {
                        if (expectedCount > 0) {
                            varCount = dist.nextPoisson(expectedCount);
                        }
                    }
                    A12.setLong(idx.set(i), varCount);
                } catch (ExpressionException ex) {
                    ex.printStackTrace(System.err);
                    throw new ExpressionException(ex.getMessage());
                }
            }
            ncfile.write("SpeciesIC", new int[1], A12);
            // Reaction_Rate_Constants(NumReactions, NumMaxDepList) ;
            ArrayDouble A13 = new ArrayDouble.D2(numReactions.getLength(), numMaxDepList.getLength());
            idx = A13.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ReactionRateLaw rl = reactionRateLaws[i];
                A13.setDouble(idx.set(i, 0), rl.getRateConstant());
            }
            ncfile.write("Reaction_Rate_Constants", A13);
        } catch (IOException ioe) {
            ioe.printStackTrace(System.err);
            throw new IOException("Error writing hybrid input file " + filename + ": " + ioe.getMessage());
        } catch (InvalidRangeException ire) {
            ire.printStackTrace(System.err);
            throw new InvalidRangeException("Error writing hybrid input file " + filename + ": " + ire.getMessage());
        }
        try {
            ncfile.close();
        } catch (IOException ioe) {
            throw new IOException("Error closing file " + filename + ". " + ioe.getMessage());
        }
    }
}

Example 25 with SimulationSymbolTable

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

the class HybridSolver method getHybridSolverResultSet.

/**
 * Get data columns and function columns to be ready for plot.
 * Creation date: (8/15/2006 11:36:43 AM)
 * @return cbit.vcell.solver.stoch.StochSolverResultSet
 */
private ODESolverResultSet getHybridSolverResultSet() {
    // read .stoch file, this funciton here equals to getODESolverRestultSet()+getStateVariableResultSet()  in ODE.
    ODESolverResultSet stSolverResultSet = new ODESolverResultSet();
    try {
        String filename = getBaseName() + NETCDF_DATA_EXTENSION;
        NetCDFEvaluator ncEva = new NetCDFEvaluator();
        NetCDFReader ncReader = null;
        try {
            ncEva.setNetCDFTarget(filename);
            ncReader = ncEva.getNetCDFReader();
        } catch (Exception e) {
            e.printStackTrace(System.err);
            throw new RuntimeException("Cannot open simulation result file: " + filename + "!");
        }
        // Read result according to trial number
        if (ncReader.getNumTrials() == 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription("t"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no time points
            ArrayDouble data = (ArrayDouble) ncEva.getTimeSeriesData(1);
            double[] timePoints = ncReader.getTimePoints();
            System.out.println("time points length is " + timePoints.length);
            // shape[0]:num of timepoints, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                System.out.println("one species in time series data and size is " + temData.getSize());
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                System.out.println("multiple species in time series, the length of time series is :" + data.getShape()[0] + ", and the total number of speceis is: " + data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
        } else if (ncReader.getNumTrials() > 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription("TrialNo"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no trial numbers
            ArrayDouble data = (ArrayDouble) ncEva.getDataOverTrials(ncReader.getTimePoints().length - 1);
            int[] trialNum = ncEva.getNetCDFReader().getTrialNumbers();
            // System.out.println("total trials are "+trialNum.length);
            // shape[0]:number of trials, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                // System.out.println("one species over trials, size is: "+temData.getSize());
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                // System.out.println("multiple species in multiple trials, the length of trials is :"+data.getShape()[0]+", and the total number of speceis is: "+data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
        } else {
            throw new RuntimeException("Number of trials should be a countable positive value, from 1 to N.");
        }
    } catch (Exception e) {
        e.printStackTrace(System.err);
        throw new RuntimeException("Problem encountered in parsing hybrid simulation results.\n" + e.getMessage());
    }
    /*
	 *Add appropriate Function columns to result set if the stochastic simulation is to display the trajectory.
	 *No function columns for the results of multiple stochastic trials.
	 *In stochastic simulation the functions include probability functions and clamped variable.
	 */
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (simSymbolTable.getSimulation().getSolverTaskDescription().getStochOpt().getNumOfTrials() == 1) {
        Function[] functions = simSymbolTable.getFunctions();
        for (int i = 0; i < functions.length; i++) {
            if (SimulationSymbolTable.isFunctionSaved(functions[i])) {
                Expression exp1 = new Expression(functions[i].getExpression());
                try {
                    exp1 = simSymbolTable.substituteFunctions(exp1);
                } catch (MathException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
                }
                try {
                    FunctionColumnDescription cd = new FunctionColumnDescription(exp1.flatten(), functions[i].getName(), null, functions[i].getName(), false);
                    stSolverResultSet.addFunctionColumn(cd);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                }
            }
        }
    }
    return stSolverResultSet;
}