Examples with RationalExp cbit.vcell.matrix.RationalExp used on opensource projects

Example 11 with RationalExp

use of cbit.vcell.matrix.RationalExp 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();
        fi.getFunction().bindExpression(this);
        PseudoConstant pc = new PseudoConstant(getAvailablePseudoConstantName(), fi.getFunction(), MathMapping_4_8.nullDomain);
        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);
    }
}

Example 12 with RationalExp

use of cbit.vcell.matrix.RationalExp in project vcell by virtualcell.

the class StochMathMapping_4_8 method getProbabilityRate.

/**
 * Get probability expression for the specific elementary reaction.
 * Input: ReactionStep, the reaction. isForwardDirection, if the elementary reaction is forward from the reactionstep.
 * Output: Expression. the probability expression.
 * Creation date: (9/14/2006 3:22:58 PM)
 */
Expression getProbabilityRate(ReactionStep rs, boolean isForwardDirection) throws MappingException {
    ReactionStep reactionStep = rs;
    Expression probExp = null;
    // get kinetics of the reaction step
    Kinetics kinetics = reactionStep.getKinetics();
    // to compose the rate constant expression e.g. Kf, Kr
    Expression rateConstantExpr = null;
    // to compose the stochastic variable(species) expression, e.g. s*(s-1)*(s-2)* speciesFactor.
    Expression rxnProbabilityExpr = null;
    // to compose the factor that the probability expression multiplies with, which convert the rate expression under stochastic context
    Expression factorExpr = null;
    // the structure where reaction happens
    StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure());
    Model model = getSimulationContext().getModel();
    try {
        if (// forward reaction
        isForwardDirection) {
            // for HMMs, it's a bit complicated. Vmax/(Km+s)-->Vmax*Size_s/(Km*Size_s+Ns)
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                KineticsParameter kfp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                rateConstantExpr = new Expression(kfp, getNameScope());
            // rateConstantExpr.bindExpression(this);
            }
            // get convert factor for rate constant( membrane:rateConstant*membrane_Size (factor is membrane_size), feature : rateConstant*(feature_size/KMole)(factor is feature_size/KMOLE)) )
            if (sm.getStructure() instanceof Membrane) {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
            } else {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
                Expression kmoleExpr = new Expression(1.0 / 602.0);
                factorExpr = Expression.mult(factorExpr, kmoleExpr);
            }
            // complete the probability expression by the reactants' stoichiometries if it is Mass Action rate law
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int i = 0; i < reacPart.length; i++) {
                    int stoichiometry = 0;
                    if (reacPart[i] instanceof Reactant) {
                        stoichiometry = ((Reactant) reacPart[i]).getStoichiometry();
                        // ******the following part is to form the s*(s-1)(s-2)..(s-stoi+1).portion of the probability rate.
                        StructureMapping reactSM = getSimulationContext().getGeometryContext().getStructureMapping(reacPart[i].getStructure());
                        // factor expression for species
                        Expression speciesFactor = null;
                        // convert speceis' unit from moles/liter to molecules.
                        if (reactSM.getStructure() instanceof Membrane) {
                            speciesFactor = Expression.invert(new Expression(reactSM.getStructure().getStructureSize(), getNameScope()));
                        } else {
                            Expression exp1 = new Expression(1.0 / 602.0);
                            Expression exp2 = new Expression(reactSM.getStructure().getStructureSize(), getNameScope());
                            speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                        }
                        // s*(s-1)(s-2)..(s-stoi+1)
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[i].getSpeciesContext());
                        Expression spCount_exp = new Expression(spCountParam, getNameScope());
                        // species from uM to No. of Particles, form s*(s-1)*(s-2)
                        Expression tempExpr = new Expression(spCount_exp);
                        for (int j = 1; j < stoichiometry; j++) {
                            tempExpr = Expression.mult(tempExpr, Expression.add(spCount_exp, new Expression(-j)));
                        }
                        // update total factor with speceies factor
                        if (stoichiometry == 1) {
                            factorExpr = Expression.mult(factorExpr, speciesFactor);
                        } else if (stoichiometry > 1) {
                            // rxnProbExpr * (structSize^stoichiometry)
                            Expression powerExpr = Expression.power(speciesFactor, new Expression(stoichiometry));
                            factorExpr = Expression.mult(factorExpr, powerExpr);
                        }
                        if (rxnProbabilityExpr == null) {
                            rxnProbabilityExpr = new Expression(tempExpr);
                        } else {
                            // for more than one reactant
                            rxnProbabilityExpr = Expression.mult(rxnProbabilityExpr, tempExpr);
                        }
                    }
                }
            }
        } else // reverse reaction
        {
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                KineticsParameter krp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                rateConstantExpr = new Expression(krp, getNameScope());
            // rateConstantExpr.bindExpression(this);
            }
            // get convert factor for rate constant( membrane:rateConstant*membrane_Size (factor is membrane_size), feature : rateConstant*(feature_size/KMole)(factor is feature_size/KMOLE)) )
            if (sm.getStructure() instanceof Membrane) {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
            } else {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
                Expression exp = new Expression(1.0 / 602.0);
                factorExpr = Expression.mult(factorExpr, exp);
            }
            // complete the remaining part of the probability expression by the products' stoichiometries.
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int i = 0; i < reacPart.length; i++) {
                    int stoichiometry = 0;
                    if (reacPart[i] instanceof Product) {
                        stoichiometry = ((Product) reacPart[i]).getStoichiometry();
                        // ******the following part is to form the s*(s-1)*(s-2)...(s-stoi+1).portion of the probability rate.
                        StructureMapping reactSM = getSimulationContext().getGeometryContext().getStructureMapping(reacPart[i].getStructure());
                        // factor expression for species
                        Expression speciesFactor = null;
                        // convert speceis' unit from moles/liter to molecules.
                        if (reactSM.getStructure() instanceof Membrane) {
                            speciesFactor = Expression.invert(new Expression(reactSM.getStructure().getStructureSize(), getNameScope()));
                        } else {
                            Expression exp1 = new Expression(1.0 / 602.0);
                            Expression exp2 = new Expression(reactSM.getStructure().getStructureSize(), getNameScope());
                            speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                        }
                        // s*(s-1)*(s-2)...(s-stoi+1)
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[i].getSpeciesContext());
                        Expression spCount_exp = new Expression(spCountParam, getNameScope());
                        // species from uM to No. of Particles, form s*(s-1)*(s-2)
                        Expression tempExpr = new Expression(spCount_exp);
                        for (int j = 1; j < stoichiometry; j++) {
                            tempExpr = Expression.mult(tempExpr, Expression.add(spCount_exp, new Expression(-j)));
                        }
                        // update total factor with speceies factor
                        if (stoichiometry == 1) {
                            factorExpr = Expression.mult(factorExpr, speciesFactor);
                        } else if (stoichiometry > 1) {
                            // rxnProbExpr * (structSize^stoichiometry)
                            Expression powerExpr = Expression.power(speciesFactor, new Expression(stoichiometry));
                            factorExpr = Expression.mult(factorExpr, powerExpr);
                        }
                        if (rxnProbabilityExpr == null) {
                            rxnProbabilityExpr = new Expression(tempExpr);
                        } else {
                            rxnProbabilityExpr = Expression.mult(rxnProbabilityExpr, tempExpr);
                        }
                    }
                }
            }
        }
        // Now construct the probability expression.
        if (rateConstantExpr == null) {
            throw new MappingException("Can not find reaction rate constant in reaction: " + reactionStep.getName());
        } else if (rxnProbabilityExpr == null) {
            probExp = new Expression(rateConstantExpr);
        } else if ((rateConstantExpr != null) && (rxnProbabilityExpr != null)) {
            probExp = Expression.mult(rateConstantExpr, rxnProbabilityExpr);
        }
        // simplify the factor
        RationalExp factorRatExp = RationalExpUtils.getRationalExp(factorExpr);
        factorExpr = new Expression(factorRatExp.infixString());
        factorExpr.bindExpression(this);
        // get probability rate with converting factor
        probExp = Expression.mult(probExp, factorExpr);
        probExp = probExp.flatten();
    // //
    // // round trip to rational expression for simplifying terms like KMOLE/KMOLE ...
    // // we don't want to loose the symbol binding ... so we make a temporary symbolTable from the original binding.
    // //
    // final Expression finalExp = new Expression(probExp);
    // SymbolTable symbolTable = new SymbolTable(){
    // public void getEntries(Map<String, SymbolTableEntry> entryMap) {
    // throw new RuntimeException("should not be called");
    // }
    // public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
    // return finalExp.getSymbolBinding(identifierString);
    // }
    // };
    // cbit.vcell.matrix.RationalExp ratExp = cbit.vcell.parser.RationalExpUtils.getRationalExp(probExp);
    // probExp = new Expression(ratExp.infixString());
    // probExp.bindExpression(symbolTable);
    } catch (ExpressionException e) {
        e.printStackTrace();
    }
    return probExp;
}

