source: trunk/Compiler/ModelicaMiddleEnd/src/jastadd/structural/Symbolic.jrag @ 13636

Last change on this file since 13636 was 13636, checked in by Jonathan Kämpe, 2 months ago

#5844 Merging to trunk. Fixes related to attribute evaluation robustness.

File size: 21.2 KB
Line 
1/*
2    Copyright (C) 2009 Modelon AB
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, version 3 of the License.
7
8    This program is distributed in the hope that it will be useful,
9    but WITHOUT ANY WARRANTY; without even the implied warranty of
10    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11    GNU General Public License for more details.
12
13    You should have received a copy of the GNU General Public License
14    along with this program.  If not, see <http://www.gnu.org/licenses/>.
15*/
16
17import java.util.Set;
18
19import org.jmodelica.util.Solvability;
20
21aspect EquationSolve {
22
23        syn Solvability FAbstractEquation.isSolved(FVariable var) = isSolved(var, false);
24       
25        syn Solvability FAbstractEquation.isSolved(FVariable var, boolean duringTearing) = Solvability.UNSOLVABLE;
26       
27        eq FEquation.isSolved(FVariable var, boolean duringTearing) {
28                boolean solvable = !(solution(var, duringTearing) instanceof FNoExp);
29                if (solvable)
30                        return Solvability.ANALYTICALLY_SOLVABLE;
31                else if (duringTearing && canUseLocalIteration())
32                        return Solvability.NUMERICALLY_SOLVABLE;
33                else
34                        return Solvability.UNSOLVABLE;
35        }
36
37        eq FFunctionCallEquation.isSolved(FVariable var, boolean duringTearing) {
38                Set<FVariable> leftVars = referencedFVariablesInLHS();
39               
40                boolean inLeft = false;
41                for (FVariable fv : leftVars) {
42                        if (fv == var)
43                                inLeft = true;
44                }
45               
46                Set<FVariable> argVars = referencedFVariablesInRHS();
47
48                boolean inArgs = false;
49                for (FVariable fv : argVars) {
50                        if (fv == var)
51                                inArgs = true;
52                }
53
54                if (inLeft && !inArgs)
55                        return Solvability.ANALYTICALLY_SOLVABLE;
56                else if (leftVars.size() == 1 && duringTearing && canUseLocalIteration())
57                        return Solvability.NUMERICALLY_SOLVABLE;
58                else
59                        return Solvability.UNSOLVABLE;
60        }
61       
62        eq FAlgorithm.isSolved(FVariable var, boolean duringTearing) {
63                Set<FVariable> leftVars = referencedFVariablesInLHS();
64               
65                boolean inLeft = false;
66                for (FVariable fv : leftVars) {
67                        if (fv == var)
68                                inLeft = true;
69                }
70               
71                if (inLeft)
72                        return Solvability.ANALYTICALLY_SOLVABLE;
73                else if (leftVars.size() == 1 && duringTearing && canUseLocalIteration())
74                        return Solvability.NUMERICALLY_SOLVABLE;
75                else
76                        return Solvability.UNSOLVABLE;
77        }
78
79    eq FIfWhenElseEquation.isSolved(FVariable var, boolean duringTearing) {
80        Set<FVariable> leftVars = referencedFVariablesInLHS();
81       
82        boolean inLeft = false;
83        for (FVariable fv : leftVars) {
84            if (fv == var)
85                inLeft = true;
86        }
87       
88        Set<FVariable> argVars = referencedFVariablesInRHS();
89
90        boolean inArgs = false;
91        for (FVariable fv : argVars) {
92            if (fv == var)
93                inArgs = true;
94        }
95
96        if (inLeft && !inArgs)
97            return Solvability.ANALYTICALLY_SOLVABLE;
98        else
99            return Solvability.UNSOLVABLE;
100    }
101   
102        syn FExp FEquation.solution(FVariable var) {
103                return solution(var, false);
104        }
105
106    syn nta FExp FEquation.solution(FVariable var, boolean duringTearing) {
107        try {
108            // Get terms
109            ArrayList<FExp> t = terms();
110           
111            ArrayList<FExp> activeTerms = new ArrayList<FExp>();
112            ArrayList<FExp> inactiveTerms = new ArrayList<FExp>();
113            //Find terms
114            for (FExp e : t) {
115                if (e.nbrUses(var.name())==1) {
116                    activeTerms.add(e);
117                } else if (e.nbrUses(var.name())==0) {
118                    inactiveTerms.add(e);
119                } else {
120                    return new FNoExp();
121                }
122            }
123           
124            if (activeTerms.isEmpty()) {
125                return new FNoExp();
126            }
127               
128           
129            // Compute new expressions for the inactive and active terms
130            FExp inactiveTerm = computeInactiveSolution(inactiveTerms);
131            FExp activeTerm = computeActiveSolution(activeTerms, var,
132                    duringTearing && !myOptions().getBooleanOption("divide_by_vars_in_tearing") ||
133                    inactiveTerm.isLiteralZero(), duringTearing ? myOptions().getRealOption("tearing_division_tolerance") : 0);
134           
135            if (activeTerm.isNoExp()) {
136                return new FNoExp();
137            }
138            if (inactiveTerm.isLiteralZero() || activeTerm.isLiteralOne()) {
139                return inactiveTerm;
140            }
141            return new FDivExp(inactiveTerm, activeTerm);
142        } catch (ModelicaException e) {
143            throw e;
144        } catch (Exception e) {
145            throw new org.jmodelica.util.exceptions.InternalCompilerError("Exception caught while solving equation '" + prettyPrint("") + "' for variable '" + var.name() + "'", e);
146        }
147    }
148
149    private static FExp FEquation.computeInactiveSolution(Collection<FExp> exps) {
150        if (exps.size() == 0) {
151            return new FRealLitExp(0);
152        }
153        java.util.List<FExp> expsCopies = new ArrayList<FExp>();
154        // Compute new expressions for the inactive terms
155        for (FExp exp : exps) {
156            FExp expCopy = exp.copySymbolic();
157            if (!exp.isNegativeTerm()) {
158                expCopy = new FNegExp(expCopy);
159            }
160            expsCopies.add(expCopy);
161        }
162        return FExp.createBalancedBinaryTree(new FAddExp(), expsCopies);
163    }
164   
165    private static FExp FEquation.computeActiveSolution(Collection<FExp> exps, FVariable var,
166            boolean dontAllowDivisionByVariable, double tol) {
167        java.util.List<FExp> parts = new ArrayList<FExp>();
168        TypePrefixVariability mulCoeffVar = Variability.CONSTANT; // Keep track of variability
169        // Compute new expressions for the active terms
170        for (FExp e : exps) {
171            ArrayList<FExp> fac = e.factors();
172            // There is only one reference to the active variable in each
173            // term - this is checked above.           
174            ArrayList<FExp> activeFactors = new ArrayList<FExp>();
175            ArrayList<FExp> inactiveFactors = new ArrayList<FExp>();
176            // Find terms
177            boolean negatedFactor = false;         
178            for (FExp ee : fac) {
179                if (ee.nbrUses(var.name())==1 && !ee.isInvertedFactor() && 
180                    (ee.isIdentifier(var.name()) || // Identifier
181                    (ee instanceof FPreExp))) { // pre expression
182                    activeFactors.add(ee);
183                } else if (ee.nbrUses(var.name())==1 && !ee.isInvertedFactor() && 
184                    ((ee instanceof FNegExp) && ((FNegExp)ee).getFExp().isIdentifier(var.name()))) { 
185                    // TODO: remove this branch since it is not general enough
186                    activeFactors.add(((FNegExp)ee).getFExp());
187                    negatedFactor = true;
188                } else if (ee.nbrUses(var.name())==0) {
189                    if (!nominalAllowsDivision(ee, var, tol)) {
190                        return new FNoExp();
191                    }
192                    inactiveFactors.add(ee);
193                } else {
194                    // This equation cannot be solved
195                    return new FNoExp();
196                }
197            }
198           
199            FExp coeff = computeInactiveFactorSolution(inactiveFactors);
200            for (FExp eee : inactiveFactors) {
201                mulCoeffVar = mulCoeffVar.combine(eee.variability());
202            }
203            if (e.isNegativeTerm() || negatedFactor) {
204                coeff = new FNegExp(coeff);
205            }
206            parts.add(coeff);
207        }
208        if (dontAllowDivisionByVariable && !mulCoeffVar.lessOrEqual(Variability.CONSTANT)) {
209            return new FNoExp();
210        }
211        return FExp.createBalancedBinaryTree(new FAddExp(), parts);
212    }
213   
214    /**
215     * Check if nominal values allow the variable to be divided by the factor.
216     */
217    private static Boolean FEquation.nominalAllowsDivision(FExp factor, FVariable var, double tol) {
218        if (factor.variability().constantVariability()) {
219            FExp nominal = factor.dynamicFExp(var.nominal());
220            if (nominal.variability().constantVariability() &&
221                    Math.abs(factor.ceval().realValue()) *
222                    nominal.ceval().realValue() < tol) {
223                return false;
224            }
225        }
226        return true;
227    }
228   
229    private static FExp FEquation.computeInactiveFactorSolution(Collection<FExp> exps) {
230        if (exps.size() == 0) {
231            return new FRealLitExp(1);
232        }
233        java.util.List<FExp> expsCopies = new ArrayList<FExp>();
234        // Compute new expressions for the inactive factor terms
235        for (FExp exp : exps) {
236            FExp expCopy = exp.copySymbolic();
237            if (exp.isInvertedFactor()) {
238                expCopy = new FDivExp(new FRealLitExp(1), expCopy);
239            }
240            expsCopies.add(expCopy);
241        }
242        return FExp.createBalancedBinaryTree(new FMulExp(), expsCopies);
243    }
244
245    syn ArrayList<FExp> FAbstractEquation.terms() = new ArrayList<FExp>();
246
247    eq FEquation.terms() {
248        ArrayList<FExp> t = new ArrayList<FExp>();
249        t.addAll(getLeft().terms());
250        t.addAll(getRight().terms());
251        return t;
252    }
253
254    syn ArrayList<FExp> FExp.terms() {
255        ArrayList<FExp> t = new ArrayList<FExp>();
256        t.add(this);
257        return t;
258    }
259
260    eq FDotAddExp.terms() {
261        ArrayList<FExp> t = new ArrayList<FExp>();
262        t.addAll(getLeft().terms());
263        t.addAll(getRight().terms());
264        return t;
265    }
266
267    eq FStringAddExp.terms() {
268        ArrayList<FExp> t = new ArrayList<FExp>();
269        t.add(this);
270        return t;
271    }
272
273    eq FDotSubExp.terms() {
274        ArrayList<FExp> t = new ArrayList<FExp>();
275        t.addAll(getLeft().terms());
276        t.addAll(getRight().terms());
277        return t;
278    }
279
280        // Get factors 
281        syn ArrayList<FExp> FExp.factors() {
282                ArrayList<FExp> t = new ArrayList<FExp>();
283                t.add(this);
284                return t;
285        }
286       
287        eq FDotMulExp.factors() {
288                ArrayList<FExp> t = new ArrayList<FExp>();
289                t.addAll(getLeft().factors());
290                t.addAll(getRight().factors());
291                return t;
292        }
293
294        eq FDotDivExp.factors() {
295                ArrayList<FExp> t = new ArrayList<FExp>();
296                t.addAll(getLeft().factors());
297                t.addAll(getRight().factors());
298                return t;
299        }
300       
301        eq FNegExp.factors() {
302                ArrayList<FExp> t = new ArrayList<FExp>();
303                t.add(dynamicFExp(new FNegExp(new FRealLitExp(1.))));
304                t.addAll(getFExp().factors());
305                return t;
306        }
307
308        // Negated terms       
309        inh boolean FExp.isNegativeTerm();
310        eq FDotSubExp.getRight().isNegativeTerm() = !isNegativeTerm();
311        eq FNegExp.getFExp().isNegativeTerm() = !isNegativeTerm();
312        eq FEquation.getRight().isNegativeTerm() = true;
313        eq FEquation.getLeft().isNegativeTerm() = false;
314        eq Root.getChild().isNegativeTerm() = false;
315
316        // Inverted factors
317        inh boolean FExp.isInvertedFactor();
318        eq FDotDivExp.getRight().isInvertedFactor() = !isInvertedFactor();
319        eq FEquation.getChild().isInvertedFactor() = false;
320        eq Root.getChild().isInvertedFactor() = false;
321       
322        // Classification of terms
323        syn int FExp.isMulTerm(String name) = 0;
324        eq FDotMulExp.isMulTerm(String name) {
325                if (getLeft().isIdentifier(name)) {
326                        return 1;
327                } else if (getRight().isIdentifier(name)) {
328                        return 2;
329                } else {
330                        return 0;
331                }
332        }
333
334        syn boolean FExp.isDivTerm(String name) = false;
335        eq FDotDivExp.isDivTerm(String name) {
336                if (getLeft().isIdentifier(name)) {
337                        return true;
338                } else {
339                        return false;
340                }
341        }
342
343        syn boolean FExp.isNegTerm(String name) = false;
344        eq FNegExp.isNegTerm(String name) = getFExp().isIdentifier(name);
345       
346        syn int ASTNode.nbrUses(String name) {
347                int n = 0;
348                for (int i=0;i<getNumChild();i++) {
349                        n += getChild(i).nbrUses(name);
350                }
351                return n;
352        }
353
354    eq FEquation.nbrUses(String name) = getLeft().nbrUses(name) + getRight().nbrUses(name);
355    eq CommonAccessExp.nbrUses(String name) = name.equals(name())? 1: 0;
356
357    syn FExp FRelExp.solutionForTime() {
358        if (getLeft() instanceof FTimeExp && getRight().variability().discreteOrLess()) {
359            return getRight();
360        } else if (getRight() instanceof FTimeExp && getLeft().variability().discreteOrLess()) {
361            return getLeft();
362        }
363        return new FNoExp();
364    }
365   
366   
367    private FRelExp FRelExp.originalFRelExp                         = null;
368    private FSampleExp FSampleExp.originalFSampleExp                = null;
369    private FDelayExp FDelayExp.originalFDelayExp                   = null;
370    private FSpatialDistExp FSpatialDistExp.originalFSpatialDistExp = null;
371   
372    syn FExp FExp.originalFExp() { throw new UnsupportedOperationException(); }
373    syn FRelExp         FRelExp.originalFExp()         = (originalFRelExp == null)          ? this : originalFRelExp;
374    syn FSampleExp      FSampleExp.originalFExp()      = (originalFSampleExp == null)       ? this : originalFSampleExp;
375    syn FDelayExp       FDelayExp.originalFExp()       = (originalFDelayExp == null)        ? this : originalFDelayExp;
376    syn FSpatialDistExp FSpatialDistExp.originalFExp() = (originalFSpatialDistExp == null)  ? this : originalFSpatialDistExp;
377   
378    public void FRelExp.setOriginalFExp(FRelExp original)                 { this.originalFRelExp         = original; }
379    public void FSampleExp.setOriginalFExp(FSampleExp original)           { this.originalFSampleExp      = original; }
380    public void FDelayExp.setOriginalFExp(FDelayExp original)             { this.originalFDelayExp       = original; }
381    public void FSpatialDistExp.setOriginalFExp(FSpatialDistExp original) { this.originalFSpatialDistExp = original; }
382   
383    /**
384     * This method is used when the original node has been discarded and the
385     * orignal references in FRelExp and FSampleExp needs to be reset.
386     *
387     * It is also used when it is known that original expression doesn't
388     * produce events but the new does. For example the derivative of
389     * smooot(1, ...).
390     */
391    protected void ASTNode.resetOriginalReferences() {
392        for (int i = 0; i < getNumChildNoTransform(); i++) {
393            getChildNoTransform(i).resetOriginalReferences();
394        }
395    }
396   
397    @Override
398    protected void FRelExp.resetOriginalReferences() {
399        setOriginalFExp(this);
400        super.resetOriginalReferences();
401    }
402   
403    @Override
404    protected void FSampleExp.resetOriginalReferences() {
405        setOriginalFExp(this);
406        super.resetOriginalReferences();
407    }
408   
409    @Override
410    protected void FDelayExp.resetOriginalReferences() {
411        setOriginalFExp(this);
412        super.resetOriginalReferences();
413    }
414   
415    @Override
416    protected void FSpatialDistExp.resetOriginalReferences() {
417        setOriginalFExp(this);
418        super.resetOriginalReferences();
419    }
420
421        protected void ASTNode.traverseSymbolic(ASTNode e) {
422                for (int i = 0; i < getNumChildNoTransform(); i++) {
423                        getChildNoTransform(i).traverseSymbolic(e.getChildNoTransform(i));
424                }
425        }
426   
427    @Override
428        protected void FRelExp.traverseSymbolic(ASTNode e) {
429        setOriginalFExp(((FRelExp)e).originalFExp());
430                super.traverseSymbolic(e);
431        }
432   
433    @Override
434        protected void FSampleExp.traverseSymbolic(ASTNode e) {
435        setOriginalFExp(((FSampleExp)e).originalFExp());
436                super.traverseSymbolic(e);
437        }
438       
439    @Override
440    protected void FDelayExp.traverseSymbolic(ASTNode e) {
441        setOriginalFExp(((FDelayExp)e).originalFExp());
442        super.traverseSymbolic(e);
443    }
444   
445    @Override
446    protected void FSpatialDistExp.traverseSymbolic(ASTNode e) {
447        setOriginalFExp(((FSpatialDistExp)e).originalFExp());
448        super.traverseSymbolic(e);
449    }
450       
451    public ASTNode ASTNode.copySymbolic() {
452                ASTNode res = fullCopy();
453                res.traverseSymbolic(this);
454                return res;
455        }
456       
457    public FExp FExp.copySymbolic() { return (FExp) super.copySymbolic(); }
458    public FStatement FStatement.copySymbolic() { return (FStatement) super.copySymbolic(); }
459    public FAbstractEquation FAbstractEquation.copySymbolic() { return (FAbstractEquation) super.copySymbolic(); }
460}
461
462
463aspect ExpressionSimplification {
464    /**
465     * Is this expression a literal zero?
466     */
467    syn boolean FExp.isLiteralZero() = false;
468    eq FNumericLitExp.isLiteralZero()= type().isScalar() && ceval().realValue() == 0.0;
469    eq FDotAddExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero() && getRightNoTransform().isLiteralZero();
470    eq FDotSubExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero() && getRightNoTransform().isLiteralZero();
471    eq FDotMulExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero() || getRightNoTransform().isLiteralZero();
472    eq FDotDivExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero();
473    eq FNegExp.isLiteralZero()       = getFExpNoTransform().isLiteralZero();
474
475    /**
476     * Is this expression a literal one?
477     */
478    syn boolean FExp.isLiteralOne() = false;
479    eq FNumericLitExp.isLiteralOne()= type().isScalar() && ceval().realValue() == 1.0;
480    eq FDotMulExp.isLiteralOne()    = getLeftNoTransform().isLiteralOne() && getRightNoTransform().isLiteralOne();
481    eq FDotDivExp.isLiteralOne()    = getLeftNoTransform().isLiteralOne() && getRightNoTransform().isLiteralOne();
482    eq FNegExp.isLiteralOne()       = getFExpNoTransform().isLiteralMinusOne();
483
484    /**
485     * Is this expression a literal minus one?
486     */
487    syn boolean FExp.isLiteralMinusOne() = false;
488    eq FNumericLitExp.isLiteralMinusOne()= type().isScalar() && ceval().realValue() == -1.0;
489    eq FNegExp.isLiteralMinusOne()       = getFExpNoTransform().isLiteralOne();
490
491    /**
492     * Convert this subtraction to an addition.
493     *
494     * Does not copy any expressions.
495     */
496    syn FExp FDotSubExp.convertToAddition() = new FDotAddExp(getLeft(), getRight().makeNegated());
497    eq FSubExp.convertToAddition()          = new FAddExp(getLeft(), getRight().makeNegated());
498
499    rewrite FDotMulExp {
500        when (getLeftNoTransform().isLiteralZero())      to FExp getRight().size().createZeroFExp();
501        when (getLeftNoTransform().isLiteralOne())       to FExp getRight();
502        when (getLeftNoTransform().isLiteralMinusOne())  to FExp getRight().makeNegated();
503        when (getRightNoTransform().isLiteralZero())     to FExp getLeft().size().createZeroFExp();
504        when (getRightNoTransform().isLiteralOne())      to FExp getLeft();
505        when (getRightNoTransform().isLiteralMinusOne()) to FExp getLeft().makeNegated();
506    }
507
508    rewrite FDotDivExp {
509        when (getLeftNoTransform().isLiteralZero())       to FExp getRight().size().createZeroFExp();
510        when (getRightNoTransform().isLiteralOne())       to FExp getLeft();
511        when (getRightNoTransform().isLiteralMinusOne())  to FExp getLeft().makeNegated();
512    }
513
514    rewrite FDotSubExp {
515        when (getLeftNoTransform().isLiteralZero())     to FExp getRight().makeNegated();
516        when (getRightNoTransform() instanceof FNegExp) to FExp convertToAddition();
517        when (getRightNoTransform().isLiteralZero())    to FExp getLeft();
518    }
519
520    rewrite FDotAddExp {
521        when (getLeftNoTransform().isLiteralZero())  to FExp getRight();
522        when (getRightNoTransform().isLiteralZero()) to FExp getLeft();
523    }
524
525    rewrite FNegExp {
526        when (getFExpNoTransform().isLiteralZero())     to FExp getFExp();
527        when (getFExpNoTransform().hasSimpleNegation()) to FExp getFExp().makeNegated();
528    }
529
530    rewrite FIfExp {
531        when (getThenExpNoTransform().equalsSymbolic(getElseExpNoTransform())) to FExp getThenExp();
532    }
533   
534    syn boolean FExp.needsNoEvent() = true;
535    eq FVarRefExp.needsNoEvent()    = false;
536    eq FLitExp.needsNoEvent()       = false;
537
538    rewrite FNoEventExp {
539        when (!getFExpNoTransform().needsNoEvent()) to FExp getFExpNoTransform();
540    }
541   
542    syn boolean FExp.needsSmooth() = true;
543    eq FVarRefExp.needsSmooth()    = false;
544    eq FLitExp.needsSmooth()       = false;
545
546    rewrite FSmoothExp {
547        when (!getFExpNoTransform().needsSmooth()) to FExp getFExpNoTransform();
548    }
549
550
551    syn boolean FExp.equalsSymbolic(FExp other) {
552        return super.equals(other);
553    }
554
555    eq CommonAccessExp.equalsSymbolic(FExp other) {
556        if (!(other instanceof CommonAccessExp))
557            return false;
558        CommonAccessExp o = other.asCommonAccessExp();
559        return toString().equals(o.toString());
560    }
561
562    eq FNumericLitExp.equalsSymbolic(FExp other) {
563        if (!(other instanceof FNumericLitExp) || !type().isScalar() || !other.type().isScalar())
564            return false;
565        return ceval().equals(other.ceval());
566    }
567
568    eq FBooleanLitExpTrue.equalsSymbolic(FExp other)  = other instanceof FBooleanLitExpTrue;
569    eq FBooleanLitExpFalse.equalsSymbolic(FExp other) = other instanceof FBooleanLitExpFalse;
570
571    eq FStringLitExp.equalsSymbolic(FExp other) {
572        if (!(other instanceof FStringLitExp))
573            return false;
574        FStringLitExp o = (FStringLitExp) other;
575        return getString().equals(o.getString()); 
576    }
577
578    eq FEnumLitExp.equalsSymbolic(FExp other) {
579        if (!(other instanceof FEnumLitExp))
580            return false;
581        FEnumLitExp o = (FEnumLitExp) other;
582        return getValue() == o.getValue(); 
583    }
584   
585    eq FBinExp.equalsSymbolic(FExp other) {
586        if (getClass().equals(other.getClass())) {
587            FBinExp o = (FBinExp) other;
588            return getRightNoTransform().equalsSymbolic(o.getRightNoTransform())
589                    && getLeftNoTransform().equalsSymbolic(o.getLeftNoTransform());
590        }
591        return false;
592    }
593
594 }
Note: See TracBrowser for help on using the repository browser.