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

Example 11 with RealInterval

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

the class ExpressionTest method testEvaluateInterval.

/**
 * This method was created in VisualAge.
 */
public static void testEvaluateInterval(int numTrials, int depth, long seed, boolean bIsConstraint) {
    int numCorrect = 0;
    int numWrong = 0;
    int numFailures = 0;
    java.util.Random random = new java.util.Random(seed);
    for (int i = 0; i < numTrials; i++) {
        try {
            cbit.vcell.parser.Expression exp = cbit.vcell.parser.ExpressionUtils.generateExpression(random, depth, bIsConstraint);
            // System.out.println("exp = '"+exp+"'");
            // 
            // generate random intervals for arguments
            // 
            String[] symbols = exp.getSymbols();
            RealInterval[] intervals = null;
            if (symbols != null && symbols.length > 0) {
                SimpleSymbolTable symbolTable = new SimpleSymbolTable(symbols);
                exp.bindExpression(symbolTable);
                intervals = new RealInterval[symbolTable.getSize()];
                for (int j = 0; j < intervals.length; j++) {
                    double value1 = random.nextGaussian();
                    double value2 = random.nextGaussian();
                    intervals[j] = new RealInterval(Math.min(value1, value2), Math.max(value1, value2));
                }
            } else {
                intervals = new RealInterval[0];
            }
            // 
            // evaluate using interval arithmetic
            // 
            RealInterval intervalResult = exp.evaluateInterval(intervals);
            // 
            // generate several "feasible points" (samples within input intervals)
            // and evaluate expression (must be within bounds of intervalResult).
            // 
            double[] doubleValues = new double[intervals.length];
            int maxExperiments = 1 + 40 * intervals.length;
            for (int experimentCount = 0; experimentCount < maxExperiments; experimentCount++) {
                // 
                for (int j = 0; j < doubleValues.length; j++) {
                    double low = intervals[j].lo();
                    double high = intervals[j].hi();
                    if (!Double.isInfinite(low) && !Double.isInfinite(high)) {
                        // 
                        // finite interval, uniformly sample between high and low
                        // 
                        doubleValues[j] = low + (random.nextDouble() * (high - low));
                    } else {
                        // 
                        // infinite interval, keep trying Gaussians until within the bounds
                        // 
                        boolean bWithinRange = false;
                        while (!bWithinRange) {
                            double sample = random.nextGaussian();
                            if (sample >= low && sample <= high) {
                                doubleValues[j] = sample;
                                bWithinRange = true;
                            }
                        }
                    }
                }
                // 
                // verify that expression evaluates within intervalResult
                // 
                double doubleResult = exp.evaluateVector(doubleValues);
                if (doubleResult < intervalResult.lo() || doubleResult > intervalResult.hi()) {
                    numWrong++;
                    System.out.println(exp);
                    for (int k = 0; symbols != null && k < symbols.length; k++) {
                        System.out.println(symbols[k] + " = " + intervals[k] + " ... sampled = " + doubleValues[k]);
                    }
                    System.out.println("interval result = " + intervalResult + ", sample result = " + doubleResult + "\n\n\n");
                } else {
                    numCorrect++;
                }
            }
        } catch (Throwable e) {
            e.printStackTrace(System.out);
            numFailures++;
        }
    }
    System.out.println("test for .testEvaluateInterval(), " + numCorrect + " correct, " + numWrong + " wrong, " + numFailures + " failures, " + numTrials + " trials");
}

Example 12 with RealInterval

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

the class ConstraintContainerImpl method toECLiPSe.

/**
 * Insert the method's description here.
 * Creation date: (12/28/2004 5:47:43 PM)
 * @return java.lang.String
 */
public String toECLiPSe() {
    StringBuffer buffer = new StringBuffer();
    for (int i = 0; i < fieldSimpleBounds.length; i++) {
        if (getActive(fieldSimpleBounds[i]) == false) {
            continue;
        }
        String symbol = fieldSimpleBounds[i].getIdentifier();
        RealInterval bounds = fieldSimpleBounds[i].getBounds();
        String lowBoundsString = Double.toString(bounds.lo());
        if (bounds.lo() == Double.POSITIVE_INFINITY) {
            lowBoundsString = "1.0Inf";
        } else if (bounds.lo() == Double.NEGATIVE_INFINITY) {
            lowBoundsString = "-1.0Inf";
        }
        String hiBoundsString = Double.toString(bounds.hi());
        if (bounds.hi() == Double.POSITIVE_INFINITY) {
            hiBoundsString = "1.0Inf";
        } else if (bounds.hi() == Double.NEGATIVE_INFINITY) {
            hiBoundsString = "-1.0Inf";
        }
        buffer.append(cbit.vcell.parser.SymbolUtils.getEscapedTokenECLiPSe(symbol) + " $>= " + lowBoundsString + ",");
        buffer.append(cbit.vcell.parser.SymbolUtils.getEscapedTokenECLiPSe(symbol) + " $=< " + hiBoundsString + ",");
    }
    for (int i = 0; i < fieldGeneralConstraints.length; i++) {
        if (getActive(fieldGeneralConstraints[i]) == false) {
            continue;
        }
        buffer.append(fieldGeneralConstraints[i].getExpression().infix_ECLiPSe() + ", ");
    }
    return buffer.toString();
}

