Examples with SimpleBounds cbit.vcell.constraints.SimpleBounds used on opensource projects

Example 6 with SimpleBounds

use of cbit.vcell.constraints.SimpleBounds in project vcell by virtualcell.

the class ApplicationConstraintsGenerator method fromApplication.

/**
 * Insert the method's description here.
 * Creation date: (6/26/01 8:25:55 AM)
 * @return cbit.vcell.constraints.ConstraintContainerImpl
 */
public static ConstraintContainerImpl fromApplication(SimulationContext simContext) {
    try {
        ConstraintContainerImpl ccImpl = new ConstraintContainerImpl();
        // ====================
        // add physical limits
        // ====================
        // 
        // no negative concentrations
        // 
        cbit.vcell.model.Model model = simContext.getModel();
        cbit.vcell.model.SpeciesContext[] speciesContexts = model.getSpeciesContexts();
        for (int i = 0; i < speciesContexts.length; i++) {
            ccImpl.addSimpleBound(new SimpleBounds(speciesContexts[i].getName(), new RealInterval(0, Double.POSITIVE_INFINITY), AbstractConstraint.PHYSICAL_LIMIT, "non-negative concentration"));
        }
        for (int i = 0; i < speciesContexts.length; i++) {
            SpeciesContextSpecParameter initParam = (simContext.getReactionContext().getSpeciesContextSpec(speciesContexts[i])).getInitialConditionParameter();
            if (initParam != null) {
                double initialValue = initParam.getExpression().evaluateConstant();
                ccImpl.addSimpleBound(new SimpleBounds(speciesContexts[i].getName(), new RealInterval(initialValue), AbstractConstraint.MODELING_ASSUMPTION, "specified \"initialCondition\""));
            }
        }
        // =========================
        // add modeling assumptions
        // =========================
        // 
        // mass action forward and reverse rates should be non-negative
        // 
        cbit.vcell.model.ReactionStep[] reactionSteps = model.getReactionSteps();
        for (int i = 0; i < reactionSteps.length; i++) {
            Kinetics kinetics = reactionSteps[i].getKinetics();
            if (kinetics instanceof MassActionKinetics) {
                Expression forwardRateConstraintExp = new Expression(((MassActionKinetics) kinetics).getForwardRateParameter().getExpression().infix() + ">=0");
                forwardRateConstraintExp = getSteadyStateExpression(forwardRateConstraintExp);
                if (!forwardRateConstraintExp.compareEqual(new Expression(1.0))) {
                    ccImpl.addGeneralConstraint(new GeneralConstraint(forwardRateConstraintExp, AbstractConstraint.MODELING_ASSUMPTION, "non-negative forward rate"));
                }
                Expression reverseRateConstraintExp = new Expression(((MassActionKinetics) kinetics).getReverseRateParameter().getExpression().infix() + ">=0");
                reverseRateConstraintExp = getSteadyStateExpression(reverseRateConstraintExp);
                if (!reverseRateConstraintExp.compareEqual(new Expression(1.0))) {
                    ccImpl.addGeneralConstraint(new GeneralConstraint(reverseRateConstraintExp, AbstractConstraint.MODELING_ASSUMPTION, "non-negative reverse rate"));
                }
            }
            KineticsParameter authoritativeParameter = kinetics.getAuthoritativeParameter();
            Expression kineticRateConstraintExp = new Expression(authoritativeParameter.getName() + "==" + authoritativeParameter.getExpression().infix());
            kineticRateConstraintExp = getSteadyStateExpression(kineticRateConstraintExp);
            if (!kineticRateConstraintExp.compareEqual(new Expression(1.0))) {
                ccImpl.addGeneralConstraint(new GeneralConstraint(kineticRateConstraintExp, AbstractConstraint.MODELING_ASSUMPTION, "definition"));
            }
        }
        // 
        for (int i = 0; i < reactionSteps.length; i++) {
            Kinetics kinetics = reactionSteps[i].getKinetics();
            Kinetics.KineticsParameter[] parameters = kinetics.getKineticsParameters();
            for (int j = 0; j < parameters.length; j++) {
                Expression exp = parameters[j].getExpression();
                if (exp.getSymbols() == null || exp.getSymbols().length == 0) {
                    // 
                    try {
                        double constantValue = exp.evaluateConstant();
                        RealInterval interval = new RealInterval(constantValue);
                        ccImpl.addSimpleBound(new SimpleBounds(parameters[j].getName(), interval, AbstractConstraint.MODELING_ASSUMPTION, "model value"));
                    } catch (cbit.vcell.parser.ExpressionException e) {
                        System.out.println("error evaluating parameter " + parameters[j].getName() + " in reaction step " + reactionSteps[i].getName());
                    }
                } else {
                    Expression parameterDefinitionExp = new Expression(parameters[j].getName() + "==" + parameters[j].getExpression().infix());
                    parameterDefinitionExp = getSteadyStateExpression(parameterDefinitionExp);
                    if (!parameterDefinitionExp.compareEqual(new Expression(1.0))) {
                        ccImpl.addGeneralConstraint(new GeneralConstraint(parameterDefinitionExp, AbstractConstraint.MODELING_ASSUMPTION, "parameter definition"));
                    }
                }
            }
        }
        ccImpl.addSimpleBound(new SimpleBounds(model.getFARADAY_CONSTANT().getName(), new RealInterval(model.getFARADAY_CONSTANT().getExpression().evaluateConstant()), AbstractConstraint.PHYSICAL_LIMIT, "Faraday's constant"));
        ccImpl.addSimpleBound(new SimpleBounds(model.getTEMPERATURE().getName(), new RealInterval(300), AbstractConstraint.PHYSICAL_LIMIT, "Absolute Temperature Kelvin"));
        ccImpl.addSimpleBound(new SimpleBounds(model.getGAS_CONSTANT().getName(), new RealInterval(model.getGAS_CONSTANT().getExpression().evaluateConstant()), AbstractConstraint.PHYSICAL_LIMIT, "ideal gas constant"));
        ccImpl.addSimpleBound(new SimpleBounds(model.getKMILLIVOLTS().getName(), new RealInterval(model.getKMILLIVOLTS().getExpression().evaluateConstant()), AbstractConstraint.PHYSICAL_LIMIT, "ideal gas constant"));
        return ccImpl;
    } catch (cbit.vcell.parser.ExpressionException e) {
        e.printStackTrace(System.out);
        return null;
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace(System.out);
        return null;
    }
}

