Search in sources :

Example 16 with UniformOutputTimeSpec

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

the class StandaloneRuleBasedTest method checkNonspatialStochasticSimContext.

private static void checkNonspatialStochasticSimContext(SimulationContext srcSimContext, File baseDirectory, int numTrials, long bngTimeoutDuration) throws Exception {
    if (!srcSimContext.getApplicationType().equals(Application.NETWORK_STOCHASTIC) || srcSimContext.getGeometry().getDimension() != 0) {
        throw new RuntimeException("simContext is of type " + srcSimContext.getApplicationType() + " and geometry dimension of " + srcSimContext.getGeometry().getDimension() + ", expecting nonspatial stochastic");
    }
    BioModel origBioModel = srcSimContext.getBioModel();
    BioModel bioModel = XmlHelper.XMLToBioModel(new XMLSource(XmlHelper.bioModelToXML(origBioModel)));
    bioModel.refreshDependencies();
    // create ODE and RuleBased
    SimulationContext newODEApp = SimulationContext.copySimulationContext(srcSimContext, "aUniqueNewODEApp", false, Application.NETWORK_DETERMINISTIC);
    SimulationContext newRuleBasedApp = SimulationContext.copySimulationContext(srcSimContext, "aUniqueNewRuleBasedApp", false, Application.RULE_BASED_STOCHASTIC);
    newODEApp.setBioModel(bioModel);
    newRuleBasedApp.setBioModel(bioModel);
    ArrayList<AnnotatedFunction> outputFunctionsList = srcSimContext.getOutputFunctionContext().getOutputFunctionsList();
    // OutputContext outputContext = new OutputContext(outputFunctionsList.toArray(new AnnotatedFunction[outputFunctionsList.size()]));
    newODEApp.getOutputFunctionContext().setOutputFunctions(outputFunctionsList);
    newRuleBasedApp.getOutputFunctionContext().setOutputFunctions(outputFunctionsList);
    NetworkGenerationRequirements networkGenRequirements = NetworkGenerationRequirements.getComputeFull(bngTimeoutDuration);
    bioModel.addSimulationContext(newODEApp);
    newODEApp.refreshMathDescription(new MathMappingCallbackTaskAdapter(null), networkGenRequirements);
    bioModel.addSimulationContext(newRuleBasedApp);
    newRuleBasedApp.refreshMathDescription(new MathMappingCallbackTaskAdapter(null), networkGenRequirements);
    srcSimContext.refreshMathDescription(new MathMappingCallbackTaskAdapter(null), networkGenRequirements);
    // Create non-spatialStoch, ODE and RuleBased sims
    Simulation nonspatialStochAppNewSim = srcSimContext.addNewSimulation(STOCH_SIM_NAME, /*SimulationOwner.DEFAULT_SIM_NAME_PREFIX*/
    new MathMappingCallbackTaskAdapter(null), networkGenRequirements);
    Simulation newODEAppNewSim = newODEApp.addNewSimulation(ODE_SIM_NAME, new MathMappingCallbackTaskAdapter(null), networkGenRequirements);
    Simulation newRuleBasedAppNewSim = newRuleBasedApp.addNewSimulation(NFS_SIM_NAME, new MathMappingCallbackTaskAdapter(null), networkGenRequirements);
    nonspatialStochAppNewSim.setSimulationOwner(srcSimContext);
    newODEAppNewSim.setSimulationOwner(newODEApp);
    newRuleBasedAppNewSim.setSimulationOwner(newRuleBasedApp);
    try {
        bioModel.getModel().getSpeciesContexts();
        ArrayList<String> varNameList = new ArrayList<String>();
        for (SpeciesContextSpec scs : srcSimContext.getReactionContext().getSpeciesContextSpecs()) {
            varNameList.add(scs.getSpeciesContext().getName());
        }
        String[] varNames = varNameList.toArray(new String[0]);
        OutputTimeSpec outputTimeSpec = nonspatialStochAppNewSim.getSolverTaskDescription().getOutputTimeSpec();
        ArrayList<Double> sampleTimeList = new ArrayList<Double>();
        if (outputTimeSpec instanceof UniformOutputTimeSpec) {
            double endingTime = nonspatialStochAppNewSim.getSolverTaskDescription().getTimeBounds().getEndingTime();
            double dT = ((UniformOutputTimeSpec) outputTimeSpec).getOutputTimeStep();
            int currTimeIndex = 0;
            while (currTimeIndex * dT <= (endingTime + 1e-8)) {
                sampleTimeList.add(currTimeIndex * dT);
                currTimeIndex++;
            }
        }
        double[] sampleTimes = new double[sampleTimeList.size()];
        for (int i = 0; i < sampleTimes.length; i++) {
            sampleTimes[i] = sampleTimeList.get(i);
        }
        TimeSeriesMultitrialData sampleDataStoch1 = new TimeSeriesMultitrialData("stochastic1", varNames, sampleTimes, numTrials);
        TimeSeriesMultitrialData sampleDataStoch2 = new TimeSeriesMultitrialData("stochastic2", varNames, sampleTimes, numTrials);
        TimeSeriesMultitrialData sampleDataDeterministic = new TimeSeriesMultitrialData("determinstic", varNames, sampleTimes, 1);
        runsolver(nonspatialStochAppNewSim, baseDirectory, numTrials, sampleDataStoch1);
        runsolver(newODEAppNewSim, baseDirectory, 1, sampleDataDeterministic);
        runsolver(newRuleBasedAppNewSim, baseDirectory, numTrials, sampleDataStoch2);
        StochtestFileUtils.writeVarDiffData(new File(baseDirectory, VARDIFF_FILE), sampleDataStoch1, sampleDataStoch2);
        StochtestFileUtils.writeKolmogorovSmirnovTest(new File(baseDirectory, KS_TEST_FILE), sampleDataStoch1, sampleDataStoch2);
        StochtestFileUtils.writeChiSquareTest(new File(baseDirectory, ChiSquared_TEST_FILE), sampleDataStoch1, sampleDataStoch2);
        StochtestFileUtils.writeData(sampleDataStoch1, new File(baseDirectory, "data." + sampleDataStoch1.datasetName + ".json"));
        StochtestFileUtils.writeData(sampleDataStoch2, new File(baseDirectory, "data." + sampleDataStoch2.datasetName + ".json"));
        StochtestFileUtils.writeData(sampleDataDeterministic, new File(baseDirectory, "data." + sampleDataDeterministic.datasetName + ".json"));
    } finally {
        srcSimContext.removeSimulation(nonspatialStochAppNewSim);
        newODEApp.removeSimulation(newODEAppNewSim);
        newRuleBasedApp.removeSimulation(newRuleBasedAppNewSim);
    }
}
Also used : MathMappingCallbackTaskAdapter(cbit.vcell.mapping.MathMappingCallbackTaskAdapter) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) TimeSeriesMultitrialData(org.vcell.stochtest.TimeSeriesMultitrialData) ArrayList(java.util.ArrayList) SimulationContext(cbit.vcell.mapping.SimulationContext) SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) OutputTimeSpec(cbit.vcell.solver.OutputTimeSpec) TempSimulation(cbit.vcell.solver.TempSimulation) Simulation(cbit.vcell.solver.Simulation) BioModel(cbit.vcell.biomodel.BioModel) NetworkGenerationRequirements(cbit.vcell.mapping.SimulationContext.NetworkGenerationRequirements) XMLSource(cbit.vcell.xml.XMLSource) File(java.io.File) AnnotatedFunction(cbit.vcell.solver.AnnotatedFunction)

