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

Last change on this file since 13720 was 13720, checked in by Jonathan Kämpe, 8 weeks ago

#5844 Merging refactoring to trunk

File size: 21.3 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 varNominalExp = factor.dynamicFExp(var.nominal());
220            if (varNominalExp.variability().constantVariability()) {
221                double factorNominal = Math.abs(factor.ceval().realValue());
222                double varNominal    = varNominalExp.ceval().realValue();
223                double totalNominal  = factorNominal * varNominal;
224                if (totalNominal < tol) {
225                    return false;
226                }
227            }
228        }
229        return true;
230    }
231   
232    private static FExp FEquation.computeInactiveFactorSolution(Collection<FExp> exps) {
233        if (exps.size() == 0) {
234            return new FRealLitExp(1);
235        }
236        java.util.List<FExp> expsCopies = new ArrayList<FExp>();
237        // Compute new expressions for the inactive factor terms
238        for (FExp exp : exps) {
239            FExp expCopy = exp.copySymbolic();
240            if (exp.isInvertedFactor()) {
241                expCopy = new FDivExp(new FRealLitExp(1), expCopy);
242            }
243            expsCopies.add(expCopy);
244        }
245        return FExp.createBalancedBinaryTree(new FMulExp(), expsCopies);
246    }
247
248    syn ArrayList<FExp> FAbstractEquation.terms() = new ArrayList<FExp>();
249
250    eq FEquation.terms() {
251        ArrayList<FExp> t = new ArrayList<FExp>();
252        t.addAll(getLeft().terms());
253        t.addAll(getRight().terms());
254        return t;
255    }
256
257    syn ArrayList<FExp> FExp.terms() {
258        ArrayList<FExp> t = new ArrayList<FExp>();
259        t.add(this);
260        return t;
261    }
262
263    eq FDotAddExp.terms() {
264        ArrayList<FExp> t = new ArrayList<FExp>();
265        t.addAll(getLeft().terms());
266        t.addAll(getRight().terms());
267        return t;
268    }
269
270    eq FStringAddExp.terms() {
271        ArrayList<FExp> t = new ArrayList<FExp>();
272        t.add(this);
273        return t;
274    }
275
276    eq FDotSubExp.terms() {
277        ArrayList<FExp> t = new ArrayList<FExp>();
278        t.addAll(getLeft().terms());
279        t.addAll(getRight().terms());
280        return t;
281    }
282
283        // Get factors 
284        syn ArrayList<FExp> FExp.factors() {
285                ArrayList<FExp> t = new ArrayList<FExp>();
286                t.add(this);
287                return t;
288        }
289       
290        eq FDotMulExp.factors() {
291                ArrayList<FExp> t = new ArrayList<FExp>();
292                t.addAll(getLeft().factors());
293                t.addAll(getRight().factors());
294                return t;
295        }
296
297        eq FDotDivExp.factors() {
298                ArrayList<FExp> t = new ArrayList<FExp>();
299                t.addAll(getLeft().factors());
300                t.addAll(getRight().factors());
301                return t;
302        }
303       
304        eq FNegExp.factors() {
305                ArrayList<FExp> t = new ArrayList<FExp>();
306                t.add(dynamicFExp(new FNegExp(new FRealLitExp(1.))));
307                t.addAll(getFExp().factors());
308                return t;
309        }
310
311        // Negated terms       
312        inh boolean FExp.isNegativeTerm();
313        eq FDotSubExp.getRight().isNegativeTerm() = !isNegativeTerm();
314        eq FNegExp.getFExp().isNegativeTerm() = !isNegativeTerm();
315        eq FEquation.getRight().isNegativeTerm() = true;
316        eq FEquation.getLeft().isNegativeTerm() = false;
317        eq Root.getChild().isNegativeTerm() = false;
318
319        // Inverted factors
320        inh boolean FExp.isInvertedFactor();
321        eq FDotDivExp.getRight().isInvertedFactor() = !isInvertedFactor();
322        eq FEquation.getChild().isInvertedFactor() = false;
323        eq Root.getChild().isInvertedFactor() = false;
324       
325        // Classification of terms
326        syn int FExp.isMulTerm(String name) = 0;
327        eq FDotMulExp.isMulTerm(String name) {
328                if (getLeft().isIdentifier(name)) {
329                        return 1;
330                } else if (getRight().isIdentifier(name)) {
331                        return 2;
332                } else {
333                        return 0;
334                }
335        }
336
337        syn boolean FExp.isDivTerm(String name) = false;
338        eq FDotDivExp.isDivTerm(String name) {
339                if (getLeft().isIdentifier(name)) {
340                        return true;
341                } else {
342                        return false;
343                }
344        }
345
346        syn boolean FExp.isNegTerm(String name) = false;
347        eq FNegExp.isNegTerm(String name) = getFExp().isIdentifier(name);
348       
349        syn int ASTNode.nbrUses(String name) {
350                int n = 0;
351                for (int i=0;i<getNumChild();i++) {
352                        n += getChild(i).nbrUses(name);
353                }
354                return n;
355        }
356
357    eq FEquation.nbrUses(String name) = getLeft().nbrUses(name) + getRight().nbrUses(name);
358    eq CommonAccessExp.nbrUses(String name) = name.equals(name())? 1: 0;
359
360    syn FExp FRelExp.solutionForTime() {
361        if (getLeft() instanceof FTimeExp && getRight().variability().discreteOrLess()) {
362            return getRight();
363        } else if (getRight() instanceof FTimeExp && getLeft().variability().discreteOrLess()) {
364            return getLeft();
365        }
366        return new FNoExp();
367    }
368   
369   
370    private FRelExp FRelExp.originalFRelExp                         = null;
371    private FSampleExp FSampleExp.originalFSampleExp                = null;
372    private FDelayExp FDelayExp.originalFDelayExp                   = null;
373    private FSpatialDistExp FSpatialDistExp.originalFSpatialDistExp = null;
374   
375    syn FExp FExp.originalFExp() { throw new UnsupportedOperationException(); }
376    syn FRelExp         FRelExp.originalFExp()         = (originalFRelExp == null)          ? this : originalFRelExp;
377    syn FSampleExp      FSampleExp.originalFExp()      = (originalFSampleExp == null)       ? this : originalFSampleExp;
378    syn FDelayExp       FDelayExp.originalFExp()       = (originalFDelayExp == null)        ? this : originalFDelayExp;
379    syn FSpatialDistExp FSpatialDistExp.originalFExp() = (originalFSpatialDistExp == null)  ? this : originalFSpatialDistExp;
380   
381    public void FRelExp.setOriginalFExp(FRelExp original)                 { this.originalFRelExp         = original; }
382    public void FSampleExp.setOriginalFExp(FSampleExp original)           { this.originalFSampleExp      = original; }
383    public void FDelayExp.setOriginalFExp(FDelayExp original)             { this.originalFDelayExp       = original; }
384    public void FSpatialDistExp.setOriginalFExp(FSpatialDistExp original) { this.originalFSpatialDistExp = original; }
385   
386    /**
387     * This method is used when the original node has been discarded and the
388     * orignal references in FRelExp and FSampleExp needs to be reset.
389     *
390     * It is also used when it is known that original expression doesn't
391     * produce events but the new does. For example the derivative of
392     * smooot(1, ...).
393     */
394    protected void ASTNode.resetOriginalReferences() {
395        for (int i = 0; i < getNumChildNoTransform(); i++) {
396            getChildNoTransform(i).resetOriginalReferences();
397        }
398    }
399   
400    @Override
401    protected void FRelExp.resetOriginalReferences() {
402        setOriginalFExp(this);
403        super.resetOriginalReferences();
404    }
405   
406    @Override
407    protected void FSampleExp.resetOriginalReferences() {
408        setOriginalFExp(this);
409        super.resetOriginalReferences();
410    }
411   
412    @Override
413    protected void FDelayExp.resetOriginalReferences() {
414        setOriginalFExp(this);
415        super.resetOriginalReferences();
416    }
417   
418    @Override
419    protected void FSpatialDistExp.resetOriginalReferences() {
420        setOriginalFExp(this);
421        super.resetOriginalReferences();
422    }
423
424        protected void ASTNode.traverseSymbolic(ASTNode e) {
425                for (int i = 0; i < getNumChildNoTransform(); i++) {
426                        getChildNoTransform(i).traverseSymbolic(e.getChildNoTransform(i));
427                }
428        }
429   
430    @Override
431        protected void FRelExp.traverseSymbolic(ASTNode e) {
432        setOriginalFExp(((FRelExp)e).originalFExp());
433                super.traverseSymbolic(e);
434        }
435   
436    @Override
437        protected void FSampleExp.traverseSymbolic(ASTNode e) {
438        setOriginalFExp(((FSampleExp)e).originalFExp());
439                super.traverseSymbolic(e);
440        }
441       
442    @Override
443    protected void FDelayExp.traverseSymbolic(ASTNode e) {
444        setOriginalFExp(((FDelayExp)e).originalFExp());
445        super.traverseSymbolic(e);
446    }
447   
448    @Override
449    protected void FSpatialDistExp.traverseSymbolic(ASTNode e) {
450        setOriginalFExp(((FSpatialDistExp)e).originalFExp());
451        super.traverseSymbolic(e);
452    }
453       
454    public ASTNode ASTNode.copySymbolic() {
455                ASTNode res = fullCopy();
456                res.traverseSymbolic(this);
457                return res;
458        }
459       
460    public FExp FExp.copySymbolic() { return (FExp) super.copySymbolic(); }
461    public FStatement FStatement.copySymbolic() { return (FStatement) super.copySymbolic(); }
462    public FAbstractEquation FAbstractEquation.copySymbolic() { return (FAbstractEquation) super.copySymbolic(); }
463}
464
465
466aspect ExpressionSimplification {
467    /**
468     * Is this expression a literal zero?
469     */
470    syn boolean FExp.isLiteralZero() = false;
471    eq FNumericLitExp.isLiteralZero()= type().isScalar() && ceval().realValue() == 0.0;
472    eq FDotAddExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero() && getRightNoTransform().isLiteralZero();
473    eq FDotSubExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero() && getRightNoTransform().isLiteralZero();
474    eq FDotMulExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero() || getRightNoTransform().isLiteralZero();
475    eq FDotDivExp.isLiteralZero()    = getLeftNoTransform().isLiteralZero();
476    eq FNegExp.isLiteralZero()       = getFExpNoTransform().isLiteralZero();
477
478    /**
479     * Is this expression a literal one?
480     */
481    syn boolean FExp.isLiteralOne() = false;
482    eq FNumericLitExp.isLiteralOne()= type().isScalar() && ceval().realValue() == 1.0;
483    eq FDotMulExp.isLiteralOne()    = getLeftNoTransform().isLiteralOne() && getRightNoTransform().isLiteralOne();
484    eq FDotDivExp.isLiteralOne()    = getLeftNoTransform().isLiteralOne() && getRightNoTransform().isLiteralOne();
485    eq FNegExp.isLiteralOne()       = getFExpNoTransform().isLiteralMinusOne();
486
487    /**
488     * Is this expression a literal minus one?
489     */
490    syn boolean FExp.isLiteralMinusOne() = false;
491    eq FNumericLitExp.isLiteralMinusOne()= type().isScalar() && ceval().realValue() == -1.0;
492    eq FNegExp.isLiteralMinusOne()       = getFExpNoTransform().isLiteralOne();
493
494    /**
495     * Convert this subtraction to an addition.
496     *
497     * Does not copy any expressions.
498     */
499    syn FExp FDotSubExp.convertToAddition() = new FDotAddExp(getLeft(), getRight().makeNegated());
500    eq FSubExp.convertToAddition()          = new FAddExp(getLeft(), getRight().makeNegated());
501
502    rewrite FDotMulExp {
503        when (getLeftNoTransform().isLiteralZero())      to FExp getRight().size().createZeroFExp();
504        when (getLeftNoTransform().isLiteralOne())       to FExp getRight();
505        when (getLeftNoTransform().isLiteralMinusOne())  to FExp getRight().makeNegated();
506        when (getRightNoTransform().isLiteralZero())     to FExp getLeft().size().createZeroFExp();
507        when (getRightNoTransform().isLiteralOne())      to FExp getLeft();
508        when (getRightNoTransform().isLiteralMinusOne()) to FExp getLeft().makeNegated();
509    }
510
511    rewrite FDotDivExp {
512        when (getLeftNoTransform().isLiteralZero())       to FExp getRight().size().createZeroFExp();
513        when (getRightNoTransform().isLiteralOne())       to FExp getLeft();
514        when (getRightNoTransform().isLiteralMinusOne())  to FExp getLeft().makeNegated();
515    }
516
517    rewrite FDotSubExp {
518        when (getLeftNoTransform().isLiteralZero())     to FExp getRight().makeNegated();
519        when (getRightNoTransform() instanceof FNegExp) to FExp convertToAddition();
520        when (getRightNoTransform().isLiteralZero())    to FExp getLeft();
521    }
522
523    rewrite FDotAddExp {
524        when (getLeftNoTransform().isLiteralZero())  to FExp getRight();
525        when (getRightNoTransform().isLiteralZero()) to FExp getLeft();
526    }
527
528    rewrite FNegExp {
529        when (getFExpNoTransform().isLiteralZero())     to FExp getFExp();
530        when (getFExpNoTransform().hasSimpleNegation()) to FExp getFExp().makeNegated();
531    }
532
533    rewrite FIfExp {
534        when (getThenExpNoTransform().equalsSymbolic(getElseExpNoTransform())) to FExp getThenExp();
535    }
536   
537    syn boolean FExp.needsNoEvent() = true;
538    eq FVarRefExp.needsNoEvent()    = false;
539    eq FLitExp.needsNoEvent()       = false;
540
541    rewrite FNoEventExp {
542        when (!getFExpNoTransform().needsNoEvent()) to FExp getFExpNoTransform();
543    }
544   
545    syn boolean FExp.needsSmooth() = true;
546    eq FVarRefExp.needsSmooth()    = false;
547    eq FLitExp.needsSmooth()       = false;
548
549    rewrite FSmoothExp {
550        when (!getFExpNoTransform().needsSmooth()) to FExp getFExpNoTransform();
551    }
552
553
554    syn boolean FExp.equalsSymbolic(FExp other) {
555        return super.equals(other);
556    }
557
558    eq CommonAccessExp.equalsSymbolic(FExp other) {
559        if (!(other instanceof CommonAccessExp))
560            return false;
561        CommonAccessExp o = other.asCommonAccessExp();
562        return toString().equals(o.toString());
563    }
564
565    eq FNumericLitExp.equalsSymbolic(FExp other) {
566        if (!(other instanceof FNumericLitExp) || !type().isScalar() || !other.type().isScalar())
567            return false;
568        return ceval().equals(other.ceval());
569    }
570
571    eq FBooleanLitExpTrue.equalsSymbolic(FExp other)  = other instanceof FBooleanLitExpTrue;
572    eq FBooleanLitExpFalse.equalsSymbolic(FExp other) = other instanceof FBooleanLitExpFalse;
573
574    eq FStringLitExp.equalsSymbolic(FExp other) {
575        if (!(other instanceof FStringLitExp))
576            return false;
577        FStringLitExp o = (FStringLitExp) other;
578        return getString().equals(o.getString()); 
579    }
580
581    eq FEnumLitExp.equalsSymbolic(FExp other) {
582        if (!(other instanceof FEnumLitExp))
583            return false;
584        FEnumLitExp o = (FEnumLitExp) other;
585        return getValue() == o.getValue(); 
586    }
587   
588    eq FBinExp.equalsSymbolic(FExp other) {
589        if (getClass().equals(other.getClass())) {
590            FBinExp o = (FBinExp) other;
591            return getRightNoTransform().equalsSymbolic(o.getRightNoTransform())
592                    && getLeftNoTransform().equalsSymbolic(o.getLeftNoTransform());
593        }
594        return false;
595    }
596
597 }
Note: See TracBrowser for help on using the repository browser.