Examples with Matchable org.vcell.util.Matchable used on opensource projects

Example 1 with Matchable

use of org.vcell.util.Matchable in project vcell by virtualcell.

the class MassActionSolver method solveMassAction.

public static MassActionFunction solveMassAction(Parameter optionalForwardRateParameter, Parameter optionalReverseRateParameter, Expression orgExp, ReactionStep rs) throws ExpressionException, ModelException, DivideByZeroException {
    MassActionFunction maFunc = new MassActionFunction();
    // get reactants, products, overlaps, non-overlap reactants and non-overlap products
    ArrayList<ReactionParticipant> reactants = new ArrayList<ReactionParticipant>();
    ArrayList<ReactionParticipant> products = new ArrayList<ReactionParticipant>();
    ReactionParticipant[] rp = rs.getReactionParticipants();
    // should use this one to compare functional equavalent since this duplicatedExp has all params substituted.
    Expression duplicatedExp = substituteParameters(orgExp, false);
    // separate the reactants and products, fluxes, catalysts
    String rxnName = rs.getName();
    // reaction with membrane current can not be transformed to mass action
    if (rs.getPhysicsOptions() == ReactionStep.PHYSICS_MOLECULAR_AND_ELECTRICAL || rs.getPhysicsOptions() == ReactionStep.PHYSICS_ELECTRICAL_ONLY) {
        throw new ModelException("Kinetics of reaction " + rxnName + " has membrane current. It can not be automatically transformed to Mass Action kinetics.");
    }
    for (int i = 0; i < rp.length; i++) {
        if (rp[i] instanceof Reactant) {
            reactants.add(rp[i]);
        } else if (rp[i] instanceof Product) {
            products.add(rp[i]);
        } else if (rp[i] instanceof Catalyst) {
            String catalystName = rp[i].getSpeciesContext().getName();
            // only if duplictedExp is not a non-linear function of catalystName.
            if (duplicatedExp.hasSymbol(catalystName)) {
                // Only when catalyst appears in ReactionRate, we add catalyst to both reactant and product
                // the stoichiometry should be set to 1.
                ReactionParticipant catalystRP = new ReactionParticipant(null, rs, rp[i].getSpeciesContext(), 1) {

                    public boolean compareEqual(Matchable obj) {
                        ReactionParticipant rp = (ReactionParticipant) obj;
                        if (rp == null) {
                            return false;
                        }
                        if (!Compare.isEqual(getSpecies(), rp.getSpecies())) {
                            return false;
                        }
                        if (!Compare.isEqual(getStructure(), rp.getStructure())) {
                            return false;
                        }
                        if (getStoichiometry() != rp.getStoichiometry()) {
                            return false;
                        }
                        return true;
                    }

                    @Override
                    public void writeTokens(PrintWriter pw) {
                    }

                    @Override
                    public void fromTokens(CommentStringTokenizer tokens, Model model) throws Exception {
                    }
                };
                products.add(catalystRP);
                reactants.add(catalystRP);
            }
        }
    }
    /**
     * The code below is going to solve reaction with kinetics that are NOT Massaction. Or Massaction with catalysts involved.
     */
    // 
    // 2x2 rational matrix
    // 
    // lets define
    // J() is the substituted total rate expression
    // P() as the theoretical "product" term with Kf = 1
    // R() as the theoretical "reactant" term with Kr = 1
    // 
    // Kf * R([1 1 1])  - Kr * P([1 1 1])  = J([1 1 1])
    // Kf * R([2 3 4])  - Kr * P([2 3 4])  = J([2 3 4])
    // 
    // in matrix form,
    // 
    // |                                        |
    // | R([1 1 1])   -P([1 1 1])    J([1 1 1]) |
    // | R([2 3 4])   -P([2 3 4])    J([2 3 4]) |
    // |                                        |
    // 
    // 
    // example: Kf*A*B^2*C - Kr*C*A
    // J() = forwardRate*43/Kabc*A*B^2*C - (myC*5-2)*C*A
    // P() = C*A
    // R() = A*B^2*C
    // 
    // |                                              |
    // | R([1 1 1])  -P([1 1 1])    J([1 1 1])        |
    // | R([2 3 4])  -P([2 3 4])    J([2 3 4])        |
    // |                                              |
    // 
    // |                                                                              |
    // | R([1 1 1])*R([2 3 4])    -P([1 1 1])*R([2 3 4])     J([1 1 1])*R([2 3 4])    |
    // | R([2 3 4])*R([1 1 1])    -P([2 3 4])*R([1 1 1])     J([2 3 4])*R([1 1 1])    |
    // |                                                                              |
    // 
    // 
    // |                                                                                                                    |
    // | R([1 1 1])*R([2 3 4])  -P([1 1 1])*R([2 3 4])                         J([1 1 1])*R([2 3 4])                        |
    // | 0                      -P([2 3 4])*R([1 1 1])+P([1 1 1])*R([2 3 4])   J([2 3 4])*R([1 1 1])-J([1 1 1])*R([2 3 4])  |
    // |                                                                                                                    |
    // 
    // 
    // 
    // 
    Expression forwardExp = null;
    Expression reverseExp = null;
    Expression R_exp = new Expression(1);
    if (reactants.size() > 0) {
        R_exp = new Expression(1.0);
        for (ReactionParticipant reactant : reactants) {
            R_exp = Expression.mult(R_exp, Expression.power(new Expression(reactant.getName()), new Expression(reactant.getStoichiometry())));
        }
    }
    Expression P_exp = new Expression(1);
    if (products.size() > 0) {
        P_exp = new Expression(1.0);
        for (ReactionParticipant product : products) {
            P_exp = Expression.mult(P_exp, Expression.power(new Expression(product.getName()), new Expression(product.getStoichiometry())));
        }
    }
    HashSet<String> reactionParticipantNames = new HashSet<String>();
    for (ReactionParticipant reactionParticipant : rs.getReactionParticipants()) {
        reactionParticipantNames.add(reactionParticipant.getName());
    }
    Expression R_1 = new Expression(R_exp);
    Expression R_2 = new Expression(R_exp);
    Expression P_1 = new Expression(P_exp);
    Expression P_2 = new Expression(P_exp);
    Expression J_1 = new Expression(duplicatedExp);
    Expression J_2 = new Expression(duplicatedExp);
    int index = 0;
    for (String rpName : reactionParticipantNames) {
        R_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
        R_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
        P_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
        P_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
        J_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
        J_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
        index++;
    }
    R_1 = R_1.flatten();
    R_2 = R_2.flatten();
    P_1 = P_1.flatten();
    P_2 = P_2.flatten();
    J_1 = J_1.flatten();
    J_2 = J_2.flatten();
    // e.g. A+B <-> A+B, A+2B <-> A+2B
    if (ExpressionUtils.functionallyEquivalent(R_exp, P_exp, false, 1e-8, 1e-8)) {
        throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Identical reactants and products not supported for stochastic models."));
    }
    if (R_exp.compareEqual(new Expression(1)) && P_exp.compareEqual(new Expression(1))) {
        // no reactants, no products ... nothing to do
        forwardExp = null;
        reverseExp = null;
    } else {
        // both reactants and products
        RationalExpMatrix matrix = new RationalExpMatrix(2, 3);
        matrix.set_elem(0, 0, RationalExpUtils.getRationalExp(R_1, true));
        matrix.set_elem(0, 1, RationalExpUtils.getRationalExp(Expression.negate(P_1), true));
        matrix.set_elem(0, 2, RationalExpUtils.getRationalExp(J_1, true));
        matrix.set_elem(1, 0, RationalExpUtils.getRationalExp(R_2, true));
        matrix.set_elem(1, 1, RationalExpUtils.getRationalExp(Expression.negate(P_2), true));
        matrix.set_elem(1, 2, RationalExpUtils.getRationalExp(J_2, true));
        RationalExp[] solution;
        try {
            // matrix.show();
            solution = matrix.solveLinearExpressions();
            // solution[0] is forward rate.
            solution[0] = solution[0].simplify();
            // solution[1] is reverse rate.
            solution[1] = solution[1].simplify();
        } catch (MatrixException e) {
            e.printStackTrace();
            throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "MatrixException: " + e.getMessage()));
        }
        forwardExp = new Expression(solution[0].infixString());
        reverseExp = new Expression(solution[1].infixString());
        // for massAction, if there is no reactant in the forward rate, the forwardExp should be set to null to avoid writing jump process for the forward reaction.
        if (R_exp.compareEqual(new Expression(1)) && rs.getKinetics().getKineticsDescription().equals(KineticsDescription.MassAction)) {
            forwardExp = null;
        }
        // for massAction, if there is no product in the reverse rate, the reverseExp should be set to null to avoid writing jump process for the reverse reaction.
        if (P_exp.compareEqual(new Expression(1)) && rs.getKinetics().getKineticsDescription().equals(KineticsDescription.MassAction)) {
            reverseExp = null;
        }
    }
    if (forwardExp != null) {
        forwardExp.bindExpression(rs);
    }
    if (reverseExp != null) {
        reverseExp.bindExpression(rs);
    }
    // Reconstruct the rate based on the extracted forward rate and reverse rate. If the reconstructed rate is not equivalent to the original rate,
    // it means the original rate is not in the form of Kf*r1^n1*r2^n2-Kr*p1^m1*p2^m2.
    Expression constructedExp = reconstructedRate(forwardExp, reverseExp, reactants, products, rs.getNameScope());
    Expression orgExp_withoutCatalyst = removeCatalystFromExp(duplicatedExp, rs);
    Expression constructedExp_withoutCatalyst = removeCatalystFromExp(constructedExp, rs);
    if (!ExpressionUtils.functionallyEquivalent(orgExp_withoutCatalyst, constructedExp_withoutCatalyst, false, 1e-8, 1e-8)) {
        throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Mathmatical form incompatible with mass action."));
    }
    // check if forward rate constant and reverse rate constant both can be evaluated to constants(numbers) after substituting all parameters.
    if (forwardExp != null) {
        Expression forwardExpCopy = new Expression(forwardExp);
        try {
            substituteParameters(forwardExpCopy, true).evaluateConstant();
        } catch (ExpressionException e) {
            throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Problem in forward rate '" + forwardExp.infix() + "', " + e.getMessage()));
        }
        // 
        if (optionalForwardRateParameter != null) {
            Expression forwardRateParameterCopy = new Expression(optionalForwardRateParameter, optionalForwardRateParameter.getNameScope());
            try {
                substituteParameters(forwardRateParameterCopy, true).evaluateConstant();
                if (forwardExpCopy.compareEqual(forwardRateParameterCopy)) {
                    forwardExp = new Expression(optionalForwardRateParameter, optionalForwardRateParameter.getNameScope());
                }
            } catch (ExpressionException e) {
                // not expecting a problem because forwardExpCopy didn't have a problem, but in any case it is ok to swallow this exception
                e.printStackTrace(System.out);
            }
        }
    }
    if (reverseExp != null) {
        Expression reverseExpCopy = new Expression(reverseExp);
        try {
            substituteParameters(reverseExpCopy, true).evaluateConstant();
        } catch (ExpressionException e) {
            throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Problem in reverse rate '" + reverseExp.infix() + "', " + e.getMessage()));
        }
        // 
        if (optionalReverseRateParameter != null) {
            Expression reverseRateParameterCopy = new Expression(optionalReverseRateParameter, optionalReverseRateParameter.getNameScope());
            try {
                substituteParameters(reverseRateParameterCopy, true).evaluateConstant();
                if (reverseExpCopy.compareEqual(reverseRateParameterCopy)) {
                    reverseExp = new Expression(optionalReverseRateParameter, optionalReverseRateParameter.getNameScope());
                }
            } catch (ExpressionException e) {
                // not expecting a problem because reverseExpCopy didn't have a problem, but in any case it is ok to swallow this exception
                e.printStackTrace(System.out);
            }
        }
    }
    maFunc.setForwardRate(forwardExp);
    maFunc.setReverseRate(reverseExp);
    maFunc.setReactants(reactants);
    maFunc.setProducts(products);
    return maFunc;
}

