use of cbit.vcell.math.Variable in project vcell by virtualcell.
the class FiniteVolumeFileWriter method writeFastSystem.
/*private void writeDataProcessor() throws DataAccessException, IOException, MathException, DivideByZeroException, ExpressionException {
Simulation simulation = simTask.getSimulation();
DataProcessingInstructions dpi = simulation.getDataProcessingInstructions();
if (dpi == null) {
printWriter.println("DATA_PROCESSOR_BEGIN " + DataProcessingInstructions.ROI_TIME_SERIES);
printWriter.println("DATA_PROCESSOR_END");
printWriter.println();
return;
}
FieldDataIdentifierSpec fdis = dpi.getSampleImageFieldData(simulation.getVersion().getOwner());
if (fdis == null) {
throw new DataAccessException("Can't find sample image in data processing instructions");
}
String secondarySimDataDir = PropertyLoader.getProperty(PropertyLoader.secondarySimDataDirProperty, null);
DataSetControllerImpl dsci = new DataSetControllerImpl(new NullSessionLog(),null,userDirectory.getParentFile(),secondarySimDataDir == null ? null : new File(secondarySimDataDir));
CartesianMesh origMesh = dsci.getMesh(fdis.getExternalDataIdentifier());
SimDataBlock simDataBlock = dsci.getSimDataBlock(null,fdis.getExternalDataIdentifier(), fdis.getFieldFuncArgs().getVariableName(), fdis.getFieldFuncArgs().getTime().evaluateConstant());
VariableType varType = fdis.getFieldFuncArgs().getVariableType();
VariableType dataVarType = simDataBlock.getVariableType();
if (!varType.equals(VariableType.UNKNOWN) && !varType.equals(dataVarType)) {
throw new IllegalArgumentException("field function variable type (" + varType.getTypeName() + ") doesn't match real variable type (" + dataVarType.getTypeName() + ")");
}
double[] origData = simDataBlock.getData();
String filename = SimulationJob.createSimulationJobID(Simulation.createSimulationID(simulation.getKey()), simulationJob.getJobIndex()) + FieldDataIdentifierSpec.getDefaultFieldDataFileNameForSimulation(fdis.getFieldFuncArgs());
File fdatFile = new File(userDirectory, filename);
DataSet.writeNew(fdatFile,
new String[] {fdis.getFieldFuncArgs().getVariableName()},
new VariableType[]{simDataBlock.getVariableType()},
new ISize(origMesh.getSizeX(),origMesh.getSizeY(),origMesh.getSizeZ()),
new double[][]{origData});
printWriter.println("DATA_PROCESSOR_BEGIN " + dpi.getScriptName());
printWriter.println(dpi.getScriptInput());
printWriter.println("SampleImageFile " + fdis.getFieldFuncArgs().getVariableName() + " " + fdis.getFieldFuncArgs().getTime().infix() + " " + fdatFile);
printWriter.println("DATA_PROCESSOR_END");
printWriter.println();
}
*/
/**
*# fast system dimension num_dependents
*FAST_SYSTEM_BEGIN 2 2
*INDEPENDENT_VARIALBES rf r
*DEPENDENT_VARIALBES rB rfB
*
*PSEUDO_CONSTANT_BEGIN
*__C0 (rfB + rf);
*__C1 (r + rB);
*PSEUDO_CONSTANT_END
*
*FAST_RATE_BEGIN
* - ((0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0) * rf) - (0.1 * ( - rf + __C0)));
*((0.02 * r * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)) - (0.1 * ( - r + __C1)));
*FAST_RATE_END
*
*FAST_DEPENDENCY_BEGIN
*rB ( - r + __C1);
*rfB ( - rf + __C0);
*FAST_DEPENDENCY_END
*
*JACOBIAN_BEGIN
* - (0.1 + (0.02 * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))) * rf) + (0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)));
* - (0.02 * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))) * rf);
*(0.02 * r * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0)))));
*(0.1 + (0.02 * ( - ( - r + __C1) - ( - rf + __C0) + (20.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))) + _VCell_FieldData_0)) + (0.02 * r * (1.0 + (0.0 * ((x > -3.0) && (x < 3.0) && (y > -5.0) && (y < 5.0))))));
*JACOBIAN_END
*
*FAST_SYSTEM_END
* @throws ExpressionException
* @throws MathException
*/
private void writeFastSystem(SubDomain subDomain) throws MathException, ExpressionException {
VariableDomain variableDomain = (subDomain instanceof CompartmentSubDomain) ? VariableDomain.VARIABLEDOMAIN_VOLUME : VariableDomain.VARIABLEDOMAIN_MEMBRANE;
FastSystem fastSystem = subDomain.getFastSystem();
if (fastSystem == null) {
return;
}
FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, simTask.getSimulationJob().getSimulationSymbolTable());
int numIndep = fs_analyzer.getNumIndependentVariables();
int numDep = fs_analyzer.getNumDependentVariables();
int numPseudo = fs_analyzer.getNumPseudoConstants();
printWriter.println("# fast system dimension num_dependents");
printWriter.println("FAST_SYSTEM_BEGIN " + numIndep + " " + numDep);
if (numIndep != 0) {
printWriter.print("INDEPENDENT_VARIALBES ");
Enumeration<Variable> enum1 = fs_analyzer.getIndependentVariables();
while (enum1.hasMoreElements()) {
Variable var = enum1.nextElement();
printWriter.print(var.getName() + " ");
}
printWriter.println();
}
if (numDep != 0) {
printWriter.print("DEPENDENT_VARIALBES ");
Enumeration<Variable> enum1 = fs_analyzer.getDependentVariables();
while (enum1.hasMoreElements()) {
Variable var = enum1.nextElement();
printWriter.print(var.getName() + " ");
}
printWriter.println();
}
printWriter.println();
if (numPseudo != 0) {
printWriter.println("PSEUDO_CONSTANT_BEGIN");
Enumeration<PseudoConstant> enum1 = fs_analyzer.getPseudoConstants();
while (enum1.hasMoreElements()) {
PseudoConstant pc = enum1.nextElement();
printWriter.println(pc.getName() + " " + subsituteExpression(pc.getPseudoExpression(), fs_analyzer, variableDomain).infix() + ";");
}
printWriter.println("PSEUDO_CONSTANT_END");
printWriter.println();
}
if (numIndep != 0) {
printWriter.println("FAST_RATE_BEGIN");
Enumeration<Expression> enum1 = fs_analyzer.getFastRateExpressions();
while (enum1.hasMoreElements()) {
Expression exp = enum1.nextElement();
printWriter.println(subsituteExpression(exp, fs_analyzer, variableDomain).infix() + ";");
}
printWriter.println("FAST_RATE_END");
printWriter.println();
}
if (numDep != 0) {
printWriter.println("FAST_DEPENDENCY_BEGIN");
Enumeration<Expression> enum_exp = fs_analyzer.getDependencyExps();
Enumeration<Variable> enum_var = fs_analyzer.getDependentVariables();
while (enum_exp.hasMoreElements()) {
Expression exp = enum_exp.nextElement();
Variable depVar = enum_var.nextElement();
printWriter.println(depVar.getName() + " " + subsituteExpression(exp, fs_analyzer, variableDomain).infix() + ";");
}
printWriter.println("FAST_DEPENDENCY_END");
printWriter.println();
}
if (numIndep != 0) {
printWriter.println("JACOBIAN_BEGIN");
Enumeration<Expression> enum_fre = fs_analyzer.getFastRateExpressions();
while (enum_fre.hasMoreElements()) {
Expression fre = enum_fre.nextElement();
Enumeration<Variable> enum_var = fs_analyzer.getIndependentVariables();
while (enum_var.hasMoreElements()) {
Variable var = enum_var.nextElement();
Expression exp = subsituteExpression(fre, fs_analyzer, variableDomain).flatten();
Expression differential = exp.differentiate(var.getName());
printWriter.println(subsituteExpression(differential, fs_analyzer, variableDomain).infix() + ";");
}
}
printWriter.println("JACOBIAN_END");
printWriter.println();
}
printWriter.println("FAST_SYSTEM_END");
printWriter.println();
}
use of cbit.vcell.math.Variable in project vcell by virtualcell.
the class FiniteVolumeFileWriter method writeVariables.
/**
*# Variables : type name unit time_dependent_flag advection_flag solve_whole_mesh_flag solve_regions
*VARIABLE_BEGIN
*VOLUME_ODE rB uM
*VOLUME_PDE rf uM false false
*VOLUME_PDE r uM false false
*VOLUME_ODE rfB uM
*VARIABLE_END
* @throws MathException
* @throws ExpressionException
* @throws IOException
*/
private void writeVariables() throws MathException, ExpressionException, IOException {
SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
printWriter.println("# Variables : type name domain time_dependent_flag advection_flag grad_flag solve_whole_mesh_flag solve_regions");
printWriter.println(FVInputFileKeyword.VARIABLE_BEGIN);
MathDescription mathDesc = simSymbolTable.getSimulation().getMathDescription();
Variable[] vars = simSymbolTable.getVariables();
ArrayList<RandomVariable> rvList = new ArrayList<RandomVariable>();
for (int i = 0; i < vars.length; i++) {
String varName = vars[i].getName();
String domainName = vars[i].getDomain() == null ? null : vars[i].getDomain().getName();
if (vars[i] instanceof VolumeRandomVariable || vars[i] instanceof MembraneRandomVariable) {
rvList.add((RandomVariable) vars[i]);
} else if (vars[i] instanceof VolVariable) {
if (bChomboSolver && domainName == null) {
throw new MathException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " requires that every variable is defined in a single domain");
}
VolVariable volVar = (VolVariable) vars[i];
if (mathDesc.isPDE(volVar)) {
boolean hasTimeVaryingDiffusionOrAdvection = simSymbolTable.hasTimeVaryingDiffusionOrAdvection(volVar);
final boolean hasVelocity = mathDesc.hasVelocity(volVar);
final boolean hasGradient = mathDesc.hasGradient(volVar);
if (mathDesc.isPdeSteady(volVar)) {
printWriter.print("VOLUME_PDE_STEADY ");
} else {
printWriter.print("VOLUME_PDE ");
}
printWriter.print(varName + " " + domainName + " " + hasTimeVaryingDiffusionOrAdvection + " " + hasVelocity + " " + hasGradient);
} else {
printWriter.print("VOLUME_ODE " + varName + " " + domainName);
}
if (domainName == null) {
Vector<SubDomain> listOfSubDomains = new Vector<SubDomain>();
int totalNumCompartments = 0;
Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
while (subDomainEnum.hasMoreElements()) {
SubDomain subDomain = subDomainEnum.nextElement();
if (subDomain instanceof CompartmentSubDomain) {
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
totalNumCompartments++;
Equation varEquation = subDomain.getEquation(vars[i]);
if (varEquation != null) {
if (!(varEquation instanceof PdeEquation) || !((PdeEquation) varEquation).isDummy(simSymbolTable, compartmentSubDomain)) {
listOfSubDomains.add(compartmentSubDomain);
}
}
}
}
if ((totalNumCompartments == listOfSubDomains.size()) || (listOfSubDomains.size() == 0 && simTask.getSimulation().getSolverTaskDescription().getSolverDescription().equals(SolverDescription.SundialsPDE))) {
printWriter.print(" true");
} else {
printWriter.print(" false");
for (int j = 0; j < listOfSubDomains.size(); j++) {
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) listOfSubDomains.elementAt(j);
printWriter.print(" " + compartmentSubDomain.getName());
}
}
printWriter.println();
} else {
printWriter.println(" false " + domainName);
}
} else if (vars[i] instanceof VolumeParticleVariable) {
printWriter.println(FVInputFileKeyword.VOLUME_PARTICLE + " " + varName + " " + domainName);
} else if (vars[i] instanceof MembraneParticleVariable) {
printWriter.println(FVInputFileKeyword.MEMBRANE_PARTICLE + " " + varName + " " + domainName);
} else if (vars[i] instanceof VolumeRegionVariable) {
printWriter.println("VOLUME_REGION " + varName + " " + domainName);
} else if (vars[i] instanceof MemVariable) {
if (bChomboSolver && domainName == null) {
throw new MathException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " requires that every variable is defined in a single domain");
}
MemVariable memVar = (MemVariable) vars[i];
if (mathDesc.isPDE(memVar)) {
printWriter.println("MEMBRANE_PDE " + varName + " " + domainName + " " + simSymbolTable.hasTimeVaryingDiffusionOrAdvection(memVar));
} else {
printWriter.println("MEMBRANE_ODE " + varName + " " + domainName);
}
} else if (vars[i] instanceof MembraneRegionVariable) {
printWriter.println("MEMBRANE_REGION " + varName + " " + domainName);
} else if (vars[i] instanceof FilamentVariable) {
throw new RuntimeException("Filament application not supported yet");
}
}
int numRandomVariables = rvList.size();
if (numRandomVariables > 0) {
ISize samplingSize = simTask.getSimulation().getMeshSpecification().getSamplingSize();
String[] varNameArr = new String[numRandomVariables];
VariableType[] varTypeArr = new VariableType[numRandomVariables];
double[][] dataArr = new double[numRandomVariables][];
for (int i = 0; i < numRandomVariables; i++) {
RandomVariable rv = rvList.get(i);
varNameArr[i] = rv.getName();
int numRandomNumbers = 0;
if (rv instanceof VolumeRandomVariable) {
printWriter.print("VOLUME_RANDOM");
varTypeArr[i] = VariableType.VOLUME;
numRandomNumbers = samplingSize.getXYZ();
} else if (rv instanceof MembraneRandomVariable) {
printWriter.print("MEMBRANE_RANDOM");
varTypeArr[i] = VariableType.MEMBRANE;
numRandomNumbers = resampledGeometry.getGeometrySurfaceDescription().getSurfaceCollection().getTotalPolygonCount();
} else {
throw new RuntimeException("Unknown RandomVariable type");
}
printWriter.println(" " + varNameArr[i]);
dataArr[i] = generateRandomNumbers(rv, numRandomNumbers);
}
File rvFile = new File(workingDirectory, simTask.getSimulationJobID() + RANDOM_VARIABLE_FILE_EXTENSION);
DataSet.writeNew(rvFile, varNameArr, varTypeArr, samplingSize, dataArr);
}
printWriter.println(FVInputFileKeyword.VARIABLE_END);
printWriter.println();
}
use of cbit.vcell.math.Variable in project vcell by virtualcell.
the class OptXmlWriter method getObjectiveFunctionXML.
public static Element getObjectiveFunctionXML(OptimizationSpec optimizationSpec) {
Element objectiveFunctionElement = new Element(OptXmlTags.ObjectiveFunction_Tag);
ObjectiveFunction objectiveFunction = optimizationSpec.getObjectiveFunction();
if (objectiveFunction instanceof ExplicitObjectiveFunction) {
ExplicitObjectiveFunction explicitObjectiveFunction = (ExplicitObjectiveFunction) objectiveFunction;
objectiveFunctionElement.setAttribute(OptXmlTags.ObjectiveFunctionType_Attr, OptXmlTags.ObjectiveFunctionType_Attr_Explicit);
Element expressionElement = new Element(OptXmlTags.Expression_Tag);
expressionElement.addContent(explicitObjectiveFunction.getExpression().infix());
objectiveFunctionElement.addContent(expressionElement);
} else if (objectiveFunction instanceof ExplicitFitObjectiveFunction) {
ExplicitFitObjectiveFunction explicitFitObjectiveFunction = (ExplicitFitObjectiveFunction) objectiveFunction;
objectiveFunctionElement.setAttribute(OptXmlTags.ObjectiveFunctionType_Attr, OptXmlTags.ObjectiveFunctionType_Attr_ExplicitFit);
// add expression list
Element expressionListElement = new Element(OptXmlTags.ExpressionList_Tag);
objectiveFunctionElement.addContent(expressionListElement);
ExplicitFitObjectiveFunction.ExpressionDataPair[] expDataPairs = explicitFitObjectiveFunction.getExpressionDataPairs();
for (int i = 0; i < expDataPairs.length; i++) {
Element expressionElement = new Element(OptXmlTags.Expression_Tag);
expressionElement.setAttribute(OptXmlTags.ExpRefDataColumnIndex_Attr, expDataPairs[i].getReferenceDataIndex() + "");
expressionElement.addContent(expDataPairs[i].getFitFunction().infix());
expressionListElement.addContent(expressionElement);
}
if (explicitFitObjectiveFunction.getReferenceData() instanceof SimpleReferenceData) {
Element dataElement = getDataXML(explicitFitObjectiveFunction.getReferenceData());
objectiveFunctionElement.addContent(dataElement);
} else {
throw new OptimizationException("Optimization XML Writer exports SimpleReferenceData only.");
}
} else if (objectiveFunction instanceof PdeObjectiveFunction) {
PdeObjectiveFunction pdeObjectiveFunction = (PdeObjectiveFunction) objectiveFunction;
objectiveFunctionElement.setAttribute(OptXmlTags.ObjectiveFunctionType_Attr, OptXmlTags.ObjectiveFunctionType_Attr_PDE);
SpatialReferenceData refData = pdeObjectiveFunction.getReferenceData();
Element dataElement = getDataXML(refData);
objectiveFunctionElement.addContent(dataElement);
Element modelElement = getModelXML(pdeObjectiveFunction, optimizationSpec.getParameterNames());
objectiveFunctionElement.addContent(modelElement);
Simulation tempSim = new Simulation(pdeObjectiveFunction.getMathDescription());
SimulationSymbolTable simSymbolTable = new SimulationSymbolTable(tempSim, 0);
//
for (int i = 1; i < refData.getNumVariables(); i++) {
Element modelMappingElement = new Element(OptXmlTags.ModelMapping_Tag);
modelMappingElement.setAttribute(OptXmlTags.ModelMappingDataColumn_Attr, refData.getVariableNames()[i]);
modelMappingElement.setAttribute(OptXmlTags.ModelMappingWeight_Attr, String.valueOf(refData.getVariableWeights()[i]));
Expression mappingExpression = null;
try {
Variable var = pdeObjectiveFunction.getMathDescription().getVariable(refData.getVariableNames()[i]);
if (var instanceof Function) {
Expression exp = new Expression(var.getExpression());
exp.bindExpression(simSymbolTable);
exp = MathUtilities.substituteFunctions(exp, simSymbolTable);
mappingExpression = exp.flatten();
} else {
mappingExpression = new Expression(var.getName());
}
} catch (ExpressionException e) {
e.printStackTrace();
throw new OptimizationException(e.getMessage());
}
modelMappingElement.addContent(mappingExpression.infix());
objectiveFunctionElement.addContent(modelMappingElement);
}
}
return objectiveFunctionElement;
}
use of cbit.vcell.math.Variable 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.math.Variable in project vcell by virtualcell.
the class MathTestingUtilities method substituteWithExactSolution.
/**
* Insert the method's description here.
* Creation date: (1/24/2003 10:18:14 AM)
* @return cbit.vcell.parser.Expression
* @param origExp cbit.vcell.parser.Expression
* @param subDomain cbit.vcell.math.SubDomain
*/
private static Expression substituteWithExactSolution(Expression origExp, MembraneSubDomain subDomain, MathDescription exactMathDesc) throws ExpressionException {
Expression substitutedExp = new Expression(origExp);
substitutedExp.bindExpression(exactMathDesc);
substitutedExp = MathUtilities.substituteFunctions(substitutedExp, exactMathDesc);
substitutedExp.bindExpression(null);
substitutedExp = substitutedExp.flatten();
substitutedExp.bindExpression(exactMathDesc);
String[] symbols = substitutedExp.getSymbols();
for (int i = 0; symbols != null && i < symbols.length; i++) {
Variable var = (Variable) substitutedExp.getSymbolBinding(symbols[i]);
if (var instanceof MemVariable) {
String exactVarName = var.getName() + "_" + subDomain.getName() + "_exact";
substitutedExp.substituteInPlace(new Expression(var.getName()), new Expression(exactVarName));
} else if (var instanceof InsideVariable) {
String exactVarName = var.getName() + "_" + subDomain.getInsideCompartment().getName() + "_exact";
substitutedExp.substituteInPlace(new Expression(var.getName()), new Expression(exactVarName));
} else if (var instanceof OutsideVariable) {
String exactVarName = var.getName() + "_" + subDomain.getOutsideCompartment().getName() + "_exact";
substitutedExp.substituteInPlace(new Expression(var.getName()), new Expression(exactVarName));
} else if (var instanceof VolumeRegionVariable || var instanceof MemVariable || var instanceof MembraneRegionVariable || var instanceof FilamentVariable || var instanceof FilamentRegionVariable) {
throw new RuntimeException("variable substitution not yet implemented for Variable type " + var.getClass().getName() + "(" + var.getName() + ")");
}
}
substitutedExp.bindExpression(null);
return substitutedExp;
}
Aggregations