Example 13 with RationalExp

use of cbit.vcell.matrix.RationalExp in project vcell by virtualcell.

the class VCUnitDefinition method raiseTo.

public VCUnitDefinition raiseTo(ucar.units_vcell.RationalNumber rationalNumber) {
    VCUnitDefinition powerUnit = null;
    if (powerMap != null) {
        powerUnit = powerMap.get(rationalNumber);
        if (powerUnit != null) {
            return powerUnit;
        }
    } else {
        powerMap = new HashMap<ucar.units_vcell.RationalNumber, VCUnitDefinition>();
    }
    RationalExp origRationalExpWithoutNumericScale = getSymbolRationalExpWithoutNumericScale(fieldUnitSymbol);
    RationalExp rationalExpWithoutNumericScale = origRationalExpWithoutNumericScale;
    // 
    if (rationalNumber.doubleValue() < 0) {
        rationalExpWithoutNumericScale = rationalExpWithoutNumericScale.inverse();
        rationalNumber = rationalNumber.minus();
    }
    // 
    if (rationalNumber.intValue() == rationalNumber.doubleValue()) {
        for (int i = 1; i < Math.abs(rationalNumber.intValue()); i++) {
            rationalExpWithoutNumericScale = rationalExpWithoutNumericScale.mult(origRationalExpWithoutNumericScale);
        }
        double newNumericScale = Math.pow(fieldUnitSymbol.getNumericScale(), rationalNumber.doubleValue());
        powerUnit = getUnit(newNumericScale, rationalExpWithoutNumericScale);
        powerMap.put(rationalNumber, powerUnit);
        return powerUnit;
    } else {
        throw new RuntimeException("raiseTo( non-integer ) not yet supported");
    }
}

