Examples with RealInterval net.sourceforge.interval.ia_math.RealInterval used on opensource projects

Example 16 with RealInterval

use of net.sourceforge.interval.ia_math.RealInterval 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 17 with RealInterval

use of net.sourceforge.interval.ia_math.RealInterval in project vcell by virtualcell.

the class Kinetics method gatherIssues.

/**
 * Insert the method's description here.
 * Creation date: (5/12/2004 2:53:13 PM)
 */
public void gatherIssues(IssueContext issueContext, List<Issue> issueList) {
    issueContext = issueContext.newChildContext(ContextType.ModelProcessDynamics, this);
    // 
    for (int i = 0; fieldUnresolvedParameters != null && i < fieldUnresolvedParameters.length; i++) {
        issueList.add(new Issue(fieldUnresolvedParameters[i], issueContext, IssueCategory.UnresolvedParameter, "Unresolved parameter '" + fieldUnresolvedParameters[i].getName() + "' in reaction '" + reactionStep.getName() + "'", Issue.SEVERITY_ERROR));
    }
    // 
    for (int i = 0; fieldKineticsParameters != null && i < fieldKineticsParameters.length; i++) {
        if (fieldKineticsParameters[i].getRole() == ROLE_UserDefined) {
            try {
                if (!isReferenced(fieldKineticsParameters[i], 0)) {
                    issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.KineticsUnreferencedParameter, "Unreferenced Kinetic Parameter '" + fieldKineticsParameters[i].getName() + "' in reaction '" + reactionStep.getName() + "'", Issue.SEVERITY_WARNING));
                }
            } catch (ExpressionException e) {
                issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.KineticsExpressionError, "error resolving expression " + e.getMessage(), Issue.SEVERITY_WARNING));
            } catch (ModelException e) {
                issueList.add(new Issue(getReactionStep(), issueContext, IssueCategory.CyclicDependency, "cyclic dependency in the parameter definitions", Issue.SEVERITY_ERROR));
            }
        }
    }
    // 
    if (fieldKineticsParameters != null) {
        for (KineticsParameter kineticsParameter : fieldKineticsParameters) {
            if (kineticsParameter.getExpression() == null) {
                issueList.add(new Issue(kineticsParameter, issueContext, IssueCategory.KineticsExpressionMissing, "expression is missing", Issue.SEVERITY_INFO));
            } else {
                Expression exp = kineticsParameter.getExpression();
                String[] symbols = exp.getSymbols();
                String issueMessagePrefix = "Kinetic parameter '" + kineticsParameter.getName() + "' in reaction '" + getReactionStep().getName() + "' ";
                if (symbols != null) {
                    for (int j = 0; j < symbols.length; j++) {
                        SymbolTableEntry ste = exp.getSymbolBinding(symbols[j]);
                        if (ste instanceof KineticsProxyParameter) {
                            ste = ((KineticsProxyParameter) ste).getTarget();
                        }
                        if (ste == null) {
                            issueList.add(new Issue(kineticsParameter, issueContext, IssueCategory.KineticsExpressionUndefinedSymbol, issueMessagePrefix + "references undefined symbol '" + symbols[j] + "'", Issue.SEVERITY_ERROR));
                        } else if (ste instanceof SpeciesContext) {
                            if (!getReactionStep().getModel().contains((SpeciesContext) ste)) {
                                issueList.add(new Issue(kineticsParameter, issueContext, IssueCategory.KineticsExpressionUndefinedSymbol, issueMessagePrefix + "references undefined species '" + symbols[j] + "'", Issue.SEVERITY_ERROR));
                            }
                            if (reactionStep.countNumReactionParticipants((SpeciesContext) ste) == 0) {
                                issueList.add(new Issue(kineticsParameter, issueContext, IssueCategory.KineticsExpressionNonParticipantSymbol, issueMessagePrefix + "references species context '" + symbols[j] + "', but it is not a reactant/product/catalyst of this reaction", Issue.SEVERITY_WARNING));
                            }
                        } else if (ste instanceof ModelParameter) {
                            if (!getReactionStep().getModel().contains((ModelParameter) ste)) {
                                issueList.add(new Issue(kineticsParameter, issueContext, IssueCategory.KineticsExpressionUndefinedSymbol, issueMessagePrefix + "references undefined global parameter '" + symbols[j] + "'", Issue.SEVERITY_ERROR));
                            }
                        }
                    }
                }
            }
        }
        // looking for local param which masks a global and issueing a warning
        for (KineticsParameter kineticsParameter : fieldKineticsParameters) {
            String name = kineticsParameter.getName();
            SymbolTableEntry ste = getReactionStep().getNameScope().getExternalEntry(name, getReactionStep());
            String steName;
            if (ste != null) {
                if (ste instanceof Displayable) {
                    steName = ((Displayable) ste).getDisplayType() + " " + ste.getName();
                } else {
                    steName = ste.getClass().getSimpleName() + " " + ste.getName();
                }
                String msg = steName + " is overriden by a local parameter " + name + " in reaction " + getReactionStep().getName();
                issueList.add(new Issue(kineticsParameter, issueContext, IssueCategory.Identifiers, msg, Issue.SEVERITY_WARNING));
            }
        }
    }
    try {
        // 
        // determine unit consistency for each expression
        // 
        ModelUnitSystem modelUnitSystem = getReactionStep().getModel().getUnitSystem();
        VCUnitEvaluator unitEvaluator = new VCUnitEvaluator(modelUnitSystem);
        for (int i = 0; i < fieldKineticsParameters.length; i++) {
            try {
                VCUnitDefinition paramUnitDef = fieldKineticsParameters[i].getUnitDefinition();
                VCUnitDefinition expUnitDef = unitEvaluator.getUnitDefinition(fieldKineticsParameters[i].getExpression());
                if (paramUnitDef == null) {
                    issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.Units, "defined unit is null", Issue.SEVERITY_WARNING));
                } else if (paramUnitDef.isTBD()) {
                    issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.Units, "undefined unit " + modelUnitSystem.getInstance_TBD().getSymbol(), Issue.SEVERITY_WARNING));
                } else if (expUnitDef == null) {
                    issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.Units, "computed unit is null", Issue.SEVERITY_WARNING));
                } else if (paramUnitDef.isTBD() || (!paramUnitDef.isEquivalent(expUnitDef) && !expUnitDef.isTBD())) {
                    issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.Units, "inconsistent units, defined=[" + fieldKineticsParameters[i].getUnitDefinition().getSymbol() + "], computed=[" + expUnitDef.getSymbol() + "]", Issue.SEVERITY_WARNING));
                }
            } catch (VCUnitException e) {
                issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.Units, e.getMessage(), Issue.SEVERITY_WARNING));
            } catch (ExpressionException e) {
                issueList.add(new Issue(fieldKineticsParameters[i], issueContext, IssueCategory.Units, e.getMessage(), Issue.SEVERITY_WARNING));
            }
        }
    } catch (Throwable e) {
        issueList.add(new Issue(getReactionStep(), issueContext, IssueCategory.Units, "unexpected exception: " + e.getMessage(), Issue.SEVERITY_INFO));
    }
    // 
    for (int i = 0; i < fieldKineticsParameters.length; i++) {
        RealInterval simpleBounds = bounds[fieldKineticsParameters[i].getRole()];
        if (simpleBounds != null) {
            String parmName = reactionStep.getNameScope().getName() + "." + fieldKineticsParameters[i].getName();
            issueList.add(new SimpleBoundsIssue(fieldKineticsParameters[i], issueContext, simpleBounds, "parameter " + parmName + ": must be within " + simpleBounds.toString()));
        }
    }
}

