OMCompiler/Compiler/NBackEnd/Modules/1_Main/NBResolveSingularities.mo
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * This file is part of OpenModelica. | ||
| 3 | * | ||
| 4 | * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), | ||
| 5 | * c/o Linköpings universitet, Department of Computer and Information Science, | ||
| 6 | * SE-58183 Linköping, Sweden. | ||
| 7 | * | ||
| 8 | * All rights reserved. | ||
| 9 | * | ||
| 10 | * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF AGPL VERSION 3 LICENSE OR | ||
| 11 | * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. | ||
| 12 | * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES | ||
| 13 | * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GNU AGPL | ||
| 14 | * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. | ||
| 15 | * | ||
| 16 | * The OpenModelica software and the OSMC (Open Source Modelica Consortium) | ||
| 17 | * Public License (OSMC-PL) are obtained from OSMC, either from the above | ||
| 18 | * address, from the URLs: | ||
| 19 | * http://www.openmodelica.org or | ||
| 20 | * https://github.com/OpenModelica/ or | ||
| 21 | * http://www.ida.liu.se/projects/OpenModelica, | ||
| 22 | * and in the OpenModelica distribution. | ||
| 23 | * | ||
| 24 | * GNU AGPL version 3 is obtained from: | ||
| 25 | * https://www.gnu.org/licenses/licenses.html#GPL | ||
| 26 | * | ||
| 27 | * This program is distributed WITHOUT ANY WARRANTY; without | ||
| 28 | * even the implied warranty of MERCHANTABILITY or FITNESS | ||
| 29 | * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH | ||
| 30 | * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL. | ||
| 31 | * | ||
| 32 | * See the full OSMC Public License conditions for more details. | ||
| 33 | * | ||
| 34 | */ | ||
| 35 | |||
| 36 | encapsulated package NBResolveSingularities | ||
| 37 | "file: NBResolveSingularities.mo | ||
| 38 | package: NBResolveSingularities | ||
| 39 | description: This file contains the functions to resolve structurally singular systems. | ||
| 40 | " | ||
| 41 | public | ||
| 42 | import Module = NBModule; | ||
| 43 | |||
| 44 | protected | ||
| 45 | // NF imports | ||
| 46 | import NFBackendExtension.{BackendInfo, VariableAttributes, StateSelect}; | ||
| 47 | import ComponentRef = NFComponentRef; | ||
| 48 | import Dimension = NFDimension; | ||
| 49 | import Expression = NFExpression; | ||
| 50 | import SimplifyExp = NFSimplifyExp; | ||
| 51 | import Call = NFCall; | ||
| 52 | import NFFunction.Function; | ||
| 53 | import Subscript = NFSubscript; | ||
| 54 | import Type = NFType; | ||
| 55 | |||
| 56 | // NB imports | ||
| 57 | import Adjacency = NBAdjacency; | ||
| 58 | import NBFunctionAlias.Call_Aux; | ||
| 59 | import Differentiate = NBDifferentiate; | ||
| 60 | import NBEquation.{Equation, EqData, EquationAttributes, EquationKind, EquationPointer, EquationPointers, SlicingStatus, Iterator}; | ||
| 61 | import Initialization = NBInitialization; | ||
| 62 | import Matching = NBMatching; | ||
| 63 | import Variable = NFVariable; | ||
| 64 | import BVariable = NBVariable; | ||
| 65 | import PointerWeak; | ||
| 66 | import NBVariable.{VarData, VariablePointer, VariablePointers}; | ||
| 67 | |||
| 68 | // util imports | ||
| 69 | import BackendUtil = NBBackendUtil; | ||
| 70 | import Slice = NBSlice; | ||
| 71 | import StringUtil; | ||
| 72 | import UnorderedSet; | ||
| 73 | |||
| 74 | public | ||
| 75 | function indexReduction | ||
| 76 | "algorithm | ||
| 77 | 1. IR | ||
| 78 | - get unkowns and eqs from markings and arrays | ||
| 79 | - collect state candidates from constraint eqs | ||
| 80 | - differentiate all eqs and collect new derivatives | ||
| 81 | |||
| 82 | 2. DUMMY DERIVATIVE | ||
| 83 | - sort vars with priority (StateSelect) | ||
| 84 | - (ToDo: remove always vars) | ||
| 85 | - create adjacency matrix from original vars/eqs | ||
| 86 | - match the system with inverse matching to respect ordering | ||
| 87 | - do not kick out never variables (provided by ordering) | ||
| 88 | - (ToDo: fail if a never variable could not be chosen) | ||
| 89 | - see if any equations are unmatched | ||
| 90 | - none unmatched -> static state selection | ||
| 91 | - any unmatched -> dynamic state selection with remaining eqs and vars | ||
| 92 | |||
| 93 | 3. STATIC AND DYNAMIC | ||
| 94 | - make all matched variables DUMMY_STATES and all corresponding derivatives DUMMY_DERIVATIVES | ||
| 95 | - move DUMMY_STATES to algebraic vars | ||
| 96 | |||
| 97 | 4. STATIC | ||
| 98 | - no additional tasks | ||
| 99 | |||
| 100 | (ToDo: 5. DYNAMIC) | ||
| 101 | - make ALL variables (besides StateSelect = always) DUMMY_STATES and corresponding derivatives DUMMY_DERIVATIVES | ||
| 102 | - create state set from remaining eqs and vars | ||
| 103 | - create a state and derivative variable for each remaining eq ($SET.x, $SET.dx) | ||
| 104 | - create state selection matrix $SET.A (parameter) | ||
| 105 | - create equations $SET.x[i] = sum($SET.A[i,j]*DUMMY_STATE[j] | forall j) | ||
| 106 | - create equations $SET.dx[i] = sum($SET.A[i,j]*DUMMY_DERIVATIVE[j] | forall j) | ||
| 107 | |||
| 108 | 6. AFTER IR | ||
| 109 | - add differentiated equations | ||
| 110 | - add adjacency matrix entries | ||
| 111 | - add new variables in correct arrays | ||
| 112 | " | ||
| 113 | extends Module.resolveSingularitiesInterface; | ||
| 114 | protected | ||
| 115 | Adjacency.Mapping mapping; | ||
| 116 | array<Boolean> excluded_eqns; | ||
| 117 | array<list<Integer>> msss; | ||
| 118 | list<Integer> marked_eqns; | ||
| 119 | Pointer<Equation> constraint, diffed_eqn; | ||
| 120 | list<Slice<VariablePointer>> states, dummy_states; | ||
| 121 | list<Pointer<Variable>> sliced_states, sliced_dummy_states, state_derivatives, dummy_derivatives = {}, dummy_slice_vars; | ||
| 122 | list<Pointer<Variable>> current_candidates, rest_candidates; | ||
| 123 | list<Slice<EquationPointer>> constraint_eqns, matched_eqns, unmatched_eqns; | ||
| 124 | list<Pointer<Equation>> new_eqns = {}, der_alias_eqns; | ||
| 125 | list<Pointer<Variable>> der_aliases; | ||
| 126 | Differentiate.DifferentiationArguments diffArguments; | ||
| 127 | Pointer<Differentiate.DifferentiationArguments> diffArguments_ptr; | ||
| 128 | VariablePointers candidate_ptrs; | ||
| 129 | EquationPointers constraint_ptrs; | ||
| 130 | Adjacency.Matrix set_adj, full_local; | ||
| 131 | Matching set_matching; | ||
| 132 | UnorderedMap<ComponentRef, Integer> vo, vn, eo, en; | ||
| 133 | list<tuple<String, BVariable.checkVar>> stages; | ||
| 134 | Option<list<Pointer<Variable>>> numeric_dummies; | ||
| 135 | UnorderedSet<ComponentRef> dummy_set; | ||
| 136 | BVariable.checkVar stageFunc; | ||
| 137 | String stageStr; | ||
| 138 | |||
| 139 | // slice handling | ||
| 140 | type SliceSet = UnorderedSet<Integer>; | ||
| 141 | UnorderedMap<ComponentRef, SliceSet> slice_map = UnorderedMap.new<SliceSet>(ComponentRef.hash, ComponentRef.isEqual); | ||
| 142 | UnorderedSet<ComponentRef> dummy_slice_set = UnorderedSet.new(ComponentRef.hash, ComponentRef.isEqual) "dummy variables to fill unslicable equations"; | ||
| 143 | |||
| 144 | // sliced state candidates get materialized as a whole alias variable + linking | ||
| 145 | // equation (see resolveSlicedCandidates); aux_index reuses VarData.getUniqueIndex | ||
| 146 | // (model-wide, already used for equation naming below) instead of a fresh counter, | ||
| 147 | // to avoid colliding with an alias from an earlier indexReduction call. | ||
| 148 | UnorderedMap<ComponentRef, Expression> alias_subst = UnorderedMap.new<Expression>(ComponentRef.hash, ComponentRef.isEqual); | ||
| 149 | list<Pointer<Equation>> alias_eqns; | ||
| 150 | |||
| 151 | Boolean debug = false; | ||
| 152 | algorithm | ||
| 153 | // get the mapping and fail if there is none | ||
| 154 | mapping := match mapping_opt | ||
| 155 | case SOME(mapping) then mapping; | ||
| 156 | else algorithm | ||
| 157 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because no mapping was provided."}); | |
| 158 | ✗ | then fail(); | |
| 159 | end match; | ||
| 160 | |||
| 161 | // mark the forbidden equations (discrete and already differentiated in previous index reduction steps) | ||
| 162 |
8/8✓ Branch 1 taken 5406 times.
✓ Branch 2 taken 371 times.
✓ Branch 3 taken 5406 times.
✓ Branch 4 taken 371 times.
✓ Branch 6 taken 5068 times.
✓ Branch 7 taken 338 times.
✓ Branch 9 taken 1371 times.
✓ Branch 10 taken 3697 times.
|
5777 | excluded_eqns := listArray(list(Equation.isDiscrete(eqn) or Equation.hasDerivative(eqn) for eqn in EquationPointers.toList(equations))); |
| 163 | |||
| 164 | // get the minimally structurally singular subset | ||
| 165 | msss := match adj | ||
| 166 | 371 | case Adjacency.FINAL() then getMSSS(adj.m, adj.mT, matching, excluded_eqns, mapping); | |
| 167 | else algorithm | ||
| 168 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " expected final matrix as adj input but got :\n" | |
| 169 | + Adjacency.Matrix.toString(adj)}); | ||
| 170 | ✗ | then fail(); | |
| 171 | end match; | ||
| 172 | |||
| 173 |
2/2✓ Branch 0 taken 43 times.
✓ Branch 1 taken 328 times.
|
371 | if not arrayLength(msss) == 0 then |
| 174 | changed := true; | ||
| 175 | // msss to flat unique list (via UnorderedSet) | ||
| 176 | 43 | marked_eqns := UnorderedSet.unique_list(List.flatten(arrayList(msss)), Util.id, intEq); | |
| 177 | // -------------------------------------------------------- | ||
| 178 | // 1. BASIC INDEX REDUCTION | ||
| 179 | // -------------------------------------------------------- | ||
| 180 | |||
| 181 | // get all unmatched eqns and state candidates | ||
| 182 | // state candidates and constraint equations are lists of slices | ||
| 183 | // adjacency matrix and arrays are only full based | ||
| 184 | // slice them before matching | ||
| 185 | 43 | (constraint_ptrs, candidate_ptrs, constraint_eqns) := getConstraintsAndCandidates(equations, marked_eqns, mapping); | |
| 186 | |||
| 187 |
2/2✓ Branch 0 taken 197 times.
✓ Branch 1 taken 43 times.
|
240 | for eq in constraint_eqns loop |
| 188 | 197 | UnorderedMap.add(Equation.getEqnName(Slice.getT(eq)), UnorderedSet.fromList(eq.indices, Util.id, intEq), slice_map); | |
| 189 | end for; | ||
| 190 | |||
| 191 | // a state derivative has to stay the derivative of its state, an alias takes its place as candidate | ||
| 192 | 43 | (candidate_ptrs, der_aliases, der_alias_eqns) := aliasStateDerivatives(candidate_ptrs, constraint_ptrs, VarData.getUniqueIndex(varData)); | |
| 193 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 42 times.
|
43 | if not listEmpty(der_aliases) then |
| 194 | 1 | varData := VarData.addTypedList(varData, der_aliases, NBVariable.VarData.VarType.ALGEBRAIC); | |
| 195 | 1 | variables := VariablePointers.addList(der_aliases, variables); | |
| 196 | 1 | new_eqns := listAppend(der_alias_eqns, new_eqns); | |
| 197 | end if; | ||
| 198 | |||
| 199 |
5/6✓ Branch 0 taken 197 times.
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 197 times.
✓ Branch 3 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 43 times.
|
240 | if VariablePointers.scalarSize(candidate_ptrs) < sum(Slice.size(eq, function Equation.size(resize = true)) for eq in constraint_eqns) then |
| 200 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because there was not enough state candidates to balance out the constraint equations.\n" | |
| 201 | + EquationPointers.toString(constraint_ptrs, "Constraint") + "\n" + VariablePointers.toString(candidate_ptrs, "State Candidate")}); | ||
| 202 | ✗ | fail(); | |
| 203 | end if; | ||
| 204 | |||
| 205 | // ToDo: differ between user dumping and developer dumping | ||
| 206 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
|
43 | if Flags.isSet(Flags.DUMMY_SELECT) then |
| 207 | ✗ | print(StringUtil.headline_1("Index Reduction") + "\n" | |
| 208 | + VariablePointers.toString(candidate_ptrs, "State Candidate") | ||
| 209 | + EquationPointers.toString(constraint_ptrs, "Constraint")); | ||
| 210 | end if; | ||
| 211 | |||
| 212 | // -------------------------------------------------------- | ||
| 213 | // 2. DUMMY DERIVATIVE | ||
| 214 | // -------------------------------------------------------- | ||
| 215 | // create full adjacency matrix and prepare data | ||
| 216 | 43 | full_local := Adjacency.Matrix.createFull(candidate_ptrs, constraint_ptrs, kind); | |
| 217 | set_adj := Adjacency.Matrix.EMPTY(NBAdjacency.MatrixStrictness.LINEAR); | ||
| 218 | 43 | rest_candidates := VariablePointers.toList(candidate_ptrs); | |
| 219 | 43 | eo := constraint_ptrs.map; | |
| 220 | 43 | en := UnorderedMap.new<Integer>(ComponentRef.hash, ComponentRef.isEqual); | |
| 221 | 43 | vo := UnorderedMap.new<Integer>(ComponentRef.hash, ComponentRef.isEqual); | |
| 222 | 43 | vn := UnorderedMap.new<Integer>(ComponentRef.hash, ComponentRef.isEqual); | |
| 223 | 43 | set_matching := NBMatching.EMPTY_MATCHING; | |
| 224 | |||
| 225 | // order of importance for variables to not be states: | ||
| 226 | 43 | stages := { | |
| 227 | ("1. StateSelect.NEVER", function BVariable.isStateSelect(stateSelect = StateSelect.NEVER)), | ||
| 228 | ("2. StateSelect.AVOID", function BVariable.isStateSelect(stateSelect = StateSelect.AVOID)), | ||
| 229 | ("3. Artificial Variables", BVariable.isArtificial), | ||
| 230 | ("4. StateSelect.DEFAULT without state order", function isDefaultWithoutStateOrder(state_order = VarData.getStateOrder(varData))), | ||
| 231 | ("5. StateSelect.DEFAULT", function BVariable.isStateSelect(stateSelect = StateSelect.DEFAULT)), | ||
| 232 | ("6. StateSelect.PREFER", function BVariable.isStateSelect(stateSelect = StateSelect.PREFER)) | ||
| 233 | }; | ||
| 234 | |||
| 235 | // linear constraints with constant coefficients: choose the dummy states numerically, | ||
| 236 | // a structural matching can choose dummy states with a singular Jacobian | ||
| 237 | 43 | numeric_dummies := numericDummySelection(constraint_ptrs, candidate_ptrs, orderCandidates(rest_candidates, stages)); | |
| 238 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 41 times.
|
43 | if isSome(numeric_dummies) then |
| 239 | 2 | SOME(current_candidates) := numeric_dummies; | |
| 240 |
4/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 2 times.
|
5 | dummy_set := UnorderedSet.fromList(list(BVariable.getVarName(var) for var in current_candidates), ComponentRef.hash, ComponentRef.isEqual); |
| 241 |
6/6✓ Branch 3 taken 1 time.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3 times.
✓ Branch 8 taken 2 times.
|
6 | dummy_states := list(Slice.SLICE(var, {}) for var guard(UnorderedSet.contains(BVariable.getVarName(var), dummy_set)) in VariablePointers.toList(candidate_ptrs)); |
| 242 |
6/6✓ Branch 3 taken 3 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 2 times.
|
6 | states := list(Slice.SLICE(var, {}) for var guard(not UnorderedSet.contains(BVariable.getVarName(var), dummy_set)) in VariablePointers.toList(candidate_ptrs)); |
| 243 | 2 | unmatched_eqns := {}; | |
| 244 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
|
2 | if Flags.isSet(Flags.DUMMY_SELECT) then |
| 245 | ✗ | print("[dummyselect] numeric selection of dummy states for linear constraints with constant coefficients\n"); | |
| 246 | end if; | ||
| 247 | else | ||
| 248 |
2/2✓ Branch 0 taken 246 times.
✓ Branch 1 taken 41 times.
|
287 | for stage in stages loop |
| 249 | 246 | (stageStr, stageFunc) := stage; | |
| 250 | // split the candidates to get all currently relevant ones | ||
| 251 | 246 | (current_candidates, rest_candidates) := List.splitOnTrue(rest_candidates, stageFunc); | |
| 252 | |||
| 253 |
2/2✓ Branch 0 taken 60 times.
✓ Branch 1 taken 186 times.
|
246 | if listEmpty(current_candidates) then |
| 254 | // nothing to do, no candidates for this stage or matching is already perfect | ||
| 255 | if debug then | ||
| 256 | print(StringUtil.headline_2("Nothing done for (" + stageStr + ") Index Reduction") + "\n"); | ||
| 257 | end if; | ||
| 258 | else | ||
| 259 | // prepare the current maps | ||
| 260 | 60 | vo := UnorderedMap.merge(vo, UnorderedMap.copy(vn), sourceInfo()); | |
| 261 |
4/4✓ Branch 0 taken 110 times.
✓ Branch 1 taken 60 times.
✓ Branch 2 taken 110 times.
✓ Branch 3 taken 60 times.
|
170 | vn := UnorderedMap.subMap(candidate_ptrs.map, list(BVariable.getVarName(var) for var in current_candidates)); |
| 262 | // expand the adjacency matrix | ||
| 263 | 60 | (set_adj, full_local) := Adjacency.Matrix.expand(set_adj, full_local, vo, vn, eo, en, candidate_ptrs, constraint_ptrs, kind); | |
| 264 | // continue matching | ||
| 265 | 60 | set_matching := Matching.regular(set_matching, set_adj, false, true, false); | |
| 266 | |||
| 267 | if debug then | ||
| 268 | print(Adjacency.Matrix.toString(set_adj, "(" + stageStr + ") Index Reduction")); | ||
| 269 | print(Matching.toString(set_matching, "(" + stageStr + ") Index Reduction")); | ||
| 270 | end if; | ||
| 271 | |||
| 272 |
1/4✗ Branch 1 not taken.
✓ Branch 2 taken 60 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
60 | if Matching.isEmpty(set_matching) and Matching.isPerfect(set_matching) then |
| 273 | if debug then | ||
| 274 | print(StringUtil.headline_2("Finished with perfect matching in stage " + stageStr + ".") + "\n"); | ||
| 275 | end if; | ||
| 276 | break; | ||
| 277 | end if; | ||
| 278 | end if; | ||
| 279 | end for; | ||
| 280 | |||
| 281 | // parse the result of the matching | ||
| 282 | 41 | (dummy_states, states, matched_eqns, unmatched_eqns) := Matching.getMatches(set_matching, Adjacency.Matrix.getMappingOpt(set_adj), candidate_ptrs, constraint_ptrs); | |
| 283 | end if; | ||
| 284 | 43 | unmatched_eqns := resolveSlicedUnmatched(unmatched_eqns, slice_map); | |
| 285 | |||
| 286 | // sliced state/dummy candidates have to be resolved before differentiation, so the | ||
| 287 | // constraint equations reference the whole alias rather than a slice of the | ||
| 288 | // original by the time they get differentiated below. Dummy and state sides are | ||
| 289 | // not symmetric: only a sliced state gets its own alias (resolveSlicedCandidates); | ||
| 290 | // a sliced dummy is upgraded to the whole variable instead once its sibling state | ||
| 291 | // slice(s) cover the rest (resolveSlicedDummyStates) -- see their docstrings. | ||
| 292 | 43 | dummy_states := resolveSlicedDummyStates(dummy_states, states); | |
| 293 | 43 | (states, alias_eqns) := resolveSlicedCandidates(states, alias_subst, VarData.getUniqueIndex(varData), VarData.getUniqueIndex(varData)); | |
| 294 |
2/2✓ Branch 1 taken 21 times.
✓ Branch 2 taken 22 times.
|
43 | if not UnorderedMap.isEmpty(alias_subst) then |
| 295 |
2/2✓ Branch 1 taken 104 times.
✓ Branch 2 taken 21 times.
|
125 | for constraint in EquationPointers.toList(constraint_ptrs) loop |
| 296 | 104 | substituteSlicedDummyEqn(constraint, alias_subst); | |
| 297 | end for; | ||
| 298 | end if; | ||
| 299 | 43 | new_eqns := listAppend(alias_eqns, new_eqns); | |
| 300 | // alias linking equations must be differentiated too (their derivative side needs | ||
| 301 | // linking as well, e.g. $DER.theta), so fold them into constraint_ptrs and let the | ||
| 302 | // loop below handle them; slice_map needs a matching (empty) entry so | ||
| 303 | // removeSlicedDerivatives treats them as unsliced. | ||
| 304 |
2/2✓ Branch 0 taken 21 times.
✓ Branch 1 taken 43 times.
|
64 | for eqn in alias_eqns loop |
| 305 | 21 | UnorderedMap.add(Equation.getEqnName(eqn), UnorderedSet.new(Util.id, intEq), slice_map); | |
| 306 | end for; | ||
| 307 | 43 | constraint_ptrs := EquationPointers.addList(alias_eqns, constraint_ptrs); | |
| 308 | |||
| 309 | // Build differentiation argument structure | ||
| 310 | 43 | diffArguments := Differentiate.DifferentiationArguments.default(NBDifferentiate.DifferentiationType.TIME, funcMap); | |
| 311 | 43 | diffArguments.diff_map := SOME(VarData.getStateOrder(varData)); | |
| 312 | 43 | diffArguments_ptr := Pointer.create(diffArguments); | |
| 313 | |||
| 314 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
|
43 | if Flags.isSet(Flags.DUMMY_SELECT) then |
| 315 | ✗ | print(StringUtil.headline_3("[dummyselect] 1. Differentiate the constraint equations")); | |
| 316 | end if; | ||
| 317 | |||
| 318 | // differentiate all eqns | ||
| 319 |
2/2✓ Branch 1 taken 218 times.
✓ Branch 2 taken 43 times.
|
261 | for constraint in EquationPointers.toList(constraint_ptrs) loop |
| 320 | 218 | diffed_eqn := Differentiate.differentiateEquationPointer(constraint, diffArguments_ptr); | |
| 321 | 218 | diffed_eqn := removeSlicedDerivatives(diffed_eqn, UnorderedMap.getSafe(Equation.getEqnName(constraint), slice_map, sourceInfo()), dummy_slice_set, VarData.getUniqueIndex(varData)); | |
| 322 | new_eqns := diffed_eqn :: new_eqns; | ||
| 323 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 218 times.
|
218 | if Flags.isSet(Flags.DUMMY_SELECT) then |
| 324 | ✗ | print("[dummyselect] constraint eqn:\t\t" + Equation.toString(Pointer.access(constraint)) + "\n"); | |
| 325 | ✗ | print("[dummyselect] differentiated eqn:\t" + Equation.toString(Pointer.access(diffed_eqn)) + "\n\n"); | |
| 326 | end if; | ||
| 327 | end for; | ||
| 328 | 43 | diffArguments := Pointer.access(diffArguments_ptr); | |
| 329 | |||
| 330 | // -------------------------------------------------------- | ||
| 331 | // 3. STATIC AND DYNAMIC STATE SELECTION | ||
| 332 | // -------------------------------------------------------- | ||
| 333 | // for both static and dynamic state selection all matched states are regarded dummys | ||
| 334 | // note: sliced candidates were already upgraded to whole above (resolveSlicedDummyStates); | ||
| 335 | // the else branch is a defensive fallback, not an expected path. | ||
| 336 |
2/2✓ Branch 0 taken 107 times.
✓ Branch 1 taken 43 times.
|
150 | for dummy in dummy_states loop |
| 337 |
1/2✓ Branch 0 taken 107 times.
✗ Branch 1 not taken.
|
107 | if listEmpty(dummy.indices) then |
| 338 | 107 | dummy_derivatives := BVariable.makeDummyState(Slice.getT(dummy)) :: dummy_derivatives; | |
| 339 | else | ||
| 340 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because slicing during index reduction is not yet supported.\n" | |
| 341 | + Slice.toString(dummy, BVariable.pointerToString, 10)}); | ||
| 342 | ✗ | fail(); | |
| 343 | end if; | ||
| 344 | end for; | ||
| 345 | |||
| 346 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
|
43 | if Flags.isSet(Flags.DUMMY_SELECT) then |
| 347 | ✗ | print(StringUtil.headline_4("[dummyselect] (" + intString(listLength(states)) + ") Selected States")); | |
| 348 | ✗ | print(Slice.lstToString(states, BVariable.pointerToString) + "\n\n"); | |
| 349 | end if; | ||
| 350 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 40 times.
|
43 | if Flags.isSet(Flags.DUMP_STATESELECTION_INFO) then |
| 351 | 3 | print(StringUtil.headline_4("[stateselection] (" + intString(listLength(diffArguments.new_vars)) + ") State Derivatives Created by Differentiation")); | |
| 352 | 3 | print(List.toString(diffArguments.new_vars, BVariable.pointerToString, List.Style.NEWLINE_TAB) + "\n\n"); | |
| 353 | 3 | print(StringUtil.headline_4("[stateselection] (" + intString(listLength(dummy_states)) + ") Selected Dummy States")); | |
| 354 | 3 | print(Slice.lstToString(dummy_states, BVariable.pointerToString) + "\n\n"); | |
| 355 | end if; | ||
| 356 | |||
| 357 |
1/2✓ Branch 0 taken 43 times.
✗ Branch 1 not taken.
|
43 | if listEmpty(unmatched_eqns) then |
| 358 | // -------------------------------------------------------- | ||
| 359 | // 4. STATIC STATE SELECTION | ||
| 360 | // -------------------------------------------------------- | ||
| 361 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
|
43 | if Flags.isSet(Flags.DUMMY_SELECT) then |
| 362 | ✗ | print(StringUtil.headline_2("\t STATIC STATE SELECTION\n\t(no unmatched equations)") + "\n"); | |
| 363 | end if; | ||
| 364 | else | ||
| 365 | // -------------------------------------------------------- | ||
| 366 | // 5. DYNAMIC STATE SELECTION | ||
| 367 | // -------------------------------------------------------- | ||
| 368 | ✗ | if Flags.isSet(Flags.DUMMY_SELECT) then | |
| 369 | ✗ | print(toStringDynamicSelect(dummy_states, unmatched_eqns)); | |
| 370 | end if; | ||
| 371 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because dynamic state selection is not yet supported."}); | |
| 372 | ✗ | fail(); | |
| 373 | end if; | ||
| 374 | |||
| 375 | // -------------------------------------------------------- | ||
| 376 | // 6. UPDATE VARIABLE AND EQUATION ARRAYS | ||
| 377 | // -------------------------------------------------------- | ||
| 378 | // filter all variables that were created during differentiation for state derivatives | ||
| 379 | // ToDo: these have to be slices as well! check if new created variables are whole dim of arrays | ||
| 380 | 43 | (state_derivatives, _) := List.extractOnTrue(diffArguments.new_vars, BVariable.isStateDerivative); | |
| 381 | |||
| 382 | // cleanup varData and expand eqData | ||
| 383 | // some algebraics -> states (to states) | ||
| 384 |
4/4✓ Branch 0 taken 29 times.
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 29 times.
✓ Branch 3 taken 43 times.
|
72 | sliced_states := list(Slice.getT(slice) for slice in states); |
| 385 | 43 | varData := VarData.addTypedList(varData, sliced_states, NBVariable.VarData.VarType.STATE); | |
| 386 | // new derivatives (to derivatives) | ||
| 387 | 43 | varData := VarData.addTypedList(varData, state_derivatives, NBVariable.VarData.VarType.STATE_DER); | |
| 388 | // some states -> dummy states (to algebraics) | ||
| 389 |
4/4✓ Branch 0 taken 107 times.
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 107 times.
✓ Branch 3 taken 43 times.
|
150 | sliced_dummy_states := list(Slice.getT(slice) for slice in dummy_states); |
| 390 | 43 | varData := VarData.addTypedList(varData, sliced_dummy_states, NBVariable.VarData.VarType.ALGEBRAIC); | |
| 391 | // some derivatives -> dummy derivatives (to algebraics) | ||
| 392 | 43 | varData := VarData.addTypedList(varData, dummy_derivatives, NBVariable.VarData.VarType.ALGEBRAIC); | |
| 393 | // new equations | ||
| 394 | 43 | eqData := EqData.addTypedList(eqData, new_eqns, EqData.EqType.CONTINUOUS); | |
| 395 | |||
| 396 | // add all new differentiated variables | ||
| 397 | 43 | variables := VariablePointers.addList(diffArguments.new_vars, variables); | |
| 398 | // add all dummy states | ||
| 399 | 43 | variables := VariablePointers.addList(sliced_dummy_states, variables); | |
| 400 | // remove all states | ||
| 401 | 43 | variables := VariablePointers.removeList(sliced_states, variables); | |
| 402 | // add new equations (after cleanup because equation names are added there) | ||
| 403 | 43 | equations := EquationPointers.addList(new_eqns, equations); | |
| 404 | |||
| 405 | // add all slice dummies that were added to fill equations which cannot be split | ||
| 406 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 43 times.
|
43 | dummy_slice_vars := list(BVariable.getVarPointer(cref, sourceInfo()) for cref in UnorderedSet.toList(dummy_slice_set)); |
| 407 | 43 | varData := VarData.addTypedList(varData, dummy_slice_vars, NBVariable.VarData.VarType.ALGEBRAIC); | |
| 408 | 43 | variables := VariablePointers.addList(dummy_slice_vars, variables); | |
| 409 | else | ||
| 410 | changed := false; | ||
| 411 | end if; | ||
| 412 | end indexReduction; | ||
| 413 | |||
| 414 | function balanceInitialization | ||
| 415 | extends Module.resolveSingularitiesInterface; | ||
| 416 | protected | ||
| 417 | list<Slice<VariablePointer>> unmatched_vars; | ||
| 418 | list<Slice<EquationPointer>> unmatched_eqns; | ||
| 419 | list<Pointer<Variable>> start_vars, failed_vars = {}; | ||
| 420 | list<Pointer<Equation>> sliced_eqns, start_eqns, kept_eqns; | ||
| 421 | list<Integer> remaining; | ||
| 422 | Pointer<Variable> var_ptr; | ||
| 423 | Pointer<list<Pointer<Variable>>> ptr_start_vars = Pointer.create({}); | ||
| 424 | Pointer<list<Pointer<Equation>>> ptr_start_eqns = Pointer.create({}); | ||
| 425 | Pointer<Integer> idx; | ||
| 426 | String error_msg; | ||
| 427 | UnorderedMap<ComponentRef, Integer> vo, vn, eo, en; | ||
| 428 | algorithm | ||
| 429 | 181 | (_, unmatched_vars, _, unmatched_eqns) := Matching.getMatches(matching, mapping_opt, variables, equations); | |
| 430 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 181 times.
|
181 | if Flags.isSet(Flags.INITIALIZATION) then |
| 431 | ✗ | print(toStringUnmatched(unmatched_vars, unmatched_eqns)); | |
| 432 | end if; | ||
| 433 |
3/4✓ Branch 0 taken 153 times.
✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 153 times.
|
181 | if not (listEmpty(unmatched_vars) and listEmpty(unmatched_eqns)) then |
| 434 | changed := true; | ||
| 435 | // -------------------------------------------------------- | ||
| 436 | // 1. Resolve Overdetermination | ||
| 437 | // -------------------------------------------------------- | ||
| 438 | // ToDo: unmatched eq -> dependencies -> matched eqns -> ... until no further dependencies | ||
| 439 | // dependency found twice on one branch --> loop --> fail | ||
| 440 | // recursively replace cref with solved equations | ||
| 441 | // simplify equation and check for 0 = 0 | ||
| 442 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
|
28 | if not listEmpty(unmatched_eqns) then |
| 443 | ✗ | Error.addMessage(Error.COMPILER_WARNING, {getInstanceName() | |
| 444 | + " reports an overdetermined initialization!\nChecking for consistency is not yet supported, following equations had to be removed:\n" | ||
| 445 | + Slice.lstToString(unmatched_eqns, function Equation.pointerToString(str = ""))}); | ||
| 446 | // copy old map to update adjacency matrix correctly | ||
| 447 | ✗ | eo := UnorderedMap.copy(equations.map); | |
| 448 | // get all unmatched equations and remove them from the system and overall equations | ||
| 449 | ✗ | sliced_eqns := list(Slice.getT(eqn) for eqn in unmatched_eqns); | |
| 450 | ✗ | equations := EquationPointers.removeList(sliced_eqns, equations); | |
| 451 | // only some indices of a for equation can be redundant, keep the other ones | ||
| 452 | kept_eqns := {}; | ||
| 453 | ✗ | for eqn_slice in unmatched_eqns loop | |
| 454 | ✗ | if not listEmpty(eqn_slice.indices) and Equation.isForEquation(Slice.getT(eqn_slice)) then | |
| 455 | ✗ | remaining := list(i for i guard(not List.contains(eqn_slice.indices, i, intEq)) in 0:(Equation.size(Slice.getT(eqn_slice)) - 1)); | |
| 456 | ✗ | (sliced_eqns, _) := Equation.slice(Slice.getT(eqn_slice), remaining); | |
| 457 | ✗ | kept_eqns := listAppend(sliced_eqns, kept_eqns); | |
| 458 | end if; | ||
| 459 | end for; | ||
| 460 | // also update adjacency matrices | ||
| 461 | ✗ | if listEmpty(kept_eqns) then | |
| 462 | ✗ | (adj, full) := Adjacency.Matrix.compress(adj, full, equations, variables, eo); | |
| 463 | else | ||
| 464 | ✗ | equations := EquationPointers.addList(kept_eqns, equations); | |
| 465 | ✗ | full := Adjacency.Matrix.createFull(variables, equations, kind); | |
| 466 | ✗ | adj := Adjacency.Matrix.fullToFinal(full, variables.map, equations.map, equations, NBAdjacency.MatrixStrictness.MATCHING); | |
| 467 | end if; | ||
| 468 | end if; | ||
| 469 | |||
| 470 | // -------------------------------------------------------- | ||
| 471 | // 2. Resolve Underdetermination | ||
| 472 | // -------------------------------------------------------- | ||
| 473 | 28 | idx := EqData.getUniqueIndex(eqData); | |
| 474 |
2/2✓ Branch 0 taken 54 times.
✓ Branch 1 taken 28 times.
|
82 | for var in unmatched_vars loop |
| 475 | 54 | var_ptr := Slice.getT(var); | |
| 476 |
1/2✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
|
54 | if BVariable.isFixable(var_ptr) then |
| 477 | // var = $START.var ($PRE.d = $START.d for previous vars) | ||
| 478 | // DO NOT SET VARIABLE TO FIXED! we might have to fix it again for Lambda=0 system | ||
| 479 | 54 | Initialization.createStartEquationSlice(var, ptr_start_vars, ptr_start_eqns, idx, true); | |
| 480 | else | ||
| 481 | failed_vars := var_ptr :: failed_vars; | ||
| 482 | end if; | ||
| 483 | end for; | ||
| 484 | |||
| 485 |
1/2✓ Branch 0 taken 28 times.
✗ Branch 1 not taken.
|
28 | if listEmpty(failed_vars) then |
| 486 | 28 | start_vars := Pointer.access(ptr_start_vars); | |
| 487 | 28 | start_eqns := Pointer.access(ptr_start_eqns); | |
| 488 | |||
| 489 | // copy old equation map to update adjacency matrices correctly | ||
| 490 | 28 | vo := variables.map; | |
| 491 | 28 | eo := UnorderedMap.copy(equations.map); | |
| 492 | |||
| 493 | // add new vars and equations to overall data | ||
| 494 | 28 | varData := VarData.addTypedList(varData, start_vars, VarData.VarType.START); | |
| 495 | 28 | eqData := EqData.addTypedList(eqData, start_eqns, EqData.EqType.INITIAL); | |
| 496 | |||
| 497 | // add new equations to system pointer arrays | ||
| 498 | 28 | equations := EquationPointers.addList(start_eqns, equations); | |
| 499 | |||
| 500 | // update adjacency matrices | ||
| 501 | 28 | vn := UnorderedMap.new<Integer>(ComponentRef.hash, ComponentRef.isEqual); | |
| 502 |
4/4✓ Branch 0 taken 54 times.
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 54 times.
✓ Branch 3 taken 28 times.
|
82 | en := UnorderedMap.subMap(equations.map, list(Equation.getEqnName(eqn) for eqn in start_eqns)); |
| 503 | 28 | (adj, full) := Adjacency.Matrix.expand(adj, full, vo, vn, eo, en, variables, equations, kind); | |
| 504 | |||
| 505 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 28 times.
|
28 | if Flags.isSet(Flags.INITIALIZATION) then |
| 506 | ✗ | print(List.toStringCustom(start_eqns, function Equation.pointerToString(str = ""), | |
| 507 | StringUtil.headline_4("Created Start Equations for balancing the Initialization (" + intString(listLength(start_eqns)) + "):"), "\t", "\n\t", "", false) + "\n\n"); | ||
| 508 | end if; | ||
| 509 | else | ||
| 510 | ✗ | error_msg := getInstanceName() | |
| 511 | + " failed because following non-fixable variables could not be solved:\n" | ||
| 512 | + List.toString(failed_vars, BVariable.pointerToString, List.Style.NEWLINE_TAB) + "\n"; | ||
| 513 | ✗ | if Flags.isSet(Flags.INITIALIZATION) then | |
| 514 | ✗ | error_msg := error_msg + "\nFollowing equations were created by fixing variables:\n" | |
| 515 | + List.toString(Pointer.access(ptr_start_eqns), function Equation.pointerToString(str = "\t"), List.Style.NEWLINE_TAB) + "\n"; | ||
| 516 | else | ||
| 517 | ✗ | error_msg := error_msg + "\nUse -d=initialization for more debug output."; | |
| 518 | end if; | ||
| 519 | ✗ | if Flags.isSet(Flags.BLT_DUMP) then | |
| 520 | ✗ | error_msg := error_msg + "\n" + VariablePointers.toString(variables, "All") + EquationPointers.toString(equations, "All") | |
| 521 | + Adjacency.Mapping.toString(Util.getOptionOrDefault(mapping_opt, Adjacency.Mapping.empty())) | ||
| 522 | + Adjacency.Matrix.toString(adj) + "\n" + Matching.toString(matching); | ||
| 523 | else | ||
| 524 | ✗ | error_msg := error_msg + "\nUse -d=bltdump for more verbose debug output."; | |
| 525 | end if; | ||
| 526 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{error_msg}); | |
| 527 | ✗ | fail(); | |
| 528 | end if; | ||
| 529 | else | ||
| 530 | changed := false; | ||
| 531 | end if; | ||
| 532 | end balanceInitialization; | ||
| 533 | |||
| 534 | protected | ||
| 535 | function getMSSS | ||
| 536 | "finds the minimal structurally singular subsets" | ||
| 537 | input Adjacency.IntMatrix m "eqn -> vars"; | ||
| 538 | input Adjacency.IntMatrix mT "var -> eqns"; | ||
| 539 | input Matching matching; | ||
| 540 | input array<Boolean> excluded_eqns; | ||
| 541 | input Adjacency.Mapping mapping; | ||
| 542 | output array<list<Integer>> msss; | ||
| 543 | protected | ||
| 544 | list<Integer> eqn_candidates = {}; | ||
| 545 | array<Integer> color_clustering; | ||
| 546 | array<Integer> eqn_coloring = arrayCreate(Adjacency.IntMatrix.rows(m), -1); | ||
| 547 | array<Integer> var_coloring = arrayCreate(Adjacency.IntMatrix.rows(mT), -1); | ||
| 548 | Integer color = 0; | ||
| 549 | algorithm | ||
| 550 | // find all unmatched equation indices | ||
| 551 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 371 times.
✓ Branch 2 taken 371 times.
✗ Branch 3 not taken.
|
10594 | for eqn in 1:arrayLength(matching.eqn_to_var) loop |
| 552 |
2/2✓ Branch 1 taken 505 times.
✓ Branch 2 taken 9347 times.
|
9852 | if matching.eqn_to_var[eqn] == -1 then |
| 553 | eqn_candidates := eqn :: eqn_candidates; | ||
| 554 | end if; | ||
| 555 | end for; | ||
| 556 |
2/2✓ Branch 1 taken 505 times.
✓ Branch 2 taken 371 times.
|
876 | color_clustering := listArray(list(i for i in 1:listLength(eqn_candidates))); |
| 557 | |||
| 558 | // use a new color for each uncolored equation | ||
| 559 |
2/2✓ Branch 0 taken 505 times.
✓ Branch 1 taken 371 times.
|
876 | for eqn in eqn_candidates loop |
| 560 |
1/2✓ Branch 1 taken 505 times.
✗ Branch 2 not taken.
|
505 | if eqn_coloring[eqn] == -1 then |
| 561 | 505 | color := color + 1; | |
| 562 | 505 | fillColorEqn(eqn, color, eqn_coloring, var_coloring, color_clustering, m, mT, matching, mapping); | |
| 563 | end if; | ||
| 564 | end for; | ||
| 565 | |||
| 566 | 371 | resolveClustering(color_clustering); | |
| 567 | |||
| 568 | // fill the msss array, sorting each equation to their respective color | ||
| 569 | 371 | msss := arrayCreate(color, {}); | |
| 570 |
1/2✓ Branch 0 taken 371 times.
✗ Branch 1 not taken.
|
10223 | for eqn in 1:arrayLength(eqn_coloring) loop |
| 571 |
4/4✓ Branch 1 taken 1888 times.
✓ Branch 2 taken 7964 times.
✓ Branch 5 taken 227 times.
✓ Branch 6 taken 1661 times.
|
9852 | if eqn_coloring[eqn] <> -1 and not excluded_eqns[mapping.eqn_StA[eqn]] then |
| 572 | 227 | color := color_clustering[eqn_coloring[eqn]]; | |
| 573 | 454 | msss[color] := eqn :: msss[color]; | |
| 574 | end if; | ||
| 575 | end for; | ||
| 576 | |||
| 577 | // remove all empty colors (purely discrete) | ||
| 578 |
6/6✓ Branch 1 taken 402 times.
✓ Branch 2 taken 103 times.
✓ Branch 3 taken 505 times.
✓ Branch 4 taken 371 times.
✓ Branch 5 taken 103 times.
✓ Branch 6 taken 371 times.
|
876 | msss := listArray(list(ms for ms guard(not listEmpty(ms)) in arrayList(msss))); |
| 579 | end getMSSS; | ||
| 580 | |||
| 581 | function fillColorEqn | ||
| 582 | "finds all connected equation nodes and colors them equally | ||
| 583 | starts at an equation" | ||
| 584 | input Integer eqn; | ||
| 585 | input Integer color; | ||
| 586 | input array<Integer> eqn_coloring; | ||
| 587 | input array<Integer> var_coloring; | ||
| 588 | input array<Integer> color_clustering; | ||
| 589 | input Adjacency.IntMatrix m "eqn -> vars"; | ||
| 590 | input Adjacency.IntMatrix mT "var -> eqns"; | ||
| 591 | input Matching matching; | ||
| 592 | input Adjacency.Mapping mapping; | ||
| 593 | protected | ||
| 594 | array<Integer> data = Adjacency.IntMatrix.entries(m); | ||
| 595 | Integer first = m.start[eqn]; | ||
| 596 | algorithm | ||
| 597 | 1888 | arrayUpdate(eqn_coloring, eqn, color); | |
| 598 |
2/2✓ Branch 1 taken 37 times.
✓ Branch 2 taken 1851 times.
|
4657 | for k in first:first + m.len[eqn] - 1 loop |
| 599 | 2769 | fillColorVar(data[k], color, eqn_coloring, var_coloring, color_clustering, m, mT, matching, mapping); | |
| 600 | end for; | ||
| 601 | end fillColorEqn; | ||
| 602 | |||
| 603 | function fillColorVar | ||
| 604 | "finds all connected equation nodes and colors them equally | ||
| 605 | starts at a variable" | ||
| 606 | input Integer var; | ||
| 607 | input Integer color; | ||
| 608 | input array<Integer> eqn_coloring; | ||
| 609 | input array<Integer> var_coloring; | ||
| 610 | input array<Integer> color_clustering; | ||
| 611 | input Adjacency.IntMatrix m "eqn -> vars"; | ||
| 612 | input Adjacency.IntMatrix mT "var -> eqns"; | ||
| 613 | input Matching matching; | ||
| 614 | input Adjacency.Mapping mapping; | ||
| 615 | protected | ||
| 616 | Integer eqn = matching.var_to_eqn[var]; | ||
| 617 | algorithm | ||
| 618 |
2/2✓ Branch 1 taken 1383 times.
✓ Branch 2 taken 1386 times.
|
2769 | if var_coloring[var] == -1 then |
| 619 | 1383 | arrayUpdate(var_coloring, var, color); | |
| 620 |
1/2✓ Branch 0 taken 1383 times.
✗ Branch 1 not taken.
|
1383 | if eqn <> -1 then |
| 621 |
1/2✓ Branch 1 taken 1383 times.
✗ Branch 2 not taken.
|
1383 | if eqn_coloring[eqn] == -1 then |
| 622 | 1383 | fillColorEqn(eqn, color, eqn_coloring, var_coloring, color_clustering, m, mT, matching, mapping); | |
| 623 | end if; | ||
| 624 | end if; | ||
| 625 | else | ||
| 626 | 1386 | colorClustering(var_coloring[var], color, color_clustering); | |
| 627 | end if; | ||
| 628 | end fillColorVar; | ||
| 629 | |||
| 630 | function colorClustering | ||
| 631 | input Integer old_color; | ||
| 632 | input Integer new_color; | ||
| 633 | input array<Integer> color_clustering; | ||
| 634 | algorithm | ||
| 635 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1386 times.
|
1388 | if color_clustering[old_color] <> old_color then |
| 636 | 2 | colorClustering(color_clustering[old_color], new_color, color_clustering); | |
| 637 | end if; | ||
| 638 | 1388 | arrayUpdate(color_clustering, old_color, new_color); | |
| 639 | end colorClustering; | ||
| 640 | |||
| 641 | function resolveClustering | ||
| 642 | input array<Integer> color_clustering; | ||
| 643 | protected | ||
| 644 | Integer color; | ||
| 645 | algorithm | ||
| 646 |
2/2✓ Branch 0 taken 328 times.
✓ Branch 1 taken 43 times.
|
876 | for i in 1:arrayLength(color_clustering) loop |
| 647 | color := i; | ||
| 648 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 505 times.
|
507 | while color_clustering[color] <> color loop |
| 649 | color := color_clustering[color]; | ||
| 650 | end while; | ||
| 651 | 505 | arrayUpdate(color_clustering, i, color); | |
| 652 | end for; | ||
| 653 | end resolveClustering; | ||
| 654 | |||
| 655 | function aliasStateDerivatives | ||
| 656 | "replaces state derivative candidates, e.g. $DER.x in v = $DER.x, by an alias a = $DER.x. | ||
| 657 | Otherwise the derivative would become a (dummy) state and stop being the derivative of x." | ||
| 658 | input output VariablePointers candidates; | ||
| 659 | input EquationPointers constraints; | ||
| 660 | input Pointer<Integer> uniqueIndex; | ||
| 661 | output list<Pointer<Variable>> aliases = {}; | ||
| 662 | output list<Pointer<Equation>> alias_eqns = {}; | ||
| 663 | protected | ||
| 664 | UnorderedMap<ComponentRef, ComponentRef> subst = UnorderedMap.new<ComponentRef>(ComponentRef.hash, ComponentRef.isEqual); | ||
| 665 | list<Pointer<Variable>> ders; | ||
| 666 | Pointer<Variable> alias_var; | ||
| 667 | ComponentRef der_cref, alias_cref; | ||
| 668 | algorithm | ||
| 669 |
6/6✓ Branch 2 taken 114 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 115 times.
✓ Branch 5 taken 43 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 43 times.
|
158 | ders := list(v for v guard(BVariable.isStateDerivative(v)) in VariablePointers.toList(candidates)); |
| 670 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 42 times.
|
43 | if listEmpty(ders) then |
| 671 | 42 | return; | |
| 672 | end if; | ||
| 673 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | for der_var in ders loop |
| 674 | 1 | der_cref := BVariable.getVarName(der_var); | |
| 675 | 1 | (alias_var, alias_cref) := BVariable.makeAuxVar(NBVariable.DUMMY_ALIAS_STR, Pointer.access(uniqueIndex), Variable.typeOf(Pointer.access(der_var)), false); | |
| 676 | 1 | Pointer.update(uniqueIndex, Pointer.access(uniqueIndex) + 1); | |
| 677 | 1 | alias_eqns := Equation.makeAssignment(Expression.fromCref(alias_cref), Expression.fromCref(der_cref), uniqueIndex, "DUM", Iterator.EMPTY(), EquationAttributes.default(EquationKind.CONTINUOUS, false)) :: alias_eqns; | |
| 678 | 1 | UnorderedMap.add(der_cref, alias_cref, subst); | |
| 679 | aliases := alias_var :: aliases; | ||
| 680 | end for; | ||
| 681 | 1 | candidates := VariablePointers.compress(VariablePointers.addList(aliases, VariablePointers.removeList(ders, candidates))); | |
| 682 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
|
3 | for constraint in EquationPointers.toList(constraints) loop |
| 683 | 2 | Pointer.update(constraint, Equation.map(Pointer.access(constraint), function substituteDerivativeAlias(subst = subst))); | |
| 684 | end for; | ||
| 685 | end aliasStateDerivatives; | ||
| 686 | |||
| 687 | function substituteDerivativeAlias | ||
| 688 | input output Expression exp; | ||
| 689 | input UnorderedMap<ComponentRef, ComponentRef> subst; | ||
| 690 | algorithm | ||
| 691 | exp := match exp | ||
| 692 | local | ||
| 693 | ComponentRef alias_cref; | ||
| 694 | case Expression.CREF() guard(UnorderedMap.contains(ComponentRef.stripSubscriptsAll(exp.cref), subst)) algorithm | ||
| 695 | 2 | alias_cref := UnorderedMap.getSafe(ComponentRef.stripSubscriptsAll(exp.cref), subst, sourceInfo()); | |
| 696 | 2 | then Expression.fromCref(ComponentRef.copySubscripts(exp.cref, alias_cref)); | |
| 697 | else exp; | ||
| 698 | end match; | ||
| 699 | end substituteDerivativeAlias; | ||
| 700 | |||
| 701 | function getConstraintsAndCandidates | ||
| 702 | input EquationPointers equations; | ||
| 703 | input list<Integer> marked_eqns; | ||
| 704 | input Adjacency.Mapping mapping; | ||
| 705 | output EquationPointers constr = EquationPointers.empty(); | ||
| 706 | output VariablePointers states = VariablePointers.empty(); | ||
| 707 | output list<Slice<EquationPointer>> sliced_constr = {}; | ||
| 708 | protected | ||
| 709 | UnorderedSet<Integer> eqn_indices = UnorderedSet.new(Util.id, intEq); | ||
| 710 | array<list<Integer>> eqn_slices = arrayCreate(EquationPointers.size(equations), {}); | ||
| 711 | UnorderedSet<ComponentRef> state_candidates = UnorderedSet.new(ComponentRef.hash, ComponentRef.isEqual); | ||
| 712 | Pointer<Equation> eqn_ptr; | ||
| 713 | Pointer<Variable> var_ptr; | ||
| 714 | algorithm | ||
| 715 | // collect all relevant constraint equations | ||
| 716 |
2/2✓ Branch 0 taken 227 times.
✓ Branch 1 taken 43 times.
|
270 | for eqn in marked_eqns loop |
| 717 | 227 | UnorderedSet.add(mapping.eqn_StA[eqn], eqn_indices); | |
| 718 | 454 | eqn_slices[mapping.eqn_StA[eqn]] := eqn :: eqn_slices[mapping.eqn_StA[eqn]]; | |
| 719 | end for; | ||
| 720 | |||
| 721 | // get the constraint equation, add it to the array and add all slices of it to the slice list | ||
| 722 | // furthermore, get all contained constrained equations | ||
| 723 |
2/2✓ Branch 1 taken 197 times.
✓ Branch 2 taken 43 times.
|
240 | for eqn in UnorderedSet.toList(eqn_indices) loop |
| 724 | 197 | eqn_ptr := EquationPointers.getEqnAt(equations, eqn); | |
| 725 | 197 | constr := EquationPointers.add(eqn_ptr, constr); | |
| 726 | 197 | sliced_constr := Slice.SLICE(eqn_ptr, eqn_slices[eqn]) :: sliced_constr; | |
| 727 |
2/2✓ Branch 2 taken 279 times.
✓ Branch 3 taken 197 times.
|
476 | for candidate in Equation.collectCrefs(Pointer.access(eqn_ptr), getStateCandidate) loop |
| 728 | 279 | UnorderedSet.add(candidate, state_candidates); | |
| 729 | end for; | ||
| 730 | end for; | ||
| 731 | |||
| 732 | // add all state candidates to the array | ||
| 733 |
2/2✓ Branch 1 taken 115 times.
✓ Branch 2 taken 43 times.
|
158 | for candidate in UnorderedSet.toList(state_candidates) loop |
| 734 | 115 | var_ptr := BVariable.getVarPointer(candidate, sourceInfo()); | |
| 735 | 115 | states := VariablePointers.add(var_ptr, states); | |
| 736 | end for; | ||
| 737 | end getConstraintsAndCandidates; | ||
| 738 | |||
| 739 | function getStateCandidate | ||
| 740 | input output ComponentRef cref "the cref to check"; | ||
| 741 | input UnorderedSet<ComponentRef> acc "accumulator for relevant crefs"; | ||
| 742 | protected | ||
| 743 | Pointer<Variable> var; | ||
| 744 | function getStateCandidateVar | ||
| 745 | input Pointer<Variable> var; | ||
| 746 | input UnorderedSet<ComponentRef> acc "accumulator for relevant crefs"; | ||
| 747 | algorithm | ||
| 748 |
10/12✓ Branch 1 taken 148 times.
✓ Branch 2 taken 401 times.
✓ Branch 4 taken 17 times.
✓ Branch 5 taken 384 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 384 times.
✓ Branch 10 taken 62 times.
✓ Branch 11 taken 322 times.
✓ Branch 13 taken 10 times.
✓ Branch 14 taken 312 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 10 times.
|
549 | if (BVariable.isContinuous(var, false) and not (BVariable.isTime(var) or BVariable.isDummyVariable(var) or BVariable.isDummyState(var) or (BVariable.isForcedState(var) and not BVariable.isStateSelect(var, StateSelect.PREFER)) )) then |
| 749 | 322 | UnorderedSet.add(BVariable.getVarName(var), acc); | |
| 750 | end if; | ||
| 751 | end getStateCandidateVar; | ||
| 752 | algorithm | ||
| 753 | 549 | var := BVariable.getVarPointer(cref, sourceInfo()); | |
| 754 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 549 times.
|
549 | if BVariable.isRecord(var) then |
| 755 | ✗ | for child in BVariable.getRecordChildrenCells(var) loop | |
| 756 | ✗ | getStateCandidateVar(PointerWeak.upgrade(child), acc); | |
| 757 | end for; | ||
| 758 | else | ||
| 759 | 549 | getStateCandidateVar(var, acc); | |
| 760 | end if; | ||
| 761 | end getStateCandidate; | ||
| 762 | |||
| 763 | function isDefaultWithoutStateOrder | ||
| 764 | "StateSelect.DEFAULT candidates whose derivative is not bound to a variable by an | ||
| 765 | equation der(x) = y. States with such an explicit derivative are kept as states if | ||
| 766 | possible: choosing other dummy states can make the constraint equations numerically | ||
| 767 | singular for them, e.g. for x1 = x2 + c*x3 and c*x3 = x2 - x4 the dummy states x2, x3 | ||
| 768 | leave the dependent states x1 = x4." | ||
| 769 | extends BVariable.checkVar; | ||
| 770 | input UnorderedMap<ComponentRef, ComponentRef> state_order; | ||
| 771 | algorithm | ||
| 772 |
4/4✓ Branch 1 taken 71 times.
✓ Branch 2 taken 16 times.
✓ Branch 5 taken 18 times.
✓ Branch 6 taken 53 times.
|
87 | b := BVariable.isStateSelect(var_ptr, StateSelect.DEFAULT) |
| 773 | and not UnorderedMap.contains(BVariable.getVarName(var_ptr), state_order); | ||
| 774 | end isDefaultWithoutStateOrder; | ||
| 775 | |||
| 776 | function orderCandidates | ||
| 777 | "orders the candidates by the stages, the preferred dummy states first" | ||
| 778 | input list<Pointer<Variable>> candidates; | ||
| 779 | input list<tuple<String, BVariable.checkVar>> stages; | ||
| 780 | output list<Pointer<Variable>> ordered; | ||
| 781 | protected | ||
| 782 | list<Pointer<Variable>> rest = candidates, current; | ||
| 783 | list<list<Pointer<Variable>>> parts = {}; | ||
| 784 | BVariable.checkVar stageFunc; | ||
| 785 | algorithm | ||
| 786 |
2/2✓ Branch 0 taken 258 times.
✓ Branch 1 taken 43 times.
|
301 | for stage in stages loop |
| 787 | 258 | (_, stageFunc) := stage; | |
| 788 | 258 | (current, rest) := List.splitOnTrue(rest, stageFunc); | |
| 789 | parts := current :: parts; | ||
| 790 | end for; | ||
| 791 | 43 | ordered := List.flatten(listReverse(rest :: parts)); | |
| 792 | end orderCandidates; | ||
| 793 | |||
| 794 | function numericDummySelection | ||
| 795 | "Selects the dummy states for constraint equations that are linear in the candidates | ||
| 796 | with constant coefficients: a candidate in the given order becomes a dummy state if it | ||
| 797 | increases the numerical rank of the coefficients of the dummy states. Returns NONE() | ||
| 798 | if the constraints are not of this kind or have no full rank. The structural matching | ||
| 799 | can choose dummy states with a singular Jacobian, e.g. for x1 = x2 + c*x3 and | ||
| 800 | c*x3 = x2 - x4 the dummy states x2, x3, since the coefficients cancel." | ||
| 801 | input EquationPointers constraints; | ||
| 802 | input VariablePointers candidates; | ||
| 803 | input list<Pointer<Variable>> ordered "preferred dummy states first"; | ||
| 804 | output Option<list<Pointer<Variable>>> dummies = NONE(); | ||
| 805 | protected | ||
| 806 | type SparseVector = list<tuple<Integer, Real>>; | ||
| 807 | UnorderedMap<ComponentRef, SparseVector> cols = UnorderedMap.new<SparseVector>(ComponentRef.hash, ComponentRef.isEqual); | ||
| 808 | Equation eqn; | ||
| 809 | Expression res, diff; | ||
| 810 | Differentiate.DifferentiationArguments args = Differentiate.DifferentiationArguments.default(NBDifferentiate.DifferentiationType.SIMPLE); | ||
| 811 | Integer m = EquationPointers.size(constraints), row = 0, rank = 0, pivot; | ||
| 812 | Real value, scale, pivot_val, a; | ||
| 813 | Boolean linear = true; | ||
| 814 | SparseVector col; | ||
| 815 | array<Integer> basis_piv; | ||
| 816 | array<Real> basis_val; | ||
| 817 | array<SparseVector> basis_vec; | ||
| 818 | list<Pointer<Variable>> selected = {}; | ||
| 819 | algorithm | ||
| 820 |
3/4✓ Branch 0 taken 43 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 36 times.
✓ Branch 4 taken 7 times.
|
43 | if m == 0 or List.any(ordered, BVariable.isArray) then |
| 821 | 36 | return; | |
| 822 | end if; | ||
| 823 | |||
| 824 | // coefficients of the candidates, column wise | ||
| 825 |
2/2✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
|
12 | for eqn_ptr in EquationPointers.toList(constraints) loop |
| 826 | 10 | row := row + 1; | |
| 827 | 10 | eqn := Pointer.access(eqn_ptr); | |
| 828 | linear := match eqn case Equation.SCALAR_EQUATION() then true; else false; end match; | ||
| 829 | if not linear then | ||
| 830 | ✗ | return; | |
| 831 | end if; | ||
| 832 | 10 | res := Equation.getResidualExp(eqn); | |
| 833 | // user functions can not be differentiated here and are hardly linear | ||
| 834 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
|
10 | if Expression.contains(res, isUserFunctionCall) then |
| 835 | ✗ | return; | |
| 836 | end if; | ||
| 837 |
2/2✓ Branch 4 taken 13 times.
✓ Branch 5 taken 5 times.
|
18 | for cref in Equation.collectCrefs(eqn, function Equation.collectFromMap(check_map = candidates.map)) loop |
| 838 | 13 | args.diffCref := cref; | |
| 839 | try | ||
| 840 | 13 | diff := SimplifyExp.simplify(Differentiate.differentiateExpression(res, args)); | |
| 841 | else | ||
| 842 | ✗ | return; | |
| 843 | end try; | ||
| 844 | (value, linear) := match diff | ||
| 845 | 8 | case Expression.REAL() then (diff.value, true); | |
| 846 | ✗ | case Expression.INTEGER() then (intReal(diff.value), true); | |
| 847 | else (0.0, false); | ||
| 848 | end match; | ||
| 849 | if not linear then | ||
| 850 | 5 | return; | |
| 851 | end if; | ||
| 852 |
1/2✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
|
8 | if value <> 0.0 then |
| 853 | 16 | UnorderedMap.add(cref, (row, value) :: UnorderedMap.getOrDefault(cref, cols, {}), cols); | |
| 854 | end if; | ||
| 855 | end for; | ||
| 856 | end for; | ||
| 857 | |||
| 858 | // greedy selection by incremental elimination | ||
| 859 | 2 | basis_piv := arrayCreate(m, 0); | |
| 860 | 2 | basis_val := arrayCreate(m, 0.0); | |
| 861 | 2 | basis_vec := arrayCreate(m, {}); | |
| 862 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | for var in ordered loop |
| 863 | 3 | col := listReverse(UnorderedMap.getOrDefault(BVariable.getVarName(var), cols, {})); | |
| 864 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 3 times.
|
7 | scale := List.fold(list(abs(Util.tuple22(e)) for e in col), realMax, 0.0); |
| 865 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
4 | for k in 1:rank loop |
| 866 | 1 | a := sparseGet(col, basis_piv[k]); | |
| 867 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if a <> 0.0 then |
| 868 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
|
1 | col := sparseAxpy(col, -a / basis_val[k], basis_vec[k], basis_piv[k]); |
| 869 | end if; | ||
| 870 | end for; | ||
| 871 | // the largest remaining entry is the pivot | ||
| 872 | pivot := 0; | ||
| 873 | pivot_val := 0.0; | ||
| 874 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
|
7 | for e in col loop |
| 875 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 time.
|
4 | if abs(Util.tuple22(e)) > abs(pivot_val) then |
| 876 | 3 | (pivot, pivot_val) := e; | |
| 877 | end if; | ||
| 878 | end for; | ||
| 879 |
2/4✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
|
3 | if pivot > 0 and abs(pivot_val) > 1e-10 * scale then |
| 880 | 3 | rank := rank + 1; | |
| 881 | 3 | basis_piv[rank] := pivot; | |
| 882 | 3 | basis_val[rank] := pivot_val; | |
| 883 | 3 | basis_vec[rank] := col; | |
| 884 | 3 | selected := var :: selected; | |
| 885 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | if rank == m then |
| 886 | 2 | dummies := SOME(listReverse(selected)); | |
| 887 | 2 | return; | |
| 888 | end if; | ||
| 889 | end if; | ||
| 890 | end for; | ||
| 891 | end numericDummySelection; | ||
| 892 | |||
| 893 | function isUserFunctionCall | ||
| 894 | input Expression exp; | ||
| 895 | output Boolean b; | ||
| 896 | algorithm | ||
| 897 | b := match exp | ||
| 898 | local | ||
| 899 | Function fn; | ||
| 900 | 3 | case Expression.CALL(call = Call.TYPED_CALL(fn = fn)) then not Function.isBuiltin(fn); | |
| 901 | else false; | ||
| 902 | end match; | ||
| 903 | end isUserFunctionCall; | ||
| 904 | |||
| 905 | function sparseGet | ||
| 906 | input list<tuple<Integer, Real>> v; | ||
| 907 | input Integer index; | ||
| 908 | output Real value = 0.0; | ||
| 909 | algorithm | ||
| 910 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | for e in v loop |
| 911 |
1/2✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
|
1 | if Util.tuple21(e) == index then |
| 912 | 1 | value := Util.tuple22(e); | |
| 913 | 1 | return; | |
| 914 | elseif Util.tuple21(e) > index then | ||
| 915 | ✗ | return; | |
| 916 | end if; | ||
| 917 | end for; | ||
| 918 | end sparseGet; | ||
| 919 | |||
| 920 | function sparseAxpy | ||
| 921 | "v + f*w for sparse vectors sorted by index, the entry at the eliminated index is removed" | ||
| 922 | input list<tuple<Integer, Real>> v; | ||
| 923 | input Real f; | ||
| 924 | input list<tuple<Integer, Real>> w; | ||
| 925 | input Integer eliminated; | ||
| 926 | output list<tuple<Integer, Real>> r = {}; | ||
| 927 | protected | ||
| 928 | list<tuple<Integer, Real>> v_rest = v, w_rest = w; | ||
| 929 | Integer i, j; | ||
| 930 | Real x, y; | ||
| 931 | algorithm | ||
| 932 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
|
3 | while not (listEmpty(v_rest) and listEmpty(w_rest)) loop |
| 933 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if listEmpty(w_rest) then |
| 934 | ✗ | r := listAppend(listReverse(v_rest), r); | |
| 935 | v_rest := {}; | ||
| 936 | elseif listEmpty(v_rest) then | ||
| 937 |
4/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
|
2 | r := listAppend(listReverse(list((Util.tuple21(e), f * Util.tuple22(e)) for e in w_rest)), r); |
| 938 | w_rest := {}; | ||
| 939 | else | ||
| 940 | 1 | (i, x) := listHead(v_rest); | |
| 941 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
|
1 | (j, y) := listHead(w_rest); |
| 942 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if i < j then |
| 943 | ✗ | r := (i, x) :: r; | |
| 944 | ✗ | v_rest := listRest(v_rest); | |
| 945 | elseif j < i then | ||
| 946 | ✗ | r := (j, f * y) :: r; | |
| 947 | ✗ | w_rest := listRest(w_rest); | |
| 948 | else | ||
| 949 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1 | if i <> eliminated and x + f * y <> 0.0 then |
| 950 | ✗ | r := (i, x + f * y) :: r; | |
| 951 | end if; | ||
| 952 | 1 | v_rest := listRest(v_rest); | |
| 953 | 1 | w_rest := listRest(w_rest); | |
| 954 | end if; | ||
| 955 | end if; | ||
| 956 | end while; | ||
| 957 | 1 | r := listReverse(r); | |
| 958 | end sparseAxpy; | ||
| 959 | |||
| 960 | function candidatePriority | ||
| 961 | "returns the priority of a variable for state selection. | ||
| 962 | higher priority -> better chance of getting picked as a state." | ||
| 963 | input Pointer<Variable> candidate; | ||
| 964 | output Integer prio; | ||
| 965 | algorithm | ||
| 966 | prio := match Pointer.access(candidate) | ||
| 967 | local | ||
| 968 | VariableAttributes attributes; | ||
| 969 | case Variable.VARIABLE(backendinfo = BackendInfo.BACKEND_INFO(attributes = attributes)) | ||
| 970 | then match VariableAttributes.getStateSelect(attributes) | ||
| 971 | case NFBackendExtension.StateSelect.NEVER then -200; | ||
| 972 | case NFBackendExtension.StateSelect.AVOID then -100; | ||
| 973 | case NFBackendExtension.StateSelect.DEFAULT then 0; | ||
| 974 | case NFBackendExtension.StateSelect.PREFER then 100; | ||
| 975 | case NFBackendExtension.StateSelect.ALWAYS then 200; | ||
| 976 | else 0; | ||
| 977 | end match; | ||
| 978 | else algorithm | ||
| 979 | then fail(); | ||
| 980 | end match; | ||
| 981 | end candidatePriority; | ||
| 982 | |||
| 983 | function sortCandidates | ||
| 984 | "sorts the state candidates" | ||
| 985 | input output list<Pointer<Variable>> candidates; | ||
| 986 | protected | ||
| 987 | list<tuple<Integer,Pointer<Variable>>> priorities = {}; | ||
| 988 | algorithm | ||
| 989 | ✗ | for candidate in candidates loop | |
| 990 | ✗ | priorities := (candidatePriority(candidate), candidate) :: priorities; | |
| 991 | end for; | ||
| 992 | ✗ | priorities := List.sort(priorities, BackendUtil.indexTplGt); | |
| 993 | ✗ | candidates := List.unzipSecond(priorities); | |
| 994 | end sortCandidates; | ||
| 995 | |||
| 996 | function resolveSlicedDummyStates | ||
| 997 | "unlike a sliced state candidate (resolveSlicedCandidates), a sliced dummy candidate | ||
| 998 | is not given its own alias -- a for-loop-indexed cref can't always be statically | ||
| 999 | resolved to one slice. Instead, once the combined state+dummy indices for a variable | ||
| 1000 | cover its full extent, the candidate is upgraded to the whole variable so the | ||
| 1001 | existing whole-variable BVariable.makeDummyState/isDummyState exclusion applies. | ||
| 1002 | Fails loudly if coverage is incomplete rather than wrongly excluding the rest of the | ||
| 1003 | variable from future candidacy." | ||
| 1004 | input output list<Slice<VariablePointer>> dummy_states; | ||
| 1005 | input list<Slice<VariablePointer>> states; | ||
| 1006 | protected | ||
| 1007 | type SliceSet = UnorderedSet<Integer>; | ||
| 1008 | UnorderedMap<ComponentRef, SliceSet> covered = UnorderedMap.new<SliceSet>(ComponentRef.hash, ComponentRef.isEqual); | ||
| 1009 | ComponentRef cref; | ||
| 1010 | SliceSet cover_set; | ||
| 1011 | Integer full_size; | ||
| 1012 | list<Slice<VariablePointer>> resolved = {}; | ||
| 1013 | algorithm | ||
| 1014 | // gather, per variable, every index matched as either state or dummy in this call | ||
| 1015 |
2/2✓ Branch 1 taken 136 times.
✓ Branch 2 taken 43 times.
|
179 | for cand in listAppend(states, dummy_states) loop |
| 1016 |
2/2✓ Branch 0 taken 42 times.
✓ Branch 1 taken 94 times.
|
136 | if not listEmpty(cand.indices) then |
| 1017 | 42 | cref := BVariable.getVarName(Slice.getT(cand)); | |
| 1018 |
2/2✓ Branch 1 taken 21 times.
✓ Branch 2 taken 21 times.
|
42 | if UnorderedMap.contains(cref, covered) then |
| 1019 | 21 | cover_set := UnorderedMap.getSafe(cref, covered, sourceInfo()); | |
| 1020 |
2/2✓ Branch 0 taken 63 times.
✓ Branch 1 taken 21 times.
|
84 | for idx in cand.indices loop |
| 1021 | 63 | UnorderedSet.add(idx, cover_set); | |
| 1022 | end for; | ||
| 1023 | else | ||
| 1024 | 21 | UnorderedMap.add(cref, UnorderedSet.fromList(cand.indices, Util.id, intEq), covered); | |
| 1025 | end if; | ||
| 1026 | end if; | ||
| 1027 | end for; | ||
| 1028 | |||
| 1029 |
2/2✓ Branch 0 taken 107 times.
✓ Branch 1 taken 43 times.
|
150 | for dummy in dummy_states loop |
| 1030 |
2/2✓ Branch 0 taken 86 times.
✓ Branch 1 taken 21 times.
|
107 | if listEmpty(dummy.indices) then |
| 1031 | resolved := dummy :: resolved; | ||
| 1032 | else | ||
| 1033 | 21 | cref := BVariable.getVarName(Slice.getT(dummy)); | |
| 1034 | 21 | cover_set := UnorderedMap.getSafe(cref, covered, sourceInfo()); | |
| 1035 | 21 | full_size := BVariable.size(Slice.getT(dummy)); | |
| 1036 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
21 | if UnorderedSet.size(cover_set) == full_size then |
| 1037 | 21 | resolved := Slice.SLICE(Slice.getT(dummy), {}) :: resolved; | |
| 1038 | else | ||
| 1039 | ✗ | Error.addMessage(Error.INTERNAL_ERROR,{getInstanceName() + " failed because the partially matched array variable " | |
| 1040 | + ComponentRef.toString(cref) + " could not be fully accounted for during index reduction (" | ||
| 1041 | + intString(UnorderedSet.size(cover_set)) + " of " + intString(full_size) | ||
| 1042 | + " elements matched as state or dummy state) -- the remainder belongs to a different, currently unresolved part of the system.\n" | ||
| 1043 | + Slice.toString(dummy, BVariable.pointerToString, 10)}); | ||
| 1044 | ✗ | fail(); | |
| 1045 | end if; | ||
| 1046 | end if; | ||
| 1047 | end for; | ||
| 1048 | 43 | dummy_states := listReverse(resolved); | |
| 1049 | end resolveSlicedDummyStates; | ||
| 1050 | |||
| 1051 | function resolveSlicedCandidates | ||
| 1052 | "resolves sliced entries of a *state* candidate list into a whole alias variable | ||
| 1053 | (sized to the selected indices) plus a linking equation, same technique as | ||
| 1054 | NBFunctionAlias.introduceSlicedStateAlias. Whole entries pass through unchanged. | ||
| 1055 | Must run before differentiation so the substitution rules in subst apply first. | ||
| 1056 | Not used for sliced dummy candidates -- see resolveSlicedDummyStates." | ||
| 1057 | input output list<Slice<VariablePointer>> candidates; | ||
| 1058 | input UnorderedMap<ComponentRef, Expression> subst; | ||
| 1059 | input Pointer<Integer> aux_index; | ||
| 1060 | input Pointer<Integer> eq_index; | ||
| 1061 | output list<Pointer<Equation>> alias_eqns = {}; | ||
| 1062 | protected | ||
| 1063 | list<Slice<VariablePointer>> resolved = {}; | ||
| 1064 | Pointer<Variable> alias_var; | ||
| 1065 | Pointer<Equation> alias_eqn; | ||
| 1066 | algorithm | ||
| 1067 |
2/2✓ Branch 0 taken 29 times.
✓ Branch 1 taken 43 times.
|
72 | for cand in candidates loop |
| 1068 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 21 times.
|
29 | if listEmpty(cand.indices) then |
| 1069 | resolved := cand :: resolved; | ||
| 1070 | else | ||
| 1071 | 21 | (alias_var, alias_eqn) := resolveSlicedDummy(cand, subst, aux_index, eq_index); | |
| 1072 | 21 | alias_eqns := alias_eqn :: alias_eqns; | |
| 1073 | 21 | resolved := Slice.SLICE(alias_var, {}) :: resolved; | |
| 1074 | end if; | ||
| 1075 | end for; | ||
| 1076 | 43 | candidates := listReverse(resolved); | |
| 1077 | end resolveSlicedCandidates; | ||
| 1078 | |||
| 1079 | function resolveSlicedDummy | ||
| 1080 | "materializes one sliced candidate as its own whole alias variable (see | ||
| 1081 | resolveSlicedCandidates): creates the alias, records cref substitution rules in | ||
| 1082 | subst for each selected index, and returns the linking equation." | ||
| 1083 | input Slice<VariablePointer> dummy; | ||
| 1084 | input UnorderedMap<ComponentRef, Expression> subst; | ||
| 1085 | input Pointer<Integer> aux_index; | ||
| 1086 | input Pointer<Integer> eq_index; | ||
| 1087 | output Pointer<Variable> alias_var; | ||
| 1088 | output Pointer<Equation> alias_eqn; | ||
| 1089 | protected | ||
| 1090 | Variable var = Pointer.access(Slice.getT(dummy)); | ||
| 1091 | ComponentRef orig_cref = BVariable.getVarName(Slice.getT(dummy)); | ||
| 1092 | Type elem_ty = Type.arrayElementType(Variable.typeOf(var)); | ||
| 1093 | list<Integer> sizes = list(Dimension.size(d) for d in Type.arrayDims(Variable.typeOf(var))); | ||
| 1094 | Integer n = listLength(dummy.indices); | ||
| 1095 | Type alias_ty; | ||
| 1096 | ComponentRef alias_cref, elem_cref, alias_elem_cref; | ||
| 1097 | list<Expression> elems = {}; | ||
| 1098 | list<Integer> loc; | ||
| 1099 | list<Subscript> subs; | ||
| 1100 | Integer i = 1; | ||
| 1101 | Expression rhs; | ||
| 1102 | algorithm | ||
| 1103 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
|
21 | alias_ty := if n == 1 then elem_ty else Type.ARRAY(elem_ty, {Dimension.fromInteger(n)}); |
| 1104 | 21 | (alias_var, alias_cref) := BVariable.makeAuxVar(NBVariable.DUMMY_ALIAS_STR, Pointer.access(aux_index), alias_ty, false); | |
| 1105 | 21 | Pointer.update(aux_index, Pointer.access(aux_index) + 1); | |
| 1106 | |||
| 1107 |
2/2✓ Branch 0 taken 21 times.
✓ Branch 1 taken 21 times.
|
42 | for idx in dummy.indices loop |
| 1108 | // idx is a zero-based flat index into the (possibly multi-dimensional) original | ||
| 1109 | // array; convert it back to a per-dimension, one-based subscript. | ||
| 1110 | 21 | loc := Slice.indexToLocation(idx, sizes); | |
| 1111 |
4/4✓ Branch 0 taken 42 times.
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 21 times.
|
63 | subs := list(Subscript.INDEX(Expression.INTEGER(l + 1)) for l in loc); |
| 1112 | // the subscripts belong to the array of records if the variable is a member of one | ||
| 1113 | 21 | elem_cref := ComponentRef.mergeSubscripts(subs, orig_cref, true, true, true); | |
| 1114 | 21 | elems := Expression.fromCref(elem_cref) :: elems; | |
| 1115 | |||
| 1116 |
1/2✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
|
21 | alias_elem_cref := if n == 1 then alias_cref else ComponentRef.setSubscripts({Subscript.INDEX(Expression.INTEGER(i))}, alias_cref); |
| 1117 | 21 | UnorderedMap.add(elem_cref, Expression.fromCref(alias_elem_cref), subst); | |
| 1118 | 21 | i := i + 1; | |
| 1119 | end for; | ||
| 1120 | 21 | elems := listReverse(elems); | |
| 1121 | |||
| 1122 |
1/2✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
|
21 | rhs := if n == 1 then listHead(elems) else Expression.makeArray(alias_ty, listArray(elems)); |
| 1123 | 21 | alias_eqn := Equation.makeAssignment(Expression.fromCref(alias_cref), rhs, eq_index, "DUM", Iterator.EMPTY(), EquationAttributes.default(EquationKind.CONTINUOUS, false)); | |
| 1124 | end resolveSlicedDummy; | ||
| 1125 | |||
| 1126 | function substituteSlicedDummyEqn | ||
| 1127 | "applies the sliced-dummy alias substitution rules (see resolveSlicedCandidates) to | ||
| 1128 | one constraint equation, in place, before it is differentiated." | ||
| 1129 | input Pointer<Equation> eqn_ptr; | ||
| 1130 | input UnorderedMap<ComponentRef, Expression> subst; | ||
| 1131 | protected | ||
| 1132 | Equation eqn = Pointer.access(eqn_ptr); | ||
| 1133 | algorithm | ||
| 1134 | 104 | eqn := Equation.map(eqn, function substituteSlicedDummyExp(subst = subst)); | |
| 1135 | 104 | Pointer.update(eqn_ptr, eqn); | |
| 1136 | end substituteSlicedDummyEqn; | ||
| 1137 | |||
| 1138 | function substituteSlicedDummyExp | ||
| 1139 | input output Expression exp; | ||
| 1140 | input UnorderedMap<ComponentRef, Expression> subst; | ||
| 1141 | algorithm | ||
| 1142 | exp := match exp | ||
| 1143 | case Expression.CREF() guard(UnorderedMap.contains(exp.cref, subst)) | ||
| 1144 | ✗ | then UnorderedMap.getSafe(exp.cref, subst, sourceInfo()); | |
| 1145 | else exp; | ||
| 1146 | end match; | ||
| 1147 | end substituteSlicedDummyExp; | ||
| 1148 | |||
| 1149 | function resolveSlicedUnmatched | ||
| 1150 | "removes all the unmatched slices that are irrelevant" | ||
| 1151 | input list<Slice<EquationPointer>> old_unmatched; | ||
| 1152 | output list<Slice<EquationPointer>> filtered_unmatched = {}; | ||
| 1153 | input UnorderedMap<ComponentRef, UnorderedSet<Integer>> slice_map; | ||
| 1154 | function resolveSlicedUnmatchedSingle | ||
| 1155 | input Slice<EquationPointer> eq; | ||
| 1156 | input output list<Slice<EquationPointer>> acc; | ||
| 1157 | input UnorderedMap<ComponentRef, UnorderedSet<Integer>> slice_map; | ||
| 1158 | protected | ||
| 1159 | UnorderedSet<Integer> relevant_indices; | ||
| 1160 | algorithm | ||
| 1161 | 41 | relevant_indices := UnorderedMap.getSafe(Equation.getEqnName(Slice.getT(eq)), slice_map, sourceInfo()); | |
| 1162 | // empty set indicates everything is relevant (just like slices without indices indicate everything) | ||
| 1163 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 41 times.
|
41 | if UnorderedSet.isEmpty(relevant_indices) then |
| 1164 | // case 1. | ||
| 1165 | acc := eq :: acc; | ||
| 1166 | else | ||
| 1167 | // case 2. | ||
| 1168 |
5/8✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 41 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 41 times.
✓ Branch 8 taken 41 times.
✗ Branch 9 not taken.
|
98 | eq.indices := list(ind for ind guard(UnorderedSet.contains(ind, relevant_indices)) in eq.indices); |
| 1169 | // if the list is empty, none are relevant. full relevance does not have to be considered here, would have been case 1. | ||
| 1170 |
1/2✓ Branch 0 taken 41 times.
✗ Branch 1 not taken.
|
41 | if not listEmpty(eq.indices) then acc := eq :: acc; end if; |
| 1171 | end if; | ||
| 1172 | end resolveSlicedUnmatchedSingle; | ||
| 1173 | algorithm | ||
| 1174 |
2/2✓ Branch 0 taken 41 times.
✓ Branch 1 taken 43 times.
|
84 | for eq in old_unmatched loop |
| 1175 | 41 | filtered_unmatched := resolveSlicedUnmatchedSingle(eq, filtered_unmatched, slice_map); | |
| 1176 | end for; | ||
| 1177 | end resolveSlicedUnmatched; | ||
| 1178 | |||
| 1179 | function removeSlicedDerivatives | ||
| 1180 | input output Pointer<Equation> derivative; | ||
| 1181 | input UnorderedSet<Integer> slice_set; | ||
| 1182 | input UnorderedSet<ComponentRef> dummy_slice_set; | ||
| 1183 | input Pointer<Integer> aux_index; | ||
| 1184 | protected | ||
| 1185 | Equation eqn; | ||
| 1186 | algorithm | ||
| 1187 | // only do something if the set is not empty implying full occurence | ||
| 1188 |
2/2✓ Branch 1 taken 21 times.
✓ Branch 2 taken 197 times.
|
218 | if not UnorderedSet.isEmpty(slice_set) then |
| 1189 | 197 | eqn := removeSlicedDerivateEqn(Pointer.access(derivative), Iterator.EMPTY(), dummy_slice_set, aux_index); | |
| 1190 | 197 | Pointer.update(derivative, eqn); | |
| 1191 | end if; | ||
| 1192 | end removeSlicedDerivatives; | ||
| 1193 | |||
| 1194 | function removeSlicedDerivateEqn | ||
| 1195 | input output Equation eqn; | ||
| 1196 | input Iterator iter; | ||
| 1197 | input UnorderedSet<ComponentRef> dummy_slice_set; | ||
| 1198 | input Pointer<Integer> aux_index; | ||
| 1199 | protected | ||
| 1200 | function replaceTupleLiterals | ||
| 1201 | input output Expression exp; | ||
| 1202 | input Iterator iter; | ||
| 1203 | input UnorderedSet<ComponentRef> dummy_slice_set; | ||
| 1204 | input Pointer<Integer> aux_index; | ||
| 1205 | protected | ||
| 1206 | ComponentRef aux; | ||
| 1207 | algorithm | ||
| 1208 | ✗ | if Expression.isLiteral(exp) then | |
| 1209 | ✗ | aux := Call_Aux.createName(Expression.typeOf(exp), iter, aux_index, NBVariable.DERIVATIVE_STR, false); | |
| 1210 | ✗ | exp := Expression.CREF(ComponentRef.getSubscriptedType(aux), aux); | |
| 1211 | ✗ | UnorderedSet.add(aux, dummy_slice_set); | |
| 1212 | end if; | ||
| 1213 | end replaceTupleLiterals; | ||
| 1214 | algorithm | ||
| 1215 | eqn := match eqn | ||
| 1216 | local | ||
| 1217 | Expression lhs; | ||
| 1218 | |||
| 1219 | // apply to each body equation | ||
| 1220 | case Equation.FOR_EQUATION() algorithm | ||
| 1221 |
4/4✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 42 times.
|
126 | eqn.body := list(removeSlicedDerivateEqn(b, eqn.iter, dummy_slice_set, aux_index) for b in eqn.body); |
| 1222 | then eqn; | ||
| 1223 | |||
| 1224 | // replace everything that evaluated to a literal expression with wildcard outputs | ||
| 1225 | case Equation.RECORD_EQUATION(lhs = lhs as Expression.TUPLE()) algorithm | ||
| 1226 | ✗ | lhs.elements := list(replaceTupleLiterals(e, iter, dummy_slice_set, aux_index) for e in lhs.elements); | |
| 1227 | ✗ | eqn.lhs := lhs; | |
| 1228 | then eqn; | ||
| 1229 | |||
| 1230 | // ToDo: more cases and fail if not doable | ||
| 1231 | |||
| 1232 | else eqn; | ||
| 1233 | end match; | ||
| 1234 | end removeSlicedDerivateEqn; | ||
| 1235 | |||
| 1236 | function toStringCandidatesConstraints | ||
| 1237 | input list<Slice<VariablePointer>> state_candidates; | ||
| 1238 | input list<Slice<EquationPointer>> constraint_eqns; | ||
| 1239 | output String str; | ||
| 1240 | algorithm | ||
| 1241 | ✗ | str := StringUtil.headline_1("Index Reduction") + "\n" | |
| 1242 | + StringUtil.headline_4("(" + intString(listLength(state_candidates)) + ") Sorted State Candidates") | ||
| 1243 | + Slice.lstToString(state_candidates, BVariable.pointerToString) + "\n" | ||
| 1244 | + StringUtil.headline_4("(" + intString(listLength(constraint_eqns)) + ") Constraint Equations") | ||
| 1245 | + Slice.lstToString(constraint_eqns, function Equation.pointerToString(str = "")) + "\n"; | ||
| 1246 | end toStringCandidatesConstraints; | ||
| 1247 | |||
| 1248 | function toStringDynamicSelect | ||
| 1249 | input list<Slice<VariablePointer>> dummy_states; | ||
| 1250 | input list<Slice<EquationPointer>> unmatched_eqns; | ||
| 1251 | output String str; | ||
| 1252 | algorithm | ||
| 1253 | ✗ | str := StringUtil.headline_2("\t DYNAMIC STATE SELECTION\n\t(some unmatched equations)") | |
| 1254 | + StringUtil.headline_4("(" + intString(listLength(dummy_states)) + ") Remaining State Candidates") | ||
| 1255 | + Slice.lstToString(dummy_states, BVariable.pointerToString) + "\n" | ||
| 1256 | + StringUtil.headline_4("(" + intString(listLength(unmatched_eqns)) + ") Remaining Equations") | ||
| 1257 | + Slice.lstToString(unmatched_eqns, function Equation.pointerToString(str = "")) + "\n"; | ||
| 1258 | end toStringDynamicSelect; | ||
| 1259 | |||
| 1260 | function toStringUnmatched | ||
| 1261 | input list<Slice<VariablePointer>> unmatched_vars; | ||
| 1262 | input list<Slice<EquationPointer>> unmatched_eqns; | ||
| 1263 | output String str; | ||
| 1264 | protected | ||
| 1265 | String s1, s2, s3, s4; | ||
| 1266 | algorithm | ||
| 1267 | ✗ | if listEmpty(unmatched_vars) then | |
| 1268 | ✗ | s1 := StringUtil.headline_4("Not underdetermined."); | |
| 1269 | s3 := ""; | ||
| 1270 | else | ||
| 1271 | ✗ | s1 := "Stage " + intString(listLength(unmatched_vars)) + " underdetermined.\n"; | |
| 1272 | ✗ | s3 := "\n" + StringUtil.headline_4("(" + intString(listLength(unmatched_vars)) + ") Unmatched variables:") | |
| 1273 | + Slice.lstToString(unmatched_vars, BVariable.pointerToString) + "\n"; | ||
| 1274 | end if; | ||
| 1275 | ✗ | if listEmpty(unmatched_eqns) then | |
| 1276 | ✗ | s2 := StringUtil.headline_4("Not overdetermined."); | |
| 1277 | s4 := ""; | ||
| 1278 | else | ||
| 1279 | ✗ | s2 := "Stage " + intString(listLength(unmatched_eqns)) + " overdetermined.\n"; | |
| 1280 | ✗ | s4 := "\n" + StringUtil.headline_4("(" + intString(listLength(unmatched_eqns)) + ") Unmatched equations:") | |
| 1281 | + Slice.lstToString(unmatched_eqns, function Equation.pointerToString(str = "")) + "\n"; | ||
| 1282 | end if; | ||
| 1283 | ✗ | str := s1 + s2 + s3 + s4 + "\n"; | |
| 1284 | end toStringUnmatched; | ||
| 1285 | |||
| 1286 | annotation(__OpenModelica_Interface="nbackend"); | ||
| 1287 | end NBResolveSingularities; | ||
| 1288 |