Example 14 with RationalExp

use of cbit.vcell.matrix.RationalExp in project vcell by virtualcell.

the class VCUnitDefinition method divideBy.

public VCUnitDefinition divideBy(VCUnitDefinition childUnit) {
    VCUnitDefinition ratioUnit = null;
    if (divideByMap != null) {
        ratioUnit = divideByMap.get(childUnit);
        if (ratioUnit != null) {
            return ratioUnit;
        }
    } else {
        divideByMap = new HashMap<VCUnitDefinition, VCUnitDefinition>();
    }
    double newNumericScale = fieldUnitSymbol.getNumericScale() / childUnit.fieldUnitSymbol.getNumericScale();
    RationalExp thisNonnumericUnit = getSymbolRationalExpWithoutNumericScale(this.fieldUnitSymbol);
    RationalExp otherNonnumericUnit = getSymbolRationalExpWithoutNumericScale(childUnit.fieldUnitSymbol);
    ratioUnit = getUnit(newNumericScale, thisNonnumericUnit.div(otherNonnumericUnit));
    divideByMap.put(childUnit, ratioUnit);
    return ratioUnit;
}

Example 15 with RationalExp

use of cbit.vcell.matrix.RationalExp in project vcell by virtualcell.

the class VCUnitDefinition method getSymbolRationalExpWithoutNumericScale.

private static RationalExp getSymbolRationalExpWithoutNumericScale(UnitSymbol unitSymbol) {
    String unitString = unitSymbol.getUnitSymbolAsInfixWithoutFloatScale();
    RationalExp rationalExp;
    try {
        rationalExp = RationalExpUtils.getRationalExp(new cbit.vcell.parser.Expression(unitString));
    } catch (ExpressionException e) {
        e.printStackTrace();
        throw new RuntimeException("unit exception " + e.getMessage());
    }
    return rationalExp;
}