Example 18 with RealInterval

use of net.sourceforge.interval.ia_math.RealInterval 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());
    }
}

Example 19 with RealInterval

use of net.sourceforge.interval.ia_math.RealInterval in project vcell by virtualcell.

the class ExpressionTest method testEvaluateNarrowingOld.

/**
 * This method was created in VisualAge.
 */
public static void testEvaluateNarrowingOld() {
    try {
        Expression[] exps = { (new Expression("a*b==c;")).flatten(), (new Expression("c<2*a;")).flatten() // (new Expression("a>1;")).flatten(),
        // (new Expression("a<3;")).flatten(),
        // (new Expression("b>2;")).flatten(),
        // (new Expression("b<4;")).flatten(),
        // (new Expression("c>4;")).flatten(),
        // (new Expression("c<100;")).flatten(),
        };
        // RealInterval v[] = { RealInterval.fullInterval(), RealInterval.fullInterval(), RealInterval.fullInterval() };
        RealInterval[] v = { new RealInterval(1, 3), new RealInterval(2, 4), new RealInterval(4, 100) };
        SimpleSymbolTable simpleSymbolTable = new SimpleSymbolTable(new String[] { "a", "b", "c" });
        for (int i = 0; i < exps.length; i++) {
            exps[i].bindExpression(simpleSymbolTable);
        }
        for (int i = 0; i < 5; i++) {
            for (int j = 0; j < exps.length; j++) {
                System.out.println(i + ") narrowing '" + exps[j] + "'");
                if (!exps[j].narrow(v)) {
                    System.out.println("narrowing failed");
                    break;
                }
                for (int k = 0; k < v.length; k++) {
                    System.out.println("     " + v[k]);
                }
            }
        }
    } catch (Throwable e) {
        e.printStackTrace(System.out);
    }
}