Example 13 with RealInterval

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

the class ConstraintSolver method resetIntervals.

/**
 * Insert the method's description here.
 * Creation date: (6/25/01 4:50:43 PM)
 */
public void resetIntervals() throws ExpressionBindingException {
    // 
    // allocate RealInterval array (one for each bound symbol)
    // 
    // if (fieldIntervals==null || fieldIntervals.length!=symbolList.size()){
    RealInterval[] intervals = new RealInterval[symbolList.size()];
    // 
    for (int i = 0; i < intervals.length; i++) {
        intervals[i] = RealInterval.fullInterval();
    }
    // 
    // pre-load with simple bounds on interval symbols.
    // 
    int symbolIndex = 0;
    Iterator<ConstraintSolver.Symbol> iter = symbolList.iterator();
    while (iter.hasNext()) {
        ConstraintSolver.Symbol symbol = iter.next();
        symbol.index = symbolIndex++;
        RealInterval bounds = constraintContainerImpl.getBounds(symbol.getName());
        intervals[symbol.index] = (RealInterval) bounds.clone();
    }
    // 
    for (int i = 0; i < expressionList.size(); i++) {
        ((Expression) expressionList.elementAt(i)).bindExpression(null);
        ((Expression) expressionList.elementAt(i)).bindExpression(this);
    }
    // 
    // reset "Consistency" flag in general constraints
    // 
    GeneralConstraint[] generalConstraints = constraintContainerImpl.getGeneralConstraints();
    for (int i = 0; i < generalConstraints.length; i++) {
        constraintContainerImpl.setConsistent(generalConstraints[i], true);
    }
    setIntervals(intervals);
}

Example 14 with RealInterval

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

the class SpeciesContextSpec method gatherIssues.

/**
 * Insert the method's description here.
 * Creation date: (11/1/2005 10:03:46 AM)
 * @param issueVector java.util.Vector
 */
public void gatherIssues(IssueContext issueContext, List<Issue> issueVector) {
    issueContext = issueContext.newChildContext(ContextType.SpeciesContextSpec, this);
    // 
    for (int i = 0; i < fieldParameters.length; i++) {
        RealInterval simpleBounds = parameterBounds[fieldParameters[i].getRole()];
        if (simpleBounds != null) {
            String parmName = fieldParameters[i].getNameScope().getName() + "." + fieldParameters[i].getName();
            issueVector.add(new SimpleBoundsIssue(fieldParameters[i], issueContext, simpleBounds, "parameter " + parmName + ": must be within " + simpleBounds.toString()));
        }
    }
    if (bForceContinuous && !bConstant && getSimulationContext().isStoch() && (getSimulationContext().getGeometry().getDimension() > 0)) {
        // if it's particle or constant we're good
        SpeciesContext sc = getSpeciesContext();
        ReactionContext rc = getSimulationContext().getReactionContext();
        ReactionSpec[] rsArray = rc.getReactionSpecs();
        for (ReactionSpec rs : rsArray) {
            if (!rs.isExcluded()) {
                // we only care about reactions which are not excluded
                // true if "this" is part of current reaction
                boolean iAmParticipant = false;
                // true if current reaction has at least a particle participant
                boolean haveParticle = false;
                ReactionStep step = rs.getReactionStep();
                for (ReactionParticipant p : step.getReactionParticipants()) {
                    if (p instanceof Product || p instanceof Reactant) {
                        SpeciesContextSpec candidate = rc.getSpeciesContextSpec(p.getSpeciesContext());
                        if (candidate == this) {
                            iAmParticipant = true;
                        } else if (!candidate.isForceContinuous() && !candidate.isConstant()) {
                            haveParticle = true;
                        }
                    }
                }
                if (iAmParticipant && haveParticle) {
                    String msg = "Continuous Species won't conserve mass in particle reaction " + rs.getReactionStep().getName() + ".";
                    String tip = "Mass conservation for reactions of binding between discrete and continuous species is handled approximately. <br>" + "To avoid any algorithmic approximation, which may produce undesired results, the user is advised to indicate <br>" + "the continuous species in those reactions as modifiers (i.e. 'catalysts') in the physiology.";
                    issueVector.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, tip, Issue.SEVERITY_WARNING));
                    // we issue warning as soon as we found the first reaction which satisfies criteria
                    break;
                }
            }
        }
    }
    if (!bForceContinuous && bConstant) {
        if (getSimulationContext().isStoch() && (getSimulationContext().getGeometry().getDimension() > 0)) {
            String msg = "Clamped Species must be continuous rather than particles.";
            String tip = "If choose 'clamped', must also choose 'forceContinuous'";
            issueVector.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, tip, Issue.SEVERITY_ERROR));
        }
    }
    if (bForceContinuous && !bConstant) {
        if (getSimulationContext().isStoch() && (getSimulationContext().getGeometry().getDimension() == 0)) {
            String msg = "Non-constant species is forced continuous, not supported for nonspatial stochastic applications.";
            issueVector.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.SEVERITY_ERROR));
        }
    }
}