Example 2 with Matchable

use of org.vcell.util.Matchable in project vcell by virtualcell.

the class GeometrySurfaceDescription method propertyChange.

/**
 * This method gets called when a bound property is changed.
 * @param evt A PropertyChangeEvent object describing the event source
 *   	and the property that has changed.
 */
public void propertyChange(java.beans.PropertyChangeEvent evt) {
    if (evt.getSource() == this && evt.getPropertyName().equals("volumeSampleSize")) {
        ISize oldValue = (ISize) evt.getOldValue();
        ISize newValue = (ISize) evt.getNewValue();
        if (!oldValue.compareEqual(newValue)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() == this && evt.getPropertyName().equals("filterCutoffFrequency")) {
        Double oldValue = (Double) evt.getOldValue();
        Double newValue = (Double) evt.getNewValue();
        if (!oldValue.equals(newValue)) {
            try {
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() == this && evt.getPropertyName().equals(PROPERTY_NAME_GEOMETRIC_REGIONS)) {
        refreshSurfaceClasses();
    }
    if (evt.getSource() == getGeometry().getGeometrySpec() && (evt.getPropertyName().equals("extent") || evt.getPropertyName().equals("origin"))) {
        Matchable oldExtentOrOrigin = (Matchable) evt.getOldValue();
        Matchable newExtentOrOrigin = (Matchable) evt.getNewValue();
        if (!Compare.isEqual(oldExtentOrOrigin, newExtentOrOrigin)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() instanceof AnalyticSubVolume && evt.getPropertyName().equals("expression")) {
        Expression oldExpression = (Expression) evt.getOldValue();
        Expression newExpression = (Expression) evt.getNewValue();
        if (!Compare.isEqual(oldExpression, newExpression)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() instanceof CSGObject && evt.getPropertyName().equals(CSGObject.PROPERTY_NAME_ROOT)) {
        try {
            // nobody listens to this, updateAll() will propagate changes
            getRegionImage0().setDirty();
            getSurfaceCollection0().setDirty();
            fieldGeometricRegions.setDirty();
        } catch (Exception e) {
            e.printStackTrace(System.out);
        }
    }
    if (evt.getSource() instanceof SubVolume && evt.getPropertyName().equals("name")) {
        String oldName = (String) evt.getOldValue();
        String newName = (String) evt.getNewValue();
        if (!Compare.isEqual(oldName, newName)) {
            try {
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (evt.getSource() == getGeometry().getGeometrySpec() && evt.getPropertyName().equals("subVolumes")) {
        SubVolume[] oldValue = (SubVolume[]) evt.getOldValue();
        SubVolume[] newValue = (SubVolume[]) evt.getNewValue();
        // 
        for (int i = 0; oldValue != null && i < oldValue.length; i++) {
            oldValue[i].removePropertyChangeListener(this);
        }
        for (int i = 0; newValue != null && i < newValue.length; i++) {
            newValue[i].addPropertyChangeListener(this);
        }
        // 
        if (oldValue == null || newValue == null || !Compare.isEqualStrict(oldValue, newValue)) {
            try {
                // nobody listens to this, updateAll() will propagate changes
                getRegionImage0().setDirty();
                getSurfaceCollection0().setDirty();
                fieldGeometricRegions.setDirty();
            } catch (Exception e) {
                e.printStackTrace(System.out);
            }
        } else if (fieldGeometricRegions.getCurrentValue() != null && oldValue != newValue) {
            // 
            for (int i = 0; i < newValue.length; i++) {
                SubVolume subVolume = newValue[i];
                for (int j = 0; j < fieldGeometricRegions.getCurrentValue().length; j++) {
                    if (fieldGeometricRegions.getCurrentValue()[j] instanceof VolumeGeometricRegion) {
                        VolumeGeometricRegion volumeRegion = (VolumeGeometricRegion) fieldGeometricRegions.getCurrentValue()[j];
                        // 
                        if (volumeRegion.getSubVolume().compareEqual(subVolume) && volumeRegion.getSubVolume() != subVolume) {
                            volumeRegion.setSubVolume(subVolume);
                        }
                    }
                }
            }
        }
    }
}

Example 3 with Matchable

use of org.vcell.util.Matchable in project vcell by virtualcell.

the class MovingBoundaryFileWriter method flattenExpression.

private Expression flattenExpression(Expression ex, VariableDomain varDomain) throws ExpressionException, MathException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    Variable normalX = new Variable("normalX", null) {

        public boolean compareEqual(Matchable object, boolean bIgnoreMissingDomains) {
            return false;
        }

        public String getVCML() throws MathException {
            return null;
        }
    };
    Variable normalY = new Variable("normalY", null) {

        public boolean compareEqual(Matchable object, boolean bIgnoreMissingDomains) {
            return false;
        }

        public String getVCML() throws MathException {
            return null;
        }
    };
    SymbolTable augmentedSymbolTable = new SymbolTable() {

        @Override
        public SymbolTableEntry getEntry(String identifierString) {
            if (identifierString.equals(normalX.getName())) {
                return normalX;
            }
            if (identifierString.equals(normalY.getName())) {
                return normalY;
            }
            return simSymbolTable.getEntry(identifierString);
        }

        @Override
        public void getEntries(Map<String, SymbolTableEntry> entryMap) {
            simSymbolTable.getEntries(entryMap);
            entryMap.put(normalX.getName(), normalX);
            entryMap.put(normalY.getName(), normalY);
        }
    };
    ex = new Expression(ex);
    ex.bindExpression(augmentedSymbolTable);
    Expression flattended = MathUtilities.substituteFunctions(ex, augmentedSymbolTable).flatten();
    Expression substituted = SolverUtilities.substituteSizeAndNormalFunctions(flattended, varDomain).flatten();
    return substituted;
// return simSymbolTable.substituteFunctions(ex).flatten().infix();
}

Example 4 with Matchable

use of org.vcell.util.Matchable in project vcell by virtualcell.

the class GeometrySpec method propertyChange.

/**
 * Insert the method's description here.
 * Creation date: (6/3/00 9:58:08 AM)
 * @param event java.beans.PropertyChangeEvent
 */
public void propertyChange(java.beans.PropertyChangeEvent event) {
    if (event.getSource() == this && event.getPropertyName().equals("subVolumes")) {
        SubVolume[] oldSubVolumes = (SubVolume[]) event.getOldValue();
        SubVolume[] newSubVolumes = (SubVolume[]) event.getNewValue();
        if (!Compare.isEqualStrict(oldSubVolumes, newSubVolumes)) {
            // ignore if just a change of instances
            getSampledImage().setDirty();
            getThumbnailImage().setDirty();
        }
    }
    if (event.getSource() == this && (event.getPropertyName().equals("extent") || event.getPropertyName().equals("origin"))) {
        Matchable oldExtentOrOrigin = (Matchable) event.getOldValue();
        Matchable newExtentOrOrigin = (Matchable) event.getNewValue();
        if (!Compare.isEqual(oldExtentOrOrigin, newExtentOrOrigin)) {
            getSampledImage().setDirty();
            getThumbnailImage().setDirty();
        }
    }
    if (event.getSource() instanceof AnalyticSubVolume && event.getPropertyName().equals("expression")) {
        Expression oldExpression = (Expression) event.getOldValue();
        Expression newExpression = (Expression) event.getNewValue();
        if (!Compare.isEqual(oldExpression, newExpression)) {
            getSampledImage().setDirty();
            getThumbnailImage().setDirty();
        }
    }
    if (event.getSource() instanceof CSGObject && event.getPropertyName().equals(CSGObject.PROPERTY_NAME_ROOT)) {
        getSampledImage().setDirty();
        getThumbnailImage().setDirty();
    }
}