Example 20 with RealInterval

use of net.sourceforge.interval.ia_math.RealInterval in project vcell by virtualcell.

the class IntervalNarrowingTest method testIntervalNarrowing.

@Test
public void testIntervalNarrowing() throws DivideByZeroException, ExpressionException {
    final Random random = new Random(0);
    String[] symbols = constraint.getSymbols();
    if (symbols == null || symbols.length < 1) {
        Assert.fail("<<" + constraint.infix() + ">> doesn't have any symbols in it, can't be a proper constraint, bad test case");
    }
    SimpleSymbolTable symbolTable = new SimpleSymbolTable(symbols);
    constraint.bindExpression(symbolTable);
    RealInterval[] initialItervals = new RealInterval[symbols.length];
    RealInterval[] values = new RealInterval[symbols.length];
    RealInterval[] prevValues = new RealInterval[symbols.length];
    // 
    // try to find an interval that satisfies the constraints
    // 
    final int MAX_INITIAL_VALUE_ATTEMPTS = 1000;
    int initialValueCount = 0;
    boolean bInitialValueSatisfiesConstraint = false;
    while (initialValueCount < MAX_INITIAL_VALUE_ATTEMPTS) {
        // 
        for (int j = 0; j < initialItervals.length; j++) {
            double middle = random.nextGaussian();
            double size = random.nextGaussian() * 10;
            initialItervals[j] = new RealInterval(middle - Math.abs(size), middle + Math.abs(size));
            values[j] = new RealInterval(initialItervals[j].lo(), initialItervals[j].hi());
            prevValues[j] = new RealInterval(initialItervals[j].lo(), initialItervals[j].hi());
        }
        initialValueCount++;
        try {
            RealInterval result = constraint.evaluateInterval(values);
            // test if it includes "true"
            result.intersect(new RealInterval(1.0));
            if (!result.nonEmpty()) {
                bInitialValueSatisfiesConstraint = true;
                break;
            }
        } catch (FunctionDomainException e) {
        } catch (IAException e) {
            Assert.assertEquals("exception indicates that intersection is empty", e.getMessage().contains("intersection is empty"), true);
        }
    }
    if (!bInitialValueSatisfiesConstraint) {
        // couldn't find a random initial constraint which satisfies constraint, try full intervals
        for (int j = 0; j < values.length; j++) {
            initialItervals[j] = RealInterval.fullInterval();
            values[j] = RealInterval.fullInterval();
            prevValues[j] = RealInterval.fullInterval();
        }
    }
    try {
        boolean bSatisfied = constraint.narrow(values);
    } catch (FunctionDomainException e) {
        Assert.fail("<<" + constraint.infix() + ">> " + getPoint(symbols, values) + " initial values result in FunctionDomainException, " + e.getMessage());
    } catch (IAException e) {
        Assert.fail("<<" + constraint.infix() + ">> " + getPoint(symbols, values) + " initial value result in IAException, " + e.getMessage());
    }
    boolean bValuesChanged = true;
    boolean bValuesEverChanged = false;
    final int maxIterations = 100;
    int iteration = 0;
    // 
    while (bValuesChanged && iteration < maxIterations) {
        iteration++;
        boolean bConstraintSatisfied = constraint.narrow(values);
        Assert.assertEquals("<<" + constraint.infix() + ">> " + getPoint(symbols, values) + " constraint not satisfied", bConstraintSatisfied, true);
        bValuesChanged = false;
        System.out.print("     iteration " + iteration);
        for (int j = 0; j < values.length; j++) {
            if (!prevValues[j].equals(values[j])) {
                bValuesChanged = true;
                bValuesEverChanged = true;
                prevValues[j] = (RealInterval) values[j].clone();
            }
            System.out.print(", " + symbols[j] + " = " + values[j].toString());
        }
        System.out.println();
    }
    boolean bNarrowedIsDifferent = false;
    // 
    for (int j = 0; j < values.length; j++) {
        Assert.assertEquals("<<" + constraint.infix() + ">> narrowing bounds new.low >= orig.low", values[j].lo() >= initialItervals[j].lo(), true);
        Assert.assertEquals("<<" + constraint.infix() + ">> narrowing bounds new.low <= orig.high", values[j].lo() <= initialItervals[j].hi(), true);
        Assert.assertEquals("<<" + constraint.infix() + ">> narrowing bounds new.high >= orig.low", values[j].hi() >= initialItervals[j].lo(), true);
        Assert.assertEquals("<<" + constraint.infix() + ">> narrowing bounds new.high <= orig.high", values[j].hi() <= initialItervals[j].hi(), true);
        if (!values[j].equals(initialItervals[j])) {
            bNarrowedIsDifferent = true;
        }
    }
    if (bNarrowedIsDifferent) {
        // 
        // verify that samples taken from the removed portions of the initial interval all fail.
        // 
        final int NUM_SAMPLES = 300;
        final int MAX_EXCEPTION_COUNT = 3000;
        int successCount = 0;
        int exceptionCount = 0;
        double[] samples = new double[values.length];
        FunctionDomainException functionDomainException = null;
        while (successCount < NUM_SAMPLES && exceptionCount < MAX_EXCEPTION_COUNT) {
            boolean bSamplesOutsideNarrowedIntervals = true;
            for (int k = 0; k < initialItervals.length; k++) {
                samples[k] = sample(initialItervals[k], random);
                if (samples[k] >= values[k].lo() && samples[k] <= values[k].hi()) {
                    bSamplesOutsideNarrowedIntervals = false;
                    exceptionCount++;
                    break;
                }
            }
            try {
                if (bSamplesOutsideNarrowedIntervals) {
                    double result = constraint.evaluateVector(samples);
                    Assert.assertEquals("<<" + constraint.infix() + ">> " + getPoint(symbols, samples) + " removed point was in solution", result == 0.0, true);
                    successCount++;
                }
            } catch (FunctionDomainException e) {
                exceptionCount++;
                functionDomainException = e;
            }
        }
        Assert.assertEquals("<<" + constraint.infix() + ">> only threw FunctionDomainExceptions, exception=" + functionDomainException, successCount < NUM_SAMPLES, false);
    }
}