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);
}
}
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);
}
}
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);
}
}
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);
}
}
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());
}
}
}
Aggregations