Example 17 with UniformOutputTimeSpec

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

the class OutputOptionsPanel method actionOutputOptionButtonState.

private void actionOutputOptionButtonState(java.awt.event.ActionEvent actionEvent) {
    try {
        if (solverTaskDescription == null) {
            return;
        }
        OutputTimeSpec outputTimeSpec = solverTaskDescription.getOutputTimeSpec();
        if (actionEvent.getSource() == getDefaultOutputRadioButton() && !outputTimeSpec.isDefault()) {
            solverTaskDescription.setOutputTimeSpec(new DefaultOutputTimeSpec());
        } else if (actionEvent.getSource() == getUniformOutputRadioButton() && !outputTimeSpec.isUniform()) {
            double outputTime = 0.0;
            if (solverTaskDescription.getSolverDescription().isSemiImplicitPdeSolver()) {
                String floatStr = "" + (float) (((DefaultOutputTimeSpec) outputTimeSpec).getKeepEvery() * solverTaskDescription.getTimeStep().getDefaultTimeStep());
                outputTime = Double.parseDouble(floatStr);
            } else {
                TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
                Range outputTimeRange = NumberUtils.getDecimalRange(timeBounds.getStartingTime(), timeBounds.getEndingTime() / 100, true, true);
                outputTime = outputTimeRange.getMax();
            }
            solverTaskDescription.setOutputTimeSpec(new UniformOutputTimeSpec(outputTime));
        } else if (actionEvent.getSource() == getExplicitOutputRadioButton() && !outputTimeSpec.isExplicit()) {
            TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
            solverTaskDescription.setOutputTimeSpec(new ExplicitOutputTimeSpec(new double[] { timeBounds.getStartingTime(), timeBounds.getEndingTime() }));
        }
    } catch (java.lang.Throwable ivjExc) {
        handleException(ivjExc);
    }
}
Also used : TimeBounds(cbit.vcell.solver.TimeBounds) DefaultOutputTimeSpec(cbit.vcell.solver.DefaultOutputTimeSpec) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) ExplicitOutputTimeSpec(cbit.vcell.solver.ExplicitOutputTimeSpec) OutputTimeSpec(cbit.vcell.solver.OutputTimeSpec) ExplicitOutputTimeSpec(cbit.vcell.solver.ExplicitOutputTimeSpec) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) Range(org.vcell.util.Range) DefaultOutputTimeSpec(cbit.vcell.solver.DefaultOutputTimeSpec)