Example 7 with SimpleBounds

use of cbit.vcell.constraints.SimpleBounds in project vcell by virtualcell.

the class StructureSizeSolver method updateUnitStructureSizes.

/**
 * Insert the method's description here.
 * Creation date: (5/17/2006 10:33:38 AM)
 * @return double[]
 * @param structName java.lang.String
 * @param structSize double
 */
public static void updateUnitStructureSizes(SimulationContext simContext, GeometryClass geometryClass) {
    if (simContext.getGeometryContext().getGeometry().getDimension() == 0) {
        return;
    }
    StructureMapping[] myStructMappings = simContext.getGeometryContext().getStructureMappings(geometryClass);
    if (myStructMappings != null && myStructMappings.length == 1) {
        // if the unitSizeParameter is dimensionless, then features are mapped to SubVolumes or Membranes are mapped to surfaces (should sum to 1)
        boolean bDimensionless = myStructMappings[0].getUnitSizeParameter().getUnitDefinition().isEquivalent(simContext.getModel().getUnitSystem().getInstance_DIMENSIONLESS());
        if (bDimensionless) {
            try {
                myStructMappings[0].getUnitSizeParameter().setExpression(new Expression(1.0));
                return;
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                throw new RuntimeException(e.getMessage());
            }
        }
    }
    if (myStructMappings != null && myStructMappings.length == 0) {
        // nothing to solve, there are no mappings for this geometryClass
        return;
    }
    StructureMapping[] structMappings = simContext.getGeometryContext().getStructureMappings();
    try {
        ConstraintContainerImpl ccImpl = new ConstraintContainerImpl();
        Structure struct = null;
        Expression totalVolExpr = new Expression(0.0);
        StructureTopology structureTopology = simContext.getModel().getStructureTopology();
        for (int i = 0; i < structMappings.length; i++) {
            if (structMappings[i].getGeometryClass() != geometryClass) {
                continue;
            }
            // new model with unit sizes already
            if (structMappings[i].getUnitSizeParameter() != null && structMappings[i].getUnitSizeParameter().getExpression() != null) {
                return;
            }
            if (struct == null) {
                struct = structMappings[i].getStructure();
            }
            if (structMappings[i] instanceof MembraneMapping) {
                MembraneMapping membraneMapping = (MembraneMapping) structMappings[i];
                Membrane membrane = membraneMapping.getMembrane();
                String membraneSizeName = TokenMangler.mangleToSName(membrane.getName() + "_size");
                ccImpl.addSimpleBound(new SimpleBounds(membraneSizeName, new RealInterval(0, 100000), AbstractConstraint.PHYSICAL_LIMIT, "definition"));
                Feature insideFeature = structureTopology.getInsideFeature(membrane);
                String volFractName = TokenMangler.mangleToSName(insideFeature.getName() + "_volFract");
                String svRatioName = TokenMangler.mangleToSName(insideFeature.getName() + "_svRatio");
                StructureMapping.StructureMappingParameter volFractParameter = membraneMapping.getVolumeFractionParameter();
                double volFractValue = volFractParameter.getExpression().evaluateConstant();
                ccImpl.addSimpleBound(new SimpleBounds(volFractName, new RealInterval(volFractValue, volFractValue), AbstractConstraint.MODELING_ASSUMPTION, "from model"));
                StructureMapping.StructureMappingParameter surfToVolParameter = membraneMapping.getSurfaceToVolumeParameter();
                double svRatioValue = surfToVolParameter.getExpression().evaluateConstant();
                ccImpl.addSimpleBound(new SimpleBounds(svRatioName, new RealInterval(svRatioValue, svRatioValue), AbstractConstraint.MODELING_ASSUMPTION, "from model"));
                // membrane mapped to volume
                if (geometryClass instanceof SubVolume) {
                    // 
                    // EC eclosing cyt, which contains er and golgi
                    // "(cyt_size+ er_size + golgi_size) * cyt_svRatio - PM_size == 0"
                    // 
                    Expression sumOfInsideVolumeExp = new Expression(0.0);
                    for (int j = 0; j < structMappings.length; j++) {
                        if (structMappings[j] instanceof FeatureMapping && structureTopology.enclosedBy(structMappings[j].getStructure(), insideFeature)) {
                            Feature childFeatureOfInside = ((FeatureMapping) structMappings[j]).getFeature();
                            if (simContext.getGeometryContext().getStructureMapping(childFeatureOfInside).getGeometryClass() == geometryClass) {
                                sumOfInsideVolumeExp = Expression.add(sumOfInsideVolumeExp, new Expression(TokenMangler.mangleToSName(childFeatureOfInside.getName() + "_size")));
                            }
                        }
                    }
                    Expression tempExpr = Expression.mult(sumOfInsideVolumeExp, new Expression(svRatioName));
                    tempExpr = Expression.add(tempExpr, new Expression("-" + membraneSizeName));
                    ccImpl.addGeneralConstraint(new GeneralConstraint(new Expression(tempExpr.infix() + "==0"), AbstractConstraint.MODELING_ASSUMPTION, "svRatio definition"));
                    // 
                    // EC eclosing cyt, which contains er and golgi
                    // (EC_size + cyt_size + er_size + golgi_size) * cyt_vfRatio - (cyt_size + er_size + golgi_size) == 0
                    // 
                    Feature outsideFeature = structureTopology.getOutsideFeature(membrane);
                    Expression sumOfParentVolumeExp = new Expression(0.0);
                    for (int j = 0; j < structMappings.length; j++) {
                        if (structMappings[j] instanceof FeatureMapping && structureTopology.enclosedBy(structMappings[j].getStructure(), outsideFeature)) {
                            Feature childFeatureOfParent = ((FeatureMapping) structMappings[j]).getFeature();
                            if (simContext.getGeometryContext().getStructureMapping(childFeatureOfParent).getGeometryClass() == geometryClass) {
                                sumOfParentVolumeExp = Expression.add(sumOfParentVolumeExp, new Expression(TokenMangler.mangleToSName(childFeatureOfParent.getName() + "_size")));
                            }
                        }
                    }
                    Expression exp = Expression.mult(sumOfParentVolumeExp, new Expression(volFractName));
                    exp = Expression.add(exp, Expression.negate(sumOfInsideVolumeExp));
                    ccImpl.addGeneralConstraint(new GeneralConstraint(new Expression(exp.infix() + "==0.0"), AbstractConstraint.MODELING_ASSUMPTION, "volFract definition"));
                }
            } else if (structMappings[i] instanceof FeatureMapping) {
                FeatureMapping featureMapping = (FeatureMapping) structMappings[i];
                String featureSizeName = TokenMangler.mangleToSName(featureMapping.getFeature().getName() + "_size");
                totalVolExpr = Expression.add(totalVolExpr, new Expression(featureSizeName));
                ccImpl.addSimpleBound(new SimpleBounds(featureSizeName, new RealInterval(0, 1), AbstractConstraint.PHYSICAL_LIMIT, "definition"));
            }
        }
        if (geometryClass instanceof SubVolume) {
            ccImpl.addGeneralConstraint(new GeneralConstraint(new Expression(totalVolExpr.infix() + "==1.0"), AbstractConstraint.MODELING_ASSUMPTION, "total volume"));
        }
        // ccImpl.show();
        ConstraintSolver constraintSolver = new ConstraintSolver(ccImpl);
        constraintSolver.resetIntervals();
        int numTimesNarrowed = 0;
        RealInterval[] lastSolution = null;
        boolean bChanged = true;
        while (constraintSolver.narrow() && bChanged && numTimesNarrowed < 125) {
            numTimesNarrowed++;
            bChanged = false;
            RealInterval[] thisSolution = constraintSolver.getIntervals();
            if (lastSolution != null) {
                for (int i = 0; i < thisSolution.length; i++) {
                    if (!thisSolution[i].equals(lastSolution[i])) {
                        bChanged = true;
                    }
                }
            } else {
                bChanged = true;
            }
            lastSolution = thisSolution;
        }
        System.out.println("num of times narrowed = " + numTimesNarrowed);
        if (numTimesNarrowed > 0) {
            String[] symbols = constraintSolver.getSymbols();
            net.sourceforge.interval.ia_math.RealInterval[] solution = constraintSolver.getIntervals();
            double totalArea = 0;
            double totalVolume = 0;
            for (int i = 0; i < symbols.length; i++) {
                System.out.println("solution[" + i + "] \"" + symbols[i] + "\" = " + solution[i]);
                for (int j = 0; j < structMappings.length; j++) {
                    if (symbols[i].equals(TokenMangler.mangleToSName(structMappings[j].getStructure().getName() + "_size"))) {
                        if (!Double.isInfinite(solution[i].lo()) && !Double.isInfinite(solution[i].hi())) {
                            double value = (solution[i].lo() + solution[i].hi()) / 2;
                            Expression exp = new Expression(value);
                            if (structMappings[j] instanceof FeatureMapping) {
                                FeatureMapping fm = (FeatureMapping) structMappings[j];
                                totalVolume += value;
                                if (geometryClass instanceof SubVolume) {
                                    fm.getVolumePerUnitVolumeParameter().setExpression(exp);
                                } else if (geometryClass instanceof SurfaceClass) {
                                    fm.getVolumePerUnitAreaParameter().setExpression(exp);
                                }
                            } else if (structMappings[j] instanceof MembraneMapping) {
                                MembraneMapping mm = (MembraneMapping) structMappings[j];
                                totalArea += value;
                                if (geometryClass instanceof SubVolume) {
                                    mm.getAreaPerUnitVolumeParameter().setExpression(exp);
                                } else if (geometryClass instanceof SurfaceClass) {
                                    mm.getAreaPerUnitAreaParameter().setExpression(exp);
                                }
                            }
                        }
                    }
                }
            }
            // 
            // normalize all so that total volume is 1.0 for subVolumes or
            // total area is 1.0 for surfaceClasses
            // 
            double scaleFactor = 1;
            if (geometryClass instanceof SubVolume) {
                scaleFactor = totalVolume;
            } else if (geometryClass instanceof SurfaceClass) {
                scaleFactor = totalArea;
            } else {
                throw new RuntimeException("unexpected GeometryClass");
            }
            for (int j = 0; j < structMappings.length; j++) {
                if (structMappings[j].getGeometryClass() == geometryClass) {
                    if (structMappings[j] instanceof FeatureMapping) {
                        FeatureMapping fm = (FeatureMapping) structMappings[j];
                        if (geometryClass instanceof SubVolume) {
                            fm.getVolumePerUnitVolumeParameter().setExpression(new Expression(fm.getVolumePerUnitVolumeParameter().getExpression().evaluateConstant() / scaleFactor));
                        } else if (geometryClass instanceof SurfaceClass) {
                            fm.getVolumePerUnitAreaParameter().setExpression(new Expression(fm.getVolumePerUnitAreaParameter().getExpression().evaluateConstant() / scaleFactor));
                        }
                    } else if (structMappings[j] instanceof MembraneMapping) {
                        MembraneMapping mm = (MembraneMapping) structMappings[j];
                        if (geometryClass instanceof SubVolume) {
                            mm.getAreaPerUnitVolumeParameter().setExpression(new Expression(mm.getAreaPerUnitVolumeParameter().getExpression().evaluateConstant() / scaleFactor));
                        } else if (geometryClass instanceof SurfaceClass) {
                            mm.getAreaPerUnitAreaParameter().setExpression(new Expression(mm.getAreaPerUnitAreaParameter().getExpression().evaluateConstant() / scaleFactor));
                        }
                    }
                }
            }
        } else {
            throw new RuntimeException("cannot solve for size");
        }
    } catch (ExpressionException e) {
        e.printStackTrace(System.out);
        throw new RuntimeException(e.getMessage());
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace(System.out);
        throw new RuntimeException(e.getMessage());
    }
}