Example 11 with SubDomain
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class ParticleMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
getSimulationContext().checkValidity();
if (getSimulationContext().getGeometry().getDimension() == 0) {
throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
}
StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = getSimulationContext().getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
if (reactionSpecs[i].isFast()) {
throw new MappingException("fast reactions not supported for particle models");
}
rsList.add(reactionSpecs[i].getReactionStep());
}
}
ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
rsList.copyInto(reactionSteps);
//
for (int i = 0; i < reactionSteps.length; i++) {
Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
// //
// // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
// //
// Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
// for (int i = 0; i < structures.length; i++){
// StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
// if (sm==null || (sm.getGeometryClass() == null)){
// throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
// }
// if (sm.getUnitSizeParameter()!=null){
// Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
// if(unitSizeExp != null)
// {
// try {
// double unitSize = unitSizeExp.evaluateConstant();
// if (unitSize != 1.0){
// throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
// }
// }catch (ExpressionException e){
// e.printStackTrace(System.out);
// throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
// }
// }
// }
// }
// {
// GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
// for (int i = 0; i < geometryClass.length; i++){
// Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
// }
// }
// }
// deals with model parameters
Model model = getSimulationContext().getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
ModelParameter[] modelParameters = model.getModelParameters();
// populate in globalParameterVariants hashtable
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = getSimulationContext().getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(getSimulationContext().getName() + "_generated");
}
//
// volume particle variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof ParticleVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
GeometryClass geometryClass = membraneMapping.getGeometryClass();
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
// Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = null;
Expression initExpr = null;
if (getSimulationContext().isUsingConcentration()) {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
initExpr = new Expression(initParm.getExpression());
// if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
// initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
// }
} else {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
initExpr = new Expression(initParm.getExpression());
}
if (initExpr != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = getSimulationContext().getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// constant species (either function or constant)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof Constant) {
varHash.addVariable(scm.getVariable());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (getSimulationContext().getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (getSimulationContext().getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(getSimulationContext().getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// create subdomains (volume and surfaces)
//
GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
for (int k = 0; k < geometryClasses.length; k++) {
if (geometryClasses[k] instanceof SubVolume) {
SubVolume subVolume = (SubVolume) geometryClasses[k];
//
// get priority of subDomain
//
// now does not have to match spatial feature, *BUT* needs to be unique
int priority = k;
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
if (getSimulationContext().getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (getSimulationContext().getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
} else if (geometryClasses[k] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
// determine membrane inside and outside subvolume
// this preserves backward compatibility so that membrane subdomain
// inside and outside correspond to structure hierarchy when present
Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
SubVolume innerSubVolume = ret.one;
SubVolume outerSubVolume = ret.two;
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
}
}
//
// create Particle Contexts for all Particle Variables
//
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
//
// initial distribution of particles
//
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
ParticleInitialCondition pic = null;
if (getSimulationContext().isUsingConcentration()) {
Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
if (particleVariable instanceof VolumeParticleVariable) {
initialDistribution = Expression.mult(initialDistribution, unitFactor);
}
pic = new ParticleInitialConditionConcentration(initialDistribution);
} else {
Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
if (initialCount == null) {
throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
}
Expression locationX = new Expression("u");
Expression locationY = new Expression("u");
Expression locationZ = new Expression("u");
pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
}
particleInitialConditions.add(pic);
//
// diffusion
//
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
Expression driftXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velX_quantities.length == 1 && numRegions == 1) {
driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velY_quantities.length == 1 && numRegions == 1) {
driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
GeometryClass myGC = sm.getGeometryClass();
if (myGC == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
subDomain.addParticleProperties(particleProperties);
}
}
for (ReactionStep reactionStep : reactionSteps) {
Kinetics kinetics = reactionStep.getKinetics();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
KineticsParameter reactionRateParameter = null;
if (kinetics instanceof LumpedKinetics) {
reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
} else {
reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
}
// macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
if (radiusExp != null) {
Expression expCopy = new Expression(radiusExp);
try {
MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ": '" + radiusExp.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
}
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
if (rp instanceof Reactant) {
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
}
}
} else if (rp instanceof Product) {
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createCreateAction(particle));
}
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
} else // other type of reactions
{
/* check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
rate constants are important in calculating the probability rate.
for Mass Action, we use KForward and KReverse,
for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
*/
Expression forwardRate = null;
Expression reverseRate = null;
// Using the MassActionFunction to write out the math description
MassActionSolver.MassActionFunction maFunc = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
} else {
if (maFunc.getForwardRate() != null) {
forwardRate = maFunc.getForwardRate();
}
if (maFunc.getReverseRate() != null) {
reverseRate = maFunc.getReverseRate();
}
}
}
if (maFunc != null) {
// if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
// we process it as forward reaction
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
List<Action> reverseActions = new ArrayList<Action>();
List<ReactionParticipant> reactants = maFunc.getReactants();
List<ReactionParticipant> products = maFunc.getProducts();
for (ReactionParticipant rp : reactants) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createCreateAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
for (ReactionParticipant rp : products) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createCreateAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createDestroyAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
//
// There are two unit conversions required:
//
// 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
// (where dim is 2 for membrane reactions and 3 for volume reactions)
//
// for forward rates:
// 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane reactants and 3 for volume reactants)
//
// or
//
// for reverse rates:
// 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane products and 3 for volume products)
//
RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
if (forwardRate != null) {
VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant reactant : maFunc.getReactants()) {
VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynReactantUnit = null;
if (bForceContinuous) {
// reactant is continuous (vcell units)
smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
} else {
// reactant is a particle (smoldyn units)
RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
}
// keep track of units of all reactants
smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
// accumulate unit factors for all reactants
forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
}
forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
// create particle jump process
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
}
// end of forward rate not null
if (reverseRate != null) {
VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant product : maFunc.getProducts()) {
VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynProductUnit = null;
if (bForceContinuous) {
smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
} else {
RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
}
// keep track of units of all products
smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
// accumulate unit factors for all products
reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
}
reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (reverseActions.size() > 0) {
ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(reverseProcess);
}
}
// end of reverse rate not null
}
// end of maFunc not null
}
// end of reaction step for loop
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
lg.warn(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
if (lg.isDebugEnabled()) {
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition)
ReactionStep(cbit.vcell.model.ReactionStep)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Kinetics(cbit.vcell.model.Kinetics)
DistributedKinetics(cbit.vcell.model.DistributedKinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
GeometryClass(cbit.vcell.geometry.GeometryClass)
Action(cbit.vcell.math.Action)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Pair(org.vcell.util.Pair)
ParticleInitialConditionConcentration(cbit.vcell.math.ParticleProperties.ParticleInitialConditionConcentration)
ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)
Expression(cbit.vcell.parser.Expression)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MathException(cbit.vcell.math.MathException)
Model(cbit.vcell.model.Model)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
MassActionSolver(cbit.vcell.model.MassActionSolver)
ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)
Example 12 with SubDomain
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class RulebasedMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
@Override
protected void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// use local variable instead of using getter all the time.
SimulationContext simContext = getSimulationContext();
GeometryClass geometryClass = simContext.getGeometry().getGeometrySpec().getSubVolumes()[0];
Domain domain = new Domain(geometryClass);
// local structure mapping list
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
// We have to check if all the reactions are able to transform to stochastic jump processes before generating the math.
String stochChkMsg = simContext.getModel().isValidForStochApp();
if (!(stochChkMsg.equals(""))) {
throw new ModelException("Problem updating math description: " + simContext.getName() + "\n" + stochChkMsg);
}
simContext.checkValidity();
//
if (simContext.getGeometry().getDimension() > 0) {
throw new MappingException("rule-based particle math mapping not implemented for spatial geometry - dimension >= 1");
}
//
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = simContext.getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
//
Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
for (int i = 0; i < structures.length; i++) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
if (sm == null || (sm instanceof FeatureMapping && ((FeatureMapping) sm).getGeometryClass() == null)) {
throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subVolume");
}
if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
try {
if (volFractExp != null) {
double volFract = volFractExp.evaluateConstant();
if (volFract >= 1.0) {
throw new MappingException("model structure '" + (getSimulationContext().getModel().getStructureTopology().getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0"));
}
}
} catch (ExpressionException e) {
e.printStackTrace(System.out);
}
}
}
SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
Structure[] mappedStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolumes[i]);
if (mappedStructures == null || mappedStructures.length == 0) {
throw new MappingException("geometry subVolume '" + subVolumes[i].getName() + "' not mapped from a model structure");
}
}
//
// gather only those reactionRules that are not "excluded"
//
ArrayList<ReactionRule> rrList = new ArrayList<ReactionRule>();
for (ReactionRuleSpec reactionRuleSpec : simContext.getReactionContext().getReactionRuleSpecs()) {
if (!reactionRuleSpec.isExcluded()) {
rrList.add(reactionRuleSpec.getReactionRule());
}
}
//
for (ReactionRule reactionRule : rrList) {
UnresolvedParameter[] unresolvedParameters = reactionRule.getKineticLaw().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException("In Application '" + simContext.getName() + "', " + reactionRule.getDisplayType() + " '" + reactionRule.getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = simContext.getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(simContext.getName() + "_generated");
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates
//
VariableHash varHash = new VariableHash();
//
// conversion factors
//
Model model = simContext.getModel();
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof StochVolVariable) {
varHash.addVariable(scm.getVariable());
}
}
// deals with model parameters
ModelParameter[] modelParameters = simContext.getModel().getModelParameters();
for (int j = 0; j < modelParameters.length; j++) {
Expression expr = getSubstitutedExpr(modelParameters[j].getExpression(), true, false);
expr = getIdentifierSubstitutions(expr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), expr, geometryClass));
}
// added July 2009, ElectricalStimulusParameter electric mapping tab
ElectricalStimulus[] elecStimulus = simContext.getElectricalStimuli();
if (elecStimulus.length > 0) {
throw new MappingException("Modles with electrophysiology are not supported for stochastic applications.");
}
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
Parameter initialVoltageParm = memMapping.getInitialVoltageParameter();
try {
Expression exp = initialVoltageParm.getExpression();
exp.evaluateConstant();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
}
}
}
//
for (ReactionRule reactionRule : rrList) {
// if (reactionRule.getKineticLaw() instanceof LumpedKinetics){
// throw new RuntimeException("Lumped Kinetics not yet supported for RuleBased Modeling");
// }
LocalParameter[] parameters = reactionRule.getKineticLaw().getLocalParameters();
for (LocalParameter parameter : parameters) {
//
if ((parameter.getRole() == RbmKineticLawParameterType.RuleRate)) {
continue;
}
//
if (!reactionRule.isReversible() && parameter.getRole() == RbmKineticLawParameterType.MassActionReverseRate) {
continue;
}
Expression expr = getSubstitutedExpr(parameter.getExpression(), true, false);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameter, geometryClass), getIdentifierSubstitutions(expr, parameter.getUnitDefinition(), geometryClass), geometryClass));
}
}
// the parameter "Size" is already put into mathsymbolmapping in refreshSpeciesContextMapping()
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm.getExpression() != null) {
try {
double value = parm.getExpression().evaluateConstant();
varHash.addVariable(new Constant(getMathSymbol(parm, sm.getGeometryClass()), new Expression(value)));
} catch (ExpressionException e) {
// varHash.addVariable(new Function(getMathSymbol0(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
e.printStackTrace(System.out);
throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
}
}
}
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
addInitialConditions(domain, speciesContextSpecs, varHash);
//
if (simContext.getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("Geometry must be defined in Application " + simContext.getName());
}
//
// create subDomains
//
SubVolume subVolume = simContext.getGeometry().getGeometrySpec().getSubVolumes()[0];
SubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), 0);
mathDesc.addSubDomain(subDomain);
//
// define all molecules and unique species patterns (add molecules to mathDesc and speciesPatterns to varHash).
//
HashMap<SpeciesPattern, VolumeParticleSpeciesPattern> speciesPatternMap = addSpeciesPatterns(domain, rrList);
HashSet<VolumeParticleSpeciesPattern> uniqueParticleSpeciesPatterns = new HashSet<>(speciesPatternMap.values());
for (VolumeParticleSpeciesPattern volumeParticleSpeciesPattern : uniqueParticleSpeciesPatterns) {
varHash.addVariable(volumeParticleSpeciesPattern);
}
//
// define observables (those explicitly declared and those corresponding to seed species.
//
List<ParticleObservable> observables = addObservables(geometryClass, domain, speciesPatternMap);
for (ParticleObservable particleObservable : observables) {
varHash.addVariable(particleObservable);
}
try {
addParticleJumpProcesses(varHash, geometryClass, subDomain, speciesPatternMap);
} catch (PropertyVetoException e1) {
e1.printStackTrace();
throw new MappingException(e1.getMessage(), e1);
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter || fieldMathMappingParameters[i] instanceof ObservableConcentrationParameter) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
for (SpeciesContext sc : model.getSpeciesContexts()) {
if (!sc.hasSpeciesPattern()) {
throw new MappingException("species " + sc.getName() + " has no molecular pattern");
}
VolumeParticleSpeciesPattern volumeParticleSpeciesPattern = speciesPatternMap.get(sc.getSpeciesPattern());
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleProperties.ParticleInitialCondition>();
// initial conditions from scs
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
Parameter initialCountParameter = scs.getInitialCountParameter();
Expression e = getIdentifierSubstitutions(new Expression(initialCountParameter, getNameScope()), initialCountParameter.getUnitDefinition(), geometryClass);
particleInitialConditions.add(new ParticleInitialConditionCount(e, new Expression(0.0), new Expression(0.0), new Expression(0.0)));
ParticleProperties particleProperies = new ParticleProperties(volumeParticleSpeciesPattern, new Expression(0.0), new Expression(0.0), new Expression(0.0), new Expression(0.0), particleInitialConditions);
subDomain.addParticleProperties(particleProperies);
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
if (fieldMathMappingParameters[i] instanceof ObservableConcentrationParameter) {
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
System.out.println(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
}
Also used :
MathDescription(cbit.vcell.math.MathDescription)
ArrayList(java.util.ArrayList)
SpeciesContext(cbit.vcell.model.SpeciesContext)
ExpressionException(cbit.vcell.parser.ExpressionException)
PropertyVetoException(java.beans.PropertyVetoException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
SubVolume(cbit.vcell.geometry.SubVolume)
HashSet(java.util.HashSet)
ModelException(cbit.vcell.model.ModelException)
VolumeParticleSpeciesPattern(cbit.vcell.math.VolumeParticleSpeciesPattern)
PropertyVetoException(java.beans.PropertyVetoException)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition)
ParticleProperties(cbit.vcell.math.ParticleProperties)
VolumeParticleObservable(cbit.vcell.math.VolumeParticleObservable)
ParticleObservable(cbit.vcell.math.ParticleObservable)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
GeometryClass(cbit.vcell.geometry.GeometryClass)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
StochVolVariable(cbit.vcell.math.StochVolVariable)
Variable(cbit.vcell.math.Variable)
VariableHash(cbit.vcell.math.VariableHash)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
UnresolvedParameter(cbit.vcell.mapping.ParameterContext.UnresolvedParameter)
VolumeParticleSpeciesPattern(cbit.vcell.math.VolumeParticleSpeciesPattern)
SpeciesPattern(org.vcell.model.rbm.SpeciesPattern)
Structure(cbit.vcell.model.Structure)
StochVolVariable(cbit.vcell.math.StochVolVariable)
ReactionRule(cbit.vcell.model.ReactionRule)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
Parameter(cbit.vcell.model.Parameter)
UnresolvedParameter(cbit.vcell.mapping.ParameterContext.UnresolvedParameter)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)
Example 13 with SubDomain
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class MathTestingUtilities method getConstructedResultSet.
/**
* Insert the method's description here.
* Creation date: (1/17/2003 3:47:43 PM)
* @return cbit.vcell.solver.ode.ODESolverResultSet
* @param sim cbit.vcell.solver.Simulation
*/
public static ODESolverResultSet getConstructedResultSet(MathDescription mathDesc, double[] time) throws Exception {
if (mathDesc.getGeometry().getDimension() != 0) {
throw new RuntimeException("can only handle non-spatial simulations.");
}
Simulation sim = new Simulation(mathDesc);
SimulationSymbolTable simSymbolTable = new SimulationSymbolTable(sim, 0);
ODESolverResultSet resultSet = new ODESolverResultSet();
resultSet.addDataColumn(new ODESolverResultSetColumnDescription("t"));
for (int i = 0; i < time.length; i++) {
resultSet.addRow(new double[] { time[i] });
}
java.util.Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
String errorString = "Variable(s) : ";
while (subDomainEnum.hasMoreElements()) {
SubDomain subDomain = subDomainEnum.nextElement();
java.util.Enumeration<Equation> enumEquations = subDomain.getEquations();
while (enumEquations.hasMoreElements()) {
Equation equation = enumEquations.nextElement();
Expression constructedSolution = equation.getExactSolution();
if (constructedSolution != null) {
constructedSolution = new Expression(constructedSolution);
constructedSolution.bindExpression(simSymbolTable);
constructedSolution = simSymbolTable.substituteFunctions(constructedSolution);
constructedSolution = constructedSolution.flatten();
resultSet.addFunctionColumn(new FunctionColumnDescription(constructedSolution, equation.getVariable().getName(), null, equation.getVariable().getName(), false));
} else {
errorString = errorString + equation.getVariable().getName() + ", ";
}
}
}
if (!errorString.equals("Variable(s) : ")) {
throw new RuntimeException(errorString + " don't have a constructed solution");
}
return resultSet;
}
Also used :
SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable)
OdeEquation(cbit.vcell.math.OdeEquation)
PdeEquation(cbit.vcell.math.PdeEquation)
Equation(cbit.vcell.math.Equation)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Simulation(cbit.vcell.solver.Simulation)
Expression(cbit.vcell.parser.Expression)
ODESolverResultSet(cbit.vcell.solver.ode.ODESolverResultSet)
ODESolverResultSetColumnDescription(cbit.vcell.math.ODESolverResultSetColumnDescription)
FunctionColumnDescription(cbit.vcell.math.FunctionColumnDescription)
Example 14 with SubDomain
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class MathTestingUtilities method constructExactMath.
/**
* constructExactMath()
*
* take an equation of the form:
*
* d A
* --- = F(A,t)
* d t
*
* and create a new equation with a known exact solution 'A_exact' by adding a forcing function R(t)
*
* d A
* --- = F(A,t) + R(t)
* d t
*
* where:
*
* d A_exact
* R(t) = --------- - F(A_exact,t)
* d t
*
* solving for R(t) is done analytically.
*
* Creation date: (1/21/2003 10:47:54 AM)
* @return cbit.vcell.math.MathDescription
* @param mathDesc cbit.vcell.math.MathDescription
*/
public static MathDescription constructExactMath(MathDescription mathDesc, java.util.Random random, ConstructedSolutionTemplate constructedSolutionTemplate) throws ExpressionException, MathException, MappingException {
if (mathDesc.hasFastSystems()) {
throw new RuntimeException("SolverTest.constructExactMath() suppport for fastSystems not yet implemented.");
}
MathDescription exactMath = null;
try {
exactMath = (MathDescription) BeanUtils.cloneSerializable(mathDesc);
exactMath.setDescription("constructed exact solution from MathDescription (" + mathDesc.getName() + ")");
exactMath.setName("exact from " + mathDesc.getName());
} catch (Throwable e) {
e.printStackTrace(System.out);
throw new RuntimeException("error cloning MathDescription: " + e.getMessage());
}
//
// preload the VariableHash with existing Variables (and Constants,Functions,etc) and then sort all at once.
//
VariableHash varHash = new VariableHash();
Enumeration<Variable> enumVar = exactMath.getVariables();
while (enumVar.hasMoreElements()) {
varHash.addVariable(enumVar.nextElement());
}
java.util.Enumeration<SubDomain> subDomainEnum = exactMath.getSubDomains();
while (subDomainEnum.hasMoreElements()) {
SubDomain subDomain = subDomainEnum.nextElement();
Domain domain = new Domain(subDomain);
java.util.Enumeration<Equation> equationEnum = subDomain.getEquations();
if (subDomain instanceof MembraneSubDomain) {
MembraneSubDomain memSubDomain = (MembraneSubDomain) subDomain;
AnalyticSubVolume insideAnalyticSubVolume = (AnalyticSubVolume) exactMath.getGeometry().getGeometrySpec().getSubVolume(memSubDomain.getInsideCompartment().getName());
Function[] outwardNormalFunctions = getOutwardNormal(insideAnalyticSubVolume.getExpression(), "_" + insideAnalyticSubVolume.getName());
for (int i = 0; i < outwardNormalFunctions.length; i++) {
varHash.addVariable(outwardNormalFunctions[i]);
}
}
while (equationEnum.hasMoreElements()) {
Equation equation = equationEnum.nextElement();
if (equation.getExactSolution() != null) {
throw new RuntimeException("exact solution already exists");
}
Enumeration<Constant> origMathConstants = mathDesc.getConstants();
if (equation instanceof OdeEquation) {
OdeEquation odeEquation = (OdeEquation) equation;
Expression substitutedRateExp = substituteWithExactSolution(odeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
String varName = odeEquation.getVariable().getName();
String initName = null;
while (origMathConstants.hasMoreElements()) {
Constant constant = origMathConstants.nextElement();
if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
initName = constant.getName();
}
}
String exactName = varName + "_" + subDomain.getName() + "_exact";
String errorName = varName + "_" + subDomain.getName() + "_error";
String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
Expression exactExp = solutionTemplate.getTemplateExpression();
Expression errorExp = new Expression(exactName + " - " + varName);
Expression origRateExp = new Expression(odeEquation.getRateExpression());
Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " + " + exactTimeDerivativeName);
Constant[] constants = solutionTemplate.getConstants();
for (int i = 0; i < constants.length; i++) {
varHash.addVariable(constants[i]);
}
Expression initExp = new Expression(exactExp);
initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
varHash.addVariable(new Function(initName, initExp.flatten(), domain));
varHash.addVariable(new Function(exactName, exactExp, domain));
varHash.addVariable(new Function(errorName, errorExp, domain));
varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
varHash.addVariable(new Function(origRateName, origRateExp, domain));
varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
odeEquation.setRateExpression(newRate);
odeEquation.setInitialExpression(new Expression(initName));
odeEquation.setExactSolution(new Expression(exactName));
} else if (equation instanceof PdeEquation) {
PdeEquation pdeEquation = (PdeEquation) equation;
Expression substitutedRateExp = substituteWithExactSolution(pdeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
String varName = pdeEquation.getVariable().getName();
String initName = null;
while (origMathConstants.hasMoreElements()) {
Constant constant = origMathConstants.nextElement();
if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
initName = constant.getName();
}
}
String diffusionRateName = "_" + varName + "_" + subDomain.getName() + "_diffusionRate";
String exactName = varName + "_" + subDomain.getName() + "_exact";
String errorName = varName + "_" + subDomain.getName() + "_error";
String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
String exactDxName = "_" + varName + "_" + subDomain.getName() + "_exact_dx";
String exactDyName = "_" + varName + "_" + subDomain.getName() + "_exact_dy";
String exactDzName = "_" + varName + "_" + subDomain.getName() + "_exact_dz";
String exactDx2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dx2";
String exactDy2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dy2";
String exactDz2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dz2";
String exactLaplacianName = "_" + varName + "_" + subDomain.getName() + "_exact_laplacian";
Expression exactExp = solutionTemplate.getTemplateExpression();
Expression errorExp = new Expression(exactName + " - " + varName);
Expression origRateExp = new Expression(pdeEquation.getRateExpression());
Expression initExp = new Expression(exactExp);
initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
initExp = initExp.flatten();
Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
Expression exactDxExp = exactExp.differentiate("x").flatten();
Expression exactDx2Exp = exactDxExp.differentiate("x").flatten();
Expression exactDyExp = exactExp.differentiate("y").flatten();
Expression exactDy2Exp = exactDxExp.differentiate("y").flatten();
Expression exactDzExp = exactExp.differentiate("z").flatten();
Expression exactDz2Exp = exactDxExp.differentiate("z").flatten();
Expression exactLaplacianExp = Expression.add(Expression.add(exactDx2Exp, exactDy2Exp), exactDz2Exp).flatten();
Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " - ((" + diffusionRateName + ")*" + exactLaplacianName + ")" + " + " + exactTimeDerivativeName);
Constant[] constants = solutionTemplate.getConstants();
for (int i = 0; i < constants.length; i++) {
varHash.addVariable(constants[i]);
}
varHash.addVariable(new Function(initName, initExp, domain));
varHash.addVariable(new Function(diffusionRateName, new Expression(pdeEquation.getDiffusionExpression()), domain));
varHash.addVariable(new Function(exactName, exactExp, domain));
varHash.addVariable(new Function(errorName, errorExp, domain));
varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
varHash.addVariable(new Function(origRateName, origRateExp, domain));
varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
varHash.addVariable(new Function(exactDxName, exactDxExp, domain));
varHash.addVariable(new Function(exactDyName, exactDyExp, domain));
varHash.addVariable(new Function(exactDzName, exactDzExp, domain));
varHash.addVariable(new Function(exactDx2Name, exactDx2Exp, domain));
varHash.addVariable(new Function(exactDy2Name, exactDy2Exp, domain));
varHash.addVariable(new Function(exactDz2Name, exactDz2Exp, domain));
varHash.addVariable(new Function(exactLaplacianName, exactLaplacianExp, domain));
pdeEquation.setRateExpression(newRate);
pdeEquation.setInitialExpression(new Expression(initName));
pdeEquation.setDiffusionExpression(new Expression(diffusionRateName));
pdeEquation.setExactSolution(new Expression(exactName));
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
if (compartmentSubDomain.getBoundaryConditionXm().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryXm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryXm(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionXm().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryXm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDxName));
} else {
pdeEquation.setBoundaryXm(new Expression("-" + diffusionRateName + "*" + exactDxName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
}
if (compartmentSubDomain.getBoundaryConditionXp().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryXp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryXp(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionXp().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryXp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDxName));
} else {
pdeEquation.setBoundaryXp(new Expression(diffusionRateName + "*" + exactDxName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXp());
}
if (compartmentSubDomain.getBoundaryConditionYm().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryYm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryYm(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionYm().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryYm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDyName));
} else {
pdeEquation.setBoundaryYm(new Expression("-" + diffusionRateName + "*" + exactDyName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYm());
}
if (compartmentSubDomain.getBoundaryConditionYp().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryYp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryYp(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionYp().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryYp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDyName));
} else {
pdeEquation.setBoundaryYp(new Expression(diffusionRateName + "*" + exactDyName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYp());
}
if (compartmentSubDomain.getBoundaryConditionZm().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryZm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryZm(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionZm().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryZm();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDzName));
} else {
pdeEquation.setBoundaryZm(new Expression("-" + diffusionRateName + "*" + exactDzName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
}
if (compartmentSubDomain.getBoundaryConditionZp().isDIRICHLET()) {
Expression origExp = pdeEquation.getBoundaryZp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
} else {
pdeEquation.setBoundaryZp(exactExp);
}
} else if (compartmentSubDomain.getBoundaryConditionZp().isNEUMANN()) {
Expression origExp = pdeEquation.getBoundaryZp();
if (origExp != null) {
Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDzName));
} else {
pdeEquation.setBoundaryZp(new Expression(diffusionRateName + "*" + exactDzName));
}
} else {
throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionZp());
}
} else {
throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
}
}
if (subDomain instanceof MembraneSubDomain) {
MembraneSubDomain membraneSubDomain = (MembraneSubDomain) subDomain;
Enumeration<JumpCondition> enumJumpConditions = membraneSubDomain.getJumpConditions();
while (enumJumpConditions.hasMoreElements()) {
JumpCondition jumpCondition = enumJumpConditions.nextElement();
Expression origInfluxExp = jumpCondition.getInFluxExpression();
Expression origOutfluxExp = jumpCondition.getOutFluxExpression();
Expression substitutedInfluxExp = substituteWithExactSolution(origInfluxExp, membraneSubDomain, exactMath);
Expression substitutedOutfluxExp = substituteWithExactSolution(origOutfluxExp, membraneSubDomain, exactMath);
String varName = jumpCondition.getVariable().getName();
String origInfluxName = "_" + varName + "_" + subDomain.getName() + "_origInflux";
String origOutfluxName = "_" + varName + "_" + subDomain.getName() + "_origOutflux";
String substitutedInfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedInflux";
String substitutedOutfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedOutflux";
String diffusionRateInsideName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_diffusionRate";
String diffusionRateOutsideName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_diffusionRate";
String exactInsideDxName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dx";
String exactInsideDyName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dy";
String exactInsideDzName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dz";
String exactOutsideDxName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dx";
String exactOutsideDyName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dy";
String exactOutsideDzName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dz";
String outwardNormalXName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nx";
String outwardNormalYName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Ny";
String outwardNormalZName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nz";
String exactInfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactInflux";
String exactOutfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactOutflux";
Expression exactInfluxExp = new Expression(diffusionRateInsideName + " * (" + outwardNormalXName + "*" + exactInsideDxName + " + " + outwardNormalYName + "*" + exactInsideDyName + " + " + outwardNormalZName + "*" + exactInsideDzName + ")");
Expression exactOutfluxExp = new Expression("-" + diffusionRateOutsideName + " * (" + outwardNormalXName + "*" + exactOutsideDxName + " + " + outwardNormalYName + "*" + exactOutsideDyName + " + " + outwardNormalZName + "*" + exactOutsideDzName + ")");
Expression newInfluxExp = new Expression(origInfluxName + " - " + substitutedInfluxName + " + " + exactInfluxName);
Expression newOutfluxExp = new Expression(origOutfluxName + " - " + substitutedOutfluxName + " + " + exactOutfluxName);
varHash.addVariable(new Function(origInfluxName, origInfluxExp, domain));
varHash.addVariable(new Function(origOutfluxName, origOutfluxExp, domain));
varHash.addVariable(new Function(exactInfluxName, exactInfluxExp, domain));
varHash.addVariable(new Function(exactOutfluxName, exactOutfluxExp, domain));
varHash.addVariable(new Function(substitutedInfluxName, substitutedInfluxExp, domain));
varHash.addVariable(new Function(substitutedOutfluxName, substitutedOutfluxExp, domain));
jumpCondition.setInFlux(newInfluxExp);
jumpCondition.setOutFlux(newOutfluxExp);
}
}
}
exactMath.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
if (!exactMath.isValid()) {
throw new RuntimeException("generated Math is not valid: " + exactMath.getWarning());
}
return exactMath;
}
Also used :
JumpCondition(cbit.vcell.math.JumpCondition)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
SensVariable(cbit.vcell.solver.ode.SensVariable)
FilamentVariable(cbit.vcell.math.FilamentVariable)
VolVariable(cbit.vcell.math.VolVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
ReservedVariable(cbit.vcell.math.ReservedVariable)
MemVariable(cbit.vcell.math.MemVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
MathDescription(cbit.vcell.math.MathDescription)
VariableHash(cbit.vcell.math.VariableHash)
Constant(cbit.vcell.math.Constant)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
PdeEquation(cbit.vcell.math.PdeEquation)
Function(cbit.vcell.math.Function)
OdeEquation(cbit.vcell.math.OdeEquation)
PdeEquation(cbit.vcell.math.PdeEquation)
Equation(cbit.vcell.math.Equation)
SolutionTemplate(cbit.vcell.numericstest.SolutionTemplate)
ConstructedSolutionTemplate(cbit.vcell.numericstest.ConstructedSolutionTemplate)
OdeEquation(cbit.vcell.math.OdeEquation)
Expression(cbit.vcell.parser.Expression)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)
Example 15 with SubDomain
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class MathTestingUtilities method comparePDEResultsWithExact.
/**
* Insert the method's description here.
* Creation date: (8/20/2003 12:58:10 PM)
*/
public static SimulationComparisonSummary comparePDEResultsWithExact(SimulationSymbolTable simSymbolTable, PDEDataManager dataManager, String type, double absErrorThreshold, double relErrorThreshold) throws DataAccessException, ExpressionException {
java.util.Hashtable<String, DataErrorSummary> tempVarHash = new java.util.Hashtable<String, DataErrorSummary>();
double[] timeArray = dataManager.getDataSetTimes();
Variable[] vars = simSymbolTable.getVariables();
CartesianMesh mesh = dataManager.getMesh();
MathDescription mathDesc = simSymbolTable.getSimulation().getMathDescription();
// Get volumeSubdomains from mathDesc/mesh and store in lookupTable
int numVol = mesh.getSizeX() * mesh.getSizeY() * mesh.getSizeZ();
CompartmentSubDomain[] volSubDomainLookup = new CompartmentSubDomain[numVol];
for (int i = 0; i < numVol; i++) {
int subVolumeIndex = mesh.getSubVolumeFromVolumeIndex(i);
SubVolume subVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(subVolumeIndex);
CompartmentSubDomain compSubDomain = mathDesc.getCompartmentSubDomain(subVolume.getName());
volSubDomainLookup[i] = compSubDomain;
}
// Get membraneSubdomains from mathDesc/mesh and store in lookupTable
int numMem = mesh.getMembraneElements().length;
MembraneSubDomain[] memSubDomainLookup = new MembraneSubDomain[numMem];
for (int i = 0; i < numMem; i++) {
int insideVolIndex = mesh.getMembraneElements()[i].getInsideVolumeIndex();
int outsideVolIndex = mesh.getMembraneElements()[i].getOutsideVolumeIndex();
MembraneSubDomain memSubDomain = mathDesc.getMembraneSubDomain(volSubDomainLookup[insideVolIndex], volSubDomainLookup[outsideVolIndex]);
memSubDomainLookup[i] = memSubDomain;
}
double[] valueArray = new double[4];
SimpleSymbolTable symbolTable = new SimpleSymbolTable(new String[] { "t", "x", "y", "z" });
int tIndex = symbolTable.getEntry("t").getIndex();
int xIndex = symbolTable.getEntry("x").getIndex();
int yIndex = symbolTable.getEntry("y").getIndex();
int zIndex = symbolTable.getEntry("z").getIndex();
SimulationComparisonSummary simComparisonSummary = new SimulationComparisonSummary();
String hashKey = new String("");
long dataLength = 0;
// for each var, do the following :
for (int i = 0; i < vars.length; i++) {
if (vars[i] instanceof VolVariable || vars[i] instanceof MemVariable || vars[i] instanceof FilamentVariable || vars[i] instanceof VolumeRegionVariable || vars[i] instanceof MembraneRegionVariable || vars[i] instanceof FilamentRegionVariable) {
// for each time in timeArray,
for (int j = 0; j < timeArray.length; j++) {
if (type.equals(TestCaseNew.EXACT_STEADY)) {
if (j != (timeArray.length - 1)) {
continue;
}
}
// get data block from varName, data from datablock
SimDataBlock simDataBlock = dataManager.getSimDataBlock(vars[i].getName(), timeArray[j]);
double[] data = simDataBlock.getData();
dataLength = data.length;
SubDomain subDomain = null;
Coordinate subDomainCoord = null;
// for each point in data block ...
for (int k = 0; k < dataLength; k++) {
// Get subdomain from mesh (from the lookupTable), get coordinates (x,y,z) from mesh, evaluate EXACT SOLN at that coord
if (vars[i] instanceof VolVariable) {
subDomain = volSubDomainLookup[k];
subDomainCoord = mesh.getCoordinateFromVolumeIndex(k);
} else if (vars[i] instanceof MemVariable) {
subDomain = memSubDomainLookup[k];
subDomainCoord = mesh.getCoordinateFromMembraneIndex(k);
} else {
throw new RuntimeException("Var " + vars[i].getName() + " not supported yet!");
}
hashKey = vars[i].getName() + ":" + subDomain.getName();
DataErrorSummary tempVar = (DataErrorSummary) tempVarHash.get(hashKey);
if (tempVar == null) {
Expression exp = new Expression(subDomain.getEquation(vars[i]).getExactSolution());
exp.bindExpression(simSymbolTable);
exp = MathUtilities.substituteFunctions(exp, simSymbolTable);
exp = exp.flatten();
exp.bindExpression(symbolTable);
tempVar = new DataErrorSummary(exp);
tempVarHash.put(hashKey, tempVar);
}
// time
valueArray[tIndex] = timeArray[j];
// x
valueArray[xIndex] = subDomainCoord.getX();
// y
valueArray[yIndex] = subDomainCoord.getY();
// z
valueArray[zIndex] = subDomainCoord.getZ();
// EXACT soln at coord subDomainCoord
double value = tempVar.getExactExp().evaluateVector(valueArray);
tempVar.addDataValues(value, data[k], timeArray[j], k, absErrorThreshold, relErrorThreshold);
}
// end for (k)
}
// end for (j)
}
// end - if (var)
}
// end for (i)
Enumeration<String> enumKeys = tempVarHash.keys();
while (enumKeys.hasMoreElements()) {
String key = enumKeys.nextElement();
DataErrorSummary tempVarSummary = tempVarHash.get(key);
simComparisonSummary.addVariableComparisonSummary(new VariableComparisonSummary(key, tempVarSummary.getMinRef(), tempVarSummary.getMaxRef(), tempVarSummary.getMaxAbsoluteError(), tempVarSummary.getMaxRelativeError(), tempVarSummary.getL2Norm(), tempVarSummary.getTimeAtMaxAbsoluteError(), tempVarSummary.getIndexAtMaxAbsoluteError(), tempVarSummary.getTimeAtMaxRelativeError(), tempVarSummary.getIndexAtMaxRelativeError()));
}
return simComparisonSummary;
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
SensVariable(cbit.vcell.solver.ode.SensVariable)
FilamentVariable(cbit.vcell.math.FilamentVariable)
VolVariable(cbit.vcell.math.VolVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
ReservedVariable(cbit.vcell.math.ReservedVariable)
MemVariable(cbit.vcell.math.MemVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
MathDescription(cbit.vcell.math.MathDescription)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
MemVariable(cbit.vcell.math.MemVariable)
SimDataBlock(cbit.vcell.simdata.SimDataBlock)
SubVolume(cbit.vcell.geometry.SubVolume)
AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume)
VolVariable(cbit.vcell.math.VolVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
CartesianMesh(cbit.vcell.solvers.CartesianMesh)
SimpleSymbolTable(cbit.vcell.parser.SimpleSymbolTable)
Coordinate(org.vcell.util.Coordinate)
Expression(cbit.vcell.parser.Expression)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
FilamentVariable(cbit.vcell.math.FilamentVariable)
Aggregations
SubDomain (cbit.vcell.math.SubDomain)50
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)35
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)34
Expression (cbit.vcell.parser.Expression)27
Variable (cbit.vcell.math.Variable)20
MathDescription (cbit.vcell.math.MathDescription)19
Equation (cbit.vcell.math.Equation)17
PdeEquation (cbit.vcell.math.PdeEquation)15
ExpressionException (cbit.vcell.parser.ExpressionException)15
ArrayList (java.util.ArrayList)15
OdeEquation (cbit.vcell.math.OdeEquation)14
Constant (cbit.vcell.math.Constant)13
VolVariable (cbit.vcell.math.VolVariable)13
Vector (java.util.Vector)12
SubVolume (cbit.vcell.geometry.SubVolume)11
Function (cbit.vcell.math.Function)11
MathException (cbit.vcell.math.MathException)11
Simulation (cbit.vcell.solver.Simulation)11
MemVariable (cbit.vcell.math.MemVariable)10
MembraneRegionVariable (cbit.vcell.math.MembraneRegionVariable)9