Example 18 with UniformOutputTimeSpec

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

the class OutputOptionsPanel method setNewOutputOption.

/**
 * Comment
 */
private void setNewOutputOption() {
    try {
        OutputTimeSpec ots = null;
        if (getDefaultOutputRadioButton().isSelected()) {
            int keepEvery = Integer.parseInt(getKeepEveryTextField().getText());
            if (solverTaskDescription.getSolverDescription().isSemiImplicitPdeSolver()) {
                ots = new DefaultOutputTimeSpec(keepEvery);
            } else {
                int keepAtMost = Integer.parseInt(getKeepAtMostTextField().getText());
                ots = new DefaultOutputTimeSpec(keepEvery, keepAtMost);
            }
        } else if (getUniformOutputRadioButton().isSelected()) {
            double outputTime = Double.parseDouble(getOutputTimeStepTextField().getText());
            ots = new UniformOutputTimeSpec(outputTime);
        } else if (getExplicitOutputRadioButton().isSelected()) {
            ots = ExplicitOutputTimeSpec.fromString(getOutputTimesTextField().getText());
        }
        try {
            solverTaskDescription.setOutputTimeSpec(ots);
        } catch (java.beans.PropertyVetoException e) {
            e.printStackTrace(System.out);
            throw new RuntimeException(e.getMessage());
        }
    } catch (Exception e) {
        DialogUtils.showErrorDialog(this, e.getMessage(), e);
    }
}
Also used : PropertyVetoException(java.beans.PropertyVetoException) DefaultOutputTimeSpec(cbit.vcell.solver.DefaultOutputTimeSpec) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) ExplicitOutputTimeSpec(cbit.vcell.solver.ExplicitOutputTimeSpec) OutputTimeSpec(cbit.vcell.solver.OutputTimeSpec) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) DefaultOutputTimeSpec(cbit.vcell.solver.DefaultOutputTimeSpec) PropertyVetoException(java.beans.PropertyVetoException)

Example 19 with UniformOutputTimeSpec

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

the class XmlReader method getChomboSolverSpec.