Example 15 with RealInterval

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

the class ApplicationConstraintsGenerator method steadyStateFromApplication.

/**
 * Insert the method's description here.
 * Creation date: (6/26/01 8:25:55 AM)
 * @return cbit.vcell.constraints.ConstraintContainerImpl
 */
public static ConstraintContainerImpl steadyStateFromApplication(SimulationContext simContext, double tolerance) {
    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();
                double lowInitialValue = Math.min(initialValue / tolerance, initialValue * tolerance);
                double highInitialValue = Math.max(initialValue / tolerance, initialValue * tolerance);
                ccImpl.addSimpleBound(new SimpleBounds(speciesContexts[i].getName(), new RealInterval(lowInitialValue, highInitialValue), AbstractConstraint.MODELING_ASSUMPTION, "close to 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"));
            }
        }
        // 
        try {
            simContext.setMathDescription(simContext.createNewMathMapping().getMathDescription());
        } catch (Throwable e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("cannot create mathDescription");
        }
        MathDescription mathDesc = simContext.getMathDescription();
        if (mathDesc.getGeometry().getDimension() > 0) {
            throw new RuntimeException("spatial simulations not yet supported");
        }
        CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDesc.getSubDomains().nextElement();
        java.util.Enumeration<Equation> enumEquations = subDomain.getEquations();
        while (enumEquations.hasMoreElements()) {
            Equation equation = (Equation) enumEquations.nextElement();
            Expression rateConstraintExp = new Expression(equation.getRateExpression().infix() + "==0");
            rateConstraintExp = getSteadyStateExpression(rateConstraintExp);
            if (!rateConstraintExp.compareEqual(new Expression(1.0))) {
                // not a trivial constraint (always true)
                ccImpl.addGeneralConstraint(new GeneralConstraint(rateConstraintExp, AbstractConstraint.PHYSICAL_LIMIT, "definition of steady state"));
            }
        }
        // 
        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();
                        double lowValue = Math.min(constantValue / tolerance, constantValue * tolerance);
                        double highValue = Math.max(constantValue / tolerance, constantValue * tolerance);
                        RealInterval interval = new RealInterval(lowValue, highValue);
                        ccImpl.addSimpleBound(new SimpleBounds(parameters[j].getName(), interval, AbstractConstraint.MODELING_ASSUMPTION, "parameter close to model default"));
                    } 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());
                    ccImpl.addGeneralConstraint(new GeneralConstraint(getSteadyStateExpression(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"));
        // 
        // add K_fluxs
        // 
        java.util.Enumeration<Variable> enumVars = mathDesc.getVariables();
        while (enumVars.hasMoreElements()) {
            Variable var = (Variable) enumVars.nextElement();
            if (var.getName().startsWith("Kflux_") && var instanceof Function) {
                Expression kfluxExp = new Expression(((Function) var).getExpression());
                kfluxExp.bindExpression(mathDesc);
                kfluxExp = MathUtilities.substituteFunctions(kfluxExp, mathDesc);
                kfluxExp = kfluxExp.flatten();
                ccImpl.addSimpleBound(new SimpleBounds(var.getName(), new RealInterval(kfluxExp.evaluateConstant()), AbstractConstraint.MODELING_ASSUMPTION, "flux conversion factor"));
            }
        }
        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;
    }
}