Example 1 with RationalNumber
use of ucar.units_vcell.RationalNumber 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()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm != null && parm.getExpression() != null) {
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 2 with RationalNumber
use of ucar.units_vcell.RationalNumber in project vcell by virtualcell.
the class VCUnitTranslator method processCellMLUnit.
// derives the ucar unit equivalent to a single cellml 'unit'
// part of the recursive retrieval of CellML units.
private static Unit processCellMLUnit(Element cellUnit, Namespace sNamespace, Namespace sAttNamespace, VCUnitSystem vcUnitSystem) {
Unit unit = null;
UnitName uName;
String kind = cellUnit.getAttributeValue(CELLMLTags.units, sAttNamespace);
if (kind == null || kind.length() == 0) {
throw new RuntimeException("No kind specified for the CellML unit >>> 2");
}
String exp = cellUnit.getAttributeValue(CELLMLTags.exponent, sAttNamespace);
String scaleStr = cellUnit.getAttributeValue(CELLMLTags.prefix, sAttNamespace);
int scale;
try {
if (scaleStr == null || scaleStr.length() == 0) {
scale = 0;
} else {
scale = Integer.parseInt(scaleStr);
}
} catch (NumberFormatException e) {
// bad practice but no better way.
scale = CELLMLTags.prefixToScale(scaleStr);
}
String multiplier = cellUnit.getAttributeValue(CELLMLTags.mult, sAttNamespace);
String offset = cellUnit.getAttributeValue(CELLMLTags.offset, sAttNamespace);
try {
StandardUnitDB unitDB = StandardUnitDB.instance();
if (kind.equals(CELLMLTags.noDimUnit)) {
// special treatment. ignore all the other attributes. Not used in computing total unit.
unit = DerivedUnitImpl.DIMENSIONLESS;
uName = unit.getUnitName();
} else if (CELLMLTags.isCellBaseUnit(kind)) {
unit = unitDB.get(kind);
uName = UnitName.newUnitName(kind);
} else {
Element owner = (Element) cellUnit.getParent().getParent();
String ownerName = owner.getAttributeValue(CELLMLTags.name, sAttNamespace);
// check if already added.
VCUnitDefinition unitDef = getMatchingCellMLUnitDef(owner, sAttNamespace, kind, vcUnitSystem);
if (unitDef != null) {
unit = unitDef.getUcarUnit();
} else {
// recursive block. assumes model level units are added before component level ones,
// so no need to recurse through those as well.
@SuppressWarnings("unchecked") ArrayList<Element> siblings = new ArrayList<Element>(owner.getChildren(CELLMLTags.UNITS, sNamespace));
Element cellUnitDef = null;
for (Element cuElement : siblings) {
cellUnitDef = cuElement;
if (cellUnitDef.getAttributeValue(CELLMLTags.name, sAttNamespace).equals(kind)) {
break;
} else {
cellUnitDef = null;
}
}
if (cellUnitDef != null) {
unitDef = CellMLToVCUnit(cellUnitDef, sNamespace, sAttNamespace, vcUnitSystem);
} else {
System.err.println("Unit definition: " + kind + " is missing. >>> 3");
return null;
}
}
uName = UnitName.newUnitName(kind);
}
if (exp != null && Integer.parseInt(exp) != 1) {
unit = unit.raiseTo(new RationalNumber(Integer.parseInt(exp)));
}
if (scale != 0) {
unit = new ScaledUnit(Double.parseDouble("1.0e" + scale), unit, uName);
}
if (multiplier != null && Double.parseDouble(multiplier) != 1) {
unit = new ScaledUnit(Double.parseDouble(multiplier), unit, uName);
}
if (offset != null && Double.parseDouble(offset) != 0) {
unit = new OffsetUnit(unit, Double.parseDouble(offset), uName);
}
} catch (UnitException e) {
e.printStackTrace();
throw new cbit.vcell.units.VCUnitException("Unable to set unit value: " + e.getMessage());
}
return unit;
}
Also used :
OffsetUnit(ucar.units_vcell.OffsetUnit)
ScaledUnit(ucar.units_vcell.ScaledUnit)
StandardUnitDB(ucar.units_vcell.StandardUnitDB)
Element(org.jdom.Element)
UnitException(ucar.units_vcell.UnitException)
ArrayList(java.util.ArrayList)
UnitName(ucar.units_vcell.UnitName)
BaseUnit(ucar.units_vcell.BaseUnit)
ScaledUnit(ucar.units_vcell.ScaledUnit)
OffsetUnit(ucar.units_vcell.OffsetUnit)
Unit(ucar.units_vcell.Unit)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
RationalNumber(ucar.units_vcell.RationalNumber)
Example 3 with RationalNumber
use of ucar.units_vcell.RationalNumber in project vcell by virtualcell.
the class VCUnitEvaluator method assignAndVerify.
private void assignAndVerify(VCUnitDefinition nodeUnit, SimpleNode node, UnitsHashMap unitsHashMap) throws ExpressionException {
if (node == null) {
throw new RuntimeException("VCUnitEvaluator.assignAndVerify(): node is null");
}
if (nodeUnit == null) {
throw new RuntimeException("VCUnitEvaluator.assignAndVerify(): nodeUnit is null");
}
//
if (node instanceof ASTAndNode || node instanceof ASTOrNode || node instanceof ASTNotNode) {
for (int i = 0; i < node.jjtGetNumChildren(); i++) {
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(i), unitsHashMap);
}
} else if (node instanceof ASTRelationalNode) {
//
// if either child has known unit, then other must also be known
//
SimpleNode child0 = (SimpleNode) node.jjtGetChild(0);
SimpleNode child1 = (SimpleNode) node.jjtGetChild(1);
VCUnitDefinition unit0 = getUnitDefinition(child0, unitsHashMap);
VCUnitDefinition unit1 = getUnitDefinition(child1, unitsHashMap);
if (unit0.isTBD() && unit1.isTBD()) {
assignAndVerify(unitSystem.getInstance_TBD(), child0, unitsHashMap);
assignAndVerify(unitSystem.getInstance_TBD(), child1, unitsHashMap);
return;
} else if (unit0.isTBD()) {
assignAndVerify(unit1, child0, unitsHashMap);
return;
} else if (unit1.isTBD()) {
assignAndVerify(unit0, child1, unitsHashMap);
return;
} else if (!unit0.isEquivalent(unit1)) {
throw new RuntimeException("operands of a relational operator are not same [" + unit0.getSymbol() + "] and [" + unit1.getSymbol() + "]");
} else {
assignAndVerify(unit0, child0, unitsHashMap);
assignAndVerify(unit0, child1, unitsHashMap);
}
} else if (node instanceof ASTAddNode || node instanceof ASTMinusTermNode) {
for (int i = 0; i < node.jjtGetNumChildren(); i++) {
assignAndVerify(nodeUnit, (SimpleNode) node.jjtGetChild(i), unitsHashMap);
}
} else if (node instanceof ASTMultNode) {
if (node.jjtGetNumChildren() == 1) {
assignAndVerify(nodeUnit, (SimpleNode) node.jjtGetChild(0), unitsHashMap);
} else {
//
// product of children should be that of parent (tough one)
// ... can only resolve if only one unknown
//
SimpleNode unknownChildNode = null;
int unknownChildCount = 0;
int constantChildCount = 0;
VCUnitDefinition accumUnit = nodeUnit;
for (int i = 0; i < node.jjtGetNumChildren(); i++) {
SimpleNode child = (SimpleNode) node.jjtGetChild(i);
VCUnitDefinition childUnit = getUnitDefinition(child, unitsHashMap);
assignAndVerify(childUnit, child, unitsHashMap);
if (childUnit.isTBD()) {
//
// ignore those that evaluate to a constant
//
boolean bConstant = false;
try {
child.evaluateConstant();
bConstant = true;
} catch (ExpressionException e) {
}
if (!bConstant) {
unknownChildNode = child;
unknownChildCount++;
} else {
constantChildCount++;
}
} else {
if (!accumUnit.isTBD()) {
accumUnit = accumUnit.divideBy(childUnit);
}
}
}
if (unknownChildCount == 0 && constantChildCount == 0) {
// should cancel to dimensionless, but if constant children (numbers...) are used in a product, then can assume appropriate units.
if (!accumUnit.isEquivalent(unitSystem.getInstance_DIMENSIONLESS())) {
// accumUnit.compareEquals(unitSystem.getInstance_DIMENSIONLESS());
throw new RuntimeException("expression '" + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + "' missing factor of '" + accumUnit.getSymbol() + "'");
}
} else if (unknownChildCount == 1) {
//
if (!accumUnit.isTBD()) {
assignAndVerify(accumUnit, unknownChildNode, unitsHashMap);
}
}
}
} else if (node instanceof ASTInvertTermNode) {
SimpleNode child = (SimpleNode) node.jjtGetChild(0);
if (nodeUnit.isTBD()) {
assignAndVerify(nodeUnit, child, unitsHashMap);
} else {
assignAndVerify(nodeUnit.getInverse(), child, unitsHashMap);
}
// } else if (node instanceof DerivativeNode) {
// unit = getUnitDefinition((SimpleNode)node.jjtGetChild(1),unitsHashMap);
// return getUnitDefinition((SimpleNode)node.jjtGetChild(0),unitsHashMap).divideBy(unit);
// } else if (node instanceof ASTLaplacianNode) {
// unit = ReservedSymbol.X.getUnitDefinition();
// return getUnitDefinition((SimpleNode)node.jjtGetChild(0),unitsHashMap).divideBy(unit).divideBy(unit);
} else if (node instanceof ASTFloatNode) {
// return TBD instead of dimensionless.
// do nothing
} else if (node instanceof ASTFuncNode) {
String functionName = ((ASTFuncNode) node).getName();
if (functionName.equalsIgnoreCase("sqrt")) {
// ?
assignAndVerify(nodeUnit.raiseTo(new RationalNumber(2)), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
} else if (functionName.equalsIgnoreCase("exp")) {
// ?
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
} else if (functionName.equalsIgnoreCase("pow")) {
// later....
// exponent should always be dimensionless
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(1), unitsHashMap);
//
try {
double exponentValue = ((SimpleNode) node.jjtGetChild(1)).evaluateConstant();
if (!nodeUnit.isTBD()) {
RationalNumber rn = RationalNumber.getApproximateFraction(exponentValue);
// exponent should always be dimensionless
assignAndVerify(nodeUnit.raiseTo(rn.inverse()), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
}
} catch (ExpressionException e) {
//
// a^b where b not constant, a must be non-dimensional
//
// exponent should always be dimensionless
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
}
} else if (functionName.equalsIgnoreCase("abs") || functionName.equalsIgnoreCase("min") || functionName.equalsIgnoreCase("max")) {
//
for (int i = 0; i < node.jjtGetNumChildren(); i++) {
assignAndVerify(nodeUnit, (SimpleNode) node.jjtGetChild(i), unitsHashMap);
}
} else {
if (!nodeUnit.isEquivalent(unitSystem.getInstance_DIMENSIONLESS())) {
throw new RuntimeException("function '" + functionName + "' should be dimensionless, assignAndVerify trying to impose '" + nodeUnit.getSymbol());
}
//
// child should be dimensionless
//
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
}
} else if (node instanceof ASTPowerNode) {
// exponent should always be dimensionless
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(1), unitsHashMap);
//
try {
double exponentValue = ((SimpleNode) node.jjtGetChild(1)).evaluateConstant();
RationalNumber rn = RationalNumber.getApproximateFraction(exponentValue);
// exponent should always be dimensionless
assignAndVerify(nodeUnit.raiseTo(rn.inverse()), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
} catch (ExpressionException e) {
//
// a^b where b not constant, a must be non-dimensional
//
// exponent should always be dimensionless
assignAndVerify(unitSystem.getInstance_DIMENSIONLESS(), (SimpleNode) node.jjtGetChild(0), unitsHashMap);
}
} else if (node instanceof ASTIdNode) {
if (!nodeUnit.isTBD()) {
unitsHashMap.put(((ASTIdNode) node).symbolTableEntry, nodeUnit);
}
} else {
throw new ExpressionException("node type " + node.getClass().toString() + " not supported yet");
}
}
Also used :
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
RationalNumber(ucar.units_vcell.RationalNumber)
Example 4 with RationalNumber
use of ucar.units_vcell.RationalNumber in project vcell by virtualcell.
the class VCUnitTranslator method expandUcarCellML.
// time mult not included.
// first call: (1.0, 0, unit, new ArrayList(), transType)
// Problem: units override each other's mult/multiplier before getting to the level of derived unit.
// To avoid that, add a dimensionless unit for cases where
protected static ArrayList<Element> expandUcarCellML(double mult, double offset, Unit unit, ArrayList<Element> transUnits) {
if (unit instanceof UnitImpl) {
UnitImpl unitImpl = (UnitImpl) unit;
if (unitImpl instanceof DerivedUnitImpl) {
DerivedUnitImpl baseUnit = (DerivedUnitImpl) unitImpl;
// String symbol = baseUnit.getID();
// String name = baseUnit.getName();
Factor[] factors = baseUnit.getDimension().getFactors();
for (int i = 0; i < factors.length; i++) {
RationalNumber exponent = factors[i].getExponent();
// String baseSymbol = factors[i].getBase().getID();
String baseName = ((BaseUnit) factors[i].getBase()).getName();
// System.out.println(baseName + " " + exponent);
if (factors.length > 1) {
if (i == 0) {
transUnits.add(getCellMLTransUnit(new RationalNumber(1), mult, offset, DerivedUnitImpl.DIMENSIONLESS.getName()));
}
transUnits.add(getCellMLTransUnit(exponent, 1.0, 0.0, baseName));
} else {
transUnits.add(getCellMLTransUnit(exponent, mult, offset, baseName));
}
}
return transUnits;
} else if (unitImpl instanceof OffsetUnit) {
OffsetUnit offsetUnit = (OffsetUnit) unitImpl;
offset += offsetUnit.getOffset();
// offset = offsetUnit.getOffset();
if (offsetUnit.getUnit() != offsetUnit.getDerivedUnit()) {
return expandUcarCellML(mult, offset, offsetUnit.getUnit(), transUnits);
}
} else if (unitImpl instanceof ScaledUnit) {
ScaledUnit multdUnit = (ScaledUnit) unitImpl;
mult *= multdUnit.getScale();
System.out.println("mult: " + mult);
if (multdUnit.getUnit() != multdUnit.getDerivedUnit()) {
return expandUcarCellML(mult, offset, multdUnit.getUnit(), transUnits);
}
}
if (unitImpl.getDerivedUnit() != unit) {
// i.e. we have not reached the base unit, yet
return expandUcarCellML(mult, offset, unitImpl.getDerivedUnit(), transUnits);
}
} else {
System.err.println("Unable to process unit translation for CellML: " + " " + unit.getSymbol());
}
return null;
}
Also used :
DerivedUnitImpl(ucar.units_vcell.DerivedUnitImpl)
OffsetUnit(ucar.units_vcell.OffsetUnit)
ScaledUnit(ucar.units_vcell.ScaledUnit)
UnitImpl(ucar.units_vcell.UnitImpl)
DerivedUnitImpl(ucar.units_vcell.DerivedUnitImpl)
Factor(ucar.units_vcell.Factor)
RationalNumber(ucar.units_vcell.RationalNumber)
BaseUnit(ucar.units_vcell.BaseUnit)
Example 5 with RationalNumber
use of ucar.units_vcell.RationalNumber in project vcell by virtualcell.
the class ParticleMathMapping method combineHybrid.
private void combineHybrid() throws MappingException, ExpressionException, MatrixException, MathException, ModelException {
ArrayList<SpeciesContext> continuousSpecies = new ArrayList<SpeciesContext>();
ArrayList<ParticleVariable> continuousSpeciesParticleVars = new ArrayList<ParticleVariable>();
ArrayList<SpeciesContext> stochSpecies = new ArrayList<SpeciesContext>();
//
// categorize speciesContexts as continuous and stochastic
//
SpeciesContextSpec[] scsArray = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
continuousSpecies = new ArrayList<SpeciesContext>();
stochSpecies = new ArrayList<SpeciesContext>();
for (SpeciesContextSpec speciesContextSpec : scsArray) {
if (!getSimulationContext().isStoch() || speciesContextSpec.isForceContinuous()) {
continuousSpecies.add(speciesContextSpec.getSpeciesContext());
Variable variable = getMathSymbolMapping().getVariable(speciesContextSpec.getSpeciesContext());
if (variable instanceof ParticleVariable) {
continuousSpeciesParticleVars.add((ParticleVariable) variable);
}
} else {
stochSpecies.add(speciesContextSpec.getSpeciesContext());
}
}
if (continuousSpecies.isEmpty()) {
return;
}
//
// create continuous mathDescription ... add stochastic variables and processes to the continuous Math and use this.
//
DiffEquMathMapping mathMapping = new DiffEquMathMapping(getSimulationContext(), callback, networkGenerationRequirements);
mathMapping.refresh(null);
MathDescription contMathDesc = mathMapping.getMathDescription();
//
// get list of all continuous variables
//
HashMap<String, Variable> allContinuousVars = new HashMap<String, Variable>();
Enumeration<Variable> enumVar = contMathDesc.getVariables();
while (enumVar.hasMoreElements()) {
Variable var = enumVar.nextElement();
allContinuousVars.put(var.getName(), var);
}
//
// replace those continuous variables and equations for stochastic speciesContexts
// with the particleVariables and particleProperties
// (ParticleJumpProcesses removed later)
//
ModelUnitSystem unitSystem = getSimulationContext().getModel().getUnitSystem();
for (SpeciesContext stochSpeciesContext : stochSpecies) {
Variable contVar = mathMapping.getMathSymbolMapping().getVariable(stochSpeciesContext);
Variable stochVar = getMathSymbolMapping().getVariable(stochSpeciesContext);
allContinuousVars.put(stochVar.getName(), stochVar);
//
// replace continuous "concentration" VolVariable/MemVariable for this particle with a Function for concentration
//
allContinuousVars.remove(contVar);
VCUnitDefinition sizeUnit = unitSystem.getLengthUnit().raiseTo(new RationalNumber(stochSpeciesContext.getStructure().getDimension()));
VCUnitDefinition stochasticDensityUnit = unitSystem.getStochasticSubstanceUnit().divideBy(sizeUnit);
VCUnitDefinition continuousDensityUnit = unitSystem.getConcentrationUnit(stochSpeciesContext.getStructure());
if (stochasticDensityUnit.isEquivalent(continuousDensityUnit)) {
allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), new Expression(stochVar, getNameScope()), contVar.getDomain()));
} else {
Expression conversionFactorExp = getUnitFactor(continuousDensityUnit.divideBy(stochasticDensityUnit));
allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), Expression.mult(new Expression(stochVar, getNameScope()), conversionFactorExp), contVar.getDomain()));
}
//
// remove continuous equation
//
Enumeration<SubDomain> contSubDomains = contMathDesc.getSubDomains();
while (contSubDomains.hasMoreElements()) {
SubDomain contSubDomain = contSubDomains.nextElement();
contSubDomain.removeEquation(contVar);
if (contSubDomain instanceof MembraneSubDomain) {
((MembraneSubDomain) contSubDomain).removeJumpCondition(contVar);
}
}
//
// remove all continuous variables for speciesContextSpec parameters (e.g. initial conditions, diffusion rates, boundary conditions, velocities)
//
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(stochSpeciesContext);
Parameter[] scsParameters = scs.getParameters();
for (Parameter parameter : scsParameters) {
Variable continuousScsParmVariable = mathMapping.getMathSymbolMapping().getVariable(parameter);
allContinuousVars.remove(continuousScsParmVariable);
}
//
// copy ParticleJumpProcess and ParticleProperties to the continuous math
//
SubDomain contSubDomain = contMathDesc.getSubDomain(contVar.getDomain().getName());
SubDomain stochSubDomain = mathDesc.getSubDomain(stochVar.getDomain().getName());
ParticleProperties particleProperties = stochSubDomain.getParticleProperties(stochVar);
contSubDomain.addParticleProperties(particleProperties);
}
//
// add all ParticleJumpProcesses to the continuous model
//
Enumeration<SubDomain> enumStochSubdomains = mathDesc.getSubDomains();
while (enumStochSubdomains.hasMoreElements()) {
SubDomain stochSubdomain = enumStochSubdomains.nextElement();
SubDomain contSubdomain = contMathDesc.getSubDomain(stochSubdomain.getName());
for (ParticleJumpProcess particleJumpProcess : stochSubdomain.getParticleJumpProcesses()) {
//
// modify "selection list" (particleVariables), probability rate, and actions if referenced particleVariable is to be "forced continuous"
//
ParticleVariable[] selectedParticles = particleJumpProcess.getParticleVariables();
for (ParticleVariable particleVariable : selectedParticles) {
if (continuousSpeciesParticleVars.contains(particleVariable)) {
particleJumpProcess.remove(particleVariable);
JumpProcessRateDefinition jumpProcessRateDefinition = particleJumpProcess.getParticleRateDefinition();
if (jumpProcessRateDefinition instanceof MacroscopicRateConstant) {
MacroscopicRateConstant macroscopicRateConstant = (MacroscopicRateConstant) jumpProcessRateDefinition;
macroscopicRateConstant.setExpression(Expression.mult(macroscopicRateConstant.getExpression(), new Expression(particleVariable, null)));
} else if (jumpProcessRateDefinition instanceof InteractionRadius) {
throw new MappingException("cannot adjust interaction radius for reaction process " + particleJumpProcess.getName() + ", particle " + particleVariable.getName() + " is continuous");
} else {
throw new MappingException("rate definition type " + jumpProcessRateDefinition.getClass().getSimpleName() + " not yet implemented for hybrid PDE/Particle math generation");
}
}
Iterator<Action> iterAction = particleJumpProcess.getActions().iterator();
while (iterAction.hasNext()) {
Action action = iterAction.next();
if (continuousSpeciesParticleVars.contains(action.getVar())) {
iterAction.remove();
}
}
}
if (!particleJumpProcess.getActions().isEmpty()) {
contSubdomain.addParticleJumpProcess(particleJumpProcess);
}
}
}
//
for (MathMappingParameter mathMappingParameter : fieldMathMappingParameters) {
if (mathMappingParameter instanceof UnitFactorParameter) {
String name = mathMappingParameter.getName();
if (!allContinuousVars.containsKey(name)) {
allContinuousVars.put(name, newFunctionOrConstant(name, mathMappingParameter.getExpression(), null));
}
}
}
//
// add constants and functions from the particle math that aren't already defined in the continuous math
//
Enumeration<Variable> enumVars = mathDesc.getVariables();
while (enumVars.hasMoreElements()) {
Variable var = enumVars.nextElement();
if (var instanceof Constant || var instanceof Function) {
String name = var.getName();
if (!allContinuousVars.containsKey(name)) {
allContinuousVars.put(name, var);
}
}
}
contMathDesc.setAllVariables(allContinuousVars.values().toArray(new Variable[0]));
mathDesc = contMathDesc;
//
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()) {
System.out.println(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used :
GeometryClass(cbit.vcell.geometry.GeometryClass)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
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)
MathDescription(cbit.vcell.math.MathDescription)
HashMap(java.util.HashMap)
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)
ArrayList(java.util.ArrayList)
SpeciesContext(cbit.vcell.model.SpeciesContext)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Function(cbit.vcell.math.Function)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
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)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Aggregations
RationalNumber (ucar.units_vcell.RationalNumber)6
VCUnitDefinition (cbit.vcell.units.VCUnitDefinition)5
ArrayList (java.util.ArrayList)4
GeometryClass (cbit.vcell.geometry.GeometryClass)2
SpeciesContextSpecParameter (cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)2
SpeciesContextSpecProxyParameter (cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)2
StructureMappingParameter (cbit.vcell.mapping.StructureMapping.StructureMappingParameter)2
Action (cbit.vcell.math.Action)2
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)2
Constant (cbit.vcell.math.Constant)2
InteractionRadius (cbit.vcell.math.InteractionRadius)2
JumpProcessRateDefinition (cbit.vcell.math.JumpProcessRateDefinition)2
MacroscopicRateConstant (cbit.vcell.math.MacroscopicRateConstant)2
MathDescription (cbit.vcell.math.MathDescription)2
MembraneParticleVariable (cbit.vcell.math.MembraneParticleVariable)2
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)2
ParticleJumpProcess (cbit.vcell.math.ParticleJumpProcess)2
ParticleProperties (cbit.vcell.math.ParticleProperties)2
ParticleVariable (cbit.vcell.math.ParticleVariable)2
SubDomain (cbit.vcell.math.SubDomain)2