private ChomboSolverSpec getChomboSolverSpec(SolverTaskDescription solverTaskDesc, Element element, int dimension) throws XmlParseException {
    int maxBoxSize = parseIntWithDefault(element, XMLTags.MaxBoxSizeTag, ChomboSolverSpec.getDefaultMaxBoxSize(dimension));
    double fillRatio = parseDoubleWithDefault(element, XMLTags.FillRatioTag, ChomboSolverSpec.getDefaultFillRatio());
    boolean bSaveVCellOutput = parseBooleanWithDefault(element, XMLTags.SaveVCellOutput, true);
    boolean bSaveChomboOutput = parseBooleanWithDefault(element, XMLTags.SaveChomboOutput, false);
    Element childElement = element.getChild(XMLTags.RefineRatios, vcNamespace);
    List<Integer> refineRatioList = null;
    if (childElement != null) {
        String text = childElement.getText();
        if (text != null && !text.isEmpty()) {
            StringTokenizer st = new StringTokenizer(text, ",");
            if (st.hasMoreTokens()) {
                refineRatioList = new ArrayList<Integer>();
                while (st.hasMoreElements()) {
                    String token = st.nextToken();
                    if (token != null) {
                        int n = Integer.parseInt(token);
                        refineRatioList.add(n);
                    }
                }
            }
        }
    }
    Integer viewLevel = null;
    try {
        ChomboSolverSpec css = new ChomboSolverSpec(maxBoxSize, fillRatio, viewLevel, bSaveVCellOutput, bSaveChomboOutput, refineRatioList);
        double smallVolfracThreshold = parseDoubleWithDefault(element, XMLTags.SmallVolfracThreshold, 0);
        int blockFactor = parseIntWithDefault(element, XMLTags.BlockFactorTag, ChomboSolverSpec.DEFAULT_BLOCK_FACTOR);
        boolean bActivateFeatureUnderDevelopment = parseBooleanWithDefault(element, XMLTags.ActivateFeatureUnderDevelopment, false);
        css.setSmallVolfracThreshold(smallVolfracThreshold);
        css.setActivateFeatureUnderDevelopment(bActivateFeatureUnderDevelopment);
        css.setBlockFactor(blockFactor);
        int tagsGrow = parseIntWithDefault(element, XMLTags.TagsGrowTag, ChomboSolverSpec.defaultTagsGrow);
        css.setTagsGrow(tagsGrow);
        Element timeBoundsElement = element.getChild(XMLTags.TimeBoundTag, vcNamespace);
        List<Element> timeIntervalElementList = null;
        boolean noTimeBounds = false;
        if (timeBoundsElement == null) {
            noTimeBounds = true;
        } else {
            timeIntervalElementList = timeBoundsElement.getChildren(XMLTags.TimeIntervalTag, vcNamespace);
            if (timeIntervalElementList.size() == 0) {
                noTimeBounds = true;
            }
        }
        if (noTimeBounds) {
            // old format
            double startTime = 0;
            double endTime = solverTaskDesc.getTimeBounds().getEndingTime();
            double timeStep = solverTaskDesc.getTimeStep().getDefaultTimeStep();
            double outputTimeStep = ((UniformOutputTimeSpec) solverTaskDesc.getOutputTimeSpec()).getOutputTimeStep();
            try {
                TimeInterval ti = new TimeInterval(startTime, endTime, timeStep, outputTimeStep);
                css.addTimeInterval(ti);
            } catch (IllegalArgumentException ex) {
                css.addTimeInterval(TimeInterval.getDefaultTimeInterval());
            }
        } else {
            for (Element e : timeIntervalElementList) {
                String s = e.getAttributeValue(XMLTags.StartTimeAttrTag);
                double startTime = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.EndTimeAttrTag);
                double endTime = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.TimeStepAttrTag);
                double timeStep = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.OutputTimeStepAttrTag);
                double outputTimeStep = Double.valueOf(s);
                TimeInterval ti = new TimeInterval(startTime, endTime, timeStep, outputTimeStep);
                css.addTimeInterval(ti);
            }
        }
        Element meshRefineElement = element.getChild(XMLTags.MeshRefinementTag, vcNamespace);
        if (meshRefineElement != null) {
            if (meshRefineElement.getChildren().size() != 0) {
                // in old model, if there is no refinement, set view level to finest
                // only set viewLevel when meshRefinement has children
                Element viewLevelChild = element.getChild(XMLTags.ViewLevelTag, vcNamespace);
                if (viewLevelChild != null) {
                    viewLevel = parseIntWithDefault(element, XMLTags.ViewLevelTag, 0);
                    css.setViewLevel(viewLevel);
                }
            }
            List<Element> levelElementList = meshRefineElement.getChildren(XMLTags.RefinementRoiTag, vcNamespace);
            for (Element levelElement : levelElementList) {
                String levelStr = levelElement.getAttributeValue(XMLTags.RefineRoiLevelAttrTag);
                int level = 1;
                if (levelStr != null) {
                    level = Integer.parseInt(levelStr);
                }
                String type = levelElement.getAttributeValue(XMLTags.RefinementRoiTypeAttrTag);
                RoiType roiType = RoiType.Membrane;
                if (type != null) {
                    try {
                        roiType = RoiType.valueOf(type);
                    } catch (Exception ex) {
                    // ignore
                    }
                }
                Element expElement = levelElement.getChild(XMLTags.ROIExpressionTag, vcNamespace);
                String roiExp = null;
                if (expElement != null) {
                    roiExp = expElement.getText();
                    RefinementRoi roi = new RefinementRoi(roiType, level, roiExp);
                    css.addRefinementRoi(roi);
                }
            }
        }
        return css;
    } catch (ExpressionException e) {
        throw new XmlParseException(e);
    }
}
Also used : UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) TimeInterval(org.vcell.chombo.TimeInterval) Element(org.jdom.Element) ChomboSolverSpec(org.vcell.chombo.ChomboSolverSpec) GeometryException(cbit.vcell.geometry.GeometryException) MathFormatException(cbit.vcell.math.MathFormatException) MappingException(cbit.vcell.mapping.MappingException) PropertyVetoException(java.beans.PropertyVetoException) ImageException(cbit.image.ImageException) ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException) ModelException(cbit.vcell.model.ModelException) DataConversionException(org.jdom.DataConversionException) ExpressionException(cbit.vcell.parser.ExpressionException) MathException(cbit.vcell.math.MathException) ExpressionException(cbit.vcell.parser.ExpressionException) RefinementRoi(org.vcell.chombo.RefinementRoi) StringTokenizer(java.util.StringTokenizer) RoiType(org.vcell.chombo.RefinementRoi.RoiType)

