use of cbit.vcell.math.Variable.Domain in project vcell by virtualcell.
the class DiffEquMathMapping method addSpatialProcesses.
private void addSpatialProcesses(VariableHash varHash, ArrayList<CompartmentSubdomainContext> compartmentSubdomainContexts, ArrayList<MembraneSubdomainContext> membraneSubdomainContexts) throws MathException, MappingException, ExpressionException {
if (simContext.getGeometry().getDimension() == 0) {
return;
}
//
for (SpatialObject spatialObject : simContext.getSpatialObjects()) {
if (spatialObject instanceof PointObject) {
PointObject pointObject = (PointObject) spatialObject;
//
// if true, have to solve for this category
//
boolean bPosition = pointObject.isQuantityCategoryEnabled(QuantityCategory.PointPosition);
boolean bVelocity = pointObject.isQuantityCategoryEnabled(QuantityCategory.PointVelocity);
boolean bDirection = pointObject.isQuantityCategoryEnabled(QuantityCategory.DirectionToPoint);
boolean bDistance = pointObject.isQuantityCategoryEnabled(QuantityCategory.PointDistanceMap);
//
// either make a point subdomain, or just define functions.
//
ArrayList<PointLocation> pointLocationProcesses = new ArrayList<PointLocation>();
ArrayList<PointKinematics> pointKinematicsProcesses = new ArrayList<PointKinematics>();
for (SpatialProcess spatialProcess : simContext.getSpatialProcesses()) {
if (spatialProcess instanceof PointLocation && ((PointLocation) spatialProcess).getPointObject() == pointObject) {
pointLocationProcesses.add((PointLocation) spatialProcess);
}
if (spatialProcess instanceof PointKinematics && ((PointKinematics) spatialProcess).getPointObject() == pointObject) {
pointKinematicsProcesses.add((PointKinematics) spatialProcess);
}
}
if (pointLocationProcesses.size() == 1 && pointKinematicsProcesses.size() == 0 && !bVelocity) {
GeometryClass gc = null;
PointLocation pointLocation = pointLocationProcesses.get(0);
if (bPosition) {
if (simContext.getGeometry().getDimension() == 1) {
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
LocalParameter posXParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXParam, gc), getIdentifierSubstitutions(posXParam.getExpression(), posXParam.getUnitDefinition(), gc), gc));
Expression posXExp = new Expression(posXParam, pointLocation.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXQuantity, gc), getIdentifierSubstitutions(posXExp, posXQuantity.getUnitDefinition(), gc), gc));
Expression posX_minus_X = Expression.add(posXExp, Expression.negate(xExp));
Expression signum_posX_minus_X = Expression.function(FunctionType.MIN, new Expression(1.0), Expression.function(FunctionType.MAX, Expression.mult(new Expression(1e10), posX_minus_X), new Expression(-1)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(signum_posX_minus_X, dirXQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
Expression abs_X_minux_posX = Expression.function(FunctionType.ABS, posX_minus_X);
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(abs_X_minux_posX, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 2) {
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
LocalParameter posXParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXParam, gc), getIdentifierSubstitutions(posXParam.getExpression(), posXParam.getUnitDefinition(), gc), gc));
Expression posXExp = new Expression(posXParam, pointLocation.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXQuantity, gc), getIdentifierSubstitutions(new Expression(posXParam, pointLocation.getNameScope()), posXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity posYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Y);
LocalParameter posYParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYParam, gc), getIdentifierSubstitutions(posYParam.getExpression(), posYParam.getUnitDefinition(), gc), gc));
Expression posYExp = new Expression(posYParam, pointLocation.getNameScope());
Expression yExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Y), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYQuantity, gc), getIdentifierSubstitutions(new Expression(posYParam, pointLocation.getNameScope()), posYQuantity.getUnitDefinition(), gc), gc));
Expression posX_minux_X = Expression.add(posXExp, Expression.negate(xExp));
Expression posY_minux_Y = Expression.add(posYExp, Expression.negate(yExp));
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression sqrt_DX2_DY2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 3) {
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
LocalParameter posXParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXParam, gc), getIdentifierSubstitutions(posXParam.getExpression(), posXParam.getUnitDefinition(), gc), gc));
Expression posXExp = new Expression(posXParam, pointLocation.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posXQuantity, gc), getIdentifierSubstitutions(new Expression(posXParam, pointLocation.getNameScope()), posXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity posYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Y);
LocalParameter posYParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYParam, gc), getIdentifierSubstitutions(posYParam.getExpression(), posYParam.getUnitDefinition(), gc), gc));
Expression posYExp = new Expression(posYParam, pointLocation.getNameScope());
Expression yExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Y), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posYQuantity, gc), getIdentifierSubstitutions(new Expression(posYParam, pointLocation.getNameScope()), posYQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity posZQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Z);
LocalParameter posZParam = pointLocation.getParameter(SpatialProcessParameterType.PointPositionZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posZParam, gc), getIdentifierSubstitutions(posZParam.getExpression(), posZParam.getUnitDefinition(), gc), gc));
Expression posZExp = new Expression(posZParam, pointLocation.getNameScope());
Expression zExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Z), simContext.getModel().getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(posZQuantity, gc), getIdentifierSubstitutions(new Expression(posZParam, pointLocation.getNameScope()), posZQuantity.getUnitDefinition(), gc), gc));
Expression posX_minux_X = Expression.add(posXExp, Expression.negate(xExp));
Expression posY_minux_Y = Expression.add(posYExp, Expression.negate(yExp));
Expression posZ_minux_Z = Expression.add(posZExp, Expression.negate(zExp));
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression DZ2 = Expression.mult(posZ_minux_Z, posZ_minux_Z);
Expression sqrt_DX2_DY2_DZ2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2, DZ2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirZ = Expression.div(posZ_minux_Z, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirZQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Z);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirZQuantity, gc), getIdentifierSubstitutions(dirZ, dirZQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2_DZ2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
} else {
throw new MappingException("PointLocation process defined for pointObject '" + pointObject.getName() + "' but Position not enabled");
}
} else if (pointLocationProcesses.size() == 0 && pointKinematicsProcesses.size() == 1) {
GeometryClass gc = null;
PointKinematics pointKinematics = pointKinematicsProcesses.get(0);
if (bPosition && bVelocity) {
LocalParameter velXParam = pointKinematics.getParameter(SpatialProcessParameterType.PointVelocityX);
LocalParameter iniPosXParam = pointKinematics.getParameter(SpatialProcessParameterType.PointInitialPositionX);
SpatialQuantity posXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.X);
SpatialQuantity velXQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointVelocity, QuantityComponent.X);
Expression posXExp = new Expression(posXQuantity, simContext.getNameScope());
Expression xExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.X), simContext.getModel().getNameScope());
Expression posX_minux_X = Expression.add(posXExp, Expression.negate(xExp));
LocalParameter velYParam = pointKinematics.getParameter(SpatialProcessParameterType.PointVelocityY);
LocalParameter iniPosYParam = pointKinematics.getParameter(SpatialProcessParameterType.PointInitialPositionY);
SpatialQuantity posYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Y);
SpatialQuantity velYQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointVelocity, QuantityComponent.Y);
Expression posYExp = new Expression(posYQuantity, simContext.getNameScope());
Expression yExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Y), simContext.getModel().getNameScope());
Expression posY_minux_Y = Expression.add(posYExp, Expression.negate(yExp));
LocalParameter velZParam = pointKinematics.getParameter(SpatialProcessParameterType.PointVelocityZ);
LocalParameter iniPosZParam = pointKinematics.getParameter(SpatialProcessParameterType.PointInitialPositionZ);
SpatialQuantity posZQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointPosition, QuantityComponent.Z);
SpatialQuantity velZQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointVelocity, QuantityComponent.Z);
Expression posZExp = new Expression(posZQuantity, simContext.getNameScope());
Expression zExp = new Expression(simContext.getModel().getReservedSymbolByRole(ReservedSymbolRole.Z), simContext.getModel().getNameScope());
Expression posZ_minux_Z = Expression.add(posZExp, Expression.negate(zExp));
String pointSubdomainName = pointObject.getName();
Domain domain = new Domain(pointSubdomainName);
PointSubDomain pointSubdomain = new PointSubDomain(pointSubdomainName);
mathDesc.addSubDomain(pointSubdomain);
if (simContext.getGeometry().getDimension() >= 1) {
PointVariable posXVar = new PointVariable(getMathSymbol(posXQuantity, gc), domain);
varHash.addVariable(posXVar);
Expression initXExp = getIdentifierSubstitutions(new Expression(iniPosXParam, pointKinematics.getNameScope()), iniPosXParam.getUnitDefinition(), gc);
Expression rateXExp = getIdentifierSubstitutions(new Expression(velXParam, pointKinematics.getNameScope()), velXParam.getUnitDefinition(), gc);
OdeEquation odeX = new OdeEquation(posXVar, initXExp, rateXExp);
pointSubdomain.addEquation(odeX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(iniPosXParam, gc), getIdentifierSubstitutions(new Expression(iniPosXParam.getExpression()), iniPosXParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXParam, gc), getIdentifierSubstitutions(new Expression(velXParam.getExpression()), velXParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXQuantity, gc), getIdentifierSubstitutions(new Expression(velXParam, pointKinematics.getNameScope()), velXQuantity.getUnitDefinition(), gc), gc));
pointSubdomain.setPositionX(getIdentifierSubstitutions(new Expression(posXQuantity, simContext.getNameScope()), posXQuantity.getUnitDefinition(), gc));
}
if (simContext.getGeometry().getDimension() >= 2) {
PointVariable posYVar = new PointVariable(getMathSymbol(posYQuantity, gc), domain);
varHash.addVariable(posYVar);
Expression initYExp = getIdentifierSubstitutions(new Expression(iniPosYParam, pointKinematics.getNameScope()), iniPosYParam.getUnitDefinition(), gc);
Expression rateYExp = getIdentifierSubstitutions(new Expression(velYParam, pointKinematics.getNameScope()), velYParam.getUnitDefinition(), gc);
OdeEquation odeY = new OdeEquation(posYVar, initYExp, rateYExp);
pointSubdomain.addEquation(odeY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(iniPosYParam, gc), getIdentifierSubstitutions(new Expression(iniPosYParam.getExpression()), iniPosYParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYParam, gc), getIdentifierSubstitutions(new Expression(velYParam.getExpression()), velYParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYQuantity, gc), getIdentifierSubstitutions(new Expression(velYParam, pointKinematics.getNameScope()), velYQuantity.getUnitDefinition(), gc), gc));
pointSubdomain.setPositionY(getIdentifierSubstitutions(new Expression(posYQuantity, simContext.getNameScope()), posYQuantity.getUnitDefinition(), gc));
}
if (simContext.getGeometry().getDimension() == 3) {
PointVariable posZVar = new PointVariable(getMathSymbol(posZQuantity, gc), domain);
varHash.addVariable(posZVar);
Expression initZExp = getIdentifierSubstitutions(new Expression(iniPosZParam, pointKinematics.getNameScope()), iniPosZParam.getUnitDefinition(), gc);
Expression rateZExp = getIdentifierSubstitutions(new Expression(velZParam, pointKinematics.getNameScope()), velZParam.getUnitDefinition(), gc);
OdeEquation odeZ = new OdeEquation(posZVar, initZExp, rateZExp);
pointSubdomain.addEquation(odeZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(iniPosZParam, gc), getIdentifierSubstitutions(new Expression(iniPosZParam.getExpression()), iniPosZParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZParam, gc), getIdentifierSubstitutions(new Expression(velZParam.getExpression()), velZParam.getUnitDefinition(), gc), gc));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZQuantity, gc), getIdentifierSubstitutions(new Expression(velZParam, pointKinematics.getNameScope()), velZQuantity.getUnitDefinition(), gc), gc));
pointSubdomain.setPositionZ(getIdentifierSubstitutions(new Expression(posZQuantity, simContext.getNameScope()), posZQuantity.getUnitDefinition(), gc));
}
if (simContext.getGeometry().getDimension() == 1) {
Expression signum_posX_minus_X = Expression.function(FunctionType.MIN, new Expression(1.0), Expression.function(FunctionType.MAX, Expression.mult(new Expression(1e10), posX_minux_X), new Expression(-1)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(signum_posX_minus_X, dirXQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
Expression abs_X_minux_posX = Expression.function(FunctionType.ABS, posX_minux_X);
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(abs_X_minux_posX, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 2) {
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression sqrt_DX2_DY2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
if (simContext.getGeometry().getDimension() == 3) {
Expression DX2 = Expression.mult(posX_minux_X, posX_minux_X);
Expression DY2 = Expression.mult(posY_minux_Y, posY_minux_Y);
Expression DZ2 = Expression.mult(posZ_minux_Z, posZ_minux_Z);
Expression sqrt_DX2_DY2_DZ2 = Expression.function(FunctionType.SQRT, Expression.add(DX2, DY2, DZ2));
Expression dirX = Expression.div(posX_minux_X, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirY = Expression.div(posY_minux_Y, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
Expression dirZ = Expression.div(posZ_minux_Z, Expression.add(sqrt_DX2_DY2_DZ2, new Expression(1e-8)));
if (bDirection) {
SpatialQuantity dirXQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.X);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirXQuantity, gc), getIdentifierSubstitutions(dirX, dirXQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirYQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Y);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirYQuantity, gc), getIdentifierSubstitutions(dirY, dirYQuantity.getUnitDefinition(), gc), gc));
SpatialQuantity dirZQuantity = pointObject.getSpatialQuantity(QuantityCategory.DirectionToPoint, QuantityComponent.Z);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dirZQuantity, gc), getIdentifierSubstitutions(dirZ, dirZQuantity.getUnitDefinition(), gc), gc));
}
if (bDistance) {
SpatialQuantity distanceQuantity = pointObject.getSpatialQuantity(QuantityCategory.PointDistanceMap, QuantityComponent.Scalar);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(distanceQuantity, gc), getIdentifierSubstitutions(sqrt_DX2_DY2_DZ2, distanceQuantity.getUnitDefinition(), gc), gc));
}
}
} else {
throw new MappingException(pointKinematics.getDescription() + " process defined for pointObject '" + pointObject.getName() + "' but Position and Velocity not enabled");
}
} else {
throw new MappingException("expecting 1 location or kinematics process for point '" + pointObject.getName() + "'");
}
} else if (spatialObject instanceof SurfaceRegionObject) {
SurfaceRegionObject surfaceRegionObject = (SurfaceRegionObject) spatialObject;
SubVolume insideSubvolume = surfaceRegionObject.getInsideSubVolume();
SubVolume outsideSubvolume = surfaceRegionObject.getOutsideSubVolume();
SurfaceClass surfaceClass = simContext.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubvolume, outsideSubvolume);
MembraneSubdomainContext memSubdomainContext = null;
for (MembraneSubdomainContext context : membraneSubdomainContexts) {
if (context.surfaceClass == surfaceClass) {
memSubdomainContext = context;
}
}
//
// if true, have to solve for this category
//
boolean bNormal = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.Normal);
boolean bVelocity = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceVelocity);
boolean bDistance = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceDistanceMap);
boolean bDirection = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.DirectionToSurface);
boolean bSize = surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceSize);
if (bVelocity) {
ArrayList<SurfaceKinematics> surfaceKinematicsList = new ArrayList<SurfaceKinematics>();
for (SpatialProcess spatialProcess : simContext.getSpatialProcesses()) {
if (spatialProcess instanceof SurfaceKinematics && ((SurfaceKinematics) spatialProcess).getSurfaceRegionObject() == surfaceRegionObject) {
surfaceKinematicsList.add((SurfaceKinematics) spatialProcess);
}
}
if (surfaceKinematicsList.size() == 1) {
SurfaceKinematics surfaceKinematics = surfaceKinematicsList.get(0);
if (simContext.getGeometry().getDimension() >= 1) {
SpatialQuantity velXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.X);
LocalParameter velXParam = surfaceKinematics.getParameter(SpatialProcessParameterType.SurfaceVelocityX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXParam, surfaceClass), getIdentifierSubstitutions(velXParam.getExpression(), velXParam.getUnitDefinition(), surfaceClass), surfaceClass));
Expression velXExp = new Expression(velXParam, surfaceKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXQuantity, surfaceClass), getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
if (bNormal) {
SpatialQuantity normXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.Normal, QuantityComponent.X);
Expression normXExp = new Expression(MathFunctionDefinitions.FUNCTION_normalX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(normXQuantity, surfaceClass), getIdentifierSubstitutions(normXExp, normXQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
}
if (simContext.getGeometry().getDimension() >= 2) {
SpatialQuantity velYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Y);
LocalParameter velYParam = surfaceKinematics.getParameter(SpatialProcessParameterType.SurfaceVelocityY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYParam, surfaceClass), getIdentifierSubstitutions(velYParam.getExpression(), velYParam.getUnitDefinition(), surfaceClass), surfaceClass));
Expression velYExp = new Expression(velYParam, surfaceKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYQuantity, surfaceClass), getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
if (bNormal) {
SpatialQuantity normYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.Normal, QuantityComponent.Y);
Expression normYExp = new Expression(MathFunctionDefinitions.FUNCTION_normalY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(normYQuantity, surfaceClass), getIdentifierSubstitutions(normYExp, normYQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
}
if (simContext.getGeometry().getDimension() >= 3) {
SpatialQuantity velZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Z);
LocalParameter velZParam = surfaceKinematics.getParameter(SpatialProcessParameterType.SurfaceVelocityZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZParam, surfaceClass), getIdentifierSubstitutions(velZParam.getExpression(), velZParam.getUnitDefinition(), surfaceClass), surfaceClass));
Expression velYExp = new Expression(velZParam, surfaceKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZQuantity, surfaceClass), getIdentifierSubstitutions(velYExp, velZQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
if (bNormal) {
SpatialQuantity normZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.Normal, QuantityComponent.Z);
// MathFunctionDefinitions.FUNCTION_normalZ);
Expression normZExp = new Expression("normalZ_not_implemented()");
varHash.addVariable(newFunctionOrConstant(getMathSymbol(normZQuantity, surfaceClass), getIdentifierSubstitutions(normZExp, normZQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
}
} else {
throw new MappingException("expecting 1 Surface Kinematics process for Surface Object '" + surfaceRegionObject.getName() + "'");
}
}
if (bSize) {
SpatialQuantity sizeQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceSize, QuantityComponent.Scalar);
String funcName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
Expression sizeExp = Expression.function(funcName, new Expression[] { new Expression("'" + surfaceClass.getName() + "'") });
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeQuantity, surfaceClass), getIdentifierSubstitutions(sizeExp, sizeQuantity.getUnitDefinition(), surfaceClass), surfaceClass));
}
if (bDirection) {
SpatialQuantity directionXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.DirectionToSurface, QuantityComponent.X);
SpatialQuantity directionYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.DirectionToSurface, QuantityComponent.Y);
SpatialQuantity directionZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.DirectionToSurface, QuantityComponent.Z);
for (SubVolume adjacentSubvolume : surfaceClass.getAdjacentSubvolumes()) {
Domain domain = new Domain(adjacentSubvolume);
CompartmentSubDomain compSubdomain = mathDesc.getCompartmentSubDomain(adjacentSubvolume.getName());
if (simContext.getGeometry().getDimension() >= 1) {
String name = adjacentSubvolume.getName() + "_dirX_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDirectionToMembraneQuantity localDirectionXQuantity = addLocalizedDirectionToMembraneQuantity(name, surfaceClass, adjacentSubvolume, QuantityComponent.X);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDirectionXQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.directionToMembraneX);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
if (simContext.getGeometry().getDimension() >= 2) {
String name = adjacentSubvolume.getName() + "_dirY_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDirectionToMembraneQuantity localDirectionYQuantity = addLocalizedDirectionToMembraneQuantity(name, surfaceClass, adjacentSubvolume, QuantityComponent.Y);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDirectionYQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.directionToMembraneY);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
if (simContext.getGeometry().getDimension() == 3) {
String name = adjacentSubvolume.getName() + "_dirZ_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDirectionToMembraneQuantity localDirectionZQuantity = addLocalizedDirectionToMembraneQuantity(name, surfaceClass, adjacentSubvolume, QuantityComponent.Z);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDirectionZQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.directionToMembraneZ);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
}
}
if (bDistance) {
SpatialQuantity distanceQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceDistanceMap, QuantityComponent.Scalar);
for (SubVolume adjacentSubvolume : surfaceClass.getAdjacentSubvolumes()) {
String name = adjacentSubvolume.getName() + "_distanceTo_" + memSubdomainContext.membraneSubdomain.getName();
LocalizedDistanceToMembraneQuantity localDistanceQuantity = addLocalizedDistanceToMembraneQuantity(name, surfaceClass, adjacentSubvolume);
CompartmentSubDomain compSubdomain = mathDesc.getCompartmentSubDomain(adjacentSubvolume.getName());
Domain domain = new Domain(adjacentSubvolume);
VolVariable distanceVar = new VolVariable(getMathSymbol(localDistanceQuantity, adjacentSubvolume), domain);
varHash.addVariable(distanceVar);
ComputeMembraneMetricEquation membraneMetricEquation = new ComputeMembraneMetricEquation(distanceVar, MembraneMetricComponent.distanceToMembrane);
membraneMetricEquation.setTargetMembraneName(memSubdomainContext.membraneSubdomain.getName());
compSubdomain.addEquation(membraneMetricEquation);
}
}
} else if (spatialObject instanceof VolumeRegionObject) {
VolumeRegionObject volumeRegionObject = (VolumeRegionObject) spatialObject;
SubVolume subvolume = volumeRegionObject.getSubVolume();
boolean bCentroid = volumeRegionObject.isQuantityCategoryEnabled(QuantityCategory.DirectionToSurface);
boolean bSize = volumeRegionObject.isQuantityCategoryEnabled(QuantityCategory.VolumeSize);
boolean bVelocity = volumeRegionObject.isQuantityCategoryEnabled(QuantityCategory.InteriorVelocity);
if (bVelocity) {
ArrayList<VolumeKinematics> volumeKinematicsList = new ArrayList<VolumeKinematics>();
for (SpatialProcess spatialProcess : simContext.getSpatialProcesses()) {
if (spatialProcess instanceof VolumeKinematics && ((VolumeKinematics) spatialProcess).getVolumeRegionObject() == volumeRegionObject) {
volumeKinematicsList.add((VolumeKinematics) spatialProcess);
}
}
if (volumeKinematicsList.size() == 1) {
VolumeKinematics volumeKinematics = volumeKinematicsList.get(0);
if (simContext.getGeometry().getDimension() >= 1) {
SpatialQuantity velXQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.InteriorVelocity, QuantityComponent.X);
LocalParameter velXParam = volumeKinematics.getParameter(SpatialProcessParameterType.InternalVelocityX);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXParam, subvolume), getIdentifierSubstitutions(velXParam.getExpression(), velXParam.getUnitDefinition(), subvolume), subvolume));
Expression velXExp = new Expression(velXParam, volumeKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velXQuantity, subvolume), getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), subvolume), subvolume));
}
if (simContext.getGeometry().getDimension() >= 2) {
SpatialQuantity velYQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.InteriorVelocity, QuantityComponent.Y);
LocalParameter velYParam = volumeKinematics.getParameter(SpatialProcessParameterType.InternalVelocityY);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYParam, subvolume), getIdentifierSubstitutions(velYParam.getExpression(), velYParam.getUnitDefinition(), subvolume), subvolume));
Expression velYExp = new Expression(velYParam, volumeKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velYQuantity, subvolume), getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), subvolume), subvolume));
}
if (simContext.getGeometry().getDimension() == 3) {
SpatialQuantity velZQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.InteriorVelocity, QuantityComponent.Z);
LocalParameter velZParam = volumeKinematics.getParameter(SpatialProcessParameterType.InternalVelocityZ);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZParam, subvolume), getIdentifierSubstitutions(velZParam.getExpression(), velZParam.getUnitDefinition(), subvolume), subvolume));
Expression velZExp = new Expression(velZParam, volumeKinematics.getNameScope());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(velZQuantity, subvolume), getIdentifierSubstitutions(velZExp, velZQuantity.getUnitDefinition(), subvolume), subvolume));
}
} else {
throw new MappingException("expecting 1 Volume Kinematics process for Volume Object '" + volumeRegionObject.getName() + "'");
}
}
if (bSize) {
SpatialQuantity sizeQuantity = volumeRegionObject.getSpatialQuantity(QuantityCategory.VolumeSize, QuantityComponent.Scalar);
String funcName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
Expression sizeExp = Expression.function(funcName, new Expression[] { new Expression("'" + subvolume.getName() + "'") });
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeQuantity, subvolume), getIdentifierSubstitutions(sizeExp, sizeQuantity.getUnitDefinition(), subvolume), subvolume));
}
if (bCentroid) {
throw new MappingException(QuantityCategory.Centroid.description + " not yet implemented in math generation");
}
}
}
}
use of cbit.vcell.math.Variable.Domain in project vcell by virtualcell.
the class FastSystemAnalyzer method refreshSubstitutedRateExps.
/**
*/
private void refreshSubstitutedRateExps() throws MathException, ExpressionException {
//
// refresh PseudoConstants (temp constants involved with fastInvariants)
//
pseudoConstantList.removeAllElements();
Enumeration<FastInvariant> enum_fi = fastSystem.getFastInvariants();
while (enum_fi.hasMoreElements()) {
FastInvariant fi = enum_fi.nextElement();
Domain domain = new Domain(fastSystem.getSubDomain());
fi.getFunction().bindExpression(this);
PseudoConstant pc = new PseudoConstant(getAvailablePseudoConstantName(), fi.getFunction(), domain);
pseudoConstantList.addElement(pc);
// System.out.println("FastSystem.refreshSubstitutedRateExps() __C"+i+" = "+fi.getFunction());
}
//
// build expressions for dependent variables in terms of independent vars
// and pseudoConstants
//
dependencyExpList.removeAllElements();
for (int row = 0; row < dependentVarList.size(); row++) {
Expression exp = new Expression("0.0;");
//
for (int col = 0; col < independentVarList.size(); col++) {
Variable indepVar = (Variable) independentVarList.elementAt(col);
RationalExp coefExp = dependencyMatrix.get(row, col).simplify();
if (!coefExp.isZero()) {
exp = Expression.add(exp, new Expression(coefExp.infixString() + "*" + indepVar.getName()));
}
}
//
for (int col = independentVarList.size(); col < dependencyMatrix.getNumCols(); col++) {
PseudoConstant pc = (PseudoConstant) pseudoConstantList.elementAt(col - independentVarList.size());
RationalExp coefExp = dependencyMatrix.get(row, col);
if (!coefExp.isZero()) {
exp = Expression.add(exp, new Expression(coefExp.infixString() + "*" + pc.getName()));
}
}
exp.bindExpression(null);
exp = exp.flatten();
exp.bindExpression(this);
// System.out.println("FastSystem.refreshSubstitutedRateExps() "+((Variable)dependentVarList.elementAt(row)).getName()+" = "+exp.toString()+";");
dependencyExpList.addElement(exp);
}
//
// flatten functions, then substitute expressions for dependent vars into rate expressions
//
fastRateExpList.removeAllElements();
// VariableSymbolTable combinedSymbolTable = getCombinedSymbolTable();
Enumeration<FastRate> enum_fr = fastSystem.getFastRates();
while (enum_fr.hasMoreElements()) {
FastRate fr = enum_fr.nextElement();
Expression exp = new Expression(MathUtilities.substituteFunctions(new Expression(fr.getFunction()), this));
// System.out.println("FastSystem.refreshSubstitutedRateExps() fast rate before substitution = "+exp.toString());
for (int j = 0; j < dependentVarList.size(); j++) {
Variable depVar = (Variable) dependentVarList.elementAt(j);
Expression subExp = new Expression((Expression) dependencyExpList.elementAt(j));
exp.substituteInPlace(new Expression(depVar.getName()), subExp);
}
exp.bindExpression(null);
exp = exp.flatten();
// System.out.println("FastSystem.refreshSubstitutedRateExps() fast rate after substitution = "+exp.toString());
exp.bindExpression(this);
fastRateExpList.addElement(exp);
}
}
use of cbit.vcell.math.Variable.Domain in project vcell by virtualcell.
the class ParticleMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
getSimulationContext().checkValidity();
if (getSimulationContext().getGeometry().getDimension() == 0) {
throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
}
StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = getSimulationContext().getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
if (reactionSpecs[i].isFast()) {
throw new MappingException("fast reactions not supported for particle models");
}
rsList.add(reactionSpecs[i].getReactionStep());
}
}
ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
rsList.copyInto(reactionSteps);
//
for (int i = 0; i < reactionSteps.length; i++) {
Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
// //
// // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
// //
// Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
// for (int i = 0; i < structures.length; i++){
// StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
// if (sm==null || (sm.getGeometryClass() == null)){
// throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
// }
// if (sm.getUnitSizeParameter()!=null){
// Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
// if(unitSizeExp != null)
// {
// try {
// double unitSize = unitSizeExp.evaluateConstant();
// if (unitSize != 1.0){
// throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
// }
// }catch (ExpressionException e){
// e.printStackTrace(System.out);
// throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
// }
// }
// }
// }
// {
// GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
// for (int i = 0; i < geometryClass.length; i++){
// Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
// }
// }
// }
// deals with model parameters
Model model = getSimulationContext().getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
ModelParameter[] modelParameters = model.getModelParameters();
// populate in globalParameterVariants hashtable
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = getSimulationContext().getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(getSimulationContext().getName() + "_generated");
}
//
// volume particle variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof ParticleVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
GeometryClass geometryClass = membraneMapping.getGeometryClass();
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
// Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = null;
Expression initExpr = null;
if (getSimulationContext().isUsingConcentration()) {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
initExpr = new Expression(initParm.getExpression());
// if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
// initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
// }
} else {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
initExpr = new Expression(initParm.getExpression());
}
if (initExpr != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = getSimulationContext().getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// constant species (either function or constant)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof Constant) {
varHash.addVariable(scm.getVariable());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (getSimulationContext().getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (getSimulationContext().getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(getSimulationContext().getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// create subdomains (volume and surfaces)
//
GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
for (int k = 0; k < geometryClasses.length; k++) {
if (geometryClasses[k] instanceof SubVolume) {
SubVolume subVolume = (SubVolume) geometryClasses[k];
//
// get priority of subDomain
//
// now does not have to match spatial feature, *BUT* needs to be unique
int priority = k;
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
if (getSimulationContext().getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (getSimulationContext().getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
} else if (geometryClasses[k] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
// determine membrane inside and outside subvolume
// this preserves backward compatibility so that membrane subdomain
// inside and outside correspond to structure hierarchy when present
Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
SubVolume innerSubVolume = ret.one;
SubVolume outerSubVolume = ret.two;
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
}
}
//
// create Particle Contexts for all Particle Variables
//
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
//
// initial distribution of particles
//
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
ParticleInitialCondition pic = null;
if (getSimulationContext().isUsingConcentration()) {
Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
if (particleVariable instanceof VolumeParticleVariable) {
initialDistribution = Expression.mult(initialDistribution, unitFactor);
}
pic = new ParticleInitialConditionConcentration(initialDistribution);
} else {
Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
if (initialCount == null) {
throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
}
Expression locationX = new Expression("u");
Expression locationY = new Expression("u");
Expression locationZ = new Expression("u");
pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
}
particleInitialConditions.add(pic);
//
// diffusion
//
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
Expression driftXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velX_quantities.length == 1 && numRegions == 1) {
driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velY_quantities.length == 1 && numRegions == 1) {
driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
GeometryClass myGC = sm.getGeometryClass();
if (myGC == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
subDomain.addParticleProperties(particleProperties);
}
}
for (ReactionStep reactionStep : reactionSteps) {
Kinetics kinetics = reactionStep.getKinetics();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
KineticsParameter reactionRateParameter = null;
if (kinetics instanceof LumpedKinetics) {
reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
} else {
reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
}
// macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
if (radiusExp != null) {
Expression expCopy = new Expression(radiusExp);
try {
MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ": '" + radiusExp.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
}
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
if (rp instanceof Reactant) {
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
}
}
} else if (rp instanceof Product) {
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createCreateAction(particle));
}
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
} else // other type of reactions
{
/* check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
rate constants are important in calculating the probability rate.
for Mass Action, we use KForward and KReverse,
for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
*/
Expression forwardRate = null;
Expression reverseRate = null;
// Using the MassActionFunction to write out the math description
MassActionSolver.MassActionFunction maFunc = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
} else {
if (maFunc.getForwardRate() != null) {
forwardRate = maFunc.getForwardRate();
}
if (maFunc.getReverseRate() != null) {
reverseRate = maFunc.getReverseRate();
}
}
}
if (maFunc != null) {
// if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
// we process it as forward reaction
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
List<Action> reverseActions = new ArrayList<Action>();
List<ReactionParticipant> reactants = maFunc.getReactants();
List<ReactionParticipant> products = maFunc.getProducts();
for (ReactionParticipant rp : reactants) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createCreateAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
for (ReactionParticipant rp : products) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createCreateAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createDestroyAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
//
// There are two unit conversions required:
//
// 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
// (where dim is 2 for membrane reactions and 3 for volume reactions)
//
// for forward rates:
// 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane reactants and 3 for volume reactants)
//
// or
//
// for reverse rates:
// 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane products and 3 for volume products)
//
RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
if (forwardRate != null) {
VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant reactant : maFunc.getReactants()) {
VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynReactantUnit = null;
if (bForceContinuous) {
// reactant is continuous (vcell units)
smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
} else {
// reactant is a particle (smoldyn units)
RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
}
// keep track of units of all reactants
smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
// accumulate unit factors for all reactants
forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
}
forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
// create particle jump process
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
}
// end of forward rate not null
if (reverseRate != null) {
VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant product : maFunc.getProducts()) {
VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynProductUnit = null;
if (bForceContinuous) {
smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
} else {
RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
}
// keep track of units of all products
smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
// accumulate unit factors for all products
reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
}
reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (reverseActions.size() > 0) {
ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(reverseProcess);
}
}
// end of reverse rate not null
}
// end of maFunc not null
}
// end of reaction step for loop
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
lg.warn(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
if (lg.isDebugEnabled()) {
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
use of cbit.vcell.math.Variable.Domain in project vcell by virtualcell.
the class RulebasedMathMapping method refreshVariables.
/**
* This method was created in VisualAge.
* @Override
*/
@Override
protected void refreshVariables() throws MappingException {
Domain defaultDomain = new Domain(this.getSimulationContext().getGeometry().getGeometryClasses()[0]);
//
// non-constant independent variables require either a membrane or volume variable
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(scm.getSpeciesContext());
// stochastic variable is always a function of size.
SpeciesCountParameter spCountParm = null;
try {
String countName = scs.getSpeciesContext().getName() + BIO_PARAM_SUFFIX_SPECIES_COUNT;
Expression countExp = new Expression(0.0);
spCountParm = addSpeciesCountParameter(countName, countExp, PARAMETER_ROLE_SPECIES_COUNT, scs.getInitialCountParameter().getUnitDefinition(), scs.getSpeciesContext());
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add concentration of species as MathMappingParameter - this will map to species concentration function
try {
String concName = scs.getSpeciesContext().getName() + BIO_PARAM_SUFFIX_SPECIES_CONCENTRATION;
Expression concExp = getExpressionAmtToConc(new Expression(spCountParm, getNameScope()), scs.getSpeciesContext().getStructure());
concExp.bindExpression(this);
addSpeciesConcentrationParameter(concName, concExp, PARAMETER_ROLE_SPECIES_CONCENRATION, scs.getSpeciesContext().getUnitDefinition(), scs.getSpeciesContext());
} catch (Exception e) {
e.printStackTrace();
throw new MappingException(e.getMessage());
}
// we always add variables, all species are independent variables, no matter they are constant or not.
String countMathSymbol = getMathSymbol(spCountParm, getSimulationContext().getGeometryContext().getStructureMapping(scs.getSpeciesContext().getStructure()).getGeometryClass());
scm.setVariable(new VolumeParticleVariable(countMathSymbol, defaultDomain));
mathSymbolMapping.put(scm.getSpeciesContext(), scm.getVariable().getName());
}
for (RbmObservable observable : simContext.getModel().getRbmModelContainer().getObservableList()) {
// stochastic variable is always a function of size.
ObservableCountParameter observableCountParm = null;
try {
String countName = observable.getName() + BIO_PARAM_SUFFIX_SPECIES_COUNT;
Expression countExp = new Expression(0.0);
observableCountParm = addObservableCountParameter(countName, countExp, PARAMETER_ROLE_OBSERVABLE_COUNT, simContext.getModel().getUnitSystem().getStochasticSubstanceUnit(), observable);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add concentration of species as MathMappingParameter - this will map to species concentration function
try {
String concName = observable.getName() + BIO_PARAM_SUFFIX_SPECIES_CONCENTRATION;
Expression concExp = getExpressionAmtToConc(new Expression(observableCountParm, getNameScope()), observable.getStructure());
concExp.bindExpression(this);
addObservableConcentrationParameter(concName, concExp, PARAMETER_ROLE_OBSERVABLE_CONCENTRATION, observable.getUnitDefinition(), observable);
} catch (Exception e) {
e.printStackTrace();
throw new MappingException(e.getMessage());
}
// //we always add variables, all species are independent variables, no matter they are constant or not.
// String countMathSymbol = getMathSymbol(observableCountParm, getSimulationContext().getGeometryContext().getStructureMapping(scs.getSpeciesContext().getStructure()).getGeometryClass());
// scm.setVariable(new VolumeParticleVariable(countMathSymbol,defaultDomain));
// mathSymbolMapping.put(observable,scm.getVariable().getName());
}
}
use of cbit.vcell.math.Variable.Domain in project vcell by virtualcell.
the class RulebasedMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
@Override
protected void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// use local variable instead of using getter all the time.
SimulationContext simContext = getSimulationContext();
GeometryClass geometryClass = simContext.getGeometry().getGeometrySpec().getSubVolumes()[0];
Domain domain = new Domain(geometryClass);
// local structure mapping list
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
// We have to check if all the reactions are able to transform to stochastic jump processes before generating the math.
String stochChkMsg = simContext.getModel().isValidForStochApp();
if (!(stochChkMsg.equals(""))) {
throw new ModelException("Problem updating math description: " + simContext.getName() + "\n" + stochChkMsg);
}
simContext.checkValidity();
//
if (simContext.getGeometry().getDimension() > 0) {
throw new MappingException("rule-based particle math mapping not implemented for spatial geometry - dimension >= 1");
}
//
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = simContext.getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
//
Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
for (int i = 0; i < structures.length; i++) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
if (sm == null || (sm instanceof FeatureMapping && ((FeatureMapping) sm).getGeometryClass() == null)) {
throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subVolume");
}
if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
try {
if (volFractExp != null) {
double volFract = volFractExp.evaluateConstant();
if (volFract >= 1.0) {
throw new MappingException("model structure '" + (getSimulationContext().getModel().getStructureTopology().getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0"));
}
}
} catch (ExpressionException e) {
e.printStackTrace(System.out);
}
}
}
SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
Structure[] mappedStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolumes[i]);
if (mappedStructures == null || mappedStructures.length == 0) {
throw new MappingException("geometry subVolume '" + subVolumes[i].getName() + "' not mapped from a model structure");
}
}
//
// gather only those reactionRules that are not "excluded"
//
ArrayList<ReactionRule> rrList = new ArrayList<ReactionRule>();
for (ReactionRuleSpec reactionRuleSpec : simContext.getReactionContext().getReactionRuleSpecs()) {
if (!reactionRuleSpec.isExcluded()) {
rrList.add(reactionRuleSpec.getReactionRule());
}
}
//
for (ReactionRule reactionRule : rrList) {
UnresolvedParameter[] unresolvedParameters = reactionRule.getKineticLaw().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException("In Application '" + simContext.getName() + "', " + reactionRule.getDisplayType() + " '" + reactionRule.getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = simContext.getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(simContext.getName() + "_generated");
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates
//
VariableHash varHash = new VariableHash();
//
// conversion factors
//
Model model = simContext.getModel();
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof StochVolVariable) {
varHash.addVariable(scm.getVariable());
}
}
// deals with model parameters
ModelParameter[] modelParameters = simContext.getModel().getModelParameters();
for (int j = 0; j < modelParameters.length; j++) {
Expression expr = getSubstitutedExpr(modelParameters[j].getExpression(), true, false);
expr = getIdentifierSubstitutions(expr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), expr, geometryClass));
}
// added July 2009, ElectricalStimulusParameter electric mapping tab
ElectricalStimulus[] elecStimulus = simContext.getElectricalStimuli();
if (elecStimulus.length > 0) {
throw new MappingException("Modles with electrophysiology are not supported for stochastic applications.");
}
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
Parameter initialVoltageParm = memMapping.getInitialVoltageParameter();
try {
Expression exp = initialVoltageParm.getExpression();
exp.evaluateConstant();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
}
}
}
//
for (ReactionRule reactionRule : rrList) {
// if (reactionRule.getKineticLaw() instanceof LumpedKinetics){
// throw new RuntimeException("Lumped Kinetics not yet supported for RuleBased Modeling");
// }
LocalParameter[] parameters = reactionRule.getKineticLaw().getLocalParameters();
for (LocalParameter parameter : parameters) {
//
if ((parameter.getRole() == RbmKineticLawParameterType.RuleRate)) {
continue;
}
//
if (!reactionRule.isReversible() && parameter.getRole() == RbmKineticLawParameterType.MassActionReverseRate) {
continue;
}
Expression expr = getSubstitutedExpr(parameter.getExpression(), true, false);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameter, geometryClass), getIdentifierSubstitutions(expr, parameter.getUnitDefinition(), geometryClass), geometryClass));
}
}
// the parameter "Size" is already put into mathsymbolmapping in refreshSpeciesContextMapping()
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm.getExpression() != null) {
try {
double value = parm.getExpression().evaluateConstant();
varHash.addVariable(new Constant(getMathSymbol(parm, sm.getGeometryClass()), new Expression(value)));
} catch (ExpressionException e) {
// varHash.addVariable(new Function(getMathSymbol0(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
e.printStackTrace(System.out);
throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
}
}
}
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
addInitialConditions(domain, speciesContextSpecs, varHash);
//
if (simContext.getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("Geometry must be defined in Application " + simContext.getName());
}
//
// create subDomains
//
SubVolume subVolume = simContext.getGeometry().getGeometrySpec().getSubVolumes()[0];
SubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), 0);
mathDesc.addSubDomain(subDomain);
//
// define all molecules and unique species patterns (add molecules to mathDesc and speciesPatterns to varHash).
//
HashMap<SpeciesPattern, VolumeParticleSpeciesPattern> speciesPatternMap = addSpeciesPatterns(domain, rrList);
HashSet<VolumeParticleSpeciesPattern> uniqueParticleSpeciesPatterns = new HashSet<>(speciesPatternMap.values());
for (VolumeParticleSpeciesPattern volumeParticleSpeciesPattern : uniqueParticleSpeciesPatterns) {
varHash.addVariable(volumeParticleSpeciesPattern);
}
//
// define observables (those explicitly declared and those corresponding to seed species.
//
List<ParticleObservable> observables = addObservables(geometryClass, domain, speciesPatternMap);
for (ParticleObservable particleObservable : observables) {
varHash.addVariable(particleObservable);
}
try {
addParticleJumpProcesses(varHash, geometryClass, subDomain, speciesPatternMap);
} catch (PropertyVetoException e1) {
e1.printStackTrace();
throw new MappingException(e1.getMessage(), e1);
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter || fieldMathMappingParameters[i] instanceof ObservableConcentrationParameter) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
for (SpeciesContext sc : model.getSpeciesContexts()) {
if (!sc.hasSpeciesPattern()) {
throw new MappingException("species " + sc.getName() + " has no molecular pattern");
}
VolumeParticleSpeciesPattern volumeParticleSpeciesPattern = speciesPatternMap.get(sc.getSpeciesPattern());
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleProperties.ParticleInitialCondition>();
// initial conditions from scs
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
Parameter initialCountParameter = scs.getInitialCountParameter();
Expression e = getIdentifierSubstitutions(new Expression(initialCountParameter, getNameScope()), initialCountParameter.getUnitDefinition(), geometryClass);
particleInitialConditions.add(new ParticleInitialConditionCount(e, new Expression(0.0), new Expression(0.0), new Expression(0.0)));
ParticleProperties particleProperies = new ParticleProperties(volumeParticleSpeciesPattern, new Expression(0.0), new Expression(0.0), new Expression(0.0), new Expression(0.0), particleInitialConditions);
subDomain.addParticleProperties(particleProperies);
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
if (fieldMathMappingParameters[i] instanceof ObservableConcentrationParameter) {
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
System.out.println(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
}
Aggregations