use of cbit.vcell.solver.TimeBounds 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.TimeBounds in project vcell by virtualcell.
the class TimeBoundsPanel method setTimeBounds.
/**
* Sets the timeBounds property (cbit.vcell.solver.TimeBounds) value.
* @param timeBounds The new value for the property.
* @see #getTimeBounds
*/
public void setTimeBounds(TimeBounds timeBounds) {
TimeBounds oldValue = fieldTimeBounds;
fieldTimeBounds = timeBounds;
if (fieldTimeBounds != null) {
getStartingTimeTextField().setText(String.valueOf(fieldTimeBounds.getStartingTime()));
getEndingTimeTextField().setText(String.valueOf(fieldTimeBounds.getEndingTime()));
}
firePropertyChange("timeBounds", oldValue, timeBounds);
}
use of cbit.vcell.solver.TimeBounds 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());
}
}
}
use of cbit.vcell.solver.TimeBounds in project vcell by virtualcell.
the class XmlReader method getTimeBounds.
/**
* This method returns a TimeBounds object from a XML Element.
* Creation date: (5/22/2001 11:41:04 AM)
* @return cbit.vcell.solver.TimeBounds
* @param param org.jdom.Element
*/
private TimeBounds getTimeBounds(Element param) {
// get Attributes
double start = Double.parseDouble(param.getAttributeValue(XMLTags.StartTimeAttrTag));
double end = Double.parseDouble(param.getAttributeValue(XMLTags.EndTimeAttrTag));
// *** create new TimeBounds object ****
TimeBounds timeBounds = new TimeBounds(start, end);
return timeBounds;
}
use of cbit.vcell.solver.TimeBounds in project vcell by virtualcell.
the class VCellSBMLSolver method solve.
public File solve(String filePrefix, File outDir, String sbmlFileName, SimSpec testSpec) throws IOException, SolverException, SbmlException {
try {
cbit.util.xml.VCLogger sbmlImportLogger = new LocalLogger();
//
// Instantiate an SBMLImporter to get the speciesUnitsHash - to compute the conversion factor from VC->SB species units.
// and import SBML (sbml->bioModel)
BioModel bioModel = importSBML(sbmlFileName, sbmlImportLogger, false);
// Hashtable<String, SBMLImporter.SBVCConcentrationUnits> speciesUnitsHash = sbmlImporter.getSpeciesUnitsHash();
// double timeFactor = sbmlImporter.getSBMLTimeUnitsFactor();
String vcml_1 = XmlHelper.bioModelToXML(bioModel);
SBMLUtils.writeStringToFile(vcml_1, new File(outDir, filePrefix + ".vcml").getAbsolutePath(), true);
if (bRoundTrip) {
// Round trip the bioModel (bioModel->sbml->bioModel).
// save imported "bioModel" as VCML
// String vcml_1 = XmlHelper.bioModelToXML(bioModel);
// SBMLUtils.writeStringToFile(vcml_1, new File(outDir,filePrefix+".vcml").getAbsolutePath());
// export bioModel as sbml and save
// String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, bioModel.getSimulationContexts(0).getName());
// SimulationJob simJob = new SimulationJob(bioModel.getSimulations(bioModel.getSimulationContexts(0))[0], null, 0);
String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, 0, false, bioModel.getSimulationContext(0), null);
SBMLUtils.writeStringToFile(vcml_sbml, new File(outDir, filePrefix + ".vcml.sbml").getAbsolutePath(), true);
// re-import bioModel from exported sbml
XMLSource vcml_sbml_Src = new XMLSource(vcml_sbml);
BioModel newBioModel = (BioModel) XmlHelper.importSBML(sbmlImportLogger, vcml_sbml_Src, false);
String vcml_sbml_vcml = XmlHelper.bioModelToXML(newBioModel);
SBMLUtils.writeStringToFile(vcml_sbml_vcml, new File(outDir, filePrefix + ".vcml.sbml.vcml").getAbsolutePath(), true);
// have rest of code use the round-tripped biomodel
bioModel = newBioModel;
}
//
// select only Application, generate math, and create a single Simulation.
//
SimulationContext simContext = bioModel.getSimulationContext(0);
MathMapping mathMapping = simContext.createNewMathMapping();
MathDescription mathDesc = mathMapping.getMathDescription();
String vcml = mathDesc.getVCML();
try (PrintWriter pw = new PrintWriter("vcmlTrace.txt")) {
pw.println(vcml);
}
simContext.setMathDescription(mathDesc);
SimulationVersion simVersion = new SimulationVersion(new KeyValue("100"), "unnamed", null, null, null, null, null, null, null, null);
Simulation sim = new Simulation(simVersion, mathDesc);
sim.setName("unnamed");
// if time factor from SBML is not 1 (i.e., it is not in secs but in minutes or hours), convert endTime to min/hr as : endTime*timeFactor
// double endTime = testSpec.getEndTime()*timeFactor;
double endTime = testSpec.getEndTime();
sim.getSolverTaskDescription().setTimeBounds(new TimeBounds(0, endTime));
TimeStep timeStep = new TimeStep();
sim.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep.getMinimumTimeStep(), timeStep.getDefaultTimeStep(), endTime / 10000));
sim.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((endTime - 0) / testSpec.getNumTimeSteps()));
sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
// sim.getSolverTaskDescription().setErrorTolerance(new cbit.vcell.solver.ErrorTolerance(1e-10, 1e-12));
// Generate .idaInput string
/* IDAFileWriter idaFileWriter = new IDAFileWriter(sim);
File idaInputFile = new File(filePathName.replace(".vcml", ".idaInput"));
PrintWriter idaPW = new java.io.PrintWriter(idaInputFile);
idaFileWriter.writeInputFile(idaPW);
idaPW.close();
// use the idastandalone solver
File idaOutputFile = new File(filePathName.replace(".vcml", ".ida"));
Executable executable = new Executable("IDAStandalone " + idaInputFile + " " + idaOutputFile);
executable.start();
*/
// Generate .cvodeInput string
File cvodeFile = new File(outDir, filePrefix + SimDataConstants.CVODEINPUT_DATA_EXTENSION);
PrintWriter cvodePW = new java.io.PrintWriter(cvodeFile);
SimulationJob simJob = new SimulationJob(sim, 0, null);
SimulationTask simTask = new SimulationTask(simJob, 0);
CVodeFileWriter cvodeFileWriter = new CVodeFileWriter(cvodePW, simTask);
cvodeFileWriter.write();
cvodePW.close();
// use the cvodeStandalone solver
File cvodeOutputFile = new File(outDir, filePrefix + SimDataConstants.IDA_DATA_EXTENSION);
String executableName = null;
try {
executableName = SolverUtilities.getExes(SolverDescription.CVODE)[0].getAbsolutePath();
} catch (IOException e) {
throw new RuntimeException("failed to get executable for solver " + SolverDescription.CVODE.getDisplayLabel() + ": " + e.getMessage(), e);
}
Executable executable = new Executable(new String[] { executableName, cvodeFile.getAbsolutePath(), cvodeOutputFile.getAbsolutePath() });
executable.start();
// get the result
ODESolverResultSet odeSolverResultSet = getODESolverResultSet(simJob, cvodeOutputFile.getPath());
//
// print header
//
File outputFile = new File(outDir, filePrefix + ".vcell.csv");
java.io.PrintStream outputStream = new java.io.PrintStream(new java.io.BufferedOutputStream(new java.io.FileOutputStream(outputFile)));
outputStream.print("time");
for (int i = 0; i < testSpec.getVarsList().length; i++) {
outputStream.print("," + testSpec.getVarsList()[i]);
}
outputStream.println();
//
// extract data for time and species
//
double[][] data = new double[testSpec.getVarsList().length + 1][];
int column = odeSolverResultSet.findColumn("t");
data[0] = odeSolverResultSet.extractColumn(column);
int origDataLength = data[0].length;
for (int i = 0; i < testSpec.getVarsList().length; i++) {
column = odeSolverResultSet.findColumn(testSpec.getVarsList()[i]);
if (column == -1) {
Variable var = simJob.getSimulationSymbolTable().getVariable(testSpec.getVarsList()[i]);
data[i + 1] = new double[data[0].length];
if (var instanceof cbit.vcell.math.Constant) {
double value = ((cbit.vcell.math.Constant) var).getExpression().evaluateConstant();
for (int j = 0; j < data[i + 1].length; j++) {
data[i + 1][j] = value;
}
} else {
throw new RuntimeException("Did not find " + testSpec.getVarsList()[i] + " in simulation");
}
} else {
data[i + 1] = odeSolverResultSet.extractColumn(column);
}
}
//
// for each time, print row
//
int index = 0;
double[] sampleTimes = new double[testSpec.getNumTimeSteps() + 1];
for (int i = 0; i <= testSpec.getNumTimeSteps(); i++) {
sampleTimes[i] = endTime * i / testSpec.getNumTimeSteps();
}
Model vcModel = bioModel.getModel();
ReservedSymbol kMole = vcModel.getKMOLE();
for (int i = 0; i < sampleTimes.length; i++) {
//
while (true) {
//
if (index == odeSolverResultSet.getRowCount() - 1) {
if (data[0][index] == sampleTimes[i]) {
break;
} else {
throw new RuntimeException("sampleTime does not match at last time point");
}
}
//
if (data[0][index + 1] > sampleTimes[i]) {
break;
}
//
// sampleTime must be later in our data list.
//
index++;
}
// if data[0][index] == sampleTime no need to interpolate
if (data[0][index] == sampleTimes[i]) {
// if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
// if (timeFactor != 1.0) {
// outputStream.print(data[0][index]/timeFactor);
// } else {
outputStream.print(data[0][index]);
// }
for (int j = 0; j < testSpec.getVarsList().length; j++) {
// SBMLImporter.SBVCConcentrationUnits spConcUnits = speciesUnitsHash.get(testSpec.getVarsList()[j]);
// if (spConcUnits != null) {
// VCUnitDefinition sbunits = spConcUnits.getSBConcentrationUnits();
// VCUnitDefinition vcunits = spConcUnits.getVCConcentrationUnits();
// SBMLUnitParameter unitFactor = SBMLUtils.getConcUnitFactor("spConcParam", vcunits, sbunits, kMole);
// outputStream.print("," + data[j + 1][index] * unitFactor.getExpression().evaluateConstant()); //earlier, hack unitfactor = 0.000001
// earlier, hack unitfactor = 0.000001
outputStream.print("," + data[j + 1][index]);
// }
}
// System.out.println("No interpolation needed!");
outputStream.println();
} else {
// if data[0][index] < sampleTime, must interpolate
double fraction = (sampleTimes[i] - data[0][index]) / (data[0][index + 1] - data[0][index]);
// if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
// if (timeFactor != 1.0) {
// outputStream.print(sampleTimes[i]/timeFactor);
// } else {
outputStream.print(sampleTimes[i]);
// }
for (int j = 0; j < testSpec.getVarsList().length; j++) {
double interpolatedValue = 0.0;
double[] speciesVals = null;
double[] times = null;
// Currently using 2nd order interpolation
if (index == 0) {
// can only do 1st order interpolation
times = new double[] { data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1] };
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
} else if (index >= 1 && index <= origDataLength - 3) {
double val_1 = Math.abs(sampleTimes[i] - data[0][index - 1]);
double val_2 = Math.abs(sampleTimes[i] - data[0][index + 2]);
if (val_1 < val_2) {
times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
} else {
times = new double[] { data[0][index], data[0][index + 1], data[0][index + 2] };
speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1], data[j + 1][index + 2] };
}
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
} else {
times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
}
// // Currently using 1st order interpolation
// times = new double[] { data[0][index], data[0][index+1] };
// speciesVals = new double[] { data[j+1][index], data[j+1][index+1] };
// interpolatedValue = taylorInterpolation(sampleTimes[i], times, speciesVals);
// interpolatedValue = interpolatedValue * unitFactor.getExpression().evaluateConstant(); //earlier, hack unitfactor = 0.000001
// System.out.println("Sample time: " + sampleTimes[i] + ", between time[" + index + "]=" + data[0][index]+" and time["+(index+1)+"]="+(data[0][index+1])+", interpolated = "+interpolatedValue);
outputStream.print("," + interpolatedValue);
}
outputStream.println();
}
}
outputStream.close();
return outputFile;
} catch (RuntimeException e) {
e.printStackTrace(System.out);
// rethrow without losing context
throw e;
} catch (Exception e) {
e.printStackTrace(System.out);
throw new SolverException(e.getMessage(), e);
}
}
Aggregations