Example 20 with UniformOutputTimeSpec

use of cbit.vcell.solver.UniformOutputTimeSpec 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());
        }
    }
}
Also used : VarIniCondition(cbit.vcell.math.VarIniCondition) NonspatialStochSimOptions(cbit.vcell.solver.NonspatialStochSimOptions) ArrayList(java.util.ArrayList) Index(ucar.ma2.Index) ExpressionException(cbit.vcell.parser.ExpressionException) SubDomain(cbit.vcell.math.SubDomain) SolverTaskDescription(cbit.vcell.solver.SolverTaskDescription) Vector(java.util.Vector) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) InvalidRangeException(ucar.ma2.InvalidRangeException) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) ArrayLong(ucar.ma2.ArrayLong) ArrayInt(ucar.ma2.ArrayInt) NetcdfFileWriteable(ucar.nc2.NetcdfFileWriteable) Action(cbit.vcell.math.Action) StochVolVariable(cbit.vcell.math.StochVolVariable) Variable(cbit.vcell.math.Variable) TimeBounds(cbit.vcell.solver.TimeBounds) ArrayChar(ucar.ma2.ArrayChar) ArrayDouble(ucar.ma2.ArrayDouble) JumpProcess(cbit.vcell.math.JumpProcess) StochVolVariable(cbit.vcell.math.StochVolVariable) VarIniCount(cbit.vcell.math.VarIniCount) Dimension(ucar.nc2.Dimension) IOException(java.io.IOException) Simulation(cbit.vcell.solver.Simulation) Expression(cbit.vcell.parser.Expression)

Aggregations

UniformOutputTimeSpec (cbit.vcell.solver.UniformOutputTimeSpec)34 Simulation (cbit.vcell.solver.Simulation)18 OutputTimeSpec (cbit.vcell.solver.OutputTimeSpec)17 TimeBounds (cbit.vcell.solver.TimeBounds)15 DefaultOutputTimeSpec (cbit.vcell.solver.DefaultOutputTimeSpec)14 BioModel (cbit.vcell.biomodel.BioModel)13 SimulationContext (cbit.vcell.mapping.SimulationContext)13 SpeciesContextSpec (cbit.vcell.mapping.SpeciesContextSpec)10 MathDescription (cbit.vcell.math.MathDescription)10 KeyValue (org.vcell.util.document.KeyValue)10 Expression (cbit.vcell.parser.Expression)9 TimeStep (cbit.vcell.solver.TimeStep)9 IOException (java.io.IOException)9 ImageException (cbit.image.ImageException)8 Model (cbit.vcell.model.Model)8 SimulationVersion (org.vcell.util.document.SimulationVersion)8 Geometry (cbit.vcell.geometry.Geometry)7 MathMapping (cbit.vcell.mapping.MathMapping)7 ErrorTolerance (cbit.vcell.solver.ErrorTolerance)7 ExplicitOutputTimeSpec (cbit.vcell.solver.ExplicitOutputTimeSpec)7