use of cbit.vcell.math.JumpProcessRateDefinition 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 ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
use of cbit.vcell.math.JumpProcessRateDefinition in project vcell by virtualcell.
the class XmlReader method getParticleJumpProcess.
private ParticleJumpProcess getParticleJumpProcess(Element param, MathDescription md) throws XmlParseException {
// name
String name = unMangle(param.getAttributeValue(XMLTags.NameAttrTag));
ProcessSymmetryFactor processSymmetryFactor = null;
Attribute symmetryFactorAttr = param.getAttribute(XMLTags.ProcessSymmetryFactorAttrTag);
if (symmetryFactorAttr != null) {
processSymmetryFactor = new ProcessSymmetryFactor(Double.parseDouble(symmetryFactorAttr.getValue()));
}
// selected particle
List<ParticleVariable> varList = new ArrayList<ParticleVariable>();
Iterator<Element> iterator = param.getChildren(XMLTags.SelectedParticleTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
String varname = unMangle(tempelement.getAttributeValue(XMLTags.NameAttrTag));
Variable var = md.getVariable(varname);
if (!(var instanceof ParticleVariable)) {
throw new XmlParseException("Not a ParticleVariable in ParticleJumpProcess.");
}
varList.add((ParticleVariable) var);
}
// probability rate
JumpProcessRateDefinition jprd = null;
// for old models
Element pb = param.getChild(XMLTags.ParticleProbabilityRateTag, vcNamespace);
if (pb != null) {
Expression exp = unMangleExpression(pb.getText());
jprd = new MacroscopicRateConstant(exp);
} else // for new models
{
pb = param.getChild(XMLTags.MacroscopicRateConstantTag, vcNamespace);
if (// jump process rate defined by macroscopic rate constant
pb != null) {
Expression exp = unMangleExpression(pb.getText());
jprd = new MacroscopicRateConstant(exp);
} else // jump process rate defined by binding radius
{
pb = param.getChild(XMLTags.InteractionRadiusTag, vcNamespace);
if (pb != null) {
Expression exp = unMangleExpression(pb.getText());
jprd = new InteractionRadius(exp);
}
}
}
// add actions
List<Action> actionList = new ArrayList<Action>();
iterator = param.getChildren(XMLTags.ActionTag, vcNamespace).iterator();
while (iterator.hasNext()) {
Element tempelement = (Element) iterator.next();
try {
actionList.add(getAction(tempelement, md));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
} catch (ExpressionException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
ParticleJumpProcess jump = new ParticleJumpProcess(name, varList, jprd, actionList, processSymmetryFactor);
return jump;
}
use of cbit.vcell.math.JumpProcessRateDefinition in project vcell by virtualcell.
the class NFsimXMLWriter method writeNFsimXML.
public static Element writeNFsimXML(SimulationTask origSimTask, long randomSeed, NFsimSimulationOptions nfsimSimulationOptions, boolean bUseLocationMarks) throws SolverException {
try {
System.out.println("VCML ORIGINAL .... START\n" + origSimTask.getSimulation().getMathDescription().getVCML_database() + "\nVCML ORIGINAL .... END\n====================\n");
} catch (MathException e1) {
// TODO Auto-generated catch block
e1.printStackTrace();
}
SimulationTask clonedSimTask = null;
try {
clonedSimTask = (SimulationTask) BeanUtils.cloneSerializable(origSimTask);
} catch (Exception eee) {
throw new SolverException("failed to clone mathDescription while preparing NFSim input: " + eee.getMessage(), eee);
}
MathDescription clonedMathDesc = clonedSimTask.getSimulation().getMathDescription();
if (bUseLocationMarks) {
try {
//
// get list of Compartment Names (stored in locations).
//
ArrayList<String> locations = new ArrayList<String>();
Enumeration<Variable> varEnum = clonedMathDesc.getVariables();
ArrayList<VolumeParticleSpeciesPattern> volumeParticleSpeciesPatterns = new ArrayList<VolumeParticleSpeciesPattern>();
while (varEnum.hasMoreElements()) {
Variable var = varEnum.nextElement();
if (var instanceof VolumeParticleSpeciesPattern) {
VolumeParticleSpeciesPattern speciesPattern = (VolumeParticleSpeciesPattern) var;
if (!locations.contains(speciesPattern.getLocationName())) {
locations.add(speciesPattern.getLocationName());
}
volumeParticleSpeciesPatterns.add(speciesPattern);
}
}
//
for (ParticleMolecularType particleMolecularType : clonedMathDesc.getParticleMolecularTypes()) {
String pmcLocationName = RbmUtils.SiteStruct;
String pmcLocationId = particleMolecularType.getName() + "_" + RbmUtils.SiteStruct;
ParticleMolecularComponent locationComponent = new ParticleMolecularComponent(pmcLocationId, pmcLocationName);
for (String location : locations) {
locationComponent.addComponentStateDefinition(new ParticleComponentStateDefinition(location));
}
particleMolecularType.insertMolecularComponent(0, locationComponent);
String pmcMarkName = RbmUtils.SiteProduct;
String pmcMarkId = particleMolecularType.getName() + "_" + RbmUtils.SiteProduct;
ParticleMolecularComponent markComponent = new ParticleMolecularComponent(pmcMarkId, pmcMarkName);
markComponent.addComponentStateDefinition(new ParticleComponentStateDefinition("0"));
markComponent.addComponentStateDefinition(new ParticleComponentStateDefinition("1"));
particleMolecularType.insertMolecularComponent(1, markComponent);
}
//
for (VolumeParticleSpeciesPattern speciesPattern : volumeParticleSpeciesPatterns) {
for (ParticleMolecularTypePattern molTypePattern : speciesPattern.getParticleMolecularTypePatterns()) {
//
// add location component to pattern ... state=<location>
//
{
final ParticleMolecularComponent locationComponentDefinition = molTypePattern.getMolecularType().getComponentList().get(0);
ParticleMolecularComponentPattern locationPattern = new ParticleMolecularComponentPattern(locationComponentDefinition);
ParticleComponentStateDefinition locationStateDefinition = null;
for (ParticleComponentStateDefinition stateDef : locationComponentDefinition.getComponentStateDefinitions()) {
if (stateDef.getName().equals(speciesPattern.getLocationName())) {
locationStateDefinition = stateDef;
}
}
ParticleComponentStatePattern locationStatePattern = new ParticleComponentStatePattern(locationStateDefinition);
locationPattern.setComponentStatePattern(locationStatePattern);
locationPattern.setBondType(ParticleBondType.None);
locationPattern.setBondId(-1);
molTypePattern.insertMolecularComponentPattern(0, locationPattern);
}
//
// add mark component to pattern ... state="0" (for observables and reactants ... later we will clone and use "1" for products).
{
final ParticleMolecularComponent markComponentDefinition = molTypePattern.getMolecularType().getComponentList().get(1);
ParticleMolecularComponentPattern markPattern = new ParticleMolecularComponentPattern(markComponentDefinition);
final int clearStateIndex = 0;
final int setStateIndex = 1;
ParticleComponentStateDefinition markStateClearedDefinition = markComponentDefinition.getComponentStateDefinitions().get(clearStateIndex);
ParticleComponentStatePattern markStatePattern = new ParticleComponentStatePattern(markStateClearedDefinition);
markPattern.setComponentStatePattern(markStatePattern);
markPattern.setBondType(ParticleBondType.None);
markPattern.setBondId(-1);
molTypePattern.insertMolecularComponentPattern(1, markPattern);
}
}
}
//
// when processing ParticleJumpProcesses, we add a new "product" species pattern (by cloning the original speciesPattern)
// and setting the mark site to "1", change name to name+"_PRODUCT", and add to math model if it doesn't already exist.
//
// cloned the "standard" reactant/observable speciesPattern, set the mark for all molecules, and add to mathDesc.
//
CompartmentSubDomain subDomain = (CompartmentSubDomain) clonedMathDesc.getSubDomains().nextElement();
for (ParticleJumpProcess particleJumpProcess : subDomain.getParticleJumpProcesses()) {
for (Action action : particleJumpProcess.getActions()) {
if (action.getOperation().equals(Action.ACTION_CREATE)) {
VolumeParticleSpeciesPattern volumeParticleSpeciesPattern = (VolumeParticleSpeciesPattern) action.getVar();
String newSpeciesPatternName = volumeParticleSpeciesPattern.getName() + "_" + particleJumpProcess.getName();
VolumeParticleSpeciesPattern productPattern = new VolumeParticleSpeciesPattern(volumeParticleSpeciesPattern, newSpeciesPatternName);
// VolumeParticleSpeciesPattern productPattern = (VolumeParticleSpeciesPattern) BeanUtils.cloneSerializable(volumeParticleSpeciesPattern);
for (ParticleMolecularTypePattern productMolTypePattern : productPattern.getParticleMolecularTypePatterns()) {
ParticleComponentStateDefinition markSet = productMolTypePattern.getMolecularType().getComponentList().get(1).getComponentStateDefinitions().get(1);
productMolTypePattern.getMolecularComponentPatternList().get(1).setComponentStatePattern(new ParticleComponentStatePattern(markSet));
}
System.out.println(productPattern.getName());
if (clonedMathDesc.getVariable(productPattern.getName()) == null) {
clonedMathDesc.addVariable(productPattern);
}
action.setVar(productPattern);
}
}
}
try {
System.out.println("===============================\n ----------- VCML HACKED .... START\n" + clonedMathDesc.getVCML_database() + "\nVCML HACKED .... END\n====================\n");
} catch (MathException e1) {
// TODO Auto-generated catch block
e1.printStackTrace();
}
} catch (Exception e) {
throw new SolverException("failed to apply location mark transformation: " + e.getMessage(), e);
}
}
Element sbmlElement = new Element("sbml");
Element modelElement = new Element("model");
modelElement.setAttribute("id", "nameless");
SimulationSymbolTable simulationSymbolTable = new SimulationSymbolTable(clonedSimTask.getSimulation(), clonedSimTask.getSimulationJob().getJobIndex());
Element listOfParametersElement = getListOfParameters(clonedMathDesc, simulationSymbolTable);
Element listOfMoleculeTypesElement = getListOfMoleculeTypes(clonedMathDesc);
Element listOfSpeciesElement = getListOfSpecies(clonedMathDesc, simulationSymbolTable);
CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) clonedMathDesc.getSubDomains().nextElement();
Element listOfReactionRules = new Element("ListOfReactionRules");
for (int reactionRuleIndex = 0; reactionRuleIndex < compartmentSubDomain.getParticleJumpProcesses().size(); reactionRuleIndex++) {
ParticleJumpProcess particleJumpProcess = compartmentSubDomain.getParticleJumpProcesses().get(reactionRuleIndex);
MathRuleFactory mathRuleFactory = new MathRuleFactory();
MathRuleEntry rule = mathRuleFactory.createRuleEntry(particleJumpProcess, reactionRuleIndex);
RuleAnalysisReport report = RuleAnalysis.analyze(rule, true);
// remember, we have to add RateLaw
Element reactionRuleElement = RuleAnalysis.getNFSimXML(rule, report);
// ArrayList<MolecularTypeOfReactionParticipant> currentReactantElementsOfReaction = new ArrayList<MolecularTypeOfReactionParticipant>();
// ArrayList<ComponentOfMolecularTypeOfReactionParticipant> currentComponentOfReactantElementsOfReaction = new ArrayList<ComponentOfMolecularTypeOfReactionParticipant>();
// ArrayList<MolecularTypeOfReactionParticipant> currentProductElementsOfReaction = new ArrayList<MolecularTypeOfReactionParticipant>();
// ArrayList<ComponentOfMolecularTypeOfReactionParticipant> currentComponentOfProductElementsOfReaction = new ArrayList<ComponentOfMolecularTypeOfReactionParticipant>();
// currentMappingOfReactionParticipants.clear();
// reactionProductBondSites.clear();
// reactionReactantBondSites.clear();
//
// Element reactionRuleElement = new Element("ReactionRule");
// String reactionRuleID = "RR" + (reactionRuleIndex + 1);
// reactionRuleElement.setAttribute("id",reactionRuleID);
// reactionRuleElement.setAttribute("name",particleJumpProcess.getName());
// reactionRuleElement.setAttribute("symmetry_factor","1");
// reactionRule.resolveBonds();
//
// ArrayList<VolumeParticleSpeciesPattern> selectedPatterns = new ArrayList<VolumeParticleSpeciesPattern>();
// for (ParticleVariable particleVariable : particleJumpProcess.getParticleVariables()){
// if (!(particleVariable instanceof VolumeParticleSpeciesPattern)){
// throw new SolverException("expecting only "+VolumeParticleSpeciesPattern.class.getSimpleName()+"s for "+ParticleJumpProcess.class.getSimpleName()+" "+particleJumpProcess.getName());
// }
// selectedPatterns.add((VolumeParticleSpeciesPattern) particleVariable);
// }
// ArrayList<VolumeParticleSpeciesPattern> createdPatterns = new ArrayList<VolumeParticleSpeciesPattern>();
// HashSet<VolumeParticleSpeciesPattern> destroyedPatterns = new HashSet<VolumeParticleSpeciesPattern>();
// for (Action action : particleJumpProcess.getActions()){
// if (!(action.getVar() instanceof VolumeParticleSpeciesPattern)){
// throw new SolverException("expecting only "+VolumeParticleSpeciesPattern.class.getSimpleName()+"s for "+ParticleJumpProcess.class.getSimpleName()+" "+particleJumpProcess.getName());
// }
// if (action.getOperation().equals(Action.ACTION_CREATE)){
// createdPatterns.add((VolumeParticleSpeciesPattern) action.getVar());
// }else if (action.getOperation().equals(Action.ACTION_DESTROY)){
// destroyedPatterns.add((VolumeParticleSpeciesPattern) action.getVar());
// }else{
// throw new RuntimeException("unexpected action operation "+action.getOperation()+" for jump process "+particleJumpProcess.getName());
// }
// }
//
// Element listOfReactantPatternsElement = new Element("ListOfReactantPatterns");
// for(int reactantPatternIndex=0; reactantPatternIndex < selectedPatterns.size(); reactantPatternIndex++) {
// VolumeParticleSpeciesPattern reactantSpeciesPattern = selectedPatterns.get(reactantPatternIndex);
// String reactantPatternID = "RP" + (reactantPatternIndex + 1);
// patternReactantBondSites.clear();
// Element reactantPatternElement = getReactionParticipantPattern1(reactionRuleID, reactantPatternID, reactantSpeciesPattern,
// currentReactantElementsOfReaction, currentComponentOfReactantElementsOfReaction, "ReactantPattern");
// listOfReactantPatternsElement.addContent(reactantPatternElement);
// reactionReactantBondSites.addAll(patternReactantBondSites);
// }
// reactionRuleElement.addContent(listOfReactantPatternsElement);
//
// Element listOfProductPatternsElement = new Element("ListOfProductPatterns");
// ArrayList<VolumeParticleSpeciesPattern> productSpeciesPatterns = new ArrayList<VolumeParticleSpeciesPattern>(selectedPatterns);
// productSpeciesPatterns.removeAll(destroyedPatterns);
// productSpeciesPatterns.addAll(createdPatterns);
// // for products, add all "created" species from Actions and all "particles" that are selected but not destroyed
// for(int productPatternIndex=0; productPatternIndex < productSpeciesPatterns.size(); productPatternIndex++) {
// VolumeParticleSpeciesPattern productSpeciesPattern = productSpeciesPatterns.get(productPatternIndex);
// String productPatternID = "PP" + (productPatternIndex + 1);
// patternProductBondSites.clear();
// Element productPatternElement = getReactionParticipantPattern1(reactionRuleID, productPatternID, productSpeciesPattern,
// currentProductElementsOfReaction, currentComponentOfProductElementsOfReaction, "ProductPattern");
// listOfProductPatternsElement.addContent(productPatternElement);
// reactionProductBondSites.addAll(patternProductBondSites);
// }
// reactionRuleElement.addContent(listOfProductPatternsElement);
// <RateLaw id="RR1_RateLaw" type="Ele" totalrate="0">
// <ListOfRateConstants>
// <RateConstant value="kon"/>
// </ListOfRateConstants>
// </RateLaw>
Element rateLawElement = new Element("RateLaw");
rateLawElement.setAttribute("id", RuleAnalysis.getID(rule));
String rateConstantValue = null;
JumpProcessRateDefinition particleProbabilityRate = particleJumpProcess.getParticleRateDefinition();
if (particleProbabilityRate.getExpressions().length > 0) {
JumpProcessRateDefinition particleRateDefinition = particleJumpProcess.getParticleRateDefinition();
Expression expression = null;
if (particleRateDefinition instanceof MacroscopicRateConstant) {
expression = ((MacroscopicRateConstant) particleProbabilityRate).getExpression();
} else {
throw new SolverException("ParticleRateDefinition type " + particleRateDefinition.getClass().getSimpleName() + " not supported");
}
rateConstantValue = expression.infixBng();
// all rates constants are being flattened and given reserved names
Expression substitutedValExpr = null;
try {
substitutedValExpr = simulationSymbolTable.substituteFunctions(expression);
} catch (MathException | ExpressionException e) {
e.printStackTrace(System.out);
throw new SolverException("ParticleJumpProcess " + particleJumpProcess.getName() + " substitution failed : exp = \"" + expression.infix() + "\": " + e.getMessage());
}
Double value = null;
try {
value = substitutedValExpr.evaluateConstant();
Element parameterElement = new Element("Parameter");
String id = "K_reserved_" + reactionRuleIndex;
parameterElement.setAttribute("id", id);
if (value != null) {
parameterElement.setAttribute("type", "Constant");
parameterElement.setAttribute("value", value.toString());
parameterElement.addContent(new Comment(rateConstantValue));
rateConstantValue = id;
listOfParametersElement.addContent(parameterElement);
}
} catch (ExpressionException e) {
System.out.println("ParticleJumpProcess " + particleJumpProcess.getName() + " = " + substitutedValExpr.infix() + " does not have a constant value");
}
}
if (isFunction(rateConstantValue, clonedMathDesc, simulationSymbolTable)) {
rateLawElement.setAttribute("type", "Function");
rateLawElement.setAttribute("totalrate", "0");
rateLawElement.setAttribute("name", rateConstantValue);
} else {
rateLawElement.setAttribute("type", "Ele");
rateLawElement.setAttribute("totalrate", "0");
Element listOfRateConstantsElement = new Element("ListOfRateConstants");
Element rateConstantElement = new Element("RateConstant");
// System.out.println(" --- " + particleJumpProcess.getParticleRateDefinition().getExpressions());
if (particleProbabilityRate.getExpressions().length > 0) {
rateConstantElement.setAttribute("value", rateConstantValue);
}
listOfRateConstantsElement.addContent(rateConstantElement);
rateLawElement.addContent(listOfRateConstantsElement);
}
reactionRuleElement.addContent(rateLawElement);
// // <Map>
// // <MapItem sourceID="RR1_RP1_M1" targetID="RR1_PP1_M1"/>
// // <MapItem sourceID="RR1_RP1_M1_C1" targetID="RR1_PP1_M1_C1"/>
// // <MapItem sourceID="RR1_RP1_M1_C2" targetID="RR1_PP1_M1_C2"/>
// // <MapItem sourceID="RR1_RP2_M1" targetID="RR1_PP1_M2"/>
// // <MapItem sourceID="RR1_RP2_M1_C1" targetID="RR1_PP1_M2_C1"/>
// // </Map>
// Element mapElement = new Element("Map");
// System.out.println("----------------------------------------------------------------------");
// for(MolecularTypeOfReactionParticipant p : currentReactantElementsOfReaction) {
// System.out.println(p.moleculeName + ", " + p.elementID);
// }
// for(ComponentOfMolecularTypeOfReactionParticipant c : currentComponentOfReactantElementsOfReaction) {
// System.out.println(c.moleculeName + ", " + c.componentName + ", " + c.elementID);
// }
// System.out.println("----------------------------------------------------------------------");
// for(MolecularTypeOfReactionParticipant p : currentProductElementsOfReaction) {
// System.out.println(p.moleculeName + ", " + p.elementID);
// }
// for(ComponentOfMolecularTypeOfReactionParticipant c : currentComponentOfProductElementsOfReaction) {
// System.out.println(c.moleculeName + ", " + c.componentName + ", " + c.elementID);
// }
// System.out.println("----------------------------------------------------------------------");
//
// List<MolecularTypeOfReactionParticipant> cloneOfReactants = new ArrayList<MolecularTypeOfReactionParticipant>(currentReactantElementsOfReaction);
// List<MolecularTypeOfReactionParticipant> cloneOfProducts = new ArrayList<MolecularTypeOfReactionParticipant>(currentProductElementsOfReaction);
// for(Iterator<MolecularTypeOfReactionParticipant> itReactant = cloneOfReactants.iterator(); itReactant.hasNext();) { // participants
// MolecularTypeOfReactionParticipant reactant = itReactant.next();
// boolean foundProduct = false;
// for(Iterator<MolecularTypeOfReactionParticipant> itProduct = cloneOfProducts.iterator(); itProduct.hasNext();) {
// MolecularTypeOfReactionParticipant product = itProduct.next();
// if(reactant.find(product)) {
// MappingOfReactionParticipants m = new MappingOfReactionParticipants(reactant.elementID, product.elementID, "");
// currentMappingOfReactionParticipants.add(m );
// itProduct.remove();
// foundProduct = true;
// break; // we exit inner loop if we find a match for current reactant
// }
// }
// if(foundProduct == false) {
// System.out.println("Did not found a match for reactant " + reactant.moleculeName + ", " + reactant.elementID);
// }
// itReactant.remove(); // found or not, we remove the reactant
// }
// if(!currentProductElementsOfReaction.isEmpty()) {
// for(MolecularTypeOfReactionParticipant p : currentProductElementsOfReaction) {
// System.out.println("Did not found a match for product " + p.moleculeName + ", " + p.elementID);
// }
// }
// for(Iterator<ComponentOfMolecularTypeOfReactionParticipant> itReactant = currentComponentOfReactantElementsOfReaction.iterator(); itReactant.hasNext();) { // components
// ComponentOfMolecularTypeOfReactionParticipant reactant = itReactant.next();
// boolean foundProduct = false;
// for(Iterator<ComponentOfMolecularTypeOfReactionParticipant> itProduct = currentComponentOfProductElementsOfReaction.iterator(); itProduct.hasNext();) {
// ComponentOfMolecularTypeOfReactionParticipant product = itProduct.next();
// String state = "";
// if(reactant.find(product)) {
// if(!reactant.state.equals(product.state)) {
// state = product.state;
// }
// MappingOfReactionParticipants m = new MappingOfReactionParticipants(reactant.elementID, product.elementID, state);
// currentMappingOfReactionParticipants.add(m );
// itProduct.remove();
// foundProduct = true;
// break; // we exit inner loop if we find a match for current reactant
// }
// }
// if(foundProduct == false) {
// System.out.println("Did not found a match for reactant " + reactant.moleculeName + ", " + reactant.elementID);
// }
// itReactant.remove(); // found or not, we remove the reactant
// }
// if(!currentComponentOfProductElementsOfReaction.isEmpty()) {
// for(ComponentOfMolecularTypeOfReactionParticipant p : currentComponentOfProductElementsOfReaction) {
// System.out.println("Did not found a match for product " + p.moleculeName + ", " + p.elementID);
// }
// }
// for(Iterator<MappingOfReactionParticipants> it = currentMappingOfReactionParticipants.iterator(); it.hasNext();) {
// MappingOfReactionParticipants m = it.next();
// Element mapItemElement = new Element("MapItem");
// mapItemElement.setAttribute("sourceID", m.reactantElementID);
// mapItemElement.setAttribute("targetID", m.productElementID);
// mapElement.addContent(mapItemElement);
// }
// reactionRuleElement.addContent(mapElement);
//
// // <ListOfOperations>
// // <AddBond site1="RR1_RP1_M1_C1" site2="RR1_RP2_M1_C1"/>
// // <StateChange site="RR0_RP0_M0_C2" finalState="Y"/>
// // </ListOfOperations>
// Element listOfOperationsElement = new Element("ListOfOperations");
//
// // AddBond elements
// // add any bond in the product which is not present in the reactant
// Iterator<BondSites> it = patternProductBondSites.iterator();
// while (it.hasNext()) {
// BondSites bs = it.next();
// String reactantS1 = MappingOfReactionParticipants.findMatchingReactant(bs.component1, currentMappingOfReactionParticipants);
// String reactantS2 = MappingOfReactionParticipants.findMatchingReactant(bs.component2, currentMappingOfReactionParticipants);
// // we check if the bonds in the product existed already in the reactant, in which case they were not "added" in this reaction
// BondSites candidate = new BondSites(reactantS1, reactantS2);
// boolean preExistent = false;
// for(BondSites bsReactant : reactionReactantBondSites) {
// if(bsReactant.equals(candidate)) {
// preExistent = true;
// break;
// }
// }
// if(preExistent == true) {
// continue; // we don't add preexisting bonds
// }
// Element addBondElement = new Element("AddBond");
// addBondElement.setAttribute("site1", reactantS1);
// addBondElement.setAttribute("site2", reactantS2);
// listOfOperationsElement.addContent(addBondElement);
// }
// // StateChange elements
// for(Iterator<MappingOfReactionParticipants> it1 = currentMappingOfReactionParticipants.iterator(); it1.hasNext();) {
// MappingOfReactionParticipants m = it1.next();
// if(!m.componentFinalState.equals("")) { // state has changed if it's different from ""
// Element stateChangeElement = new Element("StateChange");
// stateChangeElement.setAttribute("site", m.reactantElementID);
// stateChangeElement.setAttribute("finalState", m.componentFinalState);
// listOfOperationsElement.addContent(stateChangeElement);
// }
// }
// // eliminate all the common entries (molecule types) in reactants and products
// // what's left in reactants was deleted, what's left in products was added
// List<MolecularTypeOfReactionParticipant> commonParticipants = new ArrayList<MolecularTypeOfReactionParticipant>();
// for(Iterator<MolecularTypeOfReactionParticipant> itReactant = currentReactantElementsOfReaction.iterator(); itReactant.hasNext();) { // participants
// MolecularTypeOfReactionParticipant reactant = itReactant.next();
// for(Iterator<MolecularTypeOfReactionParticipant> itProduct = currentProductElementsOfReaction.iterator(); itProduct.hasNext();) {
// MolecularTypeOfReactionParticipant product = itProduct.next();
// if(reactant.find(product)) {
// // commonParticipants contains the reactant molecules with a equivalent molecule in the product (meaning they are not in the "Deleted" category)
// commonParticipants.add(reactant);
// itReactant.remove();
// itProduct.remove();
// break; // we exit inner loop if we find a match for current reactant
// }
// }
// }
// // DeleteBond element
// // there is no need to mention deletion of bond if the particleSpeciesPattern
// // or the MolecularType involved in the bond are deleted as well
// // We only keep those "Deleted" bonds which belong to the molecules (of the reactant) present in commonParticipants
// // Both components (sites) of the bond need to have their molecules in commonParticipants
// boolean foundMoleculeForComponent1 = false;
// boolean foundMoleculeForComponent2 = false;
// HashSet<BondSites> cloneOfReactantBondSites = new HashSet<BondSites>(patternReactantBondSites);
// Iterator<BondSites> itbs = cloneOfReactantBondSites.iterator();
// while (itbs.hasNext()) {
// BondSites bs = itbs.next();
// String bondComponent1MoleculeId = BondSites.extractMoleculeId(bs.component1);
// String bondComponent2MoleculeId = BondSites.extractMoleculeId(bs.component2);
// for(MolecularTypeOfReactionParticipant commonReactionMoleculeule : commonParticipants) {
// String commonReactantPatternId = commonReactionMoleculeule.elementID;
// if(bondComponent1MoleculeId.equals(commonReactantPatternId)) {
// foundMoleculeForComponent1 = true;
// }
// if(bondComponent2MoleculeId.equals(commonReactantPatternId)) {
// foundMoleculeForComponent2 = true;
// }
// }
// if(!foundMoleculeForComponent1 || !foundMoleculeForComponent2) {
// // at least one of bond's molecule is not in common, hence we don't need to report the deletion of this bond
// itbs.remove();
// }
// }
// // the clone has now all the deleted bonds whose molecules have not been deleted
// itbs = cloneOfReactantBondSites.iterator();
// while (itbs.hasNext()) {
// BondSites bs = itbs.next();
// Element addBondElement = new Element("DeleteBond");
// addBondElement.setAttribute("site1", bs.component1);
// addBondElement.setAttribute("site2", bs.component2);
// listOfOperationsElement.addContent(addBondElement);
// }
// // Add MolecularType element
// for(MolecularTypeOfReactionParticipant molecule : currentProductElementsOfReaction) {
// System.out.println("created molecule: " + molecule.elementID + "' " + molecule.moleculeName);
// Element addMolecularTypePatternElement = new Element("Add");
// addMolecularTypePatternElement.setAttribute("id", molecule.elementID);
// listOfOperationsElement.addContent(addMolecularTypePatternElement);
// }
// // Delete MolecularType element
// // if the reactant pattern of the molecule being deleted still exists as part of the common, then we only delete the molecule
// // if the reactant pattern of the molecule being deleted is not as part of the common, then it's gone completely and we delete the reactant pattern
// ArrayList<String> patternsToDelete = new ArrayList<String>();
// for(MolecularTypeOfReactionParticipant molecule : currentReactantElementsOfReaction) {
// String reactantPatternId = molecule.extractReactantPatternId();
// boolean found = false;
// for(MolecularTypeOfReactionParticipant common : commonParticipants) {
// String commonId = common.extractReactantPatternId();
// if(reactantPatternId.equals(commonId)) {
// found = true;
// break; // some other molecule of this pattern still there, we don't delete the pattern
// }
// }
// if(found == true) { // some other molecule of this pattern still there, we don't delete the pattern
// System.out.println("deleted molecule: " + molecule.elementID + "' " + molecule.moleculeName);
// Element addMolecularTypePatternElement = new Element("Delete");
// addMolecularTypePatternElement.setAttribute("id", molecule.elementID);
// addMolecularTypePatternElement.setAttribute("DeleteMolecules", "0");
// listOfOperationsElement.addContent(addMolecularTypePatternElement);
// } else { // no molecule of this pattern left, we delete the pattern
// if(patternsToDelete.contains(reactantPatternId)) {
// // nothing to do, we're already deleting this pattern
// break;
// } else {
// patternsToDelete.add(reactantPatternId);
// System.out.println("deleted pattern: " + reactantPatternId);
// Element addParticleSpeciesPatternElement = new Element("Delete");
// addParticleSpeciesPatternElement.setAttribute("id", reactantPatternId);
// addParticleSpeciesPatternElement.setAttribute("DeleteMolecules", "0");
// listOfOperationsElement.addContent(addParticleSpeciesPatternElement);
// }
// }
// }
// reactionRuleElement.addContent(listOfOperationsElement);
listOfReactionRules.addContent(reactionRuleElement);
}
Element listOfObservablesElement = getListOfObservables(clonedMathDesc);
Element listOfFunctionsElement = getListOfFunctions(clonedMathDesc, simulationSymbolTable);
modelElement.addContent(listOfParametersElement);
modelElement.addContent(listOfMoleculeTypesElement);
modelElement.addContent(listOfSpeciesElement);
modelElement.addContent(listOfReactionRules);
modelElement.addContent(listOfObservablesElement);
modelElement.addContent(listOfFunctionsElement);
sbmlElement.addContent(modelElement);
// // return e1;
return sbmlElement;
}
use of cbit.vcell.math.JumpProcessRateDefinition in project vcell by virtualcell.
the class SmoldynFileWriter method writeReactions.
private void writeReactions() throws ExpressionException, MathException {
printWriter.println("# reactions");
Enumeration<SubDomain> subdomains = mathDesc.getSubDomains();
while (subdomains.hasMoreElements()) {
SubDomain subdomain = subdomains.nextElement();
for (ParticleJumpProcess pjp : subdomain.getParticleJumpProcesses()) {
ArrayList<Variable> reactants = new ArrayList<Variable>();
ArrayList<Variable> products = new ArrayList<Variable>();
for (Action a : pjp.getActions().toArray(new Action[pjp.getActions().size()])) {
if (a.getOperation().equals(Action.ACTION_CREATE)) {
products.add(a.getVar());
} else if (a.getOperation().equals(Action.ACTION_DESTROY)) {
reactants.add(a.getVar());
}
}
Expression rateDefinition = null;
JumpProcessRateDefinition jprd = pjp.getParticleRateDefinition();
if (jprd instanceof MacroscopicRateConstant) {
rateDefinition = subsituteFlatten(((MacroscopicRateConstant) jprd).getExpression());
} else if (jprd instanceof InteractionRadius) {
rateDefinition = subsituteFlatten(((InteractionRadius) jprd).getExpression());
} else {
new RuntimeException("The jump process rate definition is not supported");
}
if (rateDefinition.isZero()) {
continue;
}
if (mathDesc.isSpatialHybrid()) {
String[] symbols = rateDefinition.getSymbols();
if (symbols != null) {
if (subdomain instanceof MembraneSubDomain) {
rateDefinition = new Expression(FiniteVolumeFileWriter.replaceVolumeVariable(getSimulationTask(), (MembraneSubDomain) subdomain, rateDefinition));
}
}
} else {
try {
rateDefinition.evaluateConstant();
} catch (ExpressionException ex) {
throw new ExpressionException("reaction rate for jump process " + pjp.getName() + " is not a constant. Constants are required for all reaction rates.");
}
}
// Smoldyn takes maximum 2nd order reaction.
if (reactants.size() > 2) {
throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 reactants.");
}
if (products.size() > 2) {
throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 products.");
}
String rateDefinitionStr = simulation.getMathDescription().isSpatialHybrid() ? rateDefinition.infix() + ";" : rateDefinition.evaluateConstant() + "";
if (subdomain instanceof CompartmentSubDomain) {
// 0th order reaction, product limited to one and we'll let the reaction know where it happens
if (reactants.size() == 0 && products.size() == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_cmpt + " " + subdomain.getName() + " " + pjp.getName() + " ");
} else {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction + " " + /* + subdomain.getName() + " "*/
pjp.getName() + " ");
}
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else if (subdomain instanceof MembraneSubDomain) {
// 0th order reaction, product limited to one and it can be on mem or in vol
if (reactants.size() == 0 && products.size() == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // consuming of a species to nothing, limited to one reactant
if (reactants.size() == 1 && products.size() == 0) {
if (// consuming a mem species in mem reaction
getMembraneVariableCount(reactants) == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // it equals to adsorption, species A from vol adsorbed to mem as again species A, and then we kill the speceis A on mem.
if (getVolumeVariableCount(reactants) == 1) {
writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
String speciesName = reactants.get(0).getName();
String killMolCmd = "cmd " + SmoldynVCellMapper.SmoldynKeyword.E + " " + SmoldynVCellMapper.SmoldynKeyword.killmol + " " + speciesName + "(" + SmoldynKeyword.up + ")";
killMolCommands.add(killMolCmd);
}
} else // Use rate command for any membrane reactions with 1 reactant and 1 product
if ((reactants.size() == 1) && (products.size() == 1)) {
// Membrane reaction (1 react to 1 product).
if (getMembraneVariableCount(products) == 1 && getMembraneVariableCount(reactants) == 1) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else // Other single molecular reactions
{
writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
}
} else // membrane reactions which are not one to one, or 0th order, or consuming species
{
if (// membrane reaction has one membrane bound reactant
(getMembraneVariableCount(reactants) == 1)) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
} else if (// bimolecular membrane reaction
getMembraneVariableCount(reactants) == 2) {
if (jprd instanceof InteractionRadius) {
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionByInteractionRadius(reactants, products, subdomain, rateDefinitionStr, pjp.getName());
} else {
// throw new MathException("Error with reaction: " + pjp.getName() + ".\nVCell Spatial stochastic modeling requires macroscopic or microscopic kinetics for bimolecular membrane reactions.");
printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
}
} else if (getMembraneVariableCount(reactants) == 0) {
throw new MathException("Error with reaction: " + pjp.getName() + ".\nIn VCell spatial stochastic modeling, the membrane reaction requires at least one membrane bound reactant.");
}
}
}
}
}
printWriter.println();
}
use of cbit.vcell.math.JumpProcessRateDefinition in project vcell by virtualcell.
the class Xmlproducer method getXML.
private org.jdom.Element getXML(ParticleJumpProcess param) {
org.jdom.Element particleJumpProcessElement = new org.jdom.Element(XMLTags.ParticleJumpProcessTag);
// name
particleJumpProcessElement.setAttribute(XMLTags.NameAttrTag, mangle(param.getName()));
if (param.getProcessSymmetryFactor() != null) {
particleJumpProcessElement.setAttribute(XMLTags.ProcessSymmetryFactorAttrTag, Double.toString(param.getProcessSymmetryFactor().getFactor()));
}
// Selected Particle
for (ParticleVariable vpv : param.getParticleVariables()) {
Element e = new Element(XMLTags.SelectedParticleTag);
e.setAttribute(XMLTags.NameAttrTag, mangle(vpv.getName()));
particleJumpProcessElement.addContent(e);
}
// probability rate
Element prob = null;
JumpProcessRateDefinition particleProbabilityRate = param.getParticleRateDefinition();
if (particleProbabilityRate instanceof MacroscopicRateConstant) {
prob = new Element(XMLTags.MacroscopicRateConstantTag);
prob.addContent(mangleExpression(((MacroscopicRateConstant) particleProbabilityRate).getExpression()));
} else if (particleProbabilityRate instanceof InteractionRadius) {
prob = new Element(XMLTags.InteractionRadiusTag);
prob.addContent(mangleExpression(((InteractionRadius) particleProbabilityRate).getExpression()));
} else {
throw new RuntimeException("ParticleRateDefinition in XmlProducer not implemented");
}
particleJumpProcessElement.addContent(prob);
// Actions
for (Action action : param.getActions()) {
particleJumpProcessElement.addContent(getXML(action));
}
return particleJumpProcessElement;
}
Aggregations