Example 51 with SubVolume
use of cbit.vcell.geometry.SubVolume in project vcell by virtualcell.
the class StochMathMapping_4_8 method refreshMathDescription.
/**
* set up a math description based on current simulationContext.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// use local variable instead of using getter all the time.
SimulationContext simContext = getSimulationContext();
// local structure mapping list
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
// We have to check if all the reactions are able to tranform to stochastic jump processes before generating the math.
String stochChkMsg = simContext.getModel().isValidForStochApp();
if (!(stochChkMsg.equals(""))) {
throw new ModelException("Problem updating math description: " + simContext.getName() + "\n" + stochChkMsg);
}
// All sizes must be set for new ODE models and ratios must be set for old ones.
simContext.checkValidity();
//
// verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
//
Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
for (int i = 0; i < structures.length; i++) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
if (sm == null || (sm instanceof FeatureMapping && getSubVolume(((FeatureMapping) sm)) == null)) {
throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subVolume");
}
if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
try {
if (volFractExp != null) {
double volFract = volFractExp.evaluateConstant();
if (volFract >= 1.0) {
throw new MappingException("model structure '" + (getSimulationContext().getModel().getStructureTopology().getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0"));
}
}
} catch (ExpressionException e) {
e.printStackTrace(System.out);
}
}
}
SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
if (getStructures(subVolumes[i]) == null || getStructures(subVolumes[i]).length == 0) {
throw new MappingException("geometry subVolume '" + subVolumes[i].getName() + "' not mapped from a model structure");
}
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
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);
}
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = simContext.getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(simContext.getName() + "_generated");
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates
//
VariableHash varHash = new VariableHash();
//
// conversion factors
//
Model model = simContext.getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
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.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof StochVolVariable) {
varHash.addVariable(scm.getVariable());
}
}
//
// add rate term for all reactions
// add current source terms for each reaction step in a membrane
//
/*for (int i = 0; i < reactionSteps.length; i++){
boolean bAllReactionParticipantsFixed = true;
ReactionParticipant rp_Array[] = reactionSteps[i].getReactionParticipants();
for (int j = 0; j < rp_Array.length; j++) {
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[j].getSpeciesContext());
if (!(rp_Array[j] instanceof Catalyst) && !scs.isConstant()){
bAllReactionParticipantsFixed = false; // found at least one reactionParticipant that is not fixed and needs this rate
}
}
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionSteps[i].getStructure());
}---don't think it's useful, isn't it?*/
// deals with model parameters
ModelParameter[] modelParameters = simContext.getModel().getModelParameters();
for (int j = 0; j < modelParameters.length; j++) {
Expression expr = getSubstitutedExpr(modelParameters[j].getExpression(), true, false);
expr = getIdentifierSubstitutions(expr, modelParameters[j].getUnitDefinition(), null);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), expr));
}
// added July 2009, ElectricalStimulusParameter electric mapping tab
ElectricalStimulus[] elecStimulus = simContext.getElectricalStimuli();
if (elecStimulus.length > 0) {
throw new MappingException("Modles with electrophysiology are not supported for stochastic applications.");
}
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
Parameter initialVoltageParm = memMapping.getInitialVoltageParameter();
try {
Expression exp = initialVoltageParm.getExpression();
exp.evaluateConstant();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
}
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
if (rs.getKinetics() instanceof LumpedKinetics) {
throw new RuntimeException("Lumped Kinetics not yet supported for Stochastic Math Generation");
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(rs.getStructure());
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
if ((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
continue;
}
// don't add rate, we'll do it later when creating the jump processes
if (parameters[i].getRole() != Kinetics.ROLE_ReactionRate) {
Expression expr = getSubstitutedExpr(parameters[i].getExpression(), true, false);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(expr, parameters[i].getUnitDefinition(), sm)));
}
}
}
}
// the parameter "Size" is already put into mathsymbolmapping in refreshSpeciesContextMapping()
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm.getExpression() != null) {
try {
double value = parm.getExpression().evaluateConstant();
varHash.addVariable(new Constant(getMathSymbol(parm, sm), new Expression(value)));
} catch (ExpressionException e) {
// varHash.addVariable(new Function(getMathSymbol0(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
e.printStackTrace(System.out);
throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
}
}
}
//
// species initial values (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
// can be concentration or amount
SpeciesContextSpec.SpeciesContextSpecParameter initParam = null;
Expression iniExp = null;
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (speciesContextSpecs[i].getInitialConcentrationParameter() != null && speciesContextSpecs[i].getInitialConcentrationParameter().getExpression() != null) {
// use concentration, need to set up amount functions
initParam = speciesContextSpecs[i].getInitialConcentrationParameter();
iniExp = initParam.getExpression();
iniExp = getSubstitutedExpr(iniExp, true, !speciesContextSpecs[i].isConstant());
// now create the appropriate function or Constant for the speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
// add function for initial amount
SpeciesContextSpec.SpeciesContextSpecParameter initAmountParam = speciesContextSpecs[i].getInitialCountParameter();
Expression iniAmountExp = getExpressionConcToAmt(new Expression(initParam, getNameScope()), speciesContextSpecs[i].getSpeciesContext());
// iniAmountExp.bindExpression(this);
varHash.addVariable(new Function(getMathSymbol(initAmountParam, sm), getIdentifierSubstitutions(iniAmountExp, initAmountParam.getUnitDefinition(), sm), nullDomain));
} else if (speciesContextSpecs[i].getInitialCountParameter() != null && speciesContextSpecs[i].getInitialCountParameter().getExpression() != null) {
// use amount
initParam = speciesContextSpecs[i].getInitialCountParameter();
iniExp = initParam.getExpression();
iniExp = getSubstitutedExpr(iniExp, false, !speciesContextSpecs[i].isConstant());
// now create the appropriate function or Constant for the speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
}
// add spConcentration (concentration of species) to varHash as function or constant
SpeciesConcentrationParameter spConcParam = getSpeciesConcentrationParameter(speciesContextSpecs[i].getSpeciesContext());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(spConcParam, sm), getIdentifierSubstitutions(spConcParam.getExpression(), spConcParam.getUnitDefinition(), sm)));
}
//
// 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());
}
}
//
if (simContext.getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// functions: species which is not a variable, but has dependency expression
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Expression exp = scm.getDependencyExpression();
exp.bindExpression(this);
SpeciesCountParameter spCountParam = getSpeciesCountParameter(scm.getSpeciesContext());
varHash.addVariable(new Function(getMathSymbol(spCountParam, sm), getIdentifierSubstitutions(exp, spCountParam.getUnitDefinition(), sm), nullDomain));
}
}
//
// create subDomains
//
SubDomain subDomain = null;
subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int j = 0; j < subVolumes.length; j++) {
SubVolume subVolume = (SubVolume) subVolumes[j];
//
// get priority of subDomain
//
int priority;
Feature spatialFeature = getResolvedFeature(subVolume);
if (spatialFeature == null) {
if (simContext.getGeometryContext().getGeometry().getDimension() > 0) {
throw new MappingException("no compartment (in Physiology) is mapped to subdomain '" + subVolume.getName() + "' (in Geometry)");
} else {
priority = CompartmentSubDomain.NON_SPATIAL_PRIORITY;
}
} else {
// now does not have to match spatial feature, *BUT* needs to be unique
priority = j;
}
subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
}
// ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded()) {
continue;
}
// get the reaction
ReactionStep reactionStep = reactionSpecs[i].getReactionStep();
Kinetics kinetics = reactionStep.getKinetics();
// the structure where reaction happens
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionStep.getStructure());
// create symbol table for jump process based on reactionStep and structure mapping
// final ReactionStep finalRS = reactionStep;
// final StructureMapping finalSM = sm;
// SymbolTable symTable = new SymbolTable(){
// public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
// SymbolTableEntry ste = finalRS.getEntry(identifierString);
// if(ste == null)
// {
// ste = finalSM.getEntry(identifierString);
// }
// return ste;
// }
// };
// Different ways to deal with simple reactions and flux reactions
// probability parameter from modelUnitSystem
VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
if (// simple reactions
reactionStep instanceof SimpleReaction) {
// 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;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward).getExpression();
reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse).getExpression();
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(null, null, 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();
}
}
}
/*else if (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_irreversible.getName())==0)
{
forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Km).getExpression();
}
else if (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_reversible.getName())==0)
{
forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmFwd).getExpression();
reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmRev).getExpression();
}*/
boolean isForwardRatePresent = false;
boolean isReverseRatePresent = false;
if (forwardRate != null) {
isForwardRatePresent = true;
}
if (reverseRate != null) {
isReverseRatePresent = true;
}
// we process it as forward reaction
if ((isForwardRatePresent)) /*|| ((forwardRate == null) && (reverseRate == null))*/
{
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// get probability
Expression exp = null;
// reactions are mass actions
exp = getProbabilityRate(reactionStep, true);
// bind symbol table before substitute identifiers in the reaction step
exp.bindExpression(this);
MathMapping_4_8.ProbabilityParameter probParm = null;
try {
probParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
// actions
ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
for (int j = 0; j < reacPart.length; j++) {
Action action = null;
SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
if (reacPart[j] instanceof Reactant) {
// check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Reactant) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(stoi)));
jp.addAction(action);
}
} else if (reacPart[j] instanceof Product) {
// check if the product is a constant. If the product is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Product) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(stoi));
jp.addAction(action);
}
}
}
// add jump process to compartment subDomain
subDomain.addJumpProcess(jp);
}
if (// one more jump process for a reversible reaction
isReverseRatePresent) {
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + "_reverse";
Expression exp = null;
// reactions are mass actions
exp = getProbabilityRate(reactionStep, false);
// bind symbol table before substitute identifiers in the reaction step
exp.bindExpression(this);
MathMapping_4_8.ProbabilityParameter probRevParm = null;
try {
probRevParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
// actions
ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
for (int j = 0; j < reacPart.length; j++) {
Action action = null;
SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
if (reacPart[j] instanceof Reactant) {
// check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Reactant) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(stoi));
jp.addAction(action);
}
} else if (reacPart[j] instanceof Product) {
// check if the product is a constant. If the product is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Product) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(stoi)));
jp.addAction(action);
}
}
}
// add jump process to compartment subDomain
subDomain.addJumpProcess(jp);
}
// end of if(isForwardRateNonZero), if(isReverseRateNonRate)
} else if (// flux reactions
reactionStep instanceof FluxReaction) {
// we could set jump processes for general flux rate in forms of p1*Sout + p2*Sin
if (kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
Expression fluxRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
// we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(null, null, fluxRate, (FluxReaction) reactionStep);
// create jump process for forward flux if it exists.
if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
// jump process name
// +"_reverse";
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// we do it here instead of fluxsolver, coz we need to use getMathSymbol0(), structuremapping...etc.
Expression rate = fluxFunc.getForwardRate();
// get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
SpeciesContext scOut = fluxFunc.getReactants().get(0).getSpeciesContext();
Expression speciesFactor = null;
if (scOut.getStructure() instanceof Feature) {
Expression exp1 = new Expression(1.0 / 602.0);
Expression exp2 = new Expression(scOut.getStructure().getStructureSize(), getNameScope());
speciesFactor = Expression.div(Expression.invert(exp1), exp2);
} else {
throw new MappingException("Species involved in a flux have to be volume species.");
}
Expression speciesExp = Expression.mult(speciesFactor, new Expression(scOut, getNameScope()));
// get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
Expression expr1 = Expression.mult(rate, speciesExp);
Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
Expression exp = new Expression(1.0 / 602.0);
Expression probExp = Expression.mult(numeratorExpr, exp);
// bind symbol table before substitute identifiers in the reaction step
probExp.bindExpression(reactionStep);
MathMapping_4_8.ProbabilityParameter probParm = null;
try {
probParm = addProbabilityParameter("P_" + jpName, probExp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(probExp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
// actions
Action action = null;
SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
jp.addAction(action);
}
sc = fluxFunc.getProducts().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
jp.addAction(action);
}
subDomain.addJumpProcess(jp);
}
if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
// jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + "_reverse";
Expression rate = fluxFunc.getReverseRate();
// get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
SpeciesContext scIn = fluxFunc.getProducts().get(0).getSpeciesContext();
Expression speciesFactor = null;
if (scIn.getStructure() instanceof Feature) {
Expression exp1 = new Expression(1.0 / 602.0);
Expression exp2 = new Expression(scIn.getStructure().getStructureSize(), getNameScope());
speciesFactor = Expression.div(Expression.invert(exp1), exp2);
} else {
throw new MappingException("Species involved in a flux have to be volume species.");
}
Expression speciesExp = Expression.mult(speciesFactor, new Expression(scIn, getNameScope()));
// get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
Expression expr1 = Expression.mult(rate, speciesExp);
Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
Expression exp = new Expression(1.0 / 602.0);
Expression probRevExp = Expression.mult(numeratorExpr, exp);
// bind symbol table before substitute identifiers in the reaction step
probRevExp.bindExpression(reactionStep);
MathMapping_4_8.ProbabilityParameter probRevParm = null;
try {
probRevParm = addProbabilityParameter("P_" + jpName, probRevExp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(probRevExp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
// actions
Action action = null;
SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
jp.addAction(action);
}
sc = fluxFunc.getProducts().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
jp.addAction(action);
}
subDomain.addJumpProcess(jp);
}
}
}
// end of if (simplereaction)...else if(fluxreaction)
}
// end of reaction step loop
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
// set up variable initial conditions in subDomain
SpeciesContextSpec[] scSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
// get stochastic variable by name
SpeciesCountParameter spCountParam = getSpeciesCountParameter(speciesContextSpecs[i].getSpeciesContext());
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String varName = getMathSymbol(spCountParam, sm);
if (scSpecs[i].isConstant()) {
continue;
}
StochVolVariable var = (StochVolVariable) mathDesc.getVariable(varName);
// stochastic use initial number of particles
SpeciesContextSpec.SpeciesContextSpecParameter initParm = scSpecs[i].getInitialCountParameter();
// stochastic variables initial expression.
if (initParm != null) {
VarIniCondition varIni = new VarIniCount(var, new Expression(getMathSymbol(initParm, sm)));
subDomain.addVarIniCondition(varIni);
}
}
if (!mathDesc.isValid()) {
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
}
Also used :
VarIniCondition(cbit.vcell.math.VarIniCondition)
MembraneMapping(cbit.vcell.mapping.MembraneMapping)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping)
Product(cbit.vcell.model.Product)
FluxReaction(cbit.vcell.model.FluxReaction)
SpeciesContext(cbit.vcell.model.SpeciesContext)
SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec)
Feature(cbit.vcell.model.Feature)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
MappingException(cbit.vcell.mapping.MappingException)
PropertyVetoException(java.beans.PropertyVetoException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
ModelException(cbit.vcell.model.ModelException)
ReactionSpec(cbit.vcell.mapping.ReactionSpec)
PropertyVetoException(java.beans.PropertyVetoException)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ReactionStep(cbit.vcell.model.ReactionStep)
Kinetics(cbit.vcell.model.Kinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
Action(cbit.vcell.math.Action)
VariableHash(cbit.vcell.math.VariableHash)
Constant(cbit.vcell.math.Constant)
StructureMapping(cbit.vcell.mapping.StructureMapping)
Function(cbit.vcell.math.Function)
FeatureMapping(cbit.vcell.mapping.FeatureMapping)
JumpProcess(cbit.vcell.math.JumpProcess)
Structure(cbit.vcell.model.Structure)
StochVolVariable(cbit.vcell.math.StochVolVariable)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
SimpleReaction(cbit.vcell.model.SimpleReaction)
VarIniCount(cbit.vcell.math.VarIniCount)
SimulationContext(cbit.vcell.mapping.SimulationContext)
ElectricalStimulus(cbit.vcell.mapping.ElectricalStimulus)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
ProxyParameter(cbit.vcell.model.ProxyParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
Example 52 with SubVolume
use of cbit.vcell.geometry.SubVolume in project vcell by virtualcell.
the class GeometryContext method setDefaultUnitSizes.
private void setDefaultUnitSizes() throws PropertyVetoException {
if (getGeometry().getDimension() > 0) {
for (StructureMapping sm : fieldStructureMappings) {
StructureMapping[] sms = getStructureMappings(sm.getGeometryClass());
StructureMappingParameter unitSizeParameter = sm.getUnitSizeParameter();
if (unitSizeParameter == null) {
continue;
}
Expression exp = unitSizeParameter.getExpression();
if (sm instanceof MembraneMapping) {
// Membrane mapped to surface or subdomain, default to 1.0
if (exp == null) {
try {
unitSizeParameter.setExpression(new Expression(1.0));
} catch (ExpressionBindingException e) {
e.printStackTrace();
}
}
} else if (sm instanceof FeatureMapping) {
// Feature mapped to subdomain
if (sm.getGeometryClass() instanceof SubVolume) {
if (sms != null && sms.length == 1) {
try {
unitSizeParameter.setExpression(new Expression(1.0));
} catch (ExpressionBindingException e) {
e.printStackTrace();
}
}
} else {
if (exp == null) {
try {
unitSizeParameter.setExpression(new Expression(1.0));
} catch (ExpressionBindingException e) {
e.printStackTrace();
}
}
}
}
}
}
}
Also used :
Expression(cbit.vcell.parser.Expression)
SubVolume(cbit.vcell.geometry.SubVolume)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException)
Example 53 with SubVolume
use of cbit.vcell.geometry.SubVolume in project vcell by virtualcell.
the class MembraneStructureAnalyzer method refreshResolvedFluxes.
/**
* This method was created in VisualAge.
*/
private void refreshResolvedFluxes() throws Exception {
// System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes()");
ModelUnitSystem unitSystem = mathMapping.getSimulationContext().getModel().getUnitSystem();
GeometryContext geoContext = mathMapping.getSimulationContext().getGeometryContext();
Vector<ResolvedFlux> resolvedFluxList = new Vector<ResolvedFlux>();
//
// for each reaction, get all fluxReactions associated with this membrane
//
Vector<FluxReaction> fluxList = new Vector<FluxReaction>();
ReactionSpec[] reactionSpecs = mathMapping.getSimulationContext().getReactionContext().getReactionSpecs();
for (int j = 0; j < reactionSpecs.length; j++) {
if (reactionSpecs[j].isExcluded()) {
continue;
}
ReactionStep rs = reactionSpecs[j].getReactionStep();
if (rs.getStructure() != null && geoContext.getStructureMapping(rs.getStructure()).getGeometryClass() == surfaceClass) {
if (rs instanceof FluxReaction) {
fluxList.addElement((FluxReaction) rs);
}
}
}
//
for (int i = 0; i < fluxList.size(); i++) {
FluxReaction fr = fluxList.elementAt(i);
ReactionParticipant[] reactionParticipants = fr.getReactionParticipants();
for (int j = 0; j < reactionParticipants.length; j++) {
if (!(reactionParticipants[j] instanceof Reactant) && !(reactionParticipants[j] instanceof Product)) {
continue;
}
ResolvedFlux rf = null;
SpeciesContext speciesContext = reactionParticipants[j].getSpeciesContext();
for (int k = 0; k < resolvedFluxList.size(); k++) {
ResolvedFlux rf_tmp = resolvedFluxList.elementAt(k);
if (rf_tmp.getSpeciesContext() == reactionParticipants[j].getSpeciesContext()) {
rf = rf_tmp;
}
}
//
// if speciesContext is not "fixed" and is mapped to a volume, add flux to ResolvedFlux
//
StructureMapping structureMapping = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionParticipants[j].getStructure());
if (structureMapping.getGeometryClass() == surfaceClass) {
// flux within surface
continue;
}
if (structureMapping.getGeometryClass() instanceof SubVolume && surfaceClass.isAdjacentTo((SubVolume) structureMapping.getGeometryClass())) {
SpeciesContextSpec speciesContextSpec = mathMapping.getSimulationContext().getReactionContext().getSpeciesContextSpec(speciesContext);
if (!speciesContextSpec.isConstant()) {
if (rf == null) {
VCUnitDefinition speciesFluxUnit = speciesContext.getUnitDefinition().multiplyBy(unitSystem.getLengthUnit()).divideBy(unitSystem.getTimeUnit());
rf = new ResolvedFlux(speciesContext, speciesFluxUnit);
resolvedFluxList.addElement(rf);
}
FeatureMapping featureMapping = (FeatureMapping) structureMapping;
Expression insideFluxCorrection = Expression.invert(new Expression(featureMapping.getVolumePerUnitVolumeParameter(), mathMapping.getNameScope()));
//
if (fr.getKinetics() instanceof DistributedKinetics) {
KineticsParameter reactionRateParameter = ((DistributedKinetics) fr.getKinetics()).getReactionRateParameter();
Expression correctedReactionRate = Expression.mult(new Expression(reactionRateParameter, mathMapping.getNameScope()), insideFluxCorrection);
if (reactionParticipants[j] instanceof Product) {
if (rf.getFluxExpression().isZero()) {
rf.setFluxExpression(correctedReactionRate.flatten());
} else {
rf.setFluxExpression(Expression.add(rf.getFluxExpression(), correctedReactionRate.flatten()));
}
} else if (reactionParticipants[j] instanceof Reactant) {
if (rf.getFluxExpression().isZero()) {
rf.setFluxExpression(Expression.negate(correctedReactionRate).flatten());
} else {
rf.setFluxExpression(Expression.add(rf.getFluxExpression(), Expression.negate(correctedReactionRate).flatten()));
}
} else {
throw new RuntimeException("expected either FluxReactant or FluxProduct");
}
} else if (fr.getKinetics() instanceof LumpedKinetics) {
throw new RuntimeException("Lumped Kinetics for fluxes not yet supported");
} else {
throw new RuntimeException("unexpected Kinetic type in MembraneStructureAnalyzer.refreshResolvedFluxes()");
}
rf.getFluxExpression().bindExpression(mathMapping);
}
}
}
}
//
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded()) {
continue;
}
ReactionStep rs = reactionSpecs[i].getReactionStep();
if (rs.getStructure() != null && geoContext.getStructureMapping(rs.getStructure()).getGeometryClass() == surfaceClass) {
if (rs instanceof SimpleReaction) {
SimpleReaction sr = (SimpleReaction) rs;
ReactionParticipant[] rp_Array = sr.getReactionParticipants();
for (int k = 0; k < rp_Array.length; k++) {
if (rp_Array[k] instanceof Reactant || rp_Array[k] instanceof Product) {
SpeciesContextSpec rpSCS = mathMapping.getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[k].getSpeciesContext());
StructureMapping rpSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(rp_Array[k].getStructure());
if (rpSM.getGeometryClass() instanceof SubVolume && !rpSCS.isConstant()) {
//
// get ResolvedFlux for this species
//
ResolvedFlux rf = null;
for (int j = 0; j < resolvedFluxList.size(); j++) {
ResolvedFlux rf_tmp = (ResolvedFlux) resolvedFluxList.elementAt(j);
if (rf_tmp.getSpeciesContext() == rp_Array[k].getSpeciesContext()) {
rf = rf_tmp;
}
}
if (rf == null) {
VCUnitDefinition speciesFluxUnit = rp_Array[k].getSpeciesContext().getUnitDefinition().multiplyBy(unitSystem.getLengthUnit()).divideBy(unitSystem.getTimeUnit());
rf = new ResolvedFlux(rp_Array[k].getSpeciesContext(), speciesFluxUnit);
resolvedFluxList.addElement(rf);
}
if (rpSM.getGeometryClass() instanceof SubVolume && surfaceClass.isAdjacentTo((SubVolume) rpSM.getGeometryClass())) {
//
// for binding on inside or outside, add to ResolvedFlux.flux
//
Expression fluxRateExpression = getCorrectedRateExpression(sr, rp_Array[k], RateType.ResolvedFluxRate).renameBoundSymbols(mathMapping.getNameScope());
if (rf.getFluxExpression().isZero()) {
rf.setFluxExpression(fluxRateExpression);
} else {
rf.setFluxExpression(Expression.add(rf.getFluxExpression(), fluxRateExpression));
}
rf.getFluxExpression().bindExpression(mathMapping);
} else {
String structureName = ((rs.getStructure() != null) ? (rs.getStructure().getName()) : ("<null>"));
throw new Exception("In Application '" + mathMapping.getSimulationContext().getName() + "', SpeciesContext '" + rp_Array[k].getSpeciesContext().getName() + "' is not mapped adjacent to structure '" + structureName + "' but reacts there");
}
}
}
}
}
}
}
//
if (resolvedFluxList.size() > 0) {
resolvedFluxes = new ResolvedFlux[resolvedFluxList.size()];
resolvedFluxList.copyInto(resolvedFluxes);
} else {
resolvedFluxes = null;
}
}
Also used :
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
Product(cbit.vcell.model.Product)
FluxReaction(cbit.vcell.model.FluxReaction)
SpeciesContext(cbit.vcell.model.SpeciesContext)
Reactant(cbit.vcell.model.Reactant)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
DistributedKinetics(cbit.vcell.model.DistributedKinetics)
SimpleReaction(cbit.vcell.model.SimpleReaction)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
ReactionStep(cbit.vcell.model.ReactionStep)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
Example 54 with SubVolume
use of cbit.vcell.geometry.SubVolume in project vcell by virtualcell.
the class AbstractMathMapping method getMathSymbol0.
/**
* Substitutes appropriate variables for speciesContext bindings
*
* @return cbit.vcell.parser.Expression
* @param origExp cbit.vcell.parser.Expression
* @param structureMapping cbit.vcell.mapping.StructureMapping
*/
protected final String getMathSymbol0(SymbolTableEntry ste, GeometryClass geometryClass) throws MappingException {
String steName = ste.getName();
if (ste instanceof Kinetics.KineticsParameter) {
Integer count = localNameCountHash.get(steName);
if (count == null) {
throw new MappingException("KineticsParameter " + steName + " not found in local name count");
}
if (count > 1 || steName.equals("J")) {
return steName + "_" + ste.getNameScope().getName();
// return getNameScope().getSymbolName(ste);
} else {
return steName;
}
}
if (ste instanceof LocalParameter && ((LocalParameter) ste).getNameScope() instanceof ReactionRule.ReactionRuleNameScope) {
Integer count = localNameCountHash.get(steName);
if (count == null) {
throw new MappingException("Reaction Rule Parameter " + steName + " not found in local name count");
}
if (count > 1 || steName.equals("J")) {
return steName + "_" + ste.getNameScope().getName();
// return getNameScope().getSymbolName(ste);
} else {
return steName;
}
}
if (ste instanceof ProbabilityParameter) {
// be careful here, to see if we need mangle the reaction name
ProbabilityParameter probParm = (ProbabilityParameter) ste;
return probParm.getName() + PARAMETER_PROBABLIITY_RATE_SUFFIX;
}
if (ste instanceof SpeciesConcentrationParameter) {
SpeciesConcentrationParameter concParm = (SpeciesConcentrationParameter) ste;
return concParm.getSpeciesContext().getName() + MATH_FUNC_SUFFIX_SPECIES_CONCENTRATION;
}
if (ste instanceof SpeciesCountParameter) {
SpeciesCountParameter countParm = (SpeciesCountParameter) ste;
return countParm.getSpeciesContext().getName() + MATH_VAR_SUFFIX_SPECIES_COUNT;
}
if (ste instanceof ObservableConcentrationParameter) {
ObservableConcentrationParameter concParm = (ObservableConcentrationParameter) ste;
return concParm.getObservable().getName() + MATH_FUNC_SUFFIX_SPECIES_CONCENTRATION;
}
if (ste instanceof ObservableCountParameter) {
ObservableCountParameter countParm = (ObservableCountParameter) ste;
return countParm.getObservable().getName() + MATH_VAR_SUFFIX_SPECIES_COUNT;
}
if (ste instanceof RbmObservable) {
RbmObservable observable = (RbmObservable) ste;
return observable.getName() + MATH_FUNC_SUFFIX_SPECIES_CONCENTRATION;
}
if (ste instanceof EventAssignmentOrRateRuleInitParameter) {
EventAssignmentOrRateRuleInitParameter eventInitParm = (EventAssignmentOrRateRuleInitParameter) ste;
// + MATH_FUNC_SUFFIX_EVENTASSIGN_OR_RATE_INIT;
return eventInitParm.getName();
}
if (ste instanceof RateRuleRateParameter) {
RateRuleRateParameter rateRuleRateParm = (RateRuleRateParameter) ste;
// + MATH_FUNC_SUFFIX_RATERULE_RATE;
return rateRuleRateParm.getName();
}
if (ste instanceof Model.ReservedSymbol) {
return steName;
}
if (ste instanceof Membrane.MembraneVoltage) {
return steName;
}
if (ste instanceof Structure.StructureSize) {
Structure structure = ((Structure.StructureSize) ste).getStructure();
StructureMapping.StructureMappingParameter sizeParameter = simContext.getGeometryContext().getStructureMapping(structure).getSizeParameter();
return getMathSymbol(sizeParameter, geometryClass);
}
if (ste instanceof ProxyParameter) {
ProxyParameter pp = (ProxyParameter) ste;
return getMathSymbol0(pp.getTarget(), geometryClass);
}
//
if (ste instanceof ModelParameter) {
ModelParameter mp = (ModelParameter) ste;
return mp.getName();
}
if (ste instanceof SpeciesContextSpec.SpeciesContextSpecParameter) {
SpeciesContextSpec.SpeciesContextSpecParameter scsParm = (SpeciesContextSpec.SpeciesContextSpecParameter) ste;
SpeciesContext speciesContext = ((SpeciesContextSpec) (scsParm.getNameScope().getScopedSymbolTable())).getSpeciesContext();
SpeciesContextMapping scm = getSpeciesContextMapping(speciesContext);
String speciesContextVarName = null;
if (scm.getVariable() != null) {
speciesContextVarName = scm.getVariable().getName();
} else {
speciesContextVarName = speciesContext.getName();
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_InitialConcentration) {
return speciesContextVarName + MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX + TokenMangler.fixTokenStrict(scsParm.getUnitDefinition().getSymbol());
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_InitialCount) {
return speciesContextVarName + MATH_FUNC_SUFFIX_SPECIES_INIT_COUNT;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_DiffusionRate) {
return speciesContextVarName + PARAMETER_DIFFUSION_RATE_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_BoundaryValueXm) {
return speciesContextVarName + PARAMETER_BOUNDARY_XM_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_BoundaryValueXp) {
return speciesContextVarName + PARAMETER_BOUNDARY_XP_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_BoundaryValueYm) {
return speciesContextVarName + PARAMETER_BOUNDARY_YM_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_BoundaryValueYp) {
return speciesContextVarName + PARAMETER_BOUNDARY_YP_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_BoundaryValueZm) {
return speciesContextVarName + PARAMETER_BOUNDARY_ZM_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_BoundaryValueZp) {
return speciesContextVarName + PARAMETER_BOUNDARY_ZP_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_VelocityX) {
return speciesContextVarName + PARAMETER_VELOCITY_X_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_VelocityY) {
return speciesContextVarName + PARAMETER_VELOCITY_Y_SUFFIX;
}
if (scsParm.getRole() == SpeciesContextSpec.ROLE_VelocityZ) {
return speciesContextVarName + PARAMETER_VELOCITY_Z_SUFFIX;
}
}
if (ste instanceof ElectricalDevice.ElectricalDeviceParameter) {
ElectricalDevice.ElectricalDeviceParameter edParm = (ElectricalDevice.ElectricalDeviceParameter) ste;
ElectricalDevice electricalDevice = (ElectricalDevice) edParm.getNameScope().getScopedSymbolTable();
if (electricalDevice instanceof MembraneElectricalDevice) {
String nameWithScope = ((MembraneElectricalDevice) electricalDevice).getMembraneMapping().getMembrane().getNameScope().getName();
if (edParm.getRole() == ElectricalDevice.ROLE_TotalCurrent) {
return PARAMETER_TOTAL_CURRENT_PREFIX + nameWithScope;
}
if (edParm.getRole() == ElectricalDevice.ROLE_TransmembraneCurrent) {
return PARAMETER_TRANSMEMBRANE_CURRENT_PREFIX + nameWithScope;
}
// }else if (electricalDevice instanceof CurrentClampElectricalDevice) {
// if (edParm.getRole()==ElectricalDevice.ROLE_TotalCurrentDensity){
// return "I_"+((CurrentClampElectricalDevice)electricalDevice).getCurrentClampStimulus().getNameScope().getName();
// }
// if (edParm.getRole()==ElectricalDevice.ROLE_TransmembraneCurrentDensity){
// return "F_"+((CurrentClampElectricalDevice)electricalDevice).getCurrentClampStimulus().getNameScope().getName();
// }
// }else if (electricalDevice instanceof VoltageClampElectricalDevice) {
// if (edParm.getRole()==ElectricalDevice.ROLE_TotalCurrentDensity){
// return "I_"+((VoltageClampElectricalDevice)electricalDevice).getVoltageClampStimulus().getNameScope().getName();
// }
// if (edParm.getRole()==ElectricalDevice.ROLE_TransmembraneCurrentDensity){
// return "F_"+((VoltageClampElectricalDevice)electricalDevice).getVoltageClampStimulus().getNameScope().getName();
// }
}
}
if (ste instanceof LocalParameter && ((LocalParameter) ste).getNameScope() instanceof ElectricalStimulus.ElectricalStimulusNameScope) {
LocalParameter esParm = (LocalParameter) ste;
String nameWithScope = esParm.getNameScope().getName();
if (esParm.getRole() == ElectricalStimulus.ElectricalStimulusParameterType.TotalCurrent) {
return PARAMETER_TOTAL_CURRENT_PREFIX + nameWithScope;
} else if (esParm.getRole() == ElectricalStimulus.ElectricalStimulusParameterType.Voltage) {
return PARAMETER_VOLTAGE_PREFIX + nameWithScope;
}
}
if (ste instanceof StructureMapping.StructureMappingParameter) {
StructureMapping.StructureMappingParameter smParm = (StructureMapping.StructureMappingParameter) ste;
Structure structure = ((StructureMapping) (smParm.getNameScope().getScopedSymbolTable())).getStructure();
String nameWithScope = structure.getNameScope().getName();
int role = smParm.getRole();
if (role == StructureMapping.ROLE_InitialVoltage) {
return smParm.getName();
} else if (role == StructureMapping.ROLE_SpecificCapacitance) {
return PARAMETER_SPECIFIC_CAPACITANCE_PREFIX + nameWithScope;
} else if (role == StructureMapping.ROLE_Size) {
if (simContext.getGeometry().getDimension() == 0) {
// if geometry is compartmental, make sure compartment sizes are set if referenced in model.
if (smParm.getExpression() == null || smParm.getExpression().isZero()) {
throw new MappingException("\nIn non-spatial application '" + getSimulationContext().getName() + "', " + "size of structure '" + structure.getName() + "' must be assigned a " + "positive value if referenced in the model.\n\nPlease go to 'Structure Mapping' tab to check the size.");
}
}
return PARAMETER_SIZE_FUNCTION_PREFIX + nameWithScope;
} else if (role == StructureMapping.ROLE_AreaPerUnitArea) {
return "AreaPerUnitArea_" + nameWithScope;
} else if (role == StructureMapping.ROLE_AreaPerUnitVolume) {
return "AreaPerUnitVolume_" + nameWithScope;
} else if (role == StructureMapping.ROLE_VolumePerUnitArea) {
return "VolumePerUnitArea_" + nameWithScope;
} else if (role == StructureMapping.ROLE_VolumePerUnitVolume) {
return "VolumePerUnitVolume_" + nameWithScope;
}
}
//
if (ste instanceof SpeciesContext) {
SpeciesContext sc = (SpeciesContext) ste;
SpeciesContextMapping scm = getSpeciesContextMapping(sc);
if (scm == null) {
throw new RuntimeException("Species '" + sc.getName() + "' is referenced in model but may have been deleted. " + "Find its references in '" + GuiConstants.DOCUMENT_EDITOR_FOLDERNAME_BIOMODEL_PARAMETERS + "'.");
}
//
if (geometryClass instanceof SubVolume) {
//
if (scm.getVariable() != null && !scm.getVariable().getName().equals(steName)) {
return scm.getVariable().getName();
}
//
// for reactions within a surface, may need "_INSIDE" or "_OUTSIDE" for jump condition
//
} else if (geometryClass instanceof SurfaceClass) {
//
// if the speciesContext is also within the surface, replace SpeciesContext name with Variable name
//
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
if (sm.getGeometryClass() == geometryClass) {
if (scm.getVariable() != null && !(scm.getVariable().getName().equals(ste.getName()))) {
return scm.getVariable().getName();
}
//
// if the speciesContext is "inside" or "outside" the membrane
//
} else if (sm.getGeometryClass() instanceof SubVolume && ((SurfaceClass) geometryClass).isAdjacentTo((SubVolume) sm.getGeometryClass())) {
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
if (!scs.isConstant()) {
if (!scs.isDiffusing() && !scs.isWellMixed()) {
throw new MappingException("Enable diffusion in Application '" + simContext.getName() + "'. This must be done for any species (e.g '" + sc.getName() + "') in flux reactions.\n\n" + "To save or run simulations, set the diffusion rate to a non-zero " + "value in Initial Conditions or disable those reactions in Specifications->Reactions.");
}
}
if (scm.getVariable() != null) {
return scm.getVariable().getName();
}
} else {
throw new MappingException("species '" + sc.getName() + "' interacts with surface '" + geometryClass.getName() + "', but is not mapped spatially adjacent");
}
}
}
return getNameScope().getSymbolName(ste);
}
Also used :
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice)
SpeciesContext(cbit.vcell.model.SpeciesContext)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
Structure(cbit.vcell.model.Structure)
ElectricalDevice(cbit.vcell.mapping.potential.ElectricalDevice)
MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice)
ReactionRule(cbit.vcell.model.ReactionRule)
RbmObservable(cbit.vcell.model.RbmObservable)
StructureSize(cbit.vcell.model.Structure.StructureSize)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
ProxyParameter(cbit.vcell.model.ProxyParameter)
MembraneVoltage(cbit.vcell.model.Membrane.MembraneVoltage)
Example 55 with SubVolume
use of cbit.vcell.geometry.SubVolume in project vcell by virtualcell.
the class RayCaster method resampleGeometry.
public static Geometry resampleGeometry(GeometryThumbnailImageFactory geometryThumbnailImageFactory, Geometry origGeometry, ISize sampleSize) throws ImageException, PropertyVetoException, GeometryException, ExpressionException {
if (origGeometry.getDimension() < 3) {
throw new GeometryException("Presently, the Raycaster resampling works only for 3d geometries.");
}
GeometrySpec origGeometrySpec = origGeometry.getGeometrySpec();
VCImage origSubvolumeImage = origGeometrySpec.getSampledImage().getCurrentValue();
if (origSubvolumeImage == null) {
throw new GeometryException("original geometry does not have a sampled image");
}
VCImage resampledSubvolumeImage = RayCaster.sampleGeometry(origGeometry, sampleSize, false);
//
// Check if resampling failed:
// if not the same number of pixelClasses (between original geometry and resampled)
// if the subvolume handles are the same
//
boolean bSameSubvolumes = true;
if (origSubvolumeImage.getNumPixelClasses() != resampledSubvolumeImage.getNumPixelClasses()) {
bSameSubvolumes = false;
}
//
for (VCPixelClass origPixelClass : origSubvolumeImage.getPixelClasses()) {
VCPixelClass resampledPixelClass = resampledSubvolumeImage.getPixelClassFromPixelValue(origPixelClass.getPixel());
if (resampledPixelClass == null) {
bSameSubvolumes = false;
break;
}
}
//
if (!bSameSubvolumes) {
StringBuffer message = new StringBuffer();
message.append("\n\nexisting geometry:\n");
for (SubVolume oldSubvolume : origGeometrySpec.getSubVolumes()) {
long count = origSubvolumeImage.countPixelsByValue((byte) oldSubvolume.getHandle());
message.append("subvolume('" + oldSubvolume.getName() + "',handle=" + oldSubvolume.getHandle() + ",numPixels=" + count + " of " + origSubvolumeImage.getNumXYZ() + "\n");
}
message.append("\n\nnew resampled handle VCImage:\n");
for (VCPixelClass newPixelClass : resampledSubvolumeImage.getPixelClasses()) {
long count = resampledSubvolumeImage.countPixelsByValue((byte) newPixelClass.getPixel());
message.append("pixelClass('" + newPixelClass.getPixelClassName() + "',pixelValue=" + newPixelClass.getPixel() + ",numPixels=" + count + " of " + resampledSubvolumeImage.getNumXYZ() + ")\n");
}
throw new GeometryException("original Geometry had " + origSubvolumeImage.getNumPixelClasses() + " subvolumes, resampled Geometry found " + resampledSubvolumeImage.getNumPixelClasses() + " subvolumes " + message.toString());
}
//
// Create new VCImage that will form the basis for a new image-based geometry.
//
VCImage newVCImage = null;
if (origGeometrySpec.getImage() != null) {
//
// was an image-based geometry - try to make new VCImage similar to the original (same pixelClass names and pixel values).
// the goal is to make identical geometries if the sample size is same as the original image size.
//
// create a new VCImage with same image pixel values (not subvolume handles) and pixel class names as original image.
//
byte[] newVCImagePixels = new byte[sampleSize.getXYZ()];
byte[] resampledSubvolumePixels = resampledSubvolumeImage.getPixels();
for (int i = 0; i < sampleSize.getXYZ(); i++) {
int subvolumeHandle = resampledSubvolumePixels[i];
ImageSubVolume imageSubvolume = (ImageSubVolume) origGeometrySpec.getSubVolume(subvolumeHandle);
newVCImagePixels[i] = (byte) imageSubvolume.getPixelValue();
}
newVCImage = new VCImageUncompressed(null, newVCImagePixels, origGeometry.getExtent(), sampleSize.getX(), sampleSize.getY(), sampleSize.getZ());
newVCImage.setName(origGeometrySpec.getImage().getName());
ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
for (VCPixelClass origPixelClass : origGeometrySpec.getImage().getPixelClasses()) {
SubVolume origSubvolume = origGeometrySpec.getImageSubVolumeFromPixelValue(origPixelClass.getPixel());
newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origPixelClass.getPixel()));
}
newVCImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
} else {
//
// was an analytic geometry - create a new image-based geometry
//
// create a new VCImage with image pixel values and pixelClass names equal to corresponding subvolumes from original geometry
// make new subvolume names equal to old subvolume names.
//
byte[] newVCImageSubvolumePixels = resampledSubvolumeImage.getPixels().clone();
newVCImage = new VCImageUncompressed(null, newVCImageSubvolumePixels, origGeometry.getExtent(), sampleSize.getX(), sampleSize.getY(), sampleSize.getZ());
ArrayList<VCPixelClass> newPixelClasses = new ArrayList<VCPixelClass>();
for (SubVolume origSubvolume : origGeometrySpec.getSubVolumes()) {
// note: newVCImage already has subvolume handle pixels.
newPixelClasses.add(new VCPixelClass(null, origSubvolume.getName(), origSubvolume.getHandle()));
}
newVCImage.setPixelClasses(newPixelClasses.toArray(new VCPixelClass[newPixelClasses.size()]));
}
//
// construct the new geometry with the sampled VCImage.
//
Geometry newGeometry = new Geometry(origGeometry.getName(), newVCImage);
newGeometry.getGeometrySpec().setExtent(origGeometry.getExtent());
newGeometry.getGeometrySpec().setOrigin(origGeometry.getOrigin());
newGeometry.setDescription(origGeometry.getDescription());
newGeometry.getGeometrySurfaceDescription().setFilterCutoffFrequency(origGeometry.getGeometrySurfaceDescription().getFilterCutoffFrequency());
newGeometry.precomputeAll(geometryThumbnailImageFactory, true, true);
return newGeometry;
}
Also used :
VCPixelClass(cbit.image.VCPixelClass)
GeometryException(cbit.vcell.geometry.GeometryException)
ArrayList(java.util.ArrayList)
VCImage(cbit.image.VCImage)
VCImageUncompressed(cbit.image.VCImageUncompressed)
GeometrySpec(cbit.vcell.geometry.GeometrySpec)
Geometry(cbit.vcell.geometry.Geometry)
SubVolume(cbit.vcell.geometry.SubVolume)
ImageSubVolume(cbit.vcell.geometry.ImageSubVolume)
ImageSubVolume(cbit.vcell.geometry.ImageSubVolume)
Aggregations
SubVolume (cbit.vcell.geometry.SubVolume)84
SurfaceClass (cbit.vcell.geometry.SurfaceClass)42
Expression (cbit.vcell.parser.Expression)33
AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)26
ImageSubVolume (cbit.vcell.geometry.ImageSubVolume)22
GeometryClass (cbit.vcell.geometry.GeometryClass)21
ArrayList (java.util.ArrayList)21
Geometry (cbit.vcell.geometry.Geometry)20
MathDescription (cbit.vcell.math.MathDescription)17
SpeciesContext (cbit.vcell.model.SpeciesContext)17
PropertyVetoException (java.beans.PropertyVetoException)17
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)16
CompartmentSubVolume (cbit.vcell.geometry.CompartmentSubVolume)15
Model (cbit.vcell.model.Model)15
ImageException (cbit.image.ImageException)14
VCImage (cbit.image.VCImage)14
Feature (cbit.vcell.model.Feature)14
Structure (cbit.vcell.model.Structure)14
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)13
ExpressionException (cbit.vcell.parser.ExpressionException)13
