1(*
2 *  Modelicac
3 *
4 *  Copyright (C) 2005 - 2007 Imagine S.A.
5 *  For more information or commercial use please contact us at www.amesim.com
6 *
7 *  This program is free software; you can redistribute it and/or
8 *  modify it under the terms of the GNU General Public License
9 *  as published by the Free Software Foundation; either version 2
10 *  of the License, or (at your option) any later version.
11 *
12 *  This program is distributed in the hope that it will be useful,
13 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 *  GNU General Public License for more details.
16 *
17 *  You should have received a copy of the GNU General Public License
18 *  along with this program; if not, write to the Free Software
19 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
20 *
21 *)
22
23open SymbolicExpression
24
25
26exception Can't_perform_inversion
27
28let union xs ys =
29  List.fold_left (fun xs y -> if List.memq y xs then xs else y :: xs) xs ys
30
31module StringMap = Map.Make (struct type t = string let compare = compare end)
32
33let map_length map = StringMap.fold (fun _ _ acc -> acc + 1) map 0
34
35module IntegerElement =
36  struct
37    type t = Infinity | Int of int
38    let zero = Int 0
39    let infinity = Infinity
40    let equal = ( = )
41    let compare x y = match x, y with
42      | Int i, Int i' -> compare i i'
43      | Int _, Infinity -> -1
44      | Infinity, Int _ -> 1
45      | Infinity, Infinity -> 0
46    let add x y = match x, y with
47      | Int i, Int i' -> Int (i + i')
48      | _, Infinity | Infinity, _ -> Infinity
49    let sub x y = match x, y with
50      | Int i, Int i' -> Int (i - i')
51      | Infinity, Int _ -> Infinity
52      | _, Infinity -> assert false
53  end
54
55module IntegerHungarianMethod =
56  HungarianMethod.Make
57    (IntegerElement)
58    (SquareSparseMatrix.Make (IntegerElement))
59    (BipartiteGraph)
60
61type identifier_maps =
62  {
63    parameters_map: parameter StringMap.t Lazy.t;
64    inputs_map: input StringMap.t Lazy.t;
65    discrete_variables_map: variable StringMap.t Lazy.t;
66    variables_map: variable StringMap.t Lazy.t
67  }
68
69and parameter =
70  {
71    parameter_kind: Instantiation.parameter_kind;
72    parameter_type: parameter_type;
73    parameter_id: int;
74    parameter_comment: string;
75    parameter_start: t Lazy.t
76  }
77
78and parameter_type = IntegerType | StringType | RealType
79
80and input =
81  {
82    input_id: int;
83    input_name: string;
84    input_comment: string;
85  }
86
87and variable =
88  {
89    variable_id: int;
90    variable_comment: string;
91    variable_start: t Lazy.t option
92  }
93
94and model =
95  {
96    mutable inputs: string array;
97    mutable parameters: parameter_description array;
98    mutable discrete_variables: discrete_variable_description array;
99    mutable variables: variable_description array;
100    mutable equations: equation_description array;
101    mutable reinitializable_variables: t list;
102    mutable when_clauses: (t * when_expression list) list;
103    mutable io_dependency: bool;
104    mutable external_functions:
105      (string list *
106      Instantiation.expression_type list *
107      Instantiation.expression_type list) list;
108    trace: string option;
109    variables_infos: (string * Compilation.variable_infos) list
110  }
111
112and parameter_description =
113  {
114    mutable main: bool;
115    mutable p_type: parameter_type;
116    mutable p_name: string;
117    mutable p_comment: string;
118    mutable value: t
119  }
120
121and discrete_variable_description =
122  {
123    mutable d_output: int option;
124    mutable d_v_name: string;
125    mutable d_v_comment: string;
126    mutable d_start_value: t option
127  }
128
129and variable_description =
130  {
131    mutable output: int option;
132    mutable state: bool;
133    mutable v_name: string;
134    mutable depth_in_hierarchy: int;
135    mutable v_comment: string;
136    mutable start_value: t option
137  }
138
139and equation_description =
140  {
141    mutable solved: bool;
142    mutable inner_variables: t list;
143    mutable inner_derivatives: t list;
144    mutable assignable_variables: t list;
145    mutable expression: t
146  }
147
148and when_expression =
149  | Assign of t * t
150  | Reinit of t * t
151
152
153let scaling_factor = Num.power_num (Num.Int 10) (Num.Int 16)
154
155let exists_array p xs =
156  let l = Array.length xs in
157  let rec exists_array_from i =
158    i < l && (p xs.(i) || exists_array_from (i + 1)) in
159  exists_array_from 0
160
161let num_of_float f =
162  let num_of_positive_float f =
163    let m, e = frexp f in
164    let sm = string_of_float m in
165    let s = String.make 16 '0' in
166    String.blit sm 2 s 0 (String.length sm - 2);
167    let e' = Num.power_num (Num.Int 2) (Num.num_of_int e) in
168    Num.div_num (Num.mult_num (Num.num_of_string s) e') scaling_factor
169  in
170  if f = 0.0 then Num.Int 0
171  else if f < 0.0 then
172    let num = num_of_positive_float (abs_float f) in
173    Num.minus_num num
174  else num_of_positive_float f
175
176let string_of_reference = function
177  | [(name, [||])] -> name
178  | _ -> failwith "string_of_reference: bad reference"
179
180let is_main_parameter id parameters_map =
181  let parameter = StringMap.find id parameters_map in
182  match parameter.parameter_kind with
183    | Instantiation.Main -> true
184    | Instantiation.Sub -> false
185
186let get_parameter_start id parameters_map =
187  let parameter = StringMap.find id parameters_map in
188  Lazy.force parameter.parameter_start
189
190let get_start_value id variables_map =
191  let variable = StringMap.find id variables_map in
192  match variable.variable_start with
193    | None -> zero
194    | Some lexpr -> Lazy.force lexpr
195
196let rec symbolic_expression_of_expression inl_par maps iexpr =
197  let rec symbolic_argument_of_argument = function
198      | Instantiation.ScalarArgument iexpr ->
199          ScalarArgument (symbolic_expression_of_expression' iexpr)
200      | Instantiation.ArrayArgument (dims, iexprs) ->
201          ArrayArgument
202            (dims, Array.map symbolic_expression_of_expression' iexprs)
203  and symbolic_expression_of_expression' iexpr =
204    match iexpr.Instantiation.tex_expression with
205      | None -> assert false
206      | Some expr -> symbolic_expression_of_expression'' expr
207  and symbolic_expression_of_expression'' = function
208      | Instantiation.Abs iexpr ->
209          let expr = symbolic_expression_of_expression' iexpr in
210          symbolic_if
211            (symbolic_ge expr zero)
212            expr
213            (symbolic_minus expr)
214      | Instantiation.Acos iexpr ->
215          symbolic_acos (symbolic_expression_of_expression' iexpr)
216      | Instantiation.Acosh iexpr ->
217          symbolic_acosh (symbolic_expression_of_expression' iexpr)
218      | Instantiation.Addition (iexpr, iexpr') ->
219          symbolic_add
220            (symbolic_expression_of_expression' iexpr)
221            (symbolic_expression_of_expression' iexpr')
222      | Instantiation.And (iexpr, iexpr') ->
223          symbolic_and
224            (symbolic_expression_of_expression' iexpr)
225            (symbolic_expression_of_expression' iexpr')
226      | Instantiation.Asin iexpr ->
227          symbolic_asin (symbolic_expression_of_expression' iexpr)
228      | Instantiation.Asinh iexpr ->
229          symbolic_asinh (symbolic_expression_of_expression' iexpr)
230      | Instantiation.Atan iexpr ->
231          symbolic_atan (symbolic_expression_of_expression' iexpr)
232      | Instantiation.Atanh iexpr ->
233          symbolic_atanh (symbolic_expression_of_expression' iexpr)
234      | Instantiation.Boolean false -> false_value
235      | Instantiation.Boolean true -> true_value
236      | Instantiation.Cardinality _ ->
237          invalid_arg "symbolic_expression_of_expression'"
238      | Instantiation.CompoundElement _ ->
239          invalid_arg "symbolic_expression_of_expression'"
240      | Instantiation.Cos iexpr ->
241          symbolic_cos (symbolic_expression_of_expression' iexpr)
242      | Instantiation.Cosh iexpr ->
243          symbolic_cosh (symbolic_expression_of_expression' iexpr)
244      | Instantiation.Der iexpr ->
245          symbolic_derivative (symbolic_expression_of_expression' iexpr)
246      | Instantiation.Division (iexpr, iexpr') ->
247          symbolic_div
248            (symbolic_expression_of_expression' iexpr)
249            (symbolic_expression_of_expression' iexpr')
250      | Instantiation.Equals (iexpr, iexpr') ->
251          symbolic_eq
252            (symbolic_expression_of_expression' iexpr)
253            (symbolic_expression_of_expression' iexpr')
254      | Instantiation.Exp iexpr ->
255          symbolic_exp (symbolic_expression_of_expression' iexpr)
256      | Instantiation.ExternalFunctionCall (path, _, _, iargs) ->
257          symbolic_blackBox
258            (function_name_of path)
259            (List.map symbolic_argument_of_argument iargs)
260      | Instantiation.Floor iexpr ->
261          symbolic_floor (symbolic_expression_of_expression' iexpr)
262      | Instantiation.GreaterEqualThan (iexpr, iexpr') ->
263          symbolic_ge
264            (symbolic_expression_of_expression' iexpr)
265            (symbolic_expression_of_expression' iexpr')
266      | Instantiation.GreaterThan (iexpr, iexpr') ->
267          symbolic_gt
268            (symbolic_expression_of_expression' iexpr)
269            (symbolic_expression_of_expression' iexpr')
270      | Instantiation.If (iif_exprs, iexpr) ->
271          List.fold_right
272            (fun (iexpr, iexpr') expr ->
273              symbolic_if
274                (symbolic_expression_of_expression' iexpr)
275                (symbolic_expression_of_expression' iexpr')
276                expr)
277            iif_exprs
278            (symbolic_expression_of_expression' iexpr)
279      | Instantiation.Integer i -> create_integer i
280      | Instantiation.Log iexpr ->
281          symbolic_log (symbolic_expression_of_expression' iexpr)
282      | Instantiation.Max (iexpr, iexpr') ->
283          symbolic_max
284            (symbolic_expression_of_expression' iexpr)
285            (symbolic_expression_of_expression' iexpr')
286      | Instantiation.Min (iexpr, iexpr') ->
287          symbolic_min
288            (symbolic_expression_of_expression' iexpr)
289            (symbolic_expression_of_expression' iexpr')
290      | Instantiation.Minus iexpr ->
291          symbolic_minus (symbolic_expression_of_expression' iexpr)
292      | Instantiation.Mod (iexpr, iexpr') ->
293          let expr = symbolic_expression_of_expression' iexpr
294          and expr' = symbolic_expression_of_expression' iexpr'
295          in
296          symbolic_sub
297            expr
298            (symbolic_mult (symbolic_floor (symbolic_div expr expr')) expr')
299      | Instantiation.Multiplication (iexpr, iexpr') ->
300          symbolic_mult
301            (symbolic_expression_of_expression' iexpr)
302            (symbolic_expression_of_expression' iexpr')
303      | Instantiation.NoEvent iexpr ->
304          create_blackBox
305            "noEvent"
306            [ScalarArgument (symbolic_expression_of_expression' iexpr)]
307      | Instantiation.Not iexpr ->
308          symbolic_not (symbolic_expression_of_expression' iexpr)
309      | Instantiation.NotEquals (iexpr, iexpr') ->
310          symbolic_neq
311            (symbolic_expression_of_expression' iexpr)
312            (symbolic_expression_of_expression' iexpr')
313      | Instantiation.Or (iexpr, iexpr') ->
314          symbolic_or
315            (symbolic_expression_of_expression' iexpr)
316            (symbolic_expression_of_expression' iexpr')
317      | Instantiation.ParameterValue (_, iref) when inl_par ->
318          let id = string_of_reference iref in
319          if is_main_parameter id (Lazy.force maps.parameters_map) then
320            create_parameter
321              (StringMap.find id (Lazy.force maps.parameters_map)).parameter_id
322          else
323            get_parameter_start id (Lazy.force maps.parameters_map)
324      | Instantiation.ParameterValue (_, iref) ->
325          let id = string_of_reference iref in
326            create_parameter
327              (StringMap.find id (Lazy.force maps.parameters_map)).parameter_id
328      | Instantiation.Power (iexpr, iexpr') ->
329          symbolic_power
330            (symbolic_expression_of_expression' iexpr)
331            (symbolic_expression_of_expression' iexpr')
332      | Instantiation.Pre iexpr ->
333          symbolic_pre (symbolic_expression_of_expression' iexpr)
334      | Instantiation.Real f -> create_number (num_of_float f)
335      | Instantiation.Sin iexpr ->
336          symbolic_sin (symbolic_expression_of_expression' iexpr)
337      | Instantiation.Sinh iexpr ->
338          symbolic_sinh (symbolic_expression_of_expression' iexpr)
339      | Instantiation.Sqrt iexpr ->
340          symbolic_sqrt (symbolic_expression_of_expression' iexpr)
341      | Instantiation.String s -> create_string s
342      | Instantiation.Subtraction (iexpr, iexpr') ->
343          symbolic_sub
344            (symbolic_expression_of_expression' iexpr)
345            (symbolic_expression_of_expression' iexpr')
346      | Instantiation.Tan iexpr ->
347          symbolic_tan (symbolic_expression_of_expression' iexpr)
348      | Instantiation.Tanh iexpr ->
349          symbolic_tanh (symbolic_expression_of_expression' iexpr)
350      | Instantiation.Time -> time
351      | Instantiation.VariableStart (_, iref) ->
352          let id = string_of_reference iref in
353          begin try get_start_value id (Lazy.force maps.variables_map) with
354            | Not_found ->
355                get_start_value id (Lazy.force maps.discrete_variables_map)
356          end
357      | Instantiation.VariableValue (_, iref) -> create_variable_value iref
358      | Instantiation.Vector iexprs ->
359          invalid_arg "symbolic_expression_of_expression'"
360  and create_variable_value iref =
361    let id = string_of_reference iref in
362    try
363      create_variable
364        (StringMap.find id (Lazy.force maps.variables_map)).variable_id
365    with
366      | _ ->
367          try
368            let dvar =
369              StringMap.find id (Lazy.force maps.discrete_variables_map)
370            in create_discrete_variable dvar.variable_id
371          with
372            | _ ->
373                try
374                  let inp = StringMap.find id (Lazy.force maps.inputs_map) in
375                  create_discrete_variable (-1 - (inp.input_id))
376                  (* Use of strictly negative values to be able to distinguish
377                     inputs from other discrete variables. *)
378                with
379                  | _ -> assert false
380  and function_name_of = function
381    | [] -> assert false
382    | [s] -> s
383    | s :: ss -> function_name_of ss
384  in symbolic_expression_of_expression' iexpr
385
386let collect_external_function_names iequs =
387  let rec add_if_not_in (name, in_types, out_types) = function
388    | [] -> [name, in_types, out_types]
389    | ((name', _, _) :: _) as name_and_types when name = name' ->
390        name_and_types
391    | name_and_type :: name_and_types ->
392        name_and_type ::
393        add_if_not_in (name, in_types, out_types) name_and_types in
394  let rec collect_in_equations funcalls = function
395    | [] -> funcalls
396    | Instantiation.Equation (iexpr, iexpr') :: iequs' ->
397        let funcalls = collect_in_expressions funcalls iexpr in
398        let funcalls = collect_in_expressions funcalls iexpr' in
399        collect_in_equations funcalls iequs'
400    | Instantiation.ConditionalEquation (iif_equs, iequs) :: iequs' ->
401        let funcalls = collect_in_if_clauses funcalls iif_equs in
402        let funcalls = collect_in_equations funcalls iequs in
403        collect_in_equations funcalls iequs'
404    | Instantiation.When iwhen_clauses :: iequs' ->
405        let funcalls = collect_in_when_clauses funcalls iwhen_clauses in
406        collect_in_equations funcalls iequs'
407    | Instantiation.FlowConnection _ :: _-> assert false
408  and collect_in_if_clauses funcalls = function
409    | [] -> funcalls
410    | (iexpr, iequs) :: iif_equs ->
411        let funcalls = collect_in_expressions funcalls iexpr in
412        let funcalls = collect_in_equations funcalls iequs in
413        collect_in_if_clauses funcalls iif_equs
414  and collect_in_when_clauses funcalls = function
415    | [] -> funcalls
416    | (iexpr, iwhen_equs) :: iwhen_clauses ->
417        let funcalls = collect_in_expressions funcalls iexpr in
418        let funcalls =
419          List.fold_left
420            (fun
421              funcalls
422              (Instantiation.Reinit (iexpr, iexpr') |
423              Instantiation.Assign (iexpr, iexpr')) ->
424              let funcalls =
425                collect_in_expressions funcalls iexpr
426              in collect_in_expressions funcalls iexpr')
427            funcalls
428            iwhen_equs
429        in collect_in_when_clauses funcalls iwhen_clauses
430  and collect_in_expressions funcalls iexpr =
431    match iexpr.Instantiation.tex_expression with
432      | None -> funcalls
433      | Some expr -> collect_in_expressions' funcalls expr
434  and collect_in_expressions' funcalls = function
435      | Instantiation.Addition (iexpr, iexpr') |
436        Instantiation.And (iexpr, iexpr') |
437        Instantiation.Division (iexpr, iexpr') |
438        Instantiation.Equals (iexpr, iexpr') |
439        Instantiation.GreaterEqualThan (iexpr, iexpr') |
440        Instantiation.GreaterThan (iexpr, iexpr') |
441        Instantiation.Max (iexpr, iexpr') | Instantiation.Min (iexpr, iexpr') |
442        Instantiation.Mod (iexpr, iexpr') |
443        Instantiation.Multiplication (iexpr, iexpr') |
444        Instantiation.NotEquals (iexpr, iexpr') |
445        Instantiation.Or (iexpr, iexpr') | Instantiation.Power (iexpr, iexpr') |
446        Instantiation.Subtraction (iexpr, iexpr') ->
447          let funcalls = collect_in_expressions funcalls iexpr in
448          collect_in_expressions funcalls iexpr'
449      | Instantiation.ExternalFunctionCall (name, in_types, out_types, iargs) ->
450          let funcalls = add_if_not_in (name, in_types, out_types) funcalls in
451          List.fold_left collect_in_arguments funcalls iargs
452      | Instantiation.If (iif_exprs, iexpr) ->
453          let funcalls =
454            List.fold_left
455              (fun funcalls (iexpr, iexpr') ->
456                let funcalls = collect_in_expressions funcalls iexpr in
457                collect_in_expressions funcalls iexpr')
458              funcalls
459              iif_exprs
460          in collect_in_expressions funcalls iexpr
461      | Instantiation.Minus iexpr | Instantiation.NoEvent iexpr |
462        Instantiation.Not iexpr | Instantiation.Abs iexpr |
463        Instantiation.Acos iexpr | Instantiation.Acosh iexpr |
464        Instantiation.Cos iexpr | Instantiation.Cosh iexpr |
465        Instantiation.Exp iexpr |
466        Instantiation.Floor iexpr | Instantiation.Log iexpr |
467        Instantiation.Asin iexpr | Instantiation.Asinh iexpr |
468        Instantiation.Sin iexpr | Instantiation.Sinh iexpr |
469        Instantiation.Sqrt iexpr |
470        Instantiation.Atan iexpr | Instantiation.Atanh iexpr |
471        Instantiation.Tan iexpr | Instantiation.Tanh iexpr |
472        Instantiation.Pre iexpr -> collect_in_expressions funcalls iexpr
473      | Instantiation.ParameterValue _ | Instantiation.Real _ |
474        Instantiation.String _ | Instantiation.Time |
475        Instantiation.VariableStart _ | Instantiation.VariableValue _ |
476        Instantiation.Boolean _ | Instantiation.Der _ | Instantiation.Integer _ ->
477          funcalls
478      | Instantiation.Cardinality _ | Instantiation.CompoundElement _ |
479        Instantiation.Vector _ -> assert false
480  and collect_in_arguments funcalls = function
481    | Instantiation.ScalarArgument iexpr ->
482        collect_in_expressions funcalls iexpr
483    | Instantiation.ArrayArgument (_, iexprs) ->
484        Array.fold_left collect_in_expressions funcalls iexprs
485  in collect_in_equations [] iequs
486
487let separate_parameters_from_variables icpnts =
488  let is_parameter = function
489    | _, Instantiation.InstantiatedParameter _ -> true
490    | _, Instantiation.InstantiatedVariable _ -> false
491  in List.partition is_parameter icpnts
492
493let separate_inputs_from_others icpnts =
494  let is_input = function
495    | Instantiation.InstantiatedIntegerVariable
496        (_, Compilation.Input, _, _) |
497      Instantiation.InstantiatedStringVariable
498        (_, Compilation.Input, _, _) |
499      Instantiation.InstantiatedDiscreteVariable
500        (_, Compilation.Input, _, _) |
501      Instantiation.InstantiatedRealVariable
502        (_, Compilation.Input, _, _, _) ->
503        true
504    | Instantiation.InstantiatedIntegerVariable _ |
505      Instantiation.InstantiatedStringVariable _ |
506      Instantiation.InstantiatedDiscreteVariable _ |
507      Instantiation.InstantiatedRealVariable _ |
508      Instantiation.InstantiatedCompoundVariable _ -> false
509  in
510  let filter_variable = function
511    | _, Instantiation.InstantiatedParameter _ -> false
512    | _, Instantiation.InstantiatedVariable ivar -> is_input ivar
513  in List.partition filter_variable icpnts
514
515let separate_discrete_variables_from_others icpnts =
516  let is_discrete = function
517    | Instantiation.InstantiatedDiscreteVariable _ -> true
518    | Instantiation.InstantiatedIntegerVariable _ |
519      Instantiation.InstantiatedStringVariable _ |
520      Instantiation.InstantiatedRealVariable _ |
521      Instantiation.InstantiatedCompoundVariable _ -> false
522  in
523  let filter_variable = function
524    | _, Instantiation.InstantiatedParameter _ -> false
525    | _, Instantiation.InstantiatedVariable ivar -> is_discrete ivar
526  in List.partition filter_variable icpnts
527
528let separate_outputs_from_others icpnts =
529  let is_output = function
530    | Instantiation.InstantiatedIntegerVariable
531        (_, Compilation.Output, _, _) |
532      Instantiation.InstantiatedStringVariable
533        (_, Compilation.Output, _, _) |
534      Instantiation.InstantiatedDiscreteVariable
535        (_, Compilation.Output, _, _) |
536      Instantiation.InstantiatedRealVariable
537        (_, Compilation.Output, _, _, _) ->
538        true
539    | Instantiation.InstantiatedIntegerVariable _ |
540      Instantiation.InstantiatedStringVariable _ |
541      Instantiation.InstantiatedDiscreteVariable _ |
542      Instantiation.InstantiatedRealVariable _ |
543      Instantiation.InstantiatedCompoundVariable _ -> false
544  in
545  let filter_variable = function
546    | _, Instantiation.InstantiatedParameter _ -> false
547    | _, Instantiation.InstantiatedVariable ivar -> is_output ivar
548  in List.partition filter_variable icpnts
549
550let create_dictionary get_contents maps icpnts =
551  let rec add_entries i map = function
552    | [] -> map
553    | (s, icpnt) :: icpnts ->
554        add_entries
555          (i + 1)
556          (StringMap.add s (get_contents maps s i icpnt) map)
557          icpnts
558  in add_entries 0 StringMap.empty icpnts
559
560let separate_whens_from_equations iequs =
561  let rec separate_whens_from_equations' whens equations = function
562    | [] -> whens, equations
563    | (Instantiation.Equation _ as iequ) :: iequs' ->
564        separate_whens_from_equations' whens (iequ :: equations) iequs'
565    | Instantiation.When iwhen_clauses :: iequs' ->
566        separate_whens_from_equations' (whens @ iwhen_clauses) equations iequs'
567    | Instantiation.ConditionalEquation _ :: _ |
568      Instantiation.FlowConnection _ :: _ -> assert false
569  in separate_whens_from_equations' [] [] iequs
570
571let rewrite_when_condition expr =
572  let rec contains_time expr = match nature expr with
573    | BooleanValue _  | Constant _ | DiscreteVariable _ | Number _ |
574      Parameter _ | Variable _ | Integer _ | String _ -> false
575    | TimeVariable  -> true
576    | ArcCosine node | ArcHyperbolicCosine node |
577      ArcHyperbolicSine node | ArcHyperbolicTangent node | ArcSine node |
578      ArcTangent node | Cosine node | Derivative (node, _) |
579      Exponential node | Floor node | HyperbolicCosine node |
580      HyperbolicSine node | HyperbolicTangent node | Logarithm node |
581      Not node | Pre node | RationalPower (node, _) | Sign node | Sine node |
582      Tangent node -> contains_time node
583    | Equality (node1, node2) | Greater (node1, node2) |
584      GreaterEqual (node1, node2) | PartialDerivative (node1, node2) ->
585        contains_time node1 || contains_time node2
586    | If (node1, node2, node3) ->
587        contains_time node1 || contains_time node2 || contains_time node3
588    | And nodes | Addition nodes | Multiplication nodes | Or nodes ->
589        List.exists contains_time nodes
590    | BlackBox (_, args) -> List.exists contains_time_argument args
591  and contains_time_argument = function
592    | ScalarArgument node -> contains_time node
593    | ArrayArgument (_, nodes) -> exists_array contains_time nodes in
594  let is_primary_condition expr =
595    assignable_variables_of expr <> [] || contains_time expr in
596  let rec rewrite expr = match nature expr with
597    | Addition nodes -> apply_addition (List.map rewrite nodes)
598    | And nodes -> apply_and (List.map rewrite nodes)
599    | ArcCosine node -> symbolic_acos (rewrite node)
600    | ArcHyperbolicCosine node -> symbolic_acosh (rewrite node)
601    | ArcHyperbolicSine node -> symbolic_asinh (rewrite node)
602    | ArcHyperbolicTangent node -> symbolic_atanh (rewrite node)
603    | ArcSine node -> symbolic_asin (rewrite node)
604    | ArcTangent node -> symbolic_atan (rewrite node)
605    | BlackBox (s, args) -> apply_blackBox s (List.map rewrite_argument args)
606    | Cosine node -> symbolic_cos (rewrite node)
607    | Derivative (node, num) -> symbolic_derive (rewrite node) num
608    | Equality (node, node') -> symbolic_eq (rewrite node) (rewrite node')
609    | Exponential node -> symbolic_exp (rewrite node)
610    | Floor node -> symbolic_floor (rewrite node)
611    | Greater (node, node') -> symbolic_gt (rewrite node) (rewrite node')
612    | GreaterEqual (node, node') -> symbolic_ge (rewrite node) (rewrite node')
613    | HyperbolicCosine node -> symbolic_cosh (rewrite node)
614    | HyperbolicSine node -> symbolic_sinh (rewrite node)
615    | HyperbolicTangent node -> symbolic_tanh (rewrite node)
616    | If (node, node', node'') ->
617        symbolic_if (rewrite node) (rewrite node') (rewrite node'')
618    | Logarithm node -> symbolic_log (rewrite node)
619    | Multiplication nodes -> apply_multiplication (List.map rewrite nodes)
620    | Not node -> symbolic_not (rewrite node)
621    | Or nodes -> apply_or (List.map rewrite nodes)
622    | PartialDerivative (node, node') ->
623        create_partialDerivative (rewrite node) (rewrite node')
624    | Pre node -> symbolic_pre (rewrite node)
625    | RationalPower (node, num) -> symbolic_rationalPower (rewrite node) num
626    | Sign node -> symbolic_sgn (rewrite node)
627    | Sine node -> symbolic_sin (rewrite node)
628    | Tangent node -> symbolic_tan (rewrite node)
629    | BooleanValue _ | Constant _ | Number _ | Parameter _ |
630      TimeVariable | Variable _ | Integer _ | String _ -> expr
631    | DiscreteVariable _ -> symbolic_pre expr
632  and rewrite_argument = function
633    | ScalarArgument expr -> ScalarArgument (rewrite expr)
634    | ArrayArgument (dims, exprs) ->
635        ArrayArgument (dims, Array.map rewrite exprs) in
636  if is_primary_condition expr then rewrite expr else expr
637
638let symbolic_equation inl_par maps = function
639  | Instantiation.Equation (iexpr, iexpr') ->
640      let expr =
641        symbolic_sub
642          (symbolic_expression_of_expression inl_par maps iexpr)
643          (symbolic_expression_of_expression inl_par maps iexpr')
644      in
645      {
646        solved = false;
647        inner_variables = variables_of expr;
648        inner_derivatives = derivatives_of expr;
649        assignable_variables = assignable_variables_of expr;
650        expression = expr
651      }
652  | _ -> assert false
653
654let symbolic_surfaces inl_par maps when_clauses =
655  List.map
656    (fun (iexpr, surfaces) ->
657      let expr = symbolic_expression_of_expression inl_par maps iexpr in
658      rewrite_when_condition expr,
659      List.map
660        (function
661          | Instantiation.Reinit (iexpr, iexpr') ->
662            let var = symbolic_expression_of_expression inl_par maps iexpr in
663            begin match nature var with
664              | Variable i ->
665                  Reinit (var, symbolic_expression_of_expression inl_par maps iexpr')
666              | _ -> assert false
667            end
668          | Instantiation.Assign (iexpr, iexpr') ->
669            let var = symbolic_expression_of_expression inl_par maps iexpr in
670            begin match nature var with
671              | DiscreteVariable i ->
672                  Assign (var, symbolic_expression_of_expression inl_par maps iexpr')
673              | Variable i ->
674                  Reinit (var, symbolic_expression_of_expression inl_par maps iexpr')
675              | _ -> assert false
676            end)
677        surfaces)
678    when_clauses
679
680let propagate_noEvent expr =
681  (* such that 'noEvent' only appears in conditions *)
682  let rec propagate_noEvent' no_event expr = match nature expr with
683    | And exprs' ->
684        create_and (sort (List.map (propagate_noEvent' no_event) exprs'))
685    | ArcCosine expr' -> create_arcCosine (propagate_noEvent' no_event expr')
686    | ArcHyperbolicCosine expr' ->
687        create_arcHyperbolicCosine (propagate_noEvent' no_event expr')
688    | ArcHyperbolicSine expr' ->
689        create_arcHyperbolicSine (propagate_noEvent' no_event expr')
690    | ArcHyperbolicTangent expr' ->
691        create_arcHyperbolicTangent (propagate_noEvent' no_event expr')
692    | ArcSine expr' -> create_arcSine (propagate_noEvent' no_event expr')
693    | ArcTangent expr' -> create_arcTangent (propagate_noEvent' no_event expr')
694    | Cosine expr' -> create_cosine (propagate_noEvent' no_event expr')
695    | Derivative (expr', num) ->
696        create_derivative (propagate_noEvent' no_event expr') num
697    | Equality (expr', expr'') ->
698        create_equality
699          (propagate_noEvent' no_event expr')
700          (propagate_noEvent' no_event expr'')
701    | Exponential expr' ->
702        create_exponential (propagate_noEvent' no_event expr')
703    | Floor expr' -> create_floor (propagate_noEvent' no_event expr')
704    | Greater (expr', expr'') ->
705        create_greater
706          (propagate_noEvent' no_event expr')
707          (propagate_noEvent' no_event expr'')
708    | GreaterEqual (expr', expr'') ->
709        create_greater_equal
710          (propagate_noEvent' no_event expr')
711          (propagate_noEvent' no_event expr'')
712    | HyperbolicCosine expr' ->
713        create_hyperbolicCosine (propagate_noEvent' no_event expr')
714    | HyperbolicSine expr' ->
715        create_hyperbolicSine (propagate_noEvent' no_event expr')
716    | HyperbolicTangent expr' ->
717        create_hyperbolicTangent (propagate_noEvent' no_event expr')
718    | Logarithm expr' -> create_logarithm (propagate_noEvent' no_event expr')
719    | Pre expr' -> create_pre (propagate_noEvent' no_event expr')
720    | RationalPower (expr', num) ->
721        create_rationalPower (propagate_noEvent' no_event expr') num
722    | Sign expr' -> create_sign (propagate_noEvent' no_event expr')
723    | Sine expr' -> create_sine (propagate_noEvent' no_event expr')
724    | Tangent expr' -> create_tangent (propagate_noEvent' no_event expr')
725    | Addition exprs' ->
726        create_addition (sort (List.map (propagate_noEvent' no_event) exprs'))
727    | BlackBox ("noEvent", [ScalarArgument expr']) -> propagate_noEvent' true expr'
728    | BlackBox (name, args) ->
729        create_blackBox
730          name
731          (List.map (propagate_noEvent_into_argument no_event) args)
732    | Multiplication exprs' ->
733        create_multiplication
734          (sort (List.map (propagate_noEvent' no_event) exprs'))
735    | Not expr' -> create_not (propagate_noEvent' no_event expr')
736    | Or exprs' ->
737        create_or (sort (List.map (propagate_noEvent' no_event) exprs'))
738    | PartialDerivative (expr', expr'') ->
739        create_partialDerivative
740          (propagate_noEvent' no_event expr')
741          (propagate_noEvent' no_event expr'')
742    | If (expr', expr'', expr''') ->
743        propagate_noEvent_into_if no_event expr' expr'' expr'''
744    | BooleanValue _ | Constant _ | DiscreteVariable _ | Number _ |
745      Parameter _ | TimeVariable | Variable _ | Integer _ | String _ -> expr
746  and propagate_noEvent_into_argument no_event = function
747    | ScalarArgument expr -> ScalarArgument (propagate_noEvent' no_event expr)
748    | ArrayArgument (dims, exprs) ->
749        ArrayArgument (dims, Array.map (propagate_noEvent' no_event) exprs)
750  and propagate_noEvent_into_if no_event expr expr' expr'' =
751    let cond =
752      if no_event then
753        create_blackBox "noEvent" [ScalarArgument (propagate_noEvent' false expr)]
754      else begin match nature expr with
755        | BlackBox ("noEvent", [ScalarArgument expr']) ->
756            create_blackBox
757              "noEvent"
758              [ScalarArgument (propagate_noEvent' false expr)]
759        | _ -> propagate_noEvent' false expr
760      end
761    in
762    create_if
763      cond
764      (propagate_noEvent' no_event expr')
765      (propagate_noEvent' no_event expr'')
766  in propagate_noEvent' false expr
767
768let create_model' trace inl_par iexpr =
769  let lazy_symbolic_expression_of_expression maps iexpr =
770    match iexpr.Instantiation.tex_expression with
771    | None -> None
772    | Some _ ->
773        Some (lazy (symbolic_expression_of_expression inl_par maps iexpr))
774  in
775  let get_parameter_info maps s i = function
776    | Instantiation.InstantiatedParameter (
777        Instantiation.InstantiatedIntegerParameter (s', kind, iexpr, _)) ->
778        {
779          parameter_kind = kind;
780          parameter_type = IntegerType;
781          parameter_id = i;
782          parameter_comment = s';
783          parameter_start = lazy
784            (symbolic_expression_of_expression inl_par maps iexpr)
785        }
786    | Instantiation.InstantiatedParameter (
787        Instantiation.InstantiatedStringParameter (s', kind, iexpr, _)) ->
788        {
789          parameter_kind = kind;
790          parameter_type = StringType;
791          parameter_id = i;
792          parameter_comment = s';
793          parameter_start = lazy
794            (symbolic_expression_of_expression inl_par maps iexpr)
795        }
796    | Instantiation.InstantiatedParameter (
797        Instantiation.InstantiatedRealParameter (s', kind, iexpr, _)) ->
798        {
799          parameter_kind = kind;
800          parameter_type = RealType;
801          parameter_id = i;
802          parameter_comment = s';
803          parameter_start = lazy
804            (symbolic_expression_of_expression inl_par maps iexpr)
805        }
806    | _ -> assert false
807  and get_input_info maps s i = function
808    | Instantiation.InstantiatedVariable (
809        Instantiation.InstantiatedDiscreteVariable (s', _, _, _)) |
810      Instantiation.InstantiatedVariable (
811        Instantiation.InstantiatedRealVariable (s', _, _, _, _)) ->
812        {
813          input_id = i;
814          input_name = s;
815          input_comment = s'
816        }
817    | _ -> assert false
818  and get_variable_info maps s i = function
819    | Instantiation.InstantiatedVariable (
820        Instantiation.InstantiatedDiscreteVariable (s', _, iexpr, _)) |
821      Instantiation.InstantiatedVariable (
822        Instantiation.InstantiatedRealVariable (s', _, _, iexpr, _)) ->
823        {
824          variable_id = i;
825          variable_comment = s';
826          variable_start = lazy_symbolic_expression_of_expression maps iexpr
827        }
828    | _ -> assert false
829  in
830  let derived_variables ders =
831    List.fold_left
832      (fun vars der ->
833        match nature der with
834          | Derivative (expr, num) when num = Num.Int 1 ->
835              begin match nature expr with
836                | Variable _ -> expr :: vars
837                | _ -> assert false
838              end
839          | _ -> assert false)
840      []
841      ders
842  in
843  let add_variable_infos acc (id, icpnt) = match icpnt with
844    | Instantiation.InstantiatedVariable
845        (Instantiation.InstantiatedDiscreteVariable (_, _, _, var_infos) |
846         Instantiation.InstantiatedRealVariable (_, _, _, _, var_infos)) ->
847        (id, var_infos) :: acc
848    | Instantiation.InstantiatedParameter
849        (Instantiation.InstantiatedIntegerParameter (_, _, _, pat_infos) |
850         Instantiation.InstantiatedRealParameter (_, _, _, pat_infos)) ->
851        (id, pat_infos) :: acc
852    | _ -> acc
853  in
854  let icpnts, iinit_equs, iequs = Instantiation.expand_class iexpr in
855  let variables_infos =
856    List.fold_left add_variable_infos [] icpnts in
857  let parameters, variables = separate_parameters_from_variables icpnts in
858  let inputs, non_inputs = separate_inputs_from_others variables in
859  let discrete_variables, others =
860    separate_discrete_variables_from_others non_inputs in
861  let outputs, _ = separate_outputs_from_others non_inputs in
862  let name_and_types =
863    collect_external_function_names iequs @
864    collect_external_function_names iinit_equs in
865  let rec maps =
866    {
867      parameters_map =
868        lazy (create_dictionary get_parameter_info maps parameters);
869      inputs_map =
870        lazy (create_dictionary get_input_info maps inputs);
871      discrete_variables_map =
872       lazy (create_dictionary get_variable_info maps discrete_variables);
873      variables_map =
874        lazy (create_dictionary get_variable_info maps others)
875    }
876  in
877  let when_clauses, equations = separate_whens_from_equations iequs in
878  let nb_parameters = map_length (Lazy.force maps.parameters_map)
879  and nb_inputs = map_length (Lazy.force maps.inputs_map)
880  and nb_discrete_vars = map_length (Lazy.force maps.discrete_variables_map)
881  and nb_vars = map_length (Lazy.force maps.variables_map)
882  and nb_equs = List.length equations in
883  if nb_equs <> nb_vars then
884    failwith
885      ("The number of equations doesn't match the number of variables: " ^
886      string_of_int nb_equs ^ " equations and " ^ string_of_int nb_vars ^
887      " variables.")
888  else
889    let parameters_array =
890      Array.init
891        nb_parameters
892        (fun _ ->
893          {
894            main = false;
895            p_type = IntegerType;
896            p_name = "";
897            p_comment = "";
898            value = zero
899          })
900    and inputs_array = Array.make nb_inputs ""
901    and discrete_variables_array =
902      Array.init
903        nb_discrete_vars
904        (fun _ ->
905          {
906            d_output = None;
907            d_v_name = "";
908            d_v_comment = "";
909            d_start_value = Some zero
910          })
911    and variables_array =
912      Array.init
913        nb_vars
914        (fun _ ->
915          {
916            output = None;
917            state = true;
918            v_name = "";
919            depth_in_hierarchy = 0;
920            v_comment = "";
921            start_value = Some zero
922          })
923    and equations_array =
924      Array.init
925        nb_equs
926        (fun _ ->
927          {
928            solved = false;
929            inner_variables = [];
930            inner_derivatives = [];
931            assignable_variables = [];
932            expression = zero
933          })
934    in
935    let output_index s outputs =
936      let rec output_index' i = function
937        | [] -> None
938        | (s', _) :: _ when s' = s -> Some i
939        | _ :: outputs' -> output_index' (i + 1) outputs'
940      in output_index' 0 outputs
941    in
942    let number_of_dots s =
943      let count = ref 0 in
944      for i = 0 to String.length s - 1 do
945        if s.[i] = '.' then incr count;
946      done;
947      !count
948    in
949    let _ =
950      List.fold_left
951        (fun i equ ->
952          assert (i < Array.length equations_array);
953          equations_array.(i) <- symbolic_equation inl_par maps equ; i + 1)
954        0
955        equations
956    in ();
957    let derived_variables =
958      Array.fold_left
959        (fun vars equation ->
960          union vars (derived_variables equation.inner_derivatives))
961        []
962        equations_array
963    in
964    StringMap.iter
965      (fun s param ->
966        assert (param.parameter_id < Array.length parameters_array);
967        let parameter = parameters_array.(param.parameter_id) in
968        parameter.main <- param.parameter_kind = Instantiation.Main;
969        parameter.p_type <- param.parameter_type;
970        parameter.p_name <- s;
971        parameter.p_comment <- param.parameter_comment;
972        parameter.value <- Lazy.force param.parameter_start)
973      (Lazy.force maps.parameters_map);
974    StringMap.iter
975      (fun _ inp ->
976        assert (inp.input_id < Array.length inputs_array);
977        inputs_array.(inp.input_id) <- inp.input_name)
978      (Lazy.force maps.inputs_map);
979    StringMap.iter
980      (fun s dvar ->
981        assert (dvar.variable_id < Array.length discrete_variables_array);
982        let variable = discrete_variables_array.(dvar.variable_id) in
983        variable.d_output <- output_index s outputs;
984        variable.d_v_name <- s;
985        variable.d_v_comment <- dvar.variable_comment;
986        variable.d_start_value <-
987          match dvar.variable_start with
988            | None -> None
989            | Some lexpr -> Some (Lazy.force lexpr))
990      (Lazy.force maps.discrete_variables_map);
991    StringMap.iter
992      (fun s var ->
993        assert (var.variable_id < Array.length variables_array);
994        let variable = variables_array.(var.variable_id) in
995        variable.output <- output_index s outputs;
996        variable.state <-
997          List.memq (create_variable var.variable_id) derived_variables;
998        variable.v_name <- s;
999        variable.depth_in_hierarchy <- number_of_dots s;
1000        variable.v_comment <- var.variable_comment;
1001        variable.start_value <-
1002          match var.variable_start with
1003            | None -> None
1004            | Some lexpr -> Some (Lazy.force lexpr))
1005      (Lazy.force maps.variables_map);
1006    let when_clauses_list = symbolic_surfaces inl_par maps when_clauses in
1007    let reinitializable_variables =
1008      let add_non_discrete_variables vars = function
1009        | Reinit (var, _) when not (List.memq var vars) -> var :: vars
1010        | _ -> vars
1011      in
1012      List.fold_left
1013        (fun vars (_, when_expr) ->
1014          List.fold_left add_non_discrete_variables vars when_expr)
1015        []
1016        when_clauses_list
1017    in
1018    Array.iter
1019      (fun equation ->
1020        equation.expression <- propagate_noEvent equation.expression)
1021      equations_array;
1022    {
1023      parameters = parameters_array;
1024      inputs = inputs_array;
1025      discrete_variables = discrete_variables_array;
1026      variables = variables_array;
1027      equations = equations_array;
1028      reinitializable_variables = reinitializable_variables;
1029      when_clauses = when_clauses_list;
1030      io_dependency = false;
1031      external_functions = name_and_types;
1032      trace = trace;
1033      variables_infos = variables_infos
1034    }
1035
1036let create_model_with_parameters trace iexpr = create_model' trace false iexpr
1037
1038let create_model trace iexpr = create_model' trace true iexpr
1039
1040let print_model oc model =
1041  Printf.fprintf
1042    oc
1043    "Number of variables before simplifications: %d\n"
1044    (Array.length model.variables);
1045  Printf.fprintf
1046    oc
1047    "Number of variables after simplifications: %d\n"
1048    (Array.fold_left
1049      (fun n equation -> if equation.solved then n else n + 1)
1050      0
1051      model.equations);
1052  Printf.fprintf
1053    oc
1054    "Direct input/ouput dependency: %s\n"
1055    (if model.io_dependency then "yes" else "no");
1056  Array.iteri
1057    (fun i variable ->
1058      Printf.fprintf
1059        oc
1060        "variable (%d) (%s): %s %s variable (%ssolved)\n"
1061        i
1062        variable.v_name
1063        (if variable.output <> None then "output" else "intermediate")
1064        (if variable.state then "state" else "algebraic")
1065        (if model.equations.(i).solved then "" else "not "))
1066    model.variables;
1067  Array.iteri
1068    (fun i equation ->
1069      Printf.fprintf oc "equation(%d): " i;
1070      if equation.solved then output oc (create_variable i)
1071      else Printf.fprintf oc "0";
1072      Printf.fprintf oc " = ";
1073      output oc equation.expression;
1074      Printf.fprintf oc "\n")
1075    model.equations
1076
1077let create_index_array a p =
1078  let size = Array.length a in
1079  let indexes = Array.make size (-1) in
1080  let j = ref 0 in
1081  Array.iteri (fun i x -> if p x then begin indexes.(i) <- !j; incr j end) a;
1082  indexes
1083
1084let final_index_of_integer_parameters model =
1085  create_index_array model.parameters (fun par -> par.p_type = IntegerType)
1086
1087let final_index_of_string_parameters model =
1088  create_index_array model.parameters (fun par -> par.p_type = StringType)
1089
1090let final_index_of_real_parameters model =
1091  create_index_array model.parameters (fun par -> par.p_type = RealType)
1092
1093let final_index_of_variables model =
1094  create_index_array model.equations (fun equation -> not equation.solved)
1095
1096let permute_equations model assocs =
1097  let equations = Array.copy model.equations in
1098  List.iter
1099    (function
1100      | i, Some j -> equations.(i) <- model.equations.(j)
1101      | _, None -> assert false)
1102    assocs;
1103  model.equations <- equations
1104
1105let perform_then_propagate_inversion model i =
1106  let update_clauses var expr' clauses =
1107    List.map
1108      (fun (expr, updates) ->
1109        replace var expr' expr,
1110        List.map
1111          (function
1112            | Assign (expr1, expr2) ->
1113                Assign (expr1, replace var expr' expr2)
1114            | Reinit (expr1, expr2) ->
1115                Reinit (expr1, replace var expr' expr2))
1116          updates)
1117      clauses
1118  in
1119  let var = create_variable i
1120  and equation = model.equations.(i) in
1121  let expr = equation.expression in
1122  try match invert_if_possible_with_respect_to var expr zero with
1123    | None -> false
1124    | Some expr' ->
1125        if not (List.memq var model.reinitializable_variables) then begin
1126          equation.expression <- expr';
1127          equation.solved <- true;
1128          let additional_variables = variables_of expr' in
1129          Array.iteri
1130            (fun j equation ->
1131              if i <> j && List.memq var equation.inner_variables then begin
1132                equation.expression <- replace var expr' equation.expression;
1133                equation.inner_variables <-
1134                  union equation.inner_variables additional_variables
1135              end)
1136            model.equations;
1137          model.when_clauses <- update_clauses var expr' model.when_clauses
1138        end;
1139        true
1140  with Invalid_argument _ -> raise Can't_perform_inversion
1141
1142let compute_io_dependency model =
1143  let rec compute_io_dependency' i =
1144    if not (i = Array.length model.variables) then
1145      if
1146        model.variables.(i).output <> None &&
1147        exists
1148          (fun node -> match nature node with
1149            | DiscreteVariable i when i < 0 -> true
1150            | _ -> false)
1151          model.equations.(i).expression
1152      then
1153        model.io_dependency <- true
1154      else
1155        compute_io_dependency' (i + 1)
1156  in compute_io_dependency' 0
1157
1158let perform_hungarian_method model =
1159  let size =
1160    Array.fold_left
1161      (fun acc equation -> if equation.solved then acc else acc + 1)
1162      0
1163      model.equations
1164  in
1165  let () =
1166    Array.iter
1167      (fun equation ->
1168        if not equation.solved then begin
1169          equation.inner_variables <- variables_of equation.expression;
1170          equation.inner_derivatives <- derivatives_of equation.expression;
1171          equation.assignable_variables <-
1172            assignable_variables_of equation.expression
1173        end)
1174      model.equations
1175  in
1176  let table = Array.make size 0 in
1177  let i = ref 0 in
1178  for j = 0 to Array.length model.equations - 1 do
1179    if not model.equations.(j).solved then begin
1180      table.(!i) <- j; incr i
1181    end
1182  done;
1183  let weight i j =
1184    let m = table.(i)
1185    and n = table.(j) in
1186    let var = create_variable m in
1187    let discontinuous_state () =
1188      let not_derivative_of_current_variable e = match nature e with
1189        | Derivative (var', _) -> var' != var
1190        | _ -> true in
1191      List.for_all
1192        not_derivative_of_current_variable
1193        model.equations.(n).inner_derivatives &&
1194      List.memq var model.reinitializable_variables in
1195    if not (List.memq var model.equations.(n).assignable_variables) then
1196      IntegerElement.Infinity
1197    else if discontinuous_state () then
1198      IntegerElement.Int (size + 1)
1199    else match inversion_difficulty var model.equations.(n).expression zero with
1200        | 0 -> IntegerElement.zero
1201        | 1 -> IntegerElement.Int 1
1202        | 2 -> IntegerElement.Int (size + 1)
1203        | _ -> IntegerElement.Infinity
1204  in
1205  let strct = IntegerHungarianMethod.init size weight in
1206  let assocs = IntegerHungarianMethod.perform strct in
1207  let assocs' =
1208    List.map
1209      (function
1210        | i, Some j -> table.(i), Some table.(j)
1211        | _, None ->
1212            failwith "perform_hungarian_method: jacobian is structurally singular")
1213      assocs
1214  in
1215  assert (
1216    let rec check_results1 = function
1217      | [] -> true
1218      | (_, x) :: xs when List.exists (fun (_, y) -> x = y) xs -> false
1219      | _ :: xs -> check_results1 xs
1220    in check_results1 assocs');
1221  assert (
1222    let check_results2 = function
1223      | (i, Some j) :: assocs ->
1224          let var = create_variable i in
1225          List.memq var model.equations.(j).inner_variables
1226      | _ -> true
1227    in check_results2 assocs');
1228  assocs'
1229
1230let eliminate_trivial_relations max_simplifs model =
1231  let max_simplifs_ref = ref max_simplifs in
1232  let choose_variable i j =
1233    let sti = model.variables.(i).state
1234    and depi = model.variables.(i).depth_in_hierarchy
1235    and stj = model.variables.(j).state
1236    and depj = model.variables.(j).depth_in_hierarchy in
1237    match sti, stj with
1238      | true, false -> j
1239      | false, true -> i
1240      | true, true | false, false ->
1241          let svi = model.variables.(i).start_value
1242          and svj = model.variables.(j).start_value in
1243          begin match svi, svj with
1244            | None, Some _ -> i
1245            | Some _, None -> j
1246            | _ when depj > depi -> i
1247            | _ -> j
1248          end
1249  in
1250  let permute_equations i j =
1251    let equation = model.equations.(i) in
1252    model.equations.(i) <- model.equations.(j);
1253    model.equations.(j) <- equation
1254  in
1255  let update_variable_attributes i j =
1256    let svi = model.variables.(i).start_value
1257    and svj = model.variables.(j).start_value in
1258    match svi, svj with
1259      | None, _ -> model.variables.(i).start_value <- svj
1260      | _, None -> model.variables.(j).start_value <- svi
1261      | _ -> ()
1262  in
1263  let simplify_trivial_relation n =
1264    match nature model.equations.(n).expression with
1265      | Addition [node; node'] when !max_simplifs_ref >= 0 ->
1266          begin match nature node, nature node' with
1267            | Variable i, Number _ | Number _, Variable i ->
1268                permute_equations i n;
1269                ignore (perform_then_propagate_inversion model i);
1270                decr max_simplifs_ref
1271            | Variable i, Multiplication [node; node'] |
1272              Multiplication [node; node'], Variable i ->
1273                begin match nature node, nature node' with
1274                  | Number (Num.Int (-1)), Variable j |
1275                    Variable j, Number (Num.Int (-1)) ->
1276                      let k = choose_variable i j in
1277                      update_variable_attributes i j;
1278                      permute_equations k n;
1279                      ignore (perform_then_propagate_inversion model k);
1280                      decr max_simplifs_ref
1281                  | _ -> ()
1282                end
1283            | _ -> ()
1284          end
1285      | _ -> ()
1286  in
1287  for i = 0 to Array.length model.equations - 1 do
1288    simplify_trivial_relation i
1289  done;
1290  !max_simplifs_ref
1291
1292let eliminate_explicit_variables max_simplifs model =
1293  let is_continuous_state var =
1294    not (List.memq var model.reinitializable_variables) in
1295  let is_state_variable expr = match nature expr with
1296    | Variable j -> model.variables.(j).state
1297    | _ -> false in
1298  let try_reducing_state var i =
1299    let vars = model.equations.(i).assignable_variables
1300    and ders = model.equations.(i).inner_derivatives in
1301    if not (List.memq (symbolic_derive var (Num.Int 1)) ders) then begin
1302      Printf.eprintf "\nTrying to reduce state...";
1303      match ders with
1304      | [] when List.for_all is_state_variable vars ->
1305          let success = perform_then_propagate_inversion model i in
1306          if success then Printf.eprintf " One state variable removed.%!"
1307          else Printf.eprintf " Failed.%!"
1308      | _ ->  Printf.eprintf " Failed.%!"
1309    end in
1310  let rec eliminate_explicit_variables' simplifs =
1311    let assocs = perform_hungarian_method model in
1312    permute_equations model assocs;
1313    let bad_variable_choice, success, simplifs =
1314      List.fold_left
1315        (fun (bad_variable_choice, success, simplifs) assoc ->
1316          match assoc with
1317            | _, None -> assert false
1318            | i, Some j when simplifs >= 0 ->
1319                begin try
1320                  if not model.variables.(i).state then
1321                    ignore (perform_then_propagate_inversion model i)
1322                  else begin
1323                    let var = create_variable i in
1324                    if is_continuous_state var then try_reducing_state var i
1325                  end;
1326                  bad_variable_choice, model.equations.(i).solved, simplifs - 1
1327                with
1328                  | Can't_perform_inversion -> true, success, simplifs
1329                end
1330            | _ -> bad_variable_choice, success, simplifs)
1331        (false, false, simplifs)
1332        assocs
1333    in
1334    if bad_variable_choice || success then
1335      eliminate_explicit_variables' simplifs
1336  in eliminate_explicit_variables' max_simplifs
1337
1338let rec rewrite_conditions_in no_event expr =
1339  let rec remove_no_event expr = match nature expr with
1340    | ArcCosine expr' -> create_arcCosine (remove_no_event expr')
1341    | ArcHyperbolicCosine expr' ->
1342        create_arcHyperbolicCosine (remove_no_event expr')
1343    | ArcHyperbolicSine expr' ->
1344        create_arcHyperbolicSine (remove_no_event expr')
1345    | ArcHyperbolicTangent expr' ->
1346        create_arcHyperbolicTangent (remove_no_event expr')
1347    | ArcSine expr' -> create_arcSine (remove_no_event expr')
1348    | ArcTangent expr' -> create_arcTangent (remove_no_event expr')
1349    | Cosine expr' -> create_cosine (remove_no_event expr')
1350    | Derivative (expr', num) ->
1351        create_derivative (remove_no_event expr') num
1352    | Exponential expr' -> create_exponential (remove_no_event expr')
1353    | Floor expr' -> create_floor (remove_no_event expr')
1354    | HyperbolicCosine expr' ->
1355        create_hyperbolicCosine (remove_no_event expr')
1356    | HyperbolicSine expr' ->
1357        create_hyperbolicSine (remove_no_event expr')
1358    | HyperbolicTangent expr' ->
1359        create_hyperbolicTangent (remove_no_event expr')
1360    | Logarithm expr' -> create_logarithm (remove_no_event expr')
1361    | Pre expr' -> create_pre (remove_no_event expr')
1362    | RationalPower (expr', num) ->
1363        create_rationalPower (remove_no_event expr') num
1364    | Sign expr' -> create_sign (remove_no_event expr')
1365    | Sine expr' -> create_sine (remove_no_event expr')
1366    | Tangent expr' -> create_tangent (remove_no_event expr')
1367    | Addition exprs' ->
1368        create_addition (sort (List.map remove_no_event exprs'))
1369    | BlackBox ("noEvent", [ScalarArgument expr']) ->
1370        rewrite_conditions_in true expr'
1371    | BlackBox (name, args) ->
1372        create_blackBox name (List.map remove_no_event_in_argument args)
1373    | Multiplication exprs' ->
1374        create_multiplication (sort (List.map remove_no_event exprs'))
1375    | PartialDerivative (expr', expr'') ->
1376        create_partialDerivative
1377          (remove_no_event expr')
1378          (remove_no_event expr'')
1379    | If (expr', expr'', expr''') ->
1380        create_if
1381          (remove_no_event expr')
1382          (remove_no_event expr'')
1383          (remove_no_event expr''')
1384    | Equality (expr', expr'') ->
1385        create_equality
1386          (remove_no_event expr')
1387          (remove_no_event expr'')
1388    | Greater (expr', expr'') ->
1389        create_greater
1390          (remove_no_event expr')
1391          (remove_no_event expr'')
1392    | GreaterEqual (expr', expr'') ->
1393        create_greater_equal
1394          (remove_no_event expr')
1395          (remove_no_event expr'')
1396    | And exprs -> create_and (List.map remove_no_event exprs)
1397    | Or exprs ->  create_or (List.map remove_no_event exprs)
1398    | Not expr' -> create_not (remove_no_event expr')
1399    | BooleanValue _ | Constant _ | DiscreteVariable _ | Number _ |
1400      Parameter _ | TimeVariable | Variable _ | Integer _ | String _ -> expr
1401  and remove_no_event_in_argument = function
1402    | ScalarArgument expr -> ScalarArgument (remove_no_event expr)
1403    | ArrayArgument (dims, exprs) ->
1404        ArrayArgument (dims, Array.map remove_no_event exprs) in
1405  match nature expr with
1406    | ArcCosine expr' -> create_arcCosine (rewrite_conditions_in no_event expr')
1407    | ArcHyperbolicCosine expr' ->
1408        create_arcHyperbolicCosine (rewrite_conditions_in no_event expr')
1409    | ArcHyperbolicSine expr' ->
1410        create_arcHyperbolicSine (rewrite_conditions_in no_event expr')
1411    | ArcHyperbolicTangent expr' ->
1412        create_arcHyperbolicTangent (rewrite_conditions_in no_event expr')
1413    | ArcSine expr' -> create_arcSine (rewrite_conditions_in no_event expr')
1414    | ArcTangent expr' -> create_arcTangent (rewrite_conditions_in no_event expr')
1415    | Cosine expr' -> create_cosine (rewrite_conditions_in no_event expr')
1416    | Derivative (expr', num) ->
1417        create_derivative (rewrite_conditions_in no_event expr') num
1418    | Exponential expr' -> create_exponential (rewrite_conditions_in no_event expr')
1419    | Floor expr' -> create_floor (rewrite_conditions_in no_event expr')
1420    | HyperbolicCosine expr' ->
1421        create_hyperbolicCosine (rewrite_conditions_in no_event expr')
1422    | HyperbolicSine expr' ->
1423        create_hyperbolicSine (rewrite_conditions_in no_event expr')
1424    | HyperbolicTangent expr' ->
1425        create_hyperbolicTangent (rewrite_conditions_in no_event expr')
1426    | Logarithm expr' -> create_logarithm (rewrite_conditions_in no_event expr')
1427    | Pre expr' -> create_pre (rewrite_conditions_in no_event expr')
1428    | RationalPower (expr', num) ->
1429        create_rationalPower (rewrite_conditions_in no_event expr') num
1430    | Sign expr' -> create_sign (rewrite_conditions_in no_event expr')
1431    | Sine expr' -> create_sine (rewrite_conditions_in no_event expr')
1432    | Tangent expr' -> create_tangent (rewrite_conditions_in no_event expr')
1433    | Addition exprs' ->
1434        create_addition (sort (List.map (rewrite_conditions_in no_event) exprs'))
1435    | BlackBox ("noEvent", [ScalarArgument expr']) ->
1436        rewrite_conditions_in true expr'
1437    | BlackBox (name, args) ->
1438        create_blackBox
1439          name
1440          (List.map (rewrite_conditions_in_argument no_event) args)
1441    | Multiplication exprs' ->
1442        create_multiplication (sort (List.map (rewrite_conditions_in no_event) exprs'))
1443    | PartialDerivative (expr', expr'') ->
1444        create_partialDerivative
1445          (rewrite_conditions_in no_event expr')
1446          (rewrite_conditions_in no_event expr'')
1447    | If (expr', expr'', expr''') ->
1448        create_if
1449          (rewrite_conditions_in no_event expr')
1450          (rewrite_conditions_in no_event expr'')
1451          (rewrite_conditions_in no_event expr''')
1452    | Equality (expr', expr'') when no_event ->
1453        create_blackBox
1454          "noEvent"
1455          [ScalarArgument
1456            (create_equality
1457              (remove_no_event expr')
1458              (remove_no_event expr''))]
1459    | Equality (expr', expr'') ->
1460        create_equality
1461          (remove_no_event expr')
1462          (remove_no_event expr'')
1463    | Greater (expr', expr'') when no_event ->
1464        create_blackBox
1465          "noEvent"
1466          [ScalarArgument
1467            (create_greater
1468              (remove_no_event expr')
1469              (remove_no_event expr''))]
1470    | Greater (expr', expr'') ->
1471        create_greater
1472          (remove_no_event expr')
1473          (remove_no_event expr'')
1474    | GreaterEqual (expr', expr'') when no_event ->
1475        create_blackBox
1476          "noEvent"
1477          [ScalarArgument
1478            (create_greater_equal
1479              (remove_no_event expr')
1480              (remove_no_event expr''))]
1481    | GreaterEqual (expr', expr'') ->
1482        create_greater_equal
1483          (remove_no_event expr')
1484          (remove_no_event expr'')
1485    | And exprs when no_event ->
1486        create_blackBox
1487          "noEvent"
1488          [ScalarArgument (create_and (List.map remove_no_event exprs))]
1489    | And exprs -> create_and (List.map remove_no_event exprs)
1490    | Or exprs when no_event ->
1491        create_blackBox
1492          "noEvent"
1493          [ScalarArgument (create_or (List.map remove_no_event exprs))]
1494    | Or exprs -> create_or (List.map remove_no_event exprs)
1495    | Not expr' when no_event ->
1496        create_blackBox
1497          "noEvent"
1498          [ScalarArgument (create_not (remove_no_event expr'))]
1499    | Not expr' -> create_not (remove_no_event expr')
1500    | BooleanValue _ | Constant _ | DiscreteVariable _ | Number _ |
1501      Parameter _ | TimeVariable | Variable _ | Integer _ | String _ -> expr
1502  and rewrite_conditions_in_argument no_event = function
1503    | ScalarArgument expr ->
1504        ScalarArgument (rewrite_conditions_in no_event expr)
1505    | ArrayArgument (dims, exprs) ->
1506        ArrayArgument
1507          (dims, Array.map (rewrite_conditions_in no_event) exprs)
1508
1509let perform_simplifications max_simplifs model =
1510  Array.iter
1511    (fun equation ->
1512      equation.expression <- rewrite_conditions_in false equation.expression)
1513    model.equations;
1514  eliminate_explicit_variables max_simplifs model;
1515  compute_io_dependency model
1516
1517let find_submodels model =
1518  let final_index_of_variables = final_index_of_variables model in
1519  let size =
1520    Array.fold_left
1521      (fun acc i -> if i >= 0 then acc + 1 else acc)
1522      0
1523      final_index_of_variables
1524  in
1525  let graph = CausalityGraph.create size in
1526  Array.iteri
1527    (fun i equation ->
1528      if not equation.solved then
1529        List.iter
1530          (fun expr ->
1531            match nature expr with
1532              | Variable j ->
1533                  if not model.equations.(j).solved then
1534                    let i' = final_index_of_variables.(i)
1535                    and j' = final_index_of_variables.(j) in
1536                    CausalityGraph.connect i' j' graph
1537              | _ -> assert false)
1538          equation.inner_variables)
1539    model.equations;
1540  CausalityGraph.strongly_connected_components graph
1541