text:
######################################################################

# A full tree on A is a tree that contains A, and also contains {a}
# for all a in A.
# (Definition defn-full-trees)

`is_element/full_trees` := (A::set) -> proc(TT)
 local a;
 global reason;

 if not(`is_element/trees`(A)(TT)) then
  reason := [convert(procname,string),"TT is not a tree",TT,reason];
  return false;
 fi;
 if not(member(A,TT)) then
  reason := [convert(procname,string),"The full set A is not in TT",A,TT];
  return false;
 fi;

 for a in A do
  if not(member({a},TT)) then
   reason := [convert(procname,string),"The singleton {a} is not in TT",a,TT];
   return false;
  fi;
 od;

 return true;
end;

`is_equal/full_trees` := eval(`is_equal/trees`);
`is_leq/full_trees` := eval(`is_leq/trees`);

# A full tree TT on A is binary if every non-minimal set in TT has
# precisely two children.  This holds iff |TT| = 2|A|-1.

`is_binary/full_trees` := (A) -> proc(TT)
 return evalb(nops(TT) = 2*nops(A)-1);
end:

`list_elements/full_trees`:= proc(A::set)
 local n,a,A1,TTT1,TTT,TT1,T1;

 # Degenerate cases for small A
 if A={} then return []; fi;
 if nops(A)=1 then return [{A}]; fi;

 # If A is not too small, we write it as {a} union A1.  Then we recursively
 # choose a random tree on A1, and attach a to it.

 n := nops(A);
 a := A[n];
 A1 := A minus {a};
 TTT1 := `list_elements/full_trees`(A1);

 TTT := [];
 for TT1 in TTT1 do
  for T1 in TT1 do
   TTT := [op(TTT),`attach_to_tree`(A1)(TT1,T1,a,true)];
   if n > 2 and nops(T1) > 1 then
    TTT := [op(TTT),`attach_to_tree`(A1)(TT1,T1,a,false)];
   fi;
  od;
 od;

 return TTT;
end:

# A000311 from OEIS

`count_elements/full_trees` := proc(A::set)
 local n,i,j,k;
 n := nops(A);
 if n = 0 then
  return 0;
 elif n  = 1 then
  return 1;
 else
  return
   add(
    (n+k-1)!*
    add(
     1/(k-j)! *
     add(
     (2^i*(-1)^(i)*Stirling2(n+j-i-1,j-i))/((n+j-i-1)!*i!),
     i=0..j),
    j=1..k),
   k=1..n-1);
 fi;
end;

`random_element/full_trees`:= (A::set) -> proc()
 local n,a,A1,T0,T1,T,above;

 # Degenerate cases for small A
 if A={} then return FAIL; fi;
 if nops(A)=1 then return {A}; fi;

 # If A is not too small, we write it as {a} union A1.  Then we recursively
 # choose a random tree on A1, and attach a to it.

 n := nops(A);
 a := A[rand(1..n)()];
 A1 := A minus {a};
 if n = 2 then
  above := true;
 else
  above := evalb(rand(0..1)() = 1);
 fi;

 T0 := `random_element/full_trees`(A1)();
 T1 := select(U -> nops(U) > 1,T0);

 if above then
  T := T0[rand(1..nops(T0))()];
 else
  T := T1[rand(1..nops(T1))()];
 fi;

 T := `attach_to_tree`(A1)(T0,T,a,above);

 return T;
end;


text:
    {  Planning.mpl  }

    {  Aggregate production planning for 12 months  }

TITLE
    Production_Planning;

INDEX
    product  :=  1..3;

    month    :=  (January,February,March,April,May,June,July,
                  August,September,October,November,December);
DATA
    price[product]              :=  (105.09, 234.00, 800.00);   [$/Boxes]
    Demand[month,product]       :=  1000 DATAFILE(Demand.dat);  [Boxes/month]
    ProductionCapacity[product] :=  1000 (10, 42, 14);          [Boxes/month]
    ProductionCost[product]     :=  (94.0, 188.10, 653.20);     [$/Boxes]
    InventoryCost               :=  8.8;                        [$/Boxes/month]

DECISION VARIABLES
    Inventory[product,month]   ->  Invt;                        [Boxes]
    Production[product,month]  ->  Prod;                        [Boxes/month]
    Sales[product,month]       ->  Sale;                        [Boxes/month]

MACRO
    Revenues   :=  SUM(product,month: price * Sales);                [$]

    TotalCost  := SUM(product,month: ProductionCost[product] * Production[product,month]
                                    + InventoryCost * Inventory[product,month]);    [$]

MODEL

   MAX   Profit =  Revenues - TotalCost;

SUBJECT TO

   InvtBal[product,month] : Inventory  =  Inventory[month-1] + Production - Sales;

BOUNDS
    Sales <  Demand ;

    Production  <  ProductionCapacity ;

    Inventory[month=January..November]  <  90000 ;
    Inventory[month=December]   =  20000 ;

END

ext: mpl

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
#@ Not autoload

s_kochman_max := 10;  # max s in H_s(BP;\pi_t(S))
t_kochman_max := 11;  # max t in H_s(BP;\pi_t(S))
v_kochman_max := 10;  # work mostly mod this power of 2

E_kochman := table():


text:
read "../IdentifiabilityODE.mpl";

sys := [
diff(x3(t), t) = -b2*x3(t) + b2*x2(t),
diff(x1(t), t) = a2*x0(t) - a2*x1(t),
diff(x2(t), t) = (ka*kc*b1*x3(t) - ka*kc*b1*x2(t) + ka*b1*x0(t)*x3(t) - ka*b1*x0(t)*x2(t) - n*kc*x2(t) + kc*b1*x3(t)*x2(t) - kc*b1*x2(t)^2) / (ka*kc + ka*x0(t) + kc*x2(t)),
diff(x0(t), t) = (-ka*n*x0(t) - ka*kc*a1*x0(t) + ka*kc*a1*x1(t) - ka*a1*x0(t)^2 + ka*a1*x0(t)*x1(t) - kc*a1*x0(t)*x2(t) + kc*a1*x1(t)*x2(t)) / (ka*kc + ka*x0(t) + kc*x2(t)),
y1(t) = x0(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_mgga_x *)

a1 := 3*Pi/X_FACTOR_C:

f := (rs, x, t, u) -> a1/(2*t - u/4):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

$define gga_c_pbe_params
$include "gga_c_regtpss.mpl"

scan_e0_g := (rs, z, t) -> (1 + 4*A(rs, z, t)*t^2)^(-1/4):
f2 := (rs, z, t) -> mbeta(rs, t)*(1 - scan_e0_g(rs, z, t))/(mgamma*A(rs, z, t)):


text:
RP2_net_complex := proc()
 local T,rt,tau,i;

 T["vertices"] := {seq(i,i=1..9)};
 T["faces"] := {{1,2,3},{1,2,5},{1,3,9},{1,4,5},{1,4,9},{2,3,7},{2,5,6},{2,6,7},{3,7,8},{3,8,9}};
 T["max_simplices"] := T["faces"];

 rt := sqrt(3);
 tau := (1+sqrt(5))/2;

 T["flat_embedding"] := table([
  1 = [ 1/2, 0   ],
  2 = [-1/4, rt/4],
  3 = [-1/4,-rt/4],
  4 = [ 1  , 0   ],
  5 = [ 1/2, rt/2],
  6 = [-1/2, rt/2],
  7 = [-1  , 0   ],
  8 = [-1/2,-rt/2],
  9 = [ 1/2,-rt/2]
 ]):

 T["icosahedron_embedding"] := table([
   1 = [   0,   1, tau],
   2 = [ tau,   0,   1],
   3 = [   1, tau,   0],
   4 = [-tau,   0,   1],
   5 = [   0,  -1, tau],
   6 = [   1,-tau,   0],
   7 = [ tau,   0,  -1],
   8 = [   0,   1,-tau],
   9 = [  -1, tau,   0]
 ]):

 T["flat_plot"] :=
   `plot/simplicial_complex`(T["vertices"])(T["faces"],2,T["flat_embedding"]);

 T["icosahedron_plot"] :=
   `plot/simplicial_complex`(T["vertices"])(T["faces"],2,T["icosahedron_embedding"]);

 return eval(T);
end:

RP2_complex := proc()
 local T0,T,i;

 T0 := RP2_net_complex();
 T["vertices"] := {seq(i,i=1..6)};
 T["faces"] := subs({7=4,8=5,9=6},T0["faces"]);
 T["max_simplices"] := T["faces"];

 return eval(T);
end:


text:
/* The shortest Min-Passy Program : sample12 */
program S;begin end.

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
######################################################################

`eta/stasheff_trees` := proc(A::set)
 if nops(A) = 1 then
  return [`eta/ord`(A),`eta/trees`(A)];
 else
  return FAIL;
 fi;
end;

`gamma/stasheff_trees` := (A::set,B::set) -> (p) -> proc(Y,XX)
 local R,TT,SS,UUUU,b;

 R,TT := op(Y);

 SS := table();
 UUUU := table();

 for b in B do
  SS[b],UUUU[b] := op(XX[b]);
 od;

 return [`gamma/ord`(A,B)(p)(R,SS),
         `gamma/trees`(A,B)(p)(TT,UUUU)];
end;


text:
type $__field_meta__ <struct {
  @fieldname i32,
  @offset i32,
  @declaringclass <* void>,
  @type <* void>}>
var $v1 <* void> public
var $v2 <* void> public
var $vs <[1] <$__field_meta__>> public = [[1= 0x35ad1b, 2= 0xabcd, 3= addrof ptr $v2, 4= addrof ptr $v1]]
var $vvs <[2][1] <$__field_meta__>> public = [[[1= 11, 2= 22, 3= addrof ptr $v1, 4= addrof ptr $v1]],
                                              [[1= 33, 2= 44, 3= addrof ptr $v2, 4= addrof ptr $v2]]]
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
make_icosahedron_complex := proc()
 global icosahedron_complex,dodecahedron_complex;
 local IC,DC,tau,vI,vD,i;

 IC := table():
 DC := table():

 IC["vertices"] := [seq(i,i=1..12)];
 DC["vertices"] := [seq(i,i=1..20)];

 IC["faces"] := [
  [ 1, 2, 9], [ 1, 2,11], [ 1, 5, 6], [ 1, 5, 9], [ 1, 6,11],
  [ 2, 7, 8], [ 2, 7, 9], [ 2, 8,11], [ 3, 4,10], [ 3, 4,12],
  [ 3, 5, 6], [ 3, 5,10], [ 3, 6,12], [ 4, 7, 8], [ 4, 7,10],
  [ 4, 8,12], [ 5, 9,10], [ 6,11,12], [ 7, 9,10], [ 8,11,12]
 ];

 IC["edges"] := {op(map(op,map(f -> [[f[1],f[2]],[f[1],f[3]],[f[2],f[3]]],IC["faces"])))};
 IC["edges"] := sort([op(IC["edges"])]);

 IC["max_simplices"] := IC["faces"];

 tau := (1+sqrt(5))/2;

 IC["embedding_dim"] := 3;

 vI := table([
   1 = [   0,   1, tau],
   2 = [   0,  -1, tau],
   3 = [   0,   1,-tau],
   4 = [   0,  -1,-tau],
   5 = [   1, tau,   0],
   6 = [  -1, tau,   0],
   7 = [   1,-tau,   0],
   8 = [  -1,-tau,   0],
   9 = [ tau,   0,   1],
  10 = [ tau,   0,  -1],
  11 = [-tau,   0,   1],
  12 = [-tau,   0,  -1]
 ]):

 for i from 1 to 12 do vI[i] := simplify(rationalize(vI[i] /~ tau^2)); od;

 IC["embedding"] := eval(vI);
 IC["normalised_embedding"] := table([
  seq(i = evalf(vI[i] /~ sqrt(5-2*sqrt(5))),i = 1..12)
 ]);

 IC["edge_centres"] :=
  map(e -> (vI[e[1]] +~ vI[e[2]]) /~ 2,IC["edges"]);

 IC["face_centres"] :=
  map(f -> (vI[f[1]] +~ vI[f[2]] +~ vI[f[3]]) /~ 3, IC["faces"]);

 `plot/simplicial_complex`(IC);

 icosahedron_complex := eval(IC);
 dodecahedron_complex := eval(DC);
 return eval(IC);
end:

make_icosahedron_complex():

icosphere_complex := proc(n::nonnegint)
 local T,T0,E,v,x,i,j;

 T := eval(icosahedron_complex);

 for i from 1 to n do
  T0 := eval(T);
  T := eval(`triangular_subdivision/simplicial_complex`(T0));
  T0 := eval(T);
  T := eval(`condense/simplicial_complex`(T0));
  E := T["embedding"];
  for v in T["vertices"] do
   x := E[v];
   x := evalf(x /~ sqrt(add(x[j]^2,j=1..3)));
   E[v] := x;
  od:
 od:

 return eval(T):
end:


text:
da1 :=

[
    [   #(1) boon ((1)-(43) are from http://homepages.math.uic.edu/~jan/demo.html)
        [s1**2+g1**2 - 1,
         s2**2+g2**2 - 1,
         C1*g1**3+C2*g2**3 - 6/5,
         C1*s1**3+C2*s2**3 - 6/5,
         C1*g1**2*s1+C2*g2**2*s2 - 7/10,
         C1*g1*s1**2+C2*g2*s2**2 - 7/10
        ],
        [s1,s2,g1,g2,C1,C2]
    ],
    [   #(2) eco6
        [(x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5)*x6 - 1,
        (x2 + x1*x3 + x2*x4 + x3*x5)*x6 - 2,
        (x3 + x1*x4 + x2*x5)*x6 - 3,
        (x4 + x1*x5)*x6 - 4,
        x5*x6 - 5,
        x1 + x2 + x3 + x4 + x5 + 1
        ],
        [x1,x2,x3,x4,x5,x6]
    ],
    [   #(3) mickey
        [ x**2 + 4*y**2 - 4,
        2*y**2 - x
        ],
        [x, y]
    ],
    [   #(4) s9_1
        [ -d2*g - 2*d1*h1,
 9*d2 + 4*b,
 -4*c*h1 - 2*d2*f - 3*d1*g,
 -7*c + 9*a - 8*f,
 -4*d1*f - 5*c*g - 6*h1 - 3*d2,
 -5*d1 - 6*c*f - 7*g + 9*b,
 9*d1 + 6*a - 5*b,
 9*c - 7*a + 8
        ],
        [a,b,c,d1,d2,f,g,h1]
    ],
    [   #(5) conform1
        [-9 - t2**2 - t3**2 - 3*t2**2*t3**2 + 8*t2*t3,
 -9 - t3**2 - t1**2 - 3*t3**2*t1**2 + 8*t3*t1,
 -9 - t1**2 - t2**2 - 3*t1**2*t2**2 + 8*t1*t2
        ],
        [t1, t2, t3]
    ],
    [   #(6) eco7
        [ (x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6)*x7 - 1,
 (x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6)*x7 - 2,
 (x3 + x1*x4 + x2*x5 + x3*x6)*x7 - 3,
 (x4 + x1*x5 + x2*x6)*x7 - 4,
 (x5 + x1*x6)*x7 - 5,
 x6*x7 - 6,
 x1 + x2 + x3 + x4 + x5 + x6 + 1
        ],
        [x1,x2,x3,x4,x5,x6,x7]
    ],
    [   #(7) geneig
        [-10*x1*x6^2+ 2*x2*x6^2-x3*x6^2+x4*x6^2+ 3*x5*x6^2+x1*x6+ 2*x2*x6+x3*x6+ 2*x4*
x6+x5*x6+ 10*x1+ 2*x2-x3+ 2*x4-2*x5,
 2*x1*x6^2-11*x2*x6^2+ 2*x3*x6^2-2*x4*x6^2+x5*x6^2+ 2*x1*x6+x2*x6+ 2*x3*x6+x4*
x6+ 3*x5*x6+ 2*x1+ 9*x2+ 3*x3-x4-2*x5,
-x1*x6^2+ 2*x2*x6^2-12*x3*x6^2-x4*x6^2+x5*x6^2+x1*x6+ 2*x2*x6-2*x4*x6-2*x5*x6-
x1+ 3*x2+ 10*x3+ 2*x4-x5,
x1*x6^2-2*x2*x6^2-x3*x6^2-10*x4*x6^2+ 2*x5*x6^2+ 2*x1*x6+x2*x6-2*x3*x6+ 2*x4*
x6+ 3*x5*x6+ 2*x1-x2+ 2*x3+ 12*x4+x5,
 3*x1*x6^2+x2*x6^2+x3*x6^2+ 2*x4*x6^2-11*x5*x6^2+x1*x6+ 3*x2*x6-2*x3*x6+ 3*x4*
x6+ 3*x5*x6-2*x1-2*x2-x3+x4+ 10*x5,
x1+x2+x3+x4+x5-1
        ],
        [x1, x2, x3, x4, x5, x6]
    ],
    [   #(8) noon3
        [x1*x2^2 + x1*x3^2 - 11/10*x1 + 1,
x2*x1^2 + x2*x3^2 - 11/10*x2 + 1,
x3*x1^2 + x3*x2^2 - 11/10*x3 + 1
        ],
        [x1, x2, x3]
    ],
    [   #(9) redcyc5
        [ 1 + y1 + y2 + y3 + y4,
 y1 + y1*y2 + y2*y3 + y3*y4 + y4,
 y1*y2 + y1*y2*y3 + y2*y3*y4 + y3*y4 + y4*y1,
 y1*y2*y3 + y1*y2*y3*y4 + y2*y3*y4 + y3*y4*y1 + y4*y1*y2,
  z0**5*y1*y2*y3*y4  - 1
        ],
        [y1,y2,y3,y4,z0]
    ],
    [   #(10) sendra
        [ -270*x**4*y**3 - 314*x*y**4 - 689*x*y**3 + 1428,
36*x**7 + 417*x**6*y - 422*x**5*y**2 - 270*x**4*y**3 + 1428*x**3*y**4
- 1475*x**2*y**5 + 510*x*y**6 - 200*x**6 - 174*x**5*y - 966*x**4*y**2
+ 529*x**3*y**3 + 269*x**2*y**4 + 49*x*y**5 - 267*y**6 + 529*x**4*y
+ 1303*x**2*y**3 - 314*x*y**4 + 262*y**5 + 36*x**4 - 788*x**2*y**2
- 689*x*y**3 + 177*y**4
        ],
        [x,y]
    ],
    [   #(11) cpdm5
        [4*x1^3+ 3*x1^2*x2+ 3*x1^2*x3+ 3*x1^2*x4+ 3*x1^2*x5+ 3*x1*x2^2+ 3*x1*x3^2+ 3*
x1*x4^2+ 3*x1*x5^2+x2^3+x3^3+x4^3+x5^3+ 2*x1^2+ 3*x1*x2+ 3*x1*x3+ 3*x1*x4+ 3*
x1*x5-x2^2-x3^2-x4^2-x5^2-6*x1,
x1^3+ 3*x1^2*x2+ 3*x1*x2^2+ 4*x2^3+ 3*x2^2*x3+ 3*x2^2*x4+ 3*x2^2*x5+ 3*x2*
x3^2+ 3*x2*x4^2+ 3*x2*x5^2+x3^3+x4^3+x5^3-x1^2+ 3*x1*x2+ 2*x2^2+ 3*x2*x3+ 3*
x2*x4+ 3*x2*x5-x3^2-x4^2-x5^2-6*x2,
x1^3+ 3*x1^2*x3+ 3*x1*x3^2+x2^3+ 3*x2^2*x3+ 3*x2*x3^2+ 4*x3^3+ 3*x3^2*x4+ 3*
x3^2*x5+ 3*x3*x4^2+ 3*x3*x5^2+x4^3+x5^3-x1^2+ 3*x1*x3-x2^2+ 3*x2*x3+ 2*x3^2+
 3*x3*x4+ 3*x3*x5-x4^2-x5^2-6*x3,
x1^3+ 3*x1^2*x4+ 3*x1*x4^2+x2^3+ 3*x2^2*x4+ 3*x2*x4^2+x3^3+ 3*x3^2*x4+ 3*x3*
x4^2+ 4*x4^3+ 3*x4^2*x5+ 3*x4*x5^2+x5^3-x1^2+ 3*x1*x4-x2^2+ 3*x2*x4-x3^2+ 3*
x3*x4+ 2*x4^2+ 3*x4*x5-x5^2-6*x4,
x1^3+ 3*x1^2*x5+ 3*x1*x5^2+x2^3+ 3*x2^2*x5+ 3*x2*x5^2+x3^3+ 3*x3^2*x5+ 3*x3*
x5^2+x4^3+ 3*x4^2*x5+ 3*x4*x5^2+ 4*x5^3-x1^2+ 3*x1*x5-x2^2+ 3*x2*x5-x3^2+ 3*
x3*x5-x4^2+ 3*x4*x5+ 2*x5^2-6*x5
        ],
        [x1, x2, x3, x4, x5]
    ],
    [   #(12) eco8
        [ (x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 + x6*x7)*x8 - 1,
         (x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6 + x5*x7)*x8 - 2,
        (x3 + x1*x4 + x2*x5 + x3*x6 + x4*x7)*x8 - 3,
        (x4 + x1*x5 + x2*x6 + x3*x7)*x8 - 4,
        (x5 + x1*x6 + x2*x7)*x8 - 5,
        (x6 + x1*x7)*x8 - 6,
        x7*x8 - 7,
        x1 + x2 + x3 + x4 + x5 + x6 + x7 + 1
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8]
    ],
    [   #(13) noon4
        [x1*x2^2 + x1*x3^2 + x1*x4^2 - 11/10*x1 + 1,
x2*x1^2 + x2*x3^2 + x2*x4^2 - 11/10*x2 + 1,
x3*x1^2 + x3*x2^2 + x3*x4^2 - 11/10*x3 + 1,
x4*x1^2 + x4*x2^2 + x4*x3^2 - 11/10*x4 + 1
        ],
        [x1,x2,x3,x4]
    ],
    [   #(14) redcyc6
        [1 + y1 + y2 + y3 + y4 + y5,
 y1 + y1*y2 + y2*y3 + y3*y4 + y4*y5 + y5,
 y1*y2 + y1*y2*y3 + y2*y3*y4 + y3*y4*y5 + y4*y5 + y5*y1,
 y1*y2*y3 + y1*y2*y3*y4 + y2*y3*y4*y5 + y3*y4*y5 + y4*y5*y1 + y5*y1*y2,
 y1*y2*y3*y4 + y1*y2*y3*y4*y5 + y2*y3*y4*y5 + y3*y4*y5*y1
 + y4*y5*y1*y2 + y5*y1*y2*y3,
 z0**6*y1*y2*y3*y4*y5 - 1
        ],
        [y1,y2,y3,y4,y5,z0]
    ],
    [   #(15) solotarev
        [ 3*x**2-2*x-a,
 x**3-x**2-x*a+a-2*b-2,
 3*y**2-2*y-a,
 y**3-y**2-y*a-a+2
        ],
        [x,y,a,b]
    ],
    [   #(16) noon5
        [x1*x2^2 + x1*x3^2 + x1*x4^2 + x1*x5^2 - 11/10*x1 + 1,
x2*x1^2 + x2*x3^2 + x2*x4^2 + x2*x5^2 - 11/10*x2 + 1,
x3*x1^2 + x3*x2^2 + x3*x4^2 + x3*x5^2 - 11/10*x3 + 1,
x4*x1^2 + x4*x2^2 + x4*x3^2 + x4*x5^2 - 11/10*x4 + 1,
x5*x1^2 + x5*x2^2 + x5*x3^2 + x5*x4^2 - 11/10*x5 + 1
        ],
        [x1,x2,x3,x4,x5]
    ],
    [   #(17) redcyc7
        [ 1 + y1 + y2 + y3 + y4 + y5 + y6,
 y1 + y1*y2 + y2*y3 + y3*y4 + y4*y5 + y5*y6 + y6,
 y1*y2 + y1*y2*y3 + y2*y3*y4 + y3*y4*y5 + y4*y5*y6 + y5*y6 + y6*y1,
 y1*y2*y3 + y1*y2*y3*y4 + y2*y3*y4*y5 + y3*y4*y5*y6 + y4*y5*y6
+ y5*y6*y1 + y6*y1*y2,
 y1*y2*y3*y4 + y1*y2*y3*y4*y5 + y2*y3*y4*y5*y6 + y3*y4*y5*y6
+ y4*y5*y6*y1 + y5*y6*y1*y2 + y6*y1*y2*y3,
 y1*y2*y3*y4*y5 + y1*y2*y3*y4*y5*y6 + y2*y3*y4*y5*y6 + y3*y4*y5*y6*y1
+ y4*y5*y6*y1*y2 + y5*y6*y1*y2*y3 + y6*y1*y2*y3*y4,
 z0**7*y1*y2*y3*y4*y5*y6 - 1
        ],
        [y1,y2,y3,y4,y5,y6,z0]
    ],
    [   #(18) sparse5
        [x1**2*x2**2*x3**2*x4**2*x5**2
    + 3*x1**2 + x2**2 + x3**2 + x4**2 + x5**2 + x1*x2*x3*x4*x5 + 5,
 x1**2*x2**2*x3**2*x4**2*x5**2
    + x1**2 + 3*x2**2 + x3**2 + x4**2 + x5**2 + x1*x2*x3*x4*x5 + 5,
 x1**2*x2**2*x3**2*x4**2*x5**2
    + x1**2 + x2**2 + 3*x3**2 + x4**2 + x5**2 + x1*x2*x3*x4*x5 + 5,
 x1**2*x2**2*x3**2*x4**2*x5**2
    + x1**2 + x2**2 + x3**2 + 3*x4**2 + x5**2 + x1*x2*x3*x4*x5 + 5,
 x1**2*x2**2*x3**2*x4**2*x5**2
    + x1**2 + x2**2 + x3**2 + x4**2 + 3*x5**2 + x1*x2*x3*x4*x5 + 5
        ],
        [x1,x2,x3,x4,x5]
    ],
    [   #(19) cyclic6
        [z0 + z1 + z2 + z3 + z4 + z5,
 z0*z1 + z1*z2 + z2*z3 + z3*z4 + z4*z5 + z5*z0,
 z0*z1*z2 + z1*z2*z3 + z2*z3*z4 + z3*z4*z5 + z4*z5*z0 + z5*z0*z1,
 z0*z1*z2*z3 + z1*z2*z3*z4 + z2*z3*z4*z5 + z3*z4*z5*z0 + z4*z5*z0*z1
 + z5*z0*z1*z2,
 z0*z1*z2*z3*z4 + z1*z2*z3*z4*z5 + z2*z3*z4*z5*z0 + z3*z4*z5*z0*z1
 + z4*z5*z0*z1*z2 + z5*z0*z1*z2*z3,
 z0*z1*z2*z3*z4*z5 - 1
        ],
        [z0,z1,z2,z3,z4,z5]
    ],
    [   #(20) extcyc6
        [ z0 + z1 + z2 + z3 + z4 + z5 - 1,
 z0*z1 + z1*z2 + z2*z3 + z3*z4 + z4*z5 + z5*z0 - 1,
 z0*z1*z2 + z1*z2*z3 + z2*z3*z4 + z3*z4*z5 + z4*z5*z0 + z5*z0*z1 - 1,
 z0*z1*z2*z3 + z1*z2*z3*z4 + z2*z3*z4*z5 + z3*z4*z5*z0 + z4*z5*z0*z1
 + z5*z0*z1*z2 - 1,
 z0*z1*z2*z3*z4 + z1*z2*z3*z4*z5 + z2*z3*z4*z5*z0 + z3*z4*z5*z0*z1
 + z4*z5*z0*z1*z2 + z5*z0*z1*z2*z3 - 1 ,
 z0*z1*z2*z3*z4*z5 - 1
        ],
        [z0,z1,z2,z3,z4,z5]
    ],
    [   #(21) cyclic7
        [z0 + z1 + z2 + z3 + z4 + z5 + z6,
 z0*z1 + z1*z2 + z2*z3 + z3*z4 + z4*z5 + z5*z6 + z6*z0,
 z0*z1*z2 + z1*z2*z3 + z2*z3*z4 + z3*z4*z5 + z4*z5*z6 + z5*z6*z0 + z6*z0*z1,
 z0*z1*z2*z3 + z1*z2*z3*z4 + z2*z3*z4*z5 + z3*z4*z5*z6 + z4*z5*z6*z0
+ z5*z6*z0*z1 + z6*z0*z1*z2,
 z0*z1*z2*z3*z4 + z1*z2*z3*z4*z5 + z2*z3*z4*z5*z6 + z3*z4*z5*z6*z0
+ z4*z5*z6*z0*z1 + z5*z6*z0*z1*z2 + z6*z0*z1*z2*z3,
 z0*z1*z2*z3*z4*z5 + z1*z2*z3*z4*z5*z6 + z2*z3*z4*z5*z6*z0 + z3*z4*z5*z6*z0*z1
+ z4*z5*z6*z0*z1*z2 + z5*z6*z0*z1*z2*z3 + z6*z0*z1*z2*z3*z4,
 z0*z1*z2*z3*z4*z5*z6 - 1
        ],
        [z0,z1,z2,z3,z4,z5,z6]
    ],
    [   #(22) puma
        [ x1^2 + x2^2 - 1,
 x3^2 + x4^2 - 1,
 x5^2 + x6^2 - 1,
 x7^2 + x8^2 - 1,
 0004731/1000000*x1*x3 - 03578/10000*x2*x3 - 01238/10000*x1 - 0001637/1000000*x2 - 09338/10000*x4 + x7
- 03571/10000,
 02238/10000*x1*x3 + 07623/10000*x2*x3 + 02638/10000*x1 - 007745/100000*x2 -06734/10000*x4 -06022/10000,
 x6*x8 + 03578/10000*x1 + 0004731/1000000*x2,
 -07623/10000*x1 + 02238/10000*x2 + 03461/10000
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8]
    ],
    [   #(23) redeco5
        [ x1 + x1*x2 + x2*x3 + x3*x4 - u5,
         x2 + x1*x3 + x2*x4 - 2*u5,
         x3 + x1*x4 - 3*u5,
         x4 - 4*u5,
         x1 + x2 + x3 + x4 + 1
        ],
        [x1,x2,x3,x4,u5]
    ],
    [   #(24) trink
        [ 45*y + 35*u - 165*v - 36,
 35*y + 25*z + 40*t - 27*u,
 25*y*u - 165*v**2 + 15*x - 18*z + 30*t,
 15*y*z + 20*t*u - 9*x,
 -11*v**3 + x*y + 2*z*t,
 -11*u*v + 3*v**2 + 99*x
        ],
        [x,y,z,u,v,t]
    ],
    [   #(25) cassou
        [ 15*b**4*c*d1**2 + 6*b**4*c**3 + 21*b**4*c**2*d1 - 144*b**2*c
 - 8*b**2*c**2*d2
 - 28*b**2*c*d1*d2 - 648*b**2*d1 + 36*b**2*d1**2*d2 + 9*b**4*d1**3 - 120,
 30*c**3*b**4*d1 - 32*d1*d2**2*c - 720*d1*b**2*c - 24*c**3*b**2*d2
 - 432*c**2*b**2 + 576*d2*c - 576*d1*d2 + 16*c*b**2*d1**2*d2
 + 16*d1**2*d2**2 + 16*d2**2*c**2 + 9*c**4*b**4 + 5184
 + 39*d1**2*b**4*c**2 + 18*d1**3*b**4*c - 432*d1**2*b**2
 + 24*d1**3*b**2*d2 - 16*c**2*b**2*d1*d2 - 240*c,
 216*d1*b**2*c - 162*d1**2*b**2 - 81*c**2*b**2 + 5184 + 1008*d2*c
 - 1008*d1*d2
 + 15*c**2*b**2*d1*d2 - 15*c**3*b**2*d2 - 80*d1*d2**2*c + 40*d1**2*d2**2
 + 40*d2**2*c**2,
 261 + 4*d1*b**2*c - 3*d1**2*b**2 - 4*c**2*b**2 + 22*d2*c - 22*d1*d2
        ],
        [b,c,d1,d2]
    ],
    [   #(26) katsura5
        [2*x**2+2*y**2+2*z**2+2*t**2+2*u**2+v**2-v,
 x*y+y*z+2*z*t+2*t*u+2*u*v-u,
 2*x*z+2*y*t+2*z*u+u**2+2*t*v-t,
 2*x*t+2*y*u+2*t*u+2*z*v-z,
 t**2+2*x*v+2*y*v+2*z*v-y,
 2*x+2*y+2*z+2*t+2*u+v-1
        ],
        [x,y,z,u,v,t]
    ],
    [   #(27) quadfor2
        [ w1 + w2 - 1,
 w1*x1 + w2*x2,
 w1*x1**2 + w2*x2**2 - 2/3,
 w1*x1**3 + w2*x2**3
        ],
        [w1,w2,x1,x2]
    ],
    [   #(28) redeco6
        [ x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 -   u6,
        x2 + x1*x3 + x2*x4 + x3*x5 - 2*u6,
        x3 + x1*x4 + x2*x5 - 3*u6,
        x4 + x1*x5 - 4*u6,
        x5 - 5*u6,
        x1 + x2 + x3 + x4 + x5 + 1
        ],
        [x1,x2,x3,x4,x5,u6]
    ],
    [   #(29) virasoro
        [8*x1^2+ 8*x1*x2+ 8*x1*x3+ 2*x1*x4+ 2*x1*x5+ 2*x1*x6+ 2*x1*x7-8*x2*x3-2*x4*x7
-2*x5*x6-x1,
 8*x1*x2-8*x1*x3+ 8*x2^2+ 8*x2*x3+ 2*x2*x4+ 2*x2*x5+ 2*x2*x6+ 2*x2*x7-2*x4*x6
-2*x5*x7-x2,
-8*x1*x2+ 8*x1*x3+ 8*x2*x3+ 8*x3^2+ 2*x3*x4+ 2*x3*x5+ 2*x3*x6+ 2*x3*x7-2*x4*
x5-2*x6*x7-x3,
 2*x1*x4-2*x1*x7+ 2*x2*x4-2*x2*x6+ 2*x3*x4-2*x3*x5+ 8*x4^2+ 8*x4*x5+ 2*x4*x6+
 2*x4*x7+ 6*x4*x8-6*x5*x8-x4,
 2*x1*x5-2*x1*x6+ 2*x2*x5-2*x2*x7-2*x3*x4+ 2*x3*x5+ 8*x4*x5-6*x4*x8+ 8*x5^2+
 2*x5*x6+ 2*x5*x7+ 6*x5*x8-x5,
-2*x1*x5+ 2*x1*x6-2*x2*x4+ 2*x2*x6+ 2*x3*x6-2*x3*x7+ 2*x4*x6+ 2*x5*x6+ 8*x6^2
+ 8*x6*x7+ 6*x6*x8-6*x7*x8-x6,
-2*x1*x4+ 2*x1*x7-2*x2*x5+ 2*x2*x7-2*x3*x6+ 2*x3*x7+ 2*x4*x7+ 2*x5*x7+ 8*x6*
x7-6*x6*x8+ 8*x7^2+ 6*x7*x8-x7,
-6*x4*x5+ 6*x4*x8+ 6*x5*x8-6*x6*x7+ 6*x6*x8+ 6*x7*x8+ 8*x8^2-x8
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8]
    ],
    [   #(30) d1
        [ x1**2  + x2**2 - 1,
 x3**2  + x4**2 - 1,
 x5**2  + x6**2 - 1,
 x7**2  + x8**2 - 1,
 x9**2  + x10**2 - 1,
 x11**2 + x12**2 - 1,
 3*x3 + 2*x5 + x7 - 39701/10000,
 3*x1*x4 + 2*x1*x6 + x1*x8 - 17172/10000,
 3*x2*x4 + 2*x2*x6 + x2*x8 - 40616/10000,
 x3*x9 + x5*x9 + x7*x9 - 19791/10000,
 x2*x4*x9 + x2*x6*x9 + x2*x8*x9 + x1*x10 - 19115/10000,
 - x3*x10*x11 - x5*x10*x11 - x7*x10*x11 + x4*x12 + x6*x12 + x8*x12 - 04077/10000
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12]
    ],
    [   #(31) kin1
        [ s1**2 + c1**2 - 1,
 s2**2 + c2**2 - 1,
 s3**2 + c3**2 - 1,
 s4**2 + c4**2 - 1,
 s5**2 + c5**2 - 1,
 s6**2 + c6**2 - 1,
 s2*c5*s6 - s3*c5*s6 - s4*c5*s6 + c2*c6 + c3*c6 + c4*c6 - 04077/10000,
 c1*c2*s5 + c1*c3*s5 + c1*c4*s5 + s1*c5 - 19115/10000,
 s2*s5 + s3*s5 + s4*s5 - 19791/10000,
 c1*c2 + c1*c3 + c1*c4 + c1*c2 + c1*c3 + c1*c2 - 40616/10000,
 s1*c2 + s1*c3 + s1*c4 + s1*c2 + s1*c3 + s1*c2 - 17172/10000,
 s2 + s3 + s4 + s2 + s3 + s2 - 39701/10000
        ],
        [s1,c1,s2,c2,s3,c3,s4,c4,s5,c5,s6,c6]
    ],
    [   #(32) redeco7
        [ x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 -   u7,
        x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6 - 2*u7,
        x3 + x1*x4 + x2*x5 + x3*x6 - 3*u7,
        x4 + x1*x5 + x2*x6 - 4*u7,
        x5 + x1*x6 - 5*u7,
        x6 - 6*u7,
        x1 + x2 + x3 + x4 + x5 + x6 + 1
        ],
        [x1,x2,x3,x4,x5,x6,u7]
    ],
    [   #(33) des18_3
        [6*a33*a10*a20 + 10*a22*a10*a31 + 8*a32*a10*a21 - 162*a10**2*a21
+ 16*a21*a30 + 14*a31*a20 + 48*a10*a30,
15*a33*a10*a21 - 162*a10**2*a22 - 312*a10*a20 + 24*a10*a30 + 27*a31*a21
+ 24*a32*a20 + 18*a22*a10*a32 + 30*a22*a30 + 84*a31*a10,
-240*a10 + 420*a33 - 64*a22 + 112*a32,
180*a33*a10 - 284*a22*a10 - 162*a10**2 + 60*a22*a32 + 50*a32*a10 + 70*a30
+ 55*a33*a21 + 260*a31 - 112*a20,
66*a33*a10 + 336*a32 + 90*a31 + 78*a22*a33 - 1056*a10 - 90*a21,
136*a33 - 136,
4*a22*a10*a30 + 2*a32*a10*a20 + 6*a20*a30 - 162*a10**2*a20 + 3*a31*a21*a10,
28*a22*a10*a33 + 192*a30 + 128*a32*a10 + 36*a31*a20 + 36*a33*a20 - 300*a10*a21 + 40*a32*a21 - 648*a10**2 + 44*a22*a31
        ],
        [a10,a20,a21,a22,a30,a31,a32,a33]
    ],
    [   #(34) kinema
        [z1**2 + z2**2 + z3**2 - 12*z1 - 68,
 z4**2 + z5**2 + z6**2 - 12*z5 - 68,
 z7**2 + z8**2 + z9**2 - 24*z8 - 12*z9 + 100,
 z1*z4 + z2*z5 + z3*z6 - 6*z1 - 6*z5 - 52,
 z1*z7 + z2*z8 + z3*z9 - 6*z1 - 12*z8 - 6*z9 + 64,
 z4*z7 + z5*z8 + z6*z9 - 6*z5 - 12*z8 - 6*z9 + 32,
 2*z2 + 2*z3 - z4 - z5 - 2*z6 - z7 - z9 + 18,
 z1 + z2 + 2*z3 + 2*z4 + 2*z6 - 2*z7 + z8 - z9 - 38,
 z1 + z3 - 2*z4 + z5 - z6 + 2*z7 - 2*z8 + 8
        ],
        [z1,z2,z3,z4,z5,z6,z7,z8,z9]
    ],
    [   #(35) rabmo
        [x1      + x3      + x5        + 2*x7                  - 1,
 x1*x2   + x3*x4   + 2*x5*x6   + 2*x7*x8   + 2*x7*x9   - 2/3,
 x1*x2^2 + x3*x4^2 + 2*x5*x6^2 + 2*x7*x8^2 + 2*x7*x9^2 - 2/5,
 x1*x2^3 + x3*x4^3 + 2*x5*x6^3 + 2*x7*x8^3 + 2*x7*x9^3 - 2/7,
 x1*x2^4 + x3*x4^4 + 2*x5*x6^4 + 2*x7*x8^4 + 2*x7*x9^4 - 2/9,
 x5*x6^2 + 2*x7*x8*x9     - 1/9,
 x5*x6^4 + 2*x7*x8^2*x9^2 - 1/25,
 x5*x6^3 + x7*x8*x9^2 + x7*x8^2*x9 - 1/15,
 x5*x6^4 + x7*x8*x9^3 + x7*x8^3*x9 - 1/21
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8,x9]
    ],
    [   #(36) redeco8
        [ x1 + x1*x2 + x2*x3 + x3*x4 + x4*x5 + x5*x6 + x6*x7 - 1*u8,
        x2 + x1*x3 + x2*x4 + x3*x5 + x4*x6 + x5*x7 - 2*u8,
        x3 + x1*x4 + x2*x5 + x3*x6 + x4*x7 - 3*u8,
        x4 + x1*x5 + x2*x6 + x3*x7 - 4*u8,
        x5 + x1*x6 + x2*x7 - 5*u8,
        x6 + x1*x7 - 6*u8,
        x7 - 7*u8,
        x1 + x2 + x3 + x4 + x5 + x6 + x7 + 1
        ],
        [x1,x2,x3,x4,x5,x6,x7,u8]
    ],
    [   #(37) wright
        [x1^2-x1+x2+x3+x4+x5-10,
x2^2+x1-x2+x3+x4+x5-10,
x3^2+x1+x2-x3+x4+x5-10,
x4^2+x1+x2+x3-x4+x5-10,
x5^2+x1+x2+x3+x4-x5-10
        ],
        [x1,x2,x3,x4,x5]
    ],
    [   #(38) des22_24
        [ 16*a20*a32 + 18*a21*a31 + 20*a22*a30,
-80*a23 + 180*a34 + 855*a35,
 7*a20*a31 + 8*a21*a30,
 210*a35 - 210,
 40*a20*a34 + 44*a21*a33 + 48*a22*a32 + 52*a23*a31 + 280*a30,
 27*a20*a33 + 30*a21*a32 + 33*a22*a31 + 36*a23*a30,
 55*a20*a35 + 60*a21*a34 + 65*a22*a33 + 70*a23*a32 + 80*a30 + 375*a31,
 78*a21*a35 + 84*a22*a34 + 90*a23*a33 - 170*a20 + 102*a31 + 480*a32,
136*a23*a35 - 114*a22 + 152*a33 + 720*a34,
105*a22*a35 + 112*a23*a34 - 144*a21 + 126*a32 + 595*a33
        ],
        [a20,a21,a22,a23,a30,a31,a32,a33,a34,a35]
    ],
    [   #(39) ku10
        [ 5*x1*x2+ 5*x1+ 3*x2+ 55,
 7*x2*x3+ 9*x2+ 9*x3+ 19,
 3*x3*x4+ 6*x3+ 5*x4-4,
 6*x4*x5+ 6*x4+ 7*x5+ 118,
x5*x6+ 3*x5+ 9*x6+ 27,
 6*x6*x7+ 7*x6+x7+ 72,
 9*x7*x8+ 7*x7+x8+ 35,
 4*x8*x9+ 4*x8+ 6*x9+ 16,
 8*x9*x10+ 4*x9+ 3*x10-51,
 3*x1*x10-6*x1+x10+ 5
        ],
        [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]
    ],
    [   #(40) lorentz
        [x1*x2-x1*x3-x4+ 1,
x2*x3-x2*x4-x1+ 1,
-x1*x3+x3*x4-x2+ 1,
x1*x4-x2*x4-x3+ 1
        ],
        [x1,x2,x3,x4]
    ],
    [   #(41) rbpl24
        [62500*x1^2 + 62500*y1^2 + 62500*z1^2 -74529,
625*x2^2 + 625*y2^2 + 625*z2^2 -1250*x2 -2624,
12500*x3^2 + 12500*y3^2 + 12500*z3^2 + 2500*x3 -44975*y3 -10982,
400000*x1*x2 + 400000*y1*y2 + 400000*z1*z2 -400000*x2 + 178837,
1000000*x1*x3 + 1000000*y1*y3 + 1000000*z1*z3 + 100000*x3 -1799000*y3 -805427,
2000000*x2*x3 + 2000000*y2*y3 + 2000000*z2*z3 -2000000*x2 + 200000*x3
-3598000*y3 -1403,
 113800000000000*x3*y2*z1
-113800000000000*x2*y3*z1 -113800000000000*x3*y1*z2 +
113800000000000*x1*y3*z2 + 113800000000000*x2*y1*z3
-113800000000000*x1*y2*z3 -206888400000000*x2*y1 +
206888400000000*x3*y1 + 206888400000000*x1*y2 -206888400000000*x3*y2
-206888400000000*x1*y3 + 206888400000000*x2*y3 -2014260000000*x2*z1 +
2014260000000*x3*z1 -61907200000000*y2*z1 + 61907200000000*y3*z1 +
2014260000000*x1*z2 -2014260000000*x3*z2 + 61907200000000*y1*z2
-61907200000000*y3*z2 -2014260000000*x1*z3 + 2014260000000*x2*z3
-61907200000000*y1*z3 + 61907200000000*y2*z3 -362960716800000*x1 +
38025201600000*x2 + 292548849600000*x3 + 11809567440000*y1 +
1475978220000*y2 -825269402280000*y3 -1212982689600000*z1
-151600474800000*z2 + 825859951200000*z3 -19295432410527,
 -777600000000*x3*y2*z1 + 777600000000*x2*y3*z1 +
777600000000*x3*y1*z2 -777600000000*x1*y3*z2 -777600000000*x2*y1*z3 +
777600000000*x1*y2*z3 -1409011200000*x2*y1 + 1409011200000*x3*y1 +
1409011200000*x1*y2 -1409011200000*x3*y2 -1409011200000*x1*y3 +
1409011200000*x2*y3 -1065312000000*x2*z1 + 1065312000000*x3*z1
-805593600000*y2*z1 + 805593600000*y3*z1 + 1065312000000*x1*z2
-1065312000000*x3*z2 + 805593600000*y1*z2 -805593600000*y3*z2
-1065312000000*x1*z3 + 1065312000000*x2*z3 -805593600000*y1*z3 +
805593600000*y2*z3 + 235685027200*x1 + 398417510400*x2 +
158626915200*x3 -311668424000*y1 -268090368000*y2 + 72704002800*y3 +
412221302400*z1 + 354583756800*z2 + 307085438400*z3 + 282499646407,
3200*x2 + 1271
        ],
        [x1,y1,z1,x2,y2,z2,x3,y3,z3]
    ],
    [   #(42) reimer5
        [-1 + 2*x**2 - 2*y**2 + 2*z**2 - 2*t**2 + 2*u**2,
-1 + 2*x**3 - 2*y**3 + 2*z**3 - 2*t**3 + 2*u**3,
-1 + 2*x**4 - 2*y**4 + 2*z**4 - 2*t**4 + 2*u**4,
-1 + 2*x**5 - 2*y**5 + 2*z**5 - 2*t**5 + 2*u**5,
-1 + 2*x**6 - 2*y**6 + 2*z**6 - 2*t**6 + 2*u**6
        ],
        [x,y,z,t,u]
    ],
    [   #(43) eco5
        [ (x1 + x1*x2 + x2*x3 + x3*x4)*x5 - 1,
        (x2 + x1*x3 + x2*x4)*x5 - 2,
        (x3 + x1*x4)*x5 - 3,
        x4*x5 - 4,
        x1 + x2 + x3 + x4 + 1
        ],
        [x1,x2,x3,x4,x5]
    ]
]:

max_stars_repo_name: lihaokun/StrongSfTriDec

max_issues_repo_name: lihaokun/StrongSfTriDec

max_forks_repo_name: lihaokun/StrongSfTriDec


text:
(*
 Copyright (C) 2017 M.A.L. Marques
               2019 Susi Lehtola

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)
(* prefix:
  gga_x_ev93_params *params;

  assert(p->params != NULL);
  params = (gga_x_ev93_params * )(p->params);
*)

ev93_f0 := s -> (1 + params_a_a1*s^2 + params_a_a2*s^4 + params_a_a3*s^6)/(1 + params_a_b1*s^2 + params_a_b2*s^4 + params_a_b3*s^6):
ev93_f  := x -> ev93_f0(X2S*x):

f := (rs, z, xt, xs0, xs1) -> gga_exchange(ev93_f, rs, z, xs0, xs1):


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_gga_c *)

$include "gga_c_tca.mpl"

msinc := x -> piecewise(x = 0, 1, sin(x)/x):
aa    := Pi*(9*Pi/4)^(1/3):

fD := (rs, z, s) -> 1 - z^4*(1 - msinc(aa*s/rs)^2):

f := (rs, z, xt, xs0, xs1) ->
  f_tcs(rs, z, xt)*fD(rs, z, X2S*2^(1/3)*xt):


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: lda_exc *)
(* prefix:
  lda_c_1d_csc_params *params;

  assert(p->params != NULL);
  params = (lda_c_1d_csc_params * )(p->params);
*)

$define xc_dimensions_1d

(* factor of 2 is the conversion from Ry to Hartree *)
f_aux := (a, rs) -> -(rs + a[5]*rs^2)*log(1+ a[8]*rs + a[9]*rs^a[10])
  / (2*(a[1] + a[2]*rs + a[3]*rs^a[6] + a[4]*rs^a[7])):

f := (rs, z) ->
  f_aux(params_a_para, rs) + (f_aux(params_a_ferro, rs) - f_aux(params_a_para, rs))*z^2:


text:
###############################################################################
# Part 1: Algorithms for computation with subfields of rational functions
###############################################################################

with(PolynomialIdeals):
#------------------------------------------------------------------------------

IdealsEq := proc(A, B)
    # Checks whether polynomial ideals are equal
    return evalb((A subset B) and (B subset A)):
end:

#------------------------------------------------------------------------------

FieldToIdeal := proc(gens)
    # Input: generators of a subfield of the field of rational functions
    # Computes the MQS ideal of the field with the new variables of the form x_aux
    # See: https://mediatum.ub.tum.de/doc/685465/document.pdf Definition 2.16
    #      https://doi.org/10.1016/j.jsc.2005.09.010          Lemma 2.1
    local all_vars, subs_dupl, all_dupl, common_denom, polys, f, gb:
    all_vars := indets(gens):
    subs_dupl := map(v -> v = cat(v, _aux), all_vars):
    all_dupl := map(v -> subs(subs_dupl, v), all_vars):
    common_denom := 1:
    polys := []:
    for f in gens do
        common_denom := lcm(common_denom, denom(f)):
        polys := [op(polys), numer(f) * subs(subs_dupl, denom(f)) - subs(subs_dupl, numer(f)) * denom(f)]:
    end do:
    #gb := Groebner[Basis]([op(polys), common_denom * t - 1], plex(t, op(all_dupl))):

    gb := Groebner[Basis]([op(polys), common_denom * t - 1], tdeg(t, op(all_dupl))):
    gb := Groebner[Walk](gb, tdeg(t, op(all_dupl)), lexdeg([t], [op(all_dupl)])):

    gb := select(p -> not (t in indets(p)), gb):
    return PolynomialIdeal(gb, variables=all_dupl):
end proc:

#------------------------------------------------------------------------------

FieldCoefficientRestriction := proc(J, msq_for_subfield)
    # Input: J - a polynomial ideals over a field of rational functions
    #        msq_for_subfield - the MQS ideal for a subfield E of coefficients (see FieldToIdeal)
    # Computes the radical of the restriction of the ideal to the subfield E
    # in the sense of https://doi.org/10.1016/j.jsc.2005.09.010 (MSQ-paper in what follows)
    #
    # NOTE: unlike the algorithm in the MQS-paper, we compute the radical, not the restriction itself
    # one can obtain the algorithm MQS-paper by replacing PrimeDecomposition with PrimaryDecomposition
    # in the code below
    local poly_vars, coeff_vars, subs_aux, coeff_aux, gens, subs_aux_msq, gens_msq, msq_ideal_aux,
    msq_components, J_ext, components, primes_to_keep, P, elim_P, comp, cleaned_ideal:

    poly_vars := IdealInfo[Variables](J):
    coeff_vars := IdealInfo[Parameters](J) union IdealInfo[Parameters](msq_for_subfield):
    subs_aux := map(v -> v = parse(cat(v, _aux_aux)), coeff_vars):
    coeff_aux := subs(subs_aux, coeff_vars):

    # list F of polynomials in the notation of the MSQ-paper, page 375
    gens := map(p -> subs(subs_aux, p * denom(p)), IdealInfo[Generators](J)):

    # Substitution to avoid names clashing between the aux variables in the msq ideal and the variables in J
    subs_aux_msq := map(v -> v = parse(cat(v, _aux)), IdealInfo[Variables](msq_for_subfield)):
    gens_msq := map(p -> subs(subs_aux_msq, p), IdealInfo[Generators](msq_for_subfield)):
    msq_ideal_aux := PolynomialIdeal(gens_msq, variables=map(s -> rhs(s), subs_aux_msq)):
    msq_components := [PrimeDecomposition(msq_ideal_aux)]:

    J_ext := PolynomialIdeal([op(gens), op(gens_msq)], variables=poly_vars union coeff_aux):
    components := [PrimeDecomposition(J_ext)]:

    # Selecting prime components as in Remark on page 377 in MSQ-paper
    primes_to_keep := []:
    for P in components do
        # Going through the prime decomposition of the MSQ-deal mimics the
        # working over the restricted field (which is what has been done in the MSQ-paper)
        elim_P := EliminationIdeal(P, coeff_aux):
        for comp in msq_components do
            if elim_P subset comp then
                primes_to_keep := [op(primes_to_keep), P]:
            end if:
        end do:
    end do:
    if nops(primes_to_keep) > 0 then
        cleaned_ideal := Intersect(op(primes_to_keep)):
    else
        cleaned_ideal := PolynomialIdeal([0], variables = poly_vars):
    end if:

    # Applying Lemma 2.2 from the MSQ-paper
    return EliminationIdeal(cleaned_ideal, poly_vars):
end proc:

#------------------------------------------------------------------------------

FilterGenerators := proc(ideal)
    # Input: ideal over a rational function field
    # Computes a simplified set of generators of the field of definition
    local gb, gens, p, cf, lc, gsorted, result, big_ideal, cur_ideal, g, new_ideal:
    gb := Groebner[Basis](ideal, tdeg(op(IdealInfo[Variables](ideal)))):
    gens := {}:
    for p in gb do
        cf := sort([coeffs(p, IdealInfo[Variables](ideal))], (a, b) -> length(convert(a, string)) < length(convert(b, string))):
        lc := cf[1]:
        cf := map(c -> c / lc, cf):
        gens := {op(gens), op(cf[2..nops(cf)])}:
    end do:
    gsorted := sort([op(gens)], (a, b) -> length(convert(a, string)) < length(convert(b, string)));
    result := {}:
    big_ideal := FieldToIdeal(gens):
    cur_ideal := FieldToIdeal(result):
    for g in gsorted do
        if big_ideal = cur_ideal then
            return result:
        end if:
        new_ideal := FieldToIdeal({op(result), g}):
        if new_ideal <> cur_ideal then
        	  # a dirty hack to transform -1/a to a
            if convert(g, string)[1] = "-" then
                g := -g:
            end:
            if convert(g, string)[1] = "1" then
                g := 1 / g:
            end:
            result := {op(result), g}:
            cur_ideal := new_ideal:
        end if:
    end do:
    return result:
end proc:

#------------------------------------------------------------------------------

FieldIntersection := proc(gens_left, gens_right)
    # Input: gens_left and gens_right - generators of a subfields E and F of a field of rational functions
    # If terminates, resturns an ideal with field of definition being the intersection of generators of E and F
    # Is guaranteed to terminate if at least one of E and F is algebraically closed in rational functions (see REF)
    # Implementation below is a version of Algorithm 2.38 from https://mediatum.ub.tum.de/doc/685465/document.pdf
    local msq_left, msq_right, poly_vars, coeff_vars, Ii, Ji, gb, result, p, cf, lc;

    msq_left := FieldToIdeal(gens_left):
    msq_right := FieldToIdeal(gens_right):
    poly_vars := IdealInfo[Variables](msq_left) union IdealInfo[Variables](msq_right):
    coeff_vars := IdealInfo[Parameters](msq_left) union IdealInfo[Parameters](msq_right):

    Ii := PolynomialIdeal([1], variables=poly_vars):
    Ji := msq_left:

    while not IdealsEq(Ii, Ji) do
        Ii := FieldCoefficientRestriction(Ji, msq_right):
        Ji := FieldCoefficientRestriction(Ii, msq_left):
    end do:

    return Ji:
end proc:

#------------------------------------------------------------------------------

CompareFields := proc(gens_l, gens_r)
    # Checks whether gens_l and gens_r generate the same subfield of rational functions
    return IdealsEq(FieldToIdeal(gens_l), FieldToIdeal(gens_r)):
end proc:

#------------------------------------------------------------------------------

###############################################################################
# Part 2: Algorithm for computing  identifiable functions
###############################################################################

with(DifferentialAlgebra):

#===============================================================================

ExtractDenominator := proc(model)
    # Input: model a list of rational functions
    # returns the function multiplied by their denominators and
    # an inequality corresponding to the LCM of the denominators
    local common_denom, r;
    common_denom := 1:
    for r in model do
        common_denom := lcm(common_denom, denom(r)):
    end do:
    return [op(map(p -> denom(p) * p, model)), common_denom <> 0]:
end proc:

#===============================================================================
CheckReducibilitySet := proc(polys, charset)
    # Input: polys - list of differential polynomials
    #        charset - a characteristic set
    # Returns true iff all the polynomials are reduced to zero wrt to the charset
    local result, e:
    result := true:
    for e in polys do
        if NormalForm(e, charset) <> 0 then
            result := false:
            break:
        end if:
    end do:
    return result:
end proc:
#===============================================================================

GetIOEquations := proc(model, states, inputs, outputs, params, use_brackets, infolevel, target)
    # Input: model - list of differential polynomials defining the model
    #        states, ios, params - list of names of state variables, input-output variables
    #                              and parameter, respectively
    # Computes a list of input-output equations of the model
    local Relim, Rorig, charsets, chset_orig, general_comps, general, c, e, gen_comp, io_eqs:
    Relim := DifferentialRing(blocks = [[op(states)], [op(outputs)], [op(inputs)]], derivations = [t], arbitrary = params):
    if not use_brackets then
        Relim := DifferentialRing(blocks = [[op(states)], op(outputs), [op(inputs)]], derivations = [t], arbitrary = params):
    end if:
    Rorig := DifferentialRing(blocks = [[op(outputs)], [op(states)], [op(inputs)]], derivations = [t], arbitrary = params):
           # DifferentialRing(blocks = [[op(inputs)], [op(outputs)], [op(states)]], derivations = [t], arbitrary = params):
    chset_orig := RosenfeldGroebner(model, Rorig)[1]:

    if infolevel > 0 then
        LogText(sprintf("    Computing the characteristic set (singsol = none)\n"), target):
    end if:
    charsets := RosenfeldGroebner(model, Relim, singsol = none):
    if CheckReducibilitySet(Equations(charsets[1]), chset_orig) then
        gen_comp := charsets[1]:
    else
        if infolevel > 0 then
            LogText(sprintf("    Did not pick the right component. Using singsol = all\n"), target):
        end if:
        charsets := RosenfeldGroebner(model, Relim):
        if infolevel > 0 then
            LogText(sprint("    Selecting the general component\n"), target):
        end if:
        general_comps := []:
        for c in charsets do
            general := true:
            for e in Equations(c) do
                if NormalForm(e, chset_orig) <> 0 then
                    general := false:
                    break:
                end if:
            end do:
            if general then
                general_comps := [op(general_comps), c]:
            end if:
        end do:
        if nops(general_comps) > 1 then
            LogText(sprintf("More than one component is considered general!", general_comps), target):
        end if:
        gen_comp := general_comps[1]:
    end if:
    io_eqs := Equations(gen_comp, leader < parse(cat(states[-1], "(t)"))):
    return io_eqs:
end proc:

#===============================================================================

# Adapted from
# https://www.mapleprimes.com/questions/36772-Extract-Specific-Coefficients-Of-A-Multivariate
# by Kitonum 15364
coefff:=proc(P, t)
    local L, H, i, k:
    L:=[coeffs(P, indets(P), 'h_aux_for_coef')]: H:=[h_aux_for_coef]: k:=0:
    for i from 1 to nops(H) do
        if H[i]=t then k:=L[i] fi:
    end do:
    return k;
end proc:

#===============================================================================

DecomposePolynomial := proc(p, vars_main, vars_coef, infolevel, target)
    # Input: p - polynomial in two groups of variables: vars_main and vars_coef
    # Computes a decomposition of minimal length of p as a linear combination
    # of products A * B, where A is a polynomial in vars_main and B
    # is a polynomial in vars_coef return two lists: list of A's and list of B's
    local cf, monoms, result_cf, result_monom, i, c, m, j, lc, lm, coeff_in_c:
    cf := [coeffs(collect(p, vars_main, 'distributed'), vars_main, 'monoms')]:
    monoms := [monoms]:
    result_cf := []:
    result_monom := Array([]):
    if infolevel > 0 then
        LogText(sprintf("        Number of monomials: %a\n", nops(cf)), target):
    end:
    for i from 1 to nops(cf) do
        c := cf[i]:
        m := monoms[i]:
        for j from 1 to nops(result_cf) do
            lc, lm := Groebner[LeadingTerm](result_cf[j], plex(op(vars_coef))):
            coeff_in_c := coefff(c, lm):
            c := c - coeff_in_c / lc * result_cf[j]:
            result_monom[j] := result_monom[j] + coeff_in_c / lc * m:
        end do:
        if c <> 0 then
            result_cf := [op(result_cf), c]:
            ArrayTools[Append](result_monom, m):
        end if:
    end do:
    if infolevel > 0 then
        LogText(sprintf("        Reduced to: %a\n", nops(result_cf)), target):
    end:
    return [result_cf, convert(result_monom, list)]:
end proc:

#===============================================================================

ConstructWronskian := proc(io_eq, model, states, inputs, outputs, params, state_space, infolevel, target)
    # Input - the same as for GetIOEquations + one IO-equation + flag subs_param
    # Computes the Wronskian for this equation using the representation
    # given by DecomposePolynomial. Return a pair of the Wronskian
    # reduced modulo the original system and a list of coefficients
    # of the compressed io_eq
    local getNormalForm, diff_to_ord, v, vt, h, v_ord, vd, p, decomp, diff_polys, Rorig, chset_orig,
    M, yus, yus_reduced, yus_list, M_sub, roll, params_sub, ios, ps_sol, ps, i:

    diff_to_ord := {}:
    ios := [op(inputs), op(outputs)]:
    for v in ios do
        vt := parse(cat(v, "(t)")):
        diff_to_ord := {op(diff_to_ord), vt = v}:
        for h from 1 to nops(states) + 1 do
            v_ord := parse(cat(v, "_", h)):
            vd := diff(vt, t$h):
            diff_to_ord := {op(diff_to_ord), vd = v_ord}:
        end do:
    end do:
    p := subs(diff_to_ord, io_eq):
    if infolevel > 0 then
        LogText(sprintf("    Combining monomials to reduce the dimension\n"), target):
    end if:
    decomp := DecomposePolynomial(p, map(e -> rhs(e), diff_to_ord), params, infolevel, target):
    diff_polys := map(p -> subs(map(e -> rhs(e) = lhs(e), diff_to_ord), p), decomp[2]):
    Rorig := DifferentialRing(blocks = [[op(outputs)], [op(states)], [op(inputs)]], derivations = [t], arbitrary = params):
           # DifferentialRing(blocks = [[op(inputs)], [op(outputs)], [op(states)]], derivations = [t], arbitrary = params):
    chset_orig := RosenfeldGroebner(model, Rorig)[1]:

    if infolevel > 0 then
        LogText(sprintf("    Computing the Wronskian\n"), target):
    end if:
    M := VectorCalculus[Wronskian](diff_polys, t):
    yus := indets(M) minus {t}:

    if infolevel > 0 then
        LogText(sprintf("    Reducing the Wronskian\n"), target):
    end if:
    if state_space <> [] then
        ps_sol := GetPSSolution(state_space, nops(yus)):
        yus_reduced := []:
        for v in ios do
            vt := parse(cat(v, "(t)")):
            ps := subs(ps_sol, v(t)):
            for i from 0 to nops(yus) - 1 do
                yus_reduced := [op(yus_reduced), vt = subs({t = 0}, ps)]:
                vt := diff(vt, t):
                ps := diff(ps, t):
            end do:
        end do:
    else
        getNormalForm := proc (p, chset_orig)
            return p=DifferentialAlgebra:-NormalForm(p, chset_orig):
        end proc:
        Rorig := DifferentialRing(blocks = [[op(outputs)], [op(states)], [op(inputs)]], derivations = [t], arbitrary = params):
               # DifferentialRing(blocks = [[op(inputs)], [op(outputs)], [op(states)]], derivations = [t], arbitrary = params):
        chset_orig := RosenfeldGroebner(model, Rorig)[1]:
        yus_list := convert(yus, list):
        yus_reduced:=Threads:-Seq[tasksize=nops(yus_list)](getNormalForm(yus_list[i], chset_orig), i =1..nops(yus_list)):
        # yus_reduced := map(p -> p = NormalForm(p, chset_orig), yus):
    end if:

    # iofun:=map(y->parse(cat(convert(y, string), "(t)")), ios):
    # derivative_orders:=table([seq(p=0, p in iofun)]):
    # for each in yus_list do
    #     indets_ := indets(each) minus {t}: # get indets of derivative to compute order
    #     function_ := op(indets_ intersect {op(iofun)}): # for which function?
    #     order:=numelems(indets_): # compute order
    #     if derivative_orders[function_]<order then
    #         derivative_orders[function_]:=order;
    #     fi:
    # end do;
    # LogExpression(derivative_orders):
    # yus_reduced:=[]:
    # count:=0:
    # # yus_reduced:=Threads:-Seq[tasksize=nops(yus_list)](getNormalForm(yus_list[i], chset_orig), i =1..nops(yus_list)):
    # for y in ios do
    #     y:=parse(cat(convert(y, string), "(t)")):
    #     y_rem:=y:
    #     fun:=y:
    #     while count<derivative_orders[fun] do
    #         yus_reduced:=[op(yus_reduced), y=DifferentialAlgebra:-NormalForm(y_rem, chset_orig)]:
    #         count := count+1:
    #         y_rem:=diff(y_rem, t):
    #         y:=diff(y, t):
    #         LogText("Done Processing %a, order: %a\n", fun, count):
    #     end do:
    #     count:=0:
    # end do:

    M_sub := subs(yus_reduced, M):
    M_sub := subs(map(x -> parse(cat(x, "(t)")) = x, states), M_sub):
    return [M_sub, decomp[1]]:
end proc:

#===============================================================================

SingleExperimentIdentifiableFunctions := proc(model, output_targets, {infolevel := 1})
    # Input: model - ODE model represented as a list of differential polynomials
    # Computes generators of the field of single-identifiable functions
    local result, start, finish, states, inputs, outputs, model_eq, ios, params, io_eqs, si_gens, eq, wrnsk, echelon_form, model_denomfree, target:
    states, inputs, outputs, params, model_eq := op(ParseInput_(model)): # states, ios, params := op(ParseInput_(model)):
    ios := [op(inputs), op(outputs)]:

    # Step 1
    if infolevel > 0 then
        LogText(sprintf("SE Step 1: Computing input-output equations\n"), output_targets[log]):
    end if:
    model_denomfree := ExtractDenominator(model_eq):
    io_eqs := GetIOEquations(model_denomfree, states, inputs, outputs, params, false, infolevel, output_targets[log]):
    if infolevel > 0 then
        LogText(sprintf("SE Total number of io-equations: %a\n", nops(io_eqs)), output_targets[log]):
    end if:

    	    si_gens := {}:
	    for eq in io_eqs do
	        # Step 2
	        if infolevel > 0 then
	            LogText(sprintf("SE Step 2: Constructing the Wronskian\n"), output_targets[log]):
	        end if:
	        wrnsk := ConstructWronskian(eq, model_denomfree, states, inputs, outputs, params, [], infolevel, output_targets[log])[1]:
	        # Step 3
	        if infolevel > 0 then
	            LogText(sprintf("SE Step 3: Computing the reduced row echelon form of the Wronskian\n"), output_targets[log]):
	        end if:
	        echelon_form := LinearAlgebra[ReducedRowEchelonForm](wrnsk):
	        si_gens := {op(si_gens), op(select(x -> not type(x, numeric), convert(echelon_form, list)))}:
	    end do:

    # Step 4
    if infolevel > 0 then
        LogText(sprintf("SE Step 4: restricting to the field of parameters"), output_targets[log]):
    end if:
    	result:=map(simplify, convert(FilterGenerators(FieldIntersection(si_gens, params)), list)):
	DocumentTools:-SetProperty(output_targets[single], expression, result, 'refresh'):
	return result:
end proc:

#===============================================================================
# Adapted from https://github.com/pogudingleb/SIAN
FunctionToVariable_ := proc(f):
    convert(convert(f, string)[1..-4], symbol):
end proc:

ParseInput_ := proc(model)
   local all_symbols, x_functions, io_functions, params, states, ios, diff_polys, lhss, out_functions, in_functions, inputs, outputs:
   diff_polys := map(eq -> lhs(eq) - rhs(eq), model):
   all_symbols := foldl(`union`, op( map(e -> indets(e), diff_polys) )) minus {t}:
   x_functions := map(f -> int(f, t), select( f -> type(int(f, t), function(name)), all_symbols )):
   io_functions := select( f -> not type(int(f, t), function(name)) and type(f, function(name)) and not f in x_functions, all_symbols ):
   lhss := map(eq -> lhs(eq), model):
   out_functions := select(f -> f in lhss, io_functions):
   in_functions := select(f -> not (f in lhss), io_functions):
   params := [op(select(f -> not type(f, function(name)) and not type(int(f, t), function(name)), all_symbols))]:
   states := [op(map(FunctionToVariable_, x_functions))]:
   inputs := [op(map(FunctionToVariable_, in_functions))]:
   outputs := [op(map(FunctionToVariable_, out_functions))]:
   return [states, inputs, outputs, params, diff_polys]:
end proc:

#===============================================================================

MultiExperimentIdentifiableFunctions := proc(model, simplified_generators, no_bound, simplify_bound, max_perms, output_targets): # {simplified_generators := false, no_bound := false})
    # Input: model - ODE model in the state-space form
    # Computes the bound from Theorem REF. Returns a list containing
    # 1. the bound
    # 2. the list of lists of coefficients of the IO-equations (after compression)
    # 3. (optional) simplified set of generators of the field of multi-experiment identifiable functions
    #
    # Optional keyword arguments:
    # 1. simplified_generators - if False, then just the coefficients of the IO-equations are returned.
    #                            if True, they are being simplified (maybe be a substantial simplification)
    #                            but this takes time
    # 2. no_bound - if True, then bound for the number of experiments is not computed (may save a lot of time)

    local states, inputs, outputs, ios, params, model_eqs, io_eqs, io_coeffs, io_coef, wrnsk, s, roll, wrnsk_sub, r, bound,
    generators, i, eq, result, model_denomfree, target, start, finish, infolevel, use_brackets, output_permutations,
    outputs_, io_coeffs_sb, io_eqs_sb, bound_sb, skip_simplify, result_sb, best_output_ordering, returned_generators, idx, old_bound:

    states, inputs, outputs, params, model_eqs := op(ParseInput_(model)):
    use_brackets:=false:
    output_permutations := combinat[permute](outputs):
    output_permutations:= output_permutations[..min(nops(output_permutations),max_perms)];

    infolevel := 1:
    model_denomfree := ExtractDenominator(model_eqs):
    skip_simplify := false:
    old_bound := 10: # this will be replaced by a smaller boud every time
    idx:=0:
    for use_brackets in [true, false] do
        for outputs in output_permutations do
        	LogText(sprintf("Permutation:\t%a\n", outputs),  output_targets[log]):
        	# the first iteration is default permutation
	    	   LogText(sprintf("Use Brackets?\t%a\n", use_brackets), output_targets[log]):
	        if infolevel > 0 then
	            LogText(sprintf("ME Computing input-output equations\n"), output_targets[log]):
	   	   end if:
	   	   io_eqs := GetIOEquations(model_denomfree, states, inputs, outputs, params, use_brackets, infolevel, output_targets[log]):
	   	   if infolevel > 0 then
	            LogText(sprintf("ME Total number of io-equations: %a\n", nops(io_eqs)), output_targets[log]):
	        end if:
	        io_coeffs := []:

	        for eq in io_eqs do
	            io_coef := DecomposePolynomial(eq, indets(eq) minus {op(params)}, params, infolevel, output_targets[log])[1]:
	            io_coeffs := [op(io_coeffs), io_coef]:
	        end do:

	        generators := {}:
    	        for io_coef in io_coeffs do
                 io_coef := sort(io_coef, (a, b) -> length(convert(a, string)) < length(convert(b, string)));
                 for i from 2 to nops(io_coef) do
                     generators := {op(generators), io_coef[i] / io_coef[1]}:
                 end do:
              end do:
              bound := 0:
		    if no_bound = false then
                  for eq in io_eqs do
                      if infolevel > 0 then
                          LogText(sprintf("ME Constructing the Wronskian\n"), output_targets[log]):
                      end if:
                      wrnsk, io_coef := op(ConstructWronskian(eq, model_denomfree, states, inputs, outputs, params, model, infolevel, output_targets[log])):
                      # in the notation of the theorem
                      s := nops(io_coef) - 1:
                      # substitution does not increase the rank, so the resulting bound will be correct
                      roll := rand(1..15):
                      wrnsk_sub := map(v -> v = roll(), indets(wrnsk)):
                      r := LinearAlgebra[Rank](subs(wrnsk_sub, wrnsk)):
                      bound := max(bound, s - r + 1):
                  end do:
               end if:
               if bound<old_bound or idx=0 then # idx makes sure we assign a new bound at the first run
               	old_bound:=bound:
               end if:
               idx := idx +1;
               LogText(sprintf("New Bound:\t%a\n", bound),  output_targets[log]):
               result := [old_bound, generators]:
               if simplified_generators then
                   if infolevel > 0 then
		             LogText(sprintf("ME WARNING: Entering simplification of generators! if this takes too long, try unchecking \"Simplified Generators\"\n"), output_targets[log]):
                   end if:
                   result := [op(result), FilterGenerators(FieldToIdeal(generators))]:
                   DocumentTools:-SetProperty(output_targets[multi], expression, map(simplify, convert(result[3], list)), 'refresh'):
                   returned_generators := result[3]:
               else
                   DocumentTools:-SetProperty(output_targets[multi], expression, map(simplify, convert(result[2], list)), 'refresh'):
                   returned_generators := result[2]:
               end if:
               if old_bound > 0 then
    	              DocumentTools:-SetProperty(output_targets[bound_], expression, old_bound, 'refresh'):
    	              if old_bound=1 then
    	  	             skip_simplify := true:
    	  	             if simplify_bound then
    	  		            LogText(sprintf("Bound on the number of experiments = 1, \"Try to refine bound\" was selected but it will be skipped"), output_targets[log]):
    	  	             end if:
	        	        if simplified_generators then
       		            DocumentTools:-SetProperty(output_targets[single], expression, map(simplify, convert(result[3], list)), 'refresh'):
   		             else
    	   		            DocumentTools:-SetProperty(output_targets[single], expression, map(simplify, convert(result[2], list)), 'refresh'):
   		             end if:
    	              end if:
               else
	              DocumentTools:-SetProperty(output_targets[bound_], expression, "",'refresh'):
               end if:
               if skip_simplify or not simplify_bound then
               	break:
              	else
              		DocumentTools:-SetProperty("being_refined", caption, "being refined", 'refresh'):
               end if:
         end do: # loop over outputs
         if skip_simplify or not simplify_bound then
             # only breaks if we do not simplify. the code runs at most once
             break:
         end if:
    end do:
    DocumentTools:-SetProperty("being_refined", caption, "", 'refresh'):
    return [bound, returned_generators];#table([bd=bound]):
end proc:

GetPSSolution := proc(model, ord)
    # Input: model and integer ord
    # Output: a truncated power series solution with random parameter values
    # and initial conditions of order ord
    local roll, states, inputs, outputs, params, model_eqs, model_sub,
    x_funcs, x_sols, x_eqs, cur_ord, i, rhs_eval, total_sub, y_funcs, y_sols,
    params_subs, input_subs:
    roll := rand(1..15):
    states, inputs, outputs, params, model_eqs := op(ParseInput_(model)):
    params_subs := map(p -> p = roll(), params):
    input_subs := map(u -> parse(cat(u, "(t)")) = add([seq(roll() * t^i, i=0..ord)]), inputs):
    model_sub := subs(params_subs, model):
    model_sub := subs(input_subs, model_sub):
    x_funcs := map(x -> parse(cat(x, "(t)")), states):
    x_sols := map(x -> x = roll(), x_funcs):
    x_eqs := map(x -> subs(model_sub, diff(x, t)), x_funcs):

    for cur_ord from 1 to ord do
        for i from 1 to nops(x_funcs) do
            rhs_eval := series(subs(x_sols, x_eqs[i]), t, ord + 1):
            x_sols[i] := (lhs(x_sols[i]) = (rhs(x_sols[i]) + t^cur_ord * coeff(rhs_eval, t, cur_ord - 1) / cur_ord)):
        end do:
    end do:

    total_sub := [op(x_sols), op(params_subs), op(input_subs)]:
    y_funcs := map(y -> parse(cat(y, "(t)")), outputs):
    y_sols := map(y -> y = subs(total_sub, subs(model_sub, y)), y_funcs):

    return [op(y_sols), op(input_subs), op(params_subs), op(x_sols)]:
end proc:

#===============================================================================

#===============================================================================
IdentifiabilityODE := proc(system_ODEs, params_to_assess, output_targets, {char:=0, p := 0.99, count_solutions:=true, infolevel := 1, method := 2, num_nodes := 6})
#IdentifiabilityODE := proc(system_ODEs, params_to_assess, p, output_targets, count_solutions, char)
	#{p := 0.99, count_solutions:=true, infolevel := 1, method := 2, num_nodes := 6})
#===============================================================================
 local i, j, k, n, m, s, all_params, all_vars, eqs, Q, X, Y, poly, d0, D1,
        sample, all_subs,alpha, beta, Et, x_theta_vars, prolongation_possible,
        eqs_i, JacX, vars, vars_to_add, ord_var, var_index, deg_variety, D2,
        y_hat, u_hat, theta_hat, Et_hat, Q_hat, theta_l, theta_g, gb, v, X_eq, Y_eq,
        poly_d, separant, leader, x_functions, y_functions, u_functions,
        all_symbols_rhs, mu, x_vars, y_vars, u_vars, theta, subst_first_order,
        subst_zero_order, x_eqs, y_eqs, param, other_params, to_add, at_node,
        prime, max_rank, R, tr, e, p_local, xy_ders, polys_to_process, new_to_process, start, finish,
        Et_x_vars, out_sian, var, G, P, solutions_table, check:

  #----------------------------------------------
  # 0. Extract inputs, outputs, states, and parameters from the system
  #----------------------------------------------
  if SearchText(".", convert(system_ODEs, string)) <> 0 then
    LogText(sprintf("WARNING: It looks like your system involves floating-point numbers. This may result into a non-meaninful result, please convert them to rationals (e.g., 0.2 -> 1/5)"), "LogAreaSIAN"):
    WARNING("WARNING: It looks like your system involves floating-point numbers. This may result into a non-meaninful result, please convert them to rationals (e.g., 0.2 -> 1/5)"):
  end if:

  if not (indets(system_ODEs, name) subset indets(system_ODEs)) then
    LogText(sprintf(cat("ERROR: you are using reserved maple symbols: ", convert(indets(system_ODEs, name) minus indets(system_ODEs), string))), output_targets[log]):
    DocumentTools:-SetProperty("GlobalParams1", expression, Error(BadName)):
    DocumentTools:-SetProperty("LocalParams1", expression,  Error(BadName)):
    DocumentTools:-SetProperty("NoIDParams1", expression,  Error(BadName)):
    error (cat("ERROR: you are using reserved maple symbols: ", convert(indets(system_ODEs, name) minus indets(system_ODEs), string))):
    return;
  end if:

  randomize():
#  infolevel:=1:
#  num_nodes:=1:
  if infolevel > 0 then
    PrintHeader("0. Extracting states, inputs, outputs, and parameters from the system", output_targets[log]):
  end if:
  x_functions := map(f -> int(f, t), select( f -> type(int(f, t), function(name)), map(lhs, system_ODEs) )):
  y_functions := select( f -> not type(int(f, t), function(name)), map(lhs, system_ODEs) ):
  all_symbols_rhs := foldl(`union`, op( map(e -> indets(rhs(e)), system_ODEs) )) minus {t}:
  xy_ders := {op(x_functions), op(y_functions), op(select(f -> (f in all_symbols_rhs), map(lhs, system_ODEs)))}:
  u_functions := select( f -> type(f, function), convert(all_symbols_rhs minus xy_ders, list)):
  mu := convert(all_symbols_rhs minus {op(xy_ders), op(u_functions)}, list):

  x_vars := map(FunctionToVariable, x_functions):
  y_vars := map(FunctionToVariable, y_functions):
  u_vars := map(FunctionToVariable, u_functions):
  theta := map(ParamToInner, params_to_assess):
  subst_first_order := {seq(diff(x_functions[i], t) = MakeDerivative(x_vars[i], 1), i = 1 .. nops(x_vars))}:
  subst_zero_order := {
    seq(x_functions[i] = MakeDerivative(x_vars[i], 0), i = 1 .. nops(x_vars)),
    seq(y_functions[i] = MakeDerivative(y_vars[i], 0), i = 1 .. nops(y_vars)),
    seq(u_functions[i] = MakeDerivative(u_vars[i], 0), i = 1 .. nops(u_vars))
  }:
  x_eqs := subs(subst_zero_order, subs(subst_first_order, select(e -> type(int(lhs(e), t), function(name)), system_ODEs))):
  y_eqs := subs(subst_zero_order, select(e -> not type(int(lhs(e), t), function(name)), system_ODEs)):

  # taking into account that fact that Groebner[Basis] is Monte Carlo with probability of error
  # at most 10^(-18) (for Maple 2017)
  p_local := p + nops(params_to_assess) * 10^(-18):
  if p_local >= 1 then
    LogText("The probability of success cannot exceed 1 - #params_to_assess 10^{-18}. We reset it to 0.99", output_targets[log]);
    p_local := 0.99:
  end if:

  if infolevel > 0 then
    LogText(sprintf("\n===== Input info =====\n"), output_targets[log]):
    LogText(sprintf("%s %a\n", `State variables:         `, x_functions), output_targets[log]):
    LogText(sprintf("%s %a\n", `Output variables:        `, y_functions), output_targets[log]):
    LogText(sprintf("%s %a\n", `Input variables:         `, u_functions), output_targets[log]):
    LogText(sprintf("%s %a\n", `Parameters in equations: `, mu), output_targets[log]):
    LogText(sprintf("===================\n\n"), output_targets[log]):
  end if:

  #----------------------------------------------
  # 1. Construct the maximal system.
  #----------------------------------------------

  if infolevel > 0 then
    PrintHeader("1. Constructing the maximal polynomial system",output_targets[log]):
  end if:

  # (a) ---------------
  n := nops(x_vars):
  m := nops(y_vars):
  s := nops(mu) + n:
  all_params := [op(mu), op(map(x -> MakeDerivative(x, 0), x_vars ))]:
  all_vars := [ op(x_vars), op(y_vars), op(u_vars) ]:
  eqs := [op(x_eqs), op(y_eqs)]:
  Q := foldl( (f, g) -> lcm(f, g), op( map(f -> denom(rhs(f)), eqs) )):

  # (b,c) ---------------
  X := []:
  X_eq := []:
  for i from 1 to n do
    X := [op(X), []]:
    poly := numer(lhs(x_eqs[i]) - rhs(x_eqs[i])):
    for j from 0 to s + 1 do
      poly_d := Differentiate(poly, all_vars, j):
      leader := MakeDerivative(x_vars[i], j + 1):
      separant := diff(poly_d, leader):
      X[i] := [op(X[i]), poly_d]:
      X_eq := [op(X_eq), leader = -(poly_d - separant * leader) / separant]:
    end do:
  end do:
  # LogExpression(X[1]);

  # (d,e) ---------------
  Y := []:
  Y_eq := []:
  for i from 1 to m do
    Y := [op(Y), []]:
    poly := numer(lhs(y_eqs[i]) - rhs(y_eqs[i])):
    for j from 0 to s + 1 do
      poly_d := Differentiate(poly, all_vars, j):
      leader := MakeDerivative(y_vars[i], j):
      separant := diff(poly_d, leader):
      Y[i] := [op(Y[i]), poly_d]:
      Y_eq := [op(Y_eq), leader = -(poly_d - separant * leader) / separant]:
    end do:
  end do:

  #----------------------------------------------
  # 2. Truncate.
  #----------------------------------------------

  if infolevel > 0 then
    PrintHeader("2. Truncating the polynomial system based on the Jacobian condition", output_targets[log]):
  end if:

  # (a) ---------------
  d0 := max(op( map(f -> degree( simplify(Q * rhs(f)) ), eqs) ), degree(Q)):

  # (b) ---------------
  # extra factor nops(theta) + 1 compared to the formula in the paper is to
  # provide probability gaurantee to the local identifiability test
  D1 := floor( (nops(theta) + 1) * 2 * d0 * s * (n + 1) * (1 + 2 * d0 * s) / (1 - p_local) ):
  # prime := nextprime(D1):
  if infolevel > 1 then
    LogText(sprintf("%s %a\n", `Bound D_1 for testing the rank of the Jacobian probabilistically: `, D1), output_targets[log]);
  end if:

  # (c, d) ---------------
  sample := SamplePoint(D1, x_vars, y_vars, u_vars, mu, X_eq, Y_eq, Q):
  all_subs := sample[4]:
  u_hat := sample[2]:
  y_hat := sample[1]:

  # (e) ------------------
  alpha := [seq(1, i = 1..n)]:
  beta := [seq(0, i = 1..m)]:
  Et := [];
  # TODO: improve for arbitrary derivatives
  x_theta_vars := all_params:
  prolongation_possible := [seq(1, i = 1..m)]:

  # (f) ------------------
  while add(prolongation_possible) > 0 do
    for i from 1 to m do
      if prolongation_possible[i] = 1 then
        eqs_i := [op(Et), Y[i][beta[i] + 1]]:
        JacX := VectorCalculus[Jacobian](subs({op(u_hat), op(y_hat)}, eqs_i), x_theta_vars = subs(all_subs, x_theta_vars)):
        if LinearAlgebra[Rank](JacX) = nops(eqs_i) then
          Et := [op(Et), Y[i][beta[i] + 1]]:
          beta[i] := beta[i] + 1:
          # adding necessary X-equations
          polys_to_process := [op(Et), seq(Y[k][beta[k] + 1], k=1..m)]:
          while nops(polys_to_process) <> 0 do
            new_to_process := []:
            vars := {};
            for poly in polys_to_process do
              vars := vars union { op(GetVars(poly, x_vars)) }:
            end do:
            vars_to_add := { op(remove(v -> evalb(v in x_theta_vars), vars)) };
            for v in vars_to_add do
              x_theta_vars := [op(x_theta_vars), v];
              ord_var := GetOrderVar(v);
              var_index := ListTools[Search](ord_var[1], x_vars):
              poly := X[ var_index ][ ord_var[2] ]:
              Et := [op(Et), poly]:
              new_to_process := [op(new_to_process), poly]:
              alpha[ var_index ] := max(alpha[ var_index ], ord_var[2] + 1):
            end do:
            polys_to_process := new_to_process:
          end do:
        else
          prolongation_possible[i] := 0;
        end if:
      end if:
    end do:
  end do:
  # is used for assessing local identifiabilty
  max_rank := nops(Et):

  # (g) --------------
  for i from 1 to m do
    for j from beta[i] + 1 to nops(Y[i]) do
      to_add := true:
      for v in GetVars(Y[i][j], x_vars) do
        if not (v in x_theta_vars) then
          to_add := false:
        end if:
      end do:
      if to_add = true then
        beta[i] := beta[i] + 1:
        Et := [op(Et), Y[i][j]]:
      end if:
    end do:
  end do:

  if infolevel > 1 then
    LogText(sprintf("%s %a\n", `Orders of prolongations of the outputs (beta) = `, beta), output_targets[log]):
    LogText(sprintf("%s %a\n", `Orders of prolongations of the state variables (alpha) = `, alpha), output_targets[log]):
  end if:
  ##############################

  if infolevel > 0 then
    PrintHeader("3. Assessing local identifiability", output_targets[log]):
  end if:

  theta_l := []:
  for param in theta do
    other_params := subs(param = NULL, x_theta_vars):
    JacX := VectorCalculus[Jacobian](
        subs( { op(u_hat), param = subs(all_subs, param), op(y_hat) }, Et),
        other_params = subs(all_subs, other_params)
    ):
    if LinearAlgebra[Rank](JacX) <> max_rank then
      theta_l := [op(theta_l), param]:
    end if:
  end do:
  DocumentTools:-SetProperty("LocalLabel" , caption, "Locally Identifiable Paramters");
  if infolevel > 1 then
    LogText(sprintf("%s %a\n", `Locally identifiable paramters: `, map(x -> ParamToOuter(x, all_vars), theta_l)), output_targets[log]);
    LogText(sprintf("%s %a\n", `Nonidentifiable parameter: `, map(x -> ParamToOuter(x, all_vars), [op({op(theta)} minus {op(theta_l)})])), output_targets[log]);
  end if:

  DocumentTools:-SetProperty(output_targets[localparams], expression, map(x -> ParamToOuter(x, all_vars), theta_l), 'refresh'): # local
  DocumentTools:-SetProperty(output_targets[noidparams], expression, map(x -> ParamToOuter(x, all_vars), [op({op(theta)} minus {op(theta_l)})]), 'refresh'): # no id

  #----------------------------------------------
  # 3. Randomize.
  #----------------------------------------------

  if infolevel > 0 then
    PrintHeader("4. Randomizing the truncated system", output_targets[log]):
  end if:

  # (a) ------------
  deg_variety := foldl(`*`, op( map(e -> degree(e), Et) )):
  D2 := floor( 6 * nops(theta_l) * deg_variety * (1 + 2 * d0 * max(op(beta))) / (1 - p_local) ):
  if infolevel > 1 then
    LogText(sprintf("%s %a\n", `Bound D_2 for assessing global identifiability: `, D2), output_targets[log]):
  end if:

  # (b, c) ---------
  sample := SamplePoint(D2, x_vars, y_vars, u_vars, mu, X_eq, Y_eq, Q):
  y_hat := sample[1]:
  u_hat := sample[2]:
  theta_hat := sample[3]:
  if infolevel > 1 then
    LogText(sprintf("%s %a\n", `Random sample for the outputs and inputs is generated from `, theta_hat), output_targets[log]):
  end if:

  # (d) ------------
  Et_hat := map(e -> subs([op(y_hat), op(u_hat)], e), Et):
  Et_x_vars := {}:
  for poly in Et_hat do
    Et_x_vars := Et_x_vars union { op(GetVars(poly, x_vars)) }:
  end do:
  if infolevel > 1 then
    LogText(sprintf("%s %a %s %a %s\n", "The polynomial system contains", nops(Et_hat), "equations in ", nops(Et_x_vars) + nops(mu), " variables"), output_targets[log]);
  end if:
  Q_hat := subs(u_hat, Q):

  vars := [
    op(sort([op(Et_x_vars)], (a, b) -> CompareDiffVar(a, b, x_vars))),
    z_aux, w_aux,
    op(sort(mu))
  ]:
  if infolevel > 1 then
    LogText(sprintf("Variable ordering to be used for Groebner basis computation %a\n", vars), output_targets[log]);
  end if:

  #----------------------------------------------
  # 4. Determine.
  #----------------------------------------------

  if infolevel > 0 then
    PrintHeader("5. Assessing global identifiability", output_targets[log]):
  end if:
  theta_g := []:
  if method = 1 then
    at_node := proc(var, args_node)
      local gb_loc, fname;
      gb_loc := Groebner[Basis](op(args_node)):
      gb_loc;
    end proc:

    if nops(theta_l) > 1 and num_nodes > 1 then
      Grid[Setup]("local", numnodes = num_nodes):
      Grid[Set](at_node):
      gb := Grid[Seq](
        at_node(theta_l[i], [
          [op(Et_hat), z_aux * Q_hat - 1, (theta_l[i] - subs(theta_hat, theta_l[i])) * w_aux - 1],
          tdeg(vars)
        ]),
        i = 1..nops(theta_l)
      ):
    else
      gb := []:
      for i from 1 to nops(theta_l) do
        gb := [
          op(gb),
          at_node(
            theta_l[i],
            [[op(Et_hat), z_aux * Q_hat - 1, (theta_l[i] - subs(theta_hat, theta_l[i])) * w_aux - 1], tdeg(op(vars))]
           )
        ]:
      end do:
    end if:
    for i from 1 to nops(theta_l) do
      if gb[i] = [1] then
        theta_g := [op(theta_g), theta_l[i]]:
      else
        if infolevel > 1 then
          LogText(sprintf("%s %a %s %a\n", `Groebner basis corresponding to the parameter `, theta_l[i], ` is `, gb[i]), output_targets[log]):
        end if:
      end if:
    end do:
  elif method = 2 then
     LogText(sprintf("%s %a\n", `Computing Groebner Basis with characteristic`, char), output_targets[log]):
    	gb := Groebner[Basis]([op(Et_hat), z_aux * Q_hat - 1], tdeg(op(vars)), characteristic=char);
     for i from 1 to nops(theta_l) do
       if char>0 then
       	check := subs(theta_hat, theta_l[i]) mod char:
       else
       	check := subs(theta_hat, theta_l[i]):
       end if:
       if Groebner[NormalForm](theta_l[i], gb, tdeg(op(vars)), characteristic=char) = check then
         theta_g := [ op(theta_g), theta_l[i] ]:
       end if:
     end do:

    if count_solutions then
     LogParameters("",output_targets[parameters]):
    	for var in theta_g do
 		LogText(sprintf("%s %a, %s %a\n",`The number of solutions for`, var, `is`, 1), output_targets[log]):
		LogParameters(sprintf("%s %a, %s %a\n",`Parameter`, ParamToOuter(var, all_vars), `number of solutions:`, 1), output_targets[parameters]):
		solutions_table[var]:=1:
    	end do:

    	for var in select(p -> not p in theta_g, theta_l) do
		G := Groebner[Walk](gb, tdeg(op(vars)), lexdeg([op({op(vars)} minus {var})], [var]), characteristic=char):
		P := select(x->evalb(indets(x)={var}), G):
		solutions_table[var]:=degree(P[1], [op(indets(P))]):
		LogText(sprintf("%s %a, %s %a\n",`The number of solutions for`, var, `is`, degree(P[1], [op(indets(P))])), output_targets[log]):
		LogParameters(sprintf("%s %a, %s %a\n",`Parameter`, ParamToOuter(var, all_vars), `number of solutions:`, degree(P[1], [op(indets(P))])),  output_targets[parameters])
    	end do:
     fi:
    	DocumentTools:-SetProperty(output_targets[globalparams], value,  [op(map(x -> ParamToOuter(x, all_vars), theta_g))], 'refresh'): # global
    	DocumentTools:-SetProperty(output_targets[localparams], value,  [op(map(x -> ParamToOuter(x, all_vars), select(p -> not p in theta_g, theta_l)))], 'refresh'): # global
  else
    LogText(sprintf(`No such method`), output_targets[log]):
    LogText(sprintf("%q\n", "Provided method: %d, allowed methods: 1, 2", method), output_targets[log]):
  end if:

  if infolevel > 0 then
    LogText(sprintf("\n=== Summary ===\n"), output_targets[log]):
    LogText(sprintf("%s %a\n", `Globally identifiable parameters:                 `, map(x -> ParamToOuter(x, all_vars), theta_g)), output_targets[log]):
    LogText(sprintf("%s %a\n", `Locally but not globally identifiable parameters: `, map(x -> ParamToOuter(x, all_vars), select(p -> not p in theta_g, theta_l))), output_targets[log]):
    LogText(sprintf("%s %a\n", `Not identifiable parameters:                      `, map(x -> ParamToOuter(x, all_vars), select(p -> not p in theta_l, theta))), output_targets[log]):
    LogText(sprintf("===============\n\n"), output_targets[log]):
  end if:
  DocumentTools:-SetProperty("LocalLabel" , caption, "Locally but not Globally Identifiable Paramters");
  out_sian := table([
    globally = [op(map(x -> ParamToOuter(x, all_vars), theta_g))],
    locally_not_globally = {op(map(x -> ParamToOuter(x, all_vars), select(p -> not p in theta_g, theta_l)))},
    non_identifiable = {op(map(x -> ParamToOuter(x, all_vars), select(p -> not p in theta_l, theta)))},
    parameters = mu,
    basis=gb,
    varordering=tdeg(op(vars)),
    allvars=vars,
    for_outer=all_vars
  ]):
  if count_solutions then
      PrintHeader("WARNING: The result of solution counting is guaranteed with high probability, however it NOT the same probability 'p' as provided in the input.", output_targets[log]):
  end if:
  out_sian[num_solutions] := solutions_table:
  return out_sian;
end proc:

#===============================================================================
PrintHeader := proc(text, output_target):
#===============================================================================
  LogText(sprintf("\n=======================================================\n"), output_target):
  LogText(sprintf(text), output_target):
  LogText(sprintf("\n=======================================================\n"), output_target):
end proc:

#===============================================================================
GetParameters := proc(system_ODEs, initial_conditions) local initial_values, all_symbols_rhs, mu:
#===============================================================================
  initial_values := map(f -> subs({t = 0}, int(f, t)), select( f -> type(int(f, t), function(name)), map(lhs, system_ODEs) )):
  all_symbols_rhs := foldl(`union`, op( map(e -> indets(rhs(e)), system_ODEs) )) minus {t}:
  mu := select(s -> not type(s, function), all_symbols_rhs):
  if initial_conditions then
    return [op(mu), op(initial_values)]:
  else
    return [op(mu)]:
  end if:
end proc:

#===============================================================================
FunctionToVariable := proc(f):
#===============================================================================
  convert(cat(convert(f, string)[1..-4], "_"), symbol):
end proc:

#===============================================================================
ParamToInner := proc(p) local s:
#===============================================================================
  s := convert(p, string):
  if length(s) > 3 and s[-3..-1] = "(0)" then
    MakeDerivative(FunctionToVariable(p), 0):
  else
    p:
  end if:
end proc:

#===============================================================================
ParamToOuter := proc(p, varnames) local s:
#===============================================================================
  s := convert(p, string):
  if length(s) > 2 and s[-2..-1] = "_0" and parse(s[1..-2] )in varnames then
    parse(cat(s[1..-3], "(0)")):
  else
    p:
  end if:
end proc:

#===============================================================================
MakeDerivative := proc(var_name, der_order):
#===============================================================================
  cat(var_name, der_order):
end proc:

#===============================================================================
DifferentiateOnce := proc(diff_poly, var_list)
#===============================================================================
  local result, aux, v, h, diff_v:
  result := 0:
  for diff_v in indets(diff_poly) do
    aux := GetOrderVar(diff_v):
    # seems that Maple does not have unpacking
    v := aux[1]:
    h := aux[2]:
    if v in var_list then
      result := result + diff(diff_poly, MakeDerivative(v, h)) * MakeDerivative(v, h + 1):
    end if:
  end do:
  simplify(result):
end proc:

#===============================================================================
Differentiate := proc(diff_poly, var_list, ord := 1)
#===============================================================================
  local result, i;
  result := diff_poly:
  for i from 1 to ord do
    result := DifferentiateOnce(result, var_list):
  end do:
  result:
end proc:

#===============================================================================
GetVars := proc(diff_poly, var_list)
#===============================================================================
  local result;
  result := select(v -> evalb(GetOrderVar(v)[1] in var_list), indets(diff_poly)):
  return result:
end proc:

#===============================================================================
GetOrderVar := proc(diff_var)
#===============================================================================
  local s, v, h;
  if not StringTools[RegMatch]("^(.*_)([0-9]+)$", diff_var, s, v, h) then
    return ["", ""]:
  end if:
  return [parse(v), parse(h)]:
end proc:

#===============================================================================
SamplePoint := proc(bound, x_vars, y_vars, u_vars, mu, X_eq, Y_eq, Q)
#===============================================================================
  local n, m, s, all_params, all_vars, theta_hat, u_variables,
        u_hat, x_hat, y_hat, eq, eq_num, eq_denom,
        v, poly, i, j, all_subs, roll, to_compute;
  n := nops(x_vars):
  m := nops(y_vars):
  s := nops(mu) + n:
  all_params := [op(mu), op(map(x -> MakeDerivative(x, 0), x_vars ))]:
  all_vars := [ op(x_vars), op(y_vars), op(u_vars) ]:

  roll := rand(0 .. bound):
  while true do
    theta_hat := map(p -> p = roll(), all_params):
    u_variables := [];
    for i from 1 to nops(u_vars) do
      u_variables := [ op(u_variables), seq(MakeDerivative(u_vars[i], j), j = 0..s + 1) ]:
    end do:
    u_hat := map(p -> p = roll(), u_variables) :

    all_subs := [op(theta_hat), op(u_hat)]:
    if subs(all_subs, Q) = 0 then
      next
    end if:
    to_compute := [op(X_eq), op(Y_eq)]:
    while nops(to_compute) <> 0 do
      to_compute := map(e -> lhs(e) = subs(all_subs, rhs(e)), to_compute);
      all_subs := [ op(all_subs), op(select(e -> type(rhs(e), numeric), to_compute)) ]:
      to_compute := remove(e -> type(rhs(e), numeric), to_compute):
    end do:
    y_hat := map(e -> lhs(e) = subs(all_subs, rhs(e)), Y_eq):
    x_hat := map(e -> lhs(e) = subs(all_subs, rhs(e)), X_eq):
    break:
  end do:

  return [y_hat, u_hat, theta_hat, all_subs];
end proc:

#===============================================================================
GenerateReplica := proc(equations, r)
#===============================================================================
  # generates a system of equations corresponding to r independent trajectories of
  # the original system. Time-dependent variabes are replicated, parameters are not
  local all_functions, zero_order, first_order, funcs, without_t, result, i, subst:
  all_functions := select(f -> type(f, function), foldl(`union`, op( map(indets, equations) ))):
  zero_order := select(f -> not type(int(f, t), function(name)), all_functions):
  first_order := map(f -> int(f, t), select(f -> type(int(f, t), function(name)), all_functions)):
  funcs := {op(zero_order), op(first_order)}:
  without_t := map(f -> convert(convert(f, string)[1..-4], symbol), funcs):
  result := []:
  for i from 1 to r do
    subst := map(f -> f = convert(cat(convert(f, string), "_r", i), symbol), without_t):
    result := [op(result), op(map(e -> subs(subst, e), equations))]:
  end do:
  return result:
end proc:

#===============================================================================
CompareDiffVar := proc(dvl, dvr, var_list)
#===============================================================================
  local vl, vr, hl, hr;
  vl, hl := op(GetOrderVar(dvl, var_list)):
  vr, hr := op(GetOrderVar(dvr, var_list)):
  if evalb(hl <> hr) then
    return evalb(hl > hr):
  end if:
  if evalb(length(vl) <> length(vr)) then
    return evalb(length(vl) > length(vr)):
  end if:
  return StringTools[Compare](vr, vl):
end proc:

LogText:= proc(text, target)
	# must pass text as sprintf(text)!!!
	UpdateLog(text, target);
end proc:

LogExpression:= proc(text, target)
	# must pass text as sprintf("%q\n", text)!!!
	UpdateLog(text, target);
end proc:

#LogTextSIAN := ()->UpdateLog(sprintf(_passed), "LogAreaSIAN");
#LogExpressionSIAN := ()->UpdateLog(sprintf("%q\n", _passed), "SIAN");

#LogTextME := ()->UpdateLog(sprintf(_passed), "LogAreaME");
#LogExpressionME := ()->UpdateLog(sprintf("%q\n", _passed), "LogAreaME");

#LogTextSE := ()->UpdateLog(sprintf(_passed), "LogAreaSE");
#LogExpressionSE := ()->UpdateLog(sprintf("%q\n", _passed), "LogAreaSE");

UpdateLog := proc(s, target)
local logsofar;

	logsofar := DocumentTools:-GetProperty(target, value);
	if logsofar <> "" then
		logsofar := logsofar, "\n";
	end if;

	DocumentTools:-SetProperty(target, value, cat(logsofar,s), 'refresh');

end proc:

examples := table([
  	"Biohydrogenation" = [ "Taken from R. Munoz-Tamayo, L. Puillet, J.B. Daniel, D. Sauvant, O. Martin, M. Taghipoor, P. Blavy\n Review: To be or not to be an identifiable model. Is this a relevant question in animal science modelling?\ndoi.org/10.1017/S1751731117002774\nSystem (3) in Supplementary Material 2, initial conditions are assumed to be unknown",
  	[
  "dx4/dt = - k5 * x4 / (k6 + x4);\n",
  "dx5/dt = k5 * x4 / (k6 + x4) - k7 * x5/(k8 + x5 + x6);\n",
  "dx6/dt = k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10;\n",
  "dx7/dt = k9 * x6 * (k10 - x6) / k10;\n",
  "y1 = x4;\n",
  "y2 = x5"]],

  	"Chemical Reaction Network" = ["Taken from  Conradi, C., Shiu, A., Dynamics of post-translational modification systems: recent progress and future directions Eq. 3.4",
	[
  "dx1/dt = -k1 * x1 * x2 + k2 * x4 + k4 * x6;\n",
  "dx2/dt = k1 * x1 * x2 + k2 * x4 + k3 * x4;\n",
  "dx3/dt = k3 * x4 + k5 * x6 - k6 * x3 * x5;\n",
  "dx4/dt = k1 * x1 * x2 - k2 * x4 - k3 * x4;\n",
  "dx5/dt = k4 * x6 + k5 * x6 - k6 * x3 * x5;\n",
  "dx6/dt = -k4 * x6 - k5 * x6 + k6 * x3 * x5;\n",
  "y1 = x3;\n"
  "y2 = x2;\n" ]],

	"DAISY Ex. 3" = ["DAISY Example 3", [
  "dx1/dt = -1 * p1 * x1 + p2 * x2 + u(t);\n",
  "dx2/dt = p3 * x1 - p4 * x2 + p5 * x3;\n",
  "dx3/dt = p6 * x1 - p7 * x3;\n",
  "y1 = x1;\n"]],

	"DAISY_mamil3" = ["DAISY mamil 3",
	[
  "dx1/dt = -(a21 + a31 + a01) * x1 + a12 * x2 + a13 * x3 + u(t);\n",
  "dx2/dt = a21 * x1 - a12 * x2;\n",
  "dx3/dt = a31 * x1 - a13 * x3;\n",
  "y = x1"]],

	"DAISY_mamil4" = ["DAISY mamil 4", [
  "dx1/dt = -k01 * x1 + k12 * x2 + k13 * x3 + k14 * x4 - k21 * x1 - k31 * x1 - k41 * x1 + u(t);\n",
  "dx2/dt = -k12 * x2 + k21 * x1;\n",
  "dx3/dt = -k13 * x3 + k31 * x1;\n",
  "dx4/dt = -k14 * x4 + k41 * x1;\n",
  "y = x1"]],

	"HIV" = ["Example (with initial conditions assumed being unknown) from Section IV of 'DAISY: an Efficient Tool to Test Global Identifiability. Some Case Studies' by G. Bellu, M.P. Saccomani",
	[
  "dx1/dt = -b * x1 * x4 - d * x1 + s;\n",
  "dx2/dt = b * q1 * x1 * x4 - k1 * x2 - mu1 * x2;\n",
  "dx3/dt = b * q2 * x1 * x4 + k1 * x2 - mu2 * x3;\n",
  "dx4/dt = -c * x4 + k2 * x3;\n",
  "y1 = x1;\n",
  "y2 = x4"]],

	"HIV2" = ["The system is taken from Wodarz, D., Nowak, M.\nSpecific therapy regimes could lead to long-term immunological control of HIV\nhttps://doi.org/10.1073/pnas.96.25.14464\nPage 1",
	[
  "dx/dt = lm - d * x - beta * x * v;\n",
  "dy/dt = beta * x * v - a * y;\n",
  "dv/dt = k * y - u * v;\n",
  "dw/dt = c * x * y * w - c * q * y * w - b * w;\n",
  "dz/dt = c * q * y * w - h * z;\n",
  "y1 = w;\n",
  "y2 = z"]],

	"Lipolysis" = ["Taken from R. Munoz-Tamayo, L. Puillet, J.B. Daniel, D. Sauvant, O. Martin, M. Taghipoor, P. Blavy\nReview: To be or not to be an identifiable model. Is this a relevant question in animal science modelling?\ndoi.org/10.1017/S1751731117002774\nSystem (1) in Supplementary Material 2, initial conditions are assumed to be unknown\nbrought to the rational function form by introducing new state variable x5 = k1 e^(-k3 t)",
	[
  "dx1/dt = -x1 * x5 / (k2 + x1);\n",
  "dx2/dt = 2 * x1 * x5 / ((k2 + x1) * 3) - k4 * x2;\n",
  "dx3/dt = k4*(x2)/2 - k4*x3;\n",
  "dx4/dt = x1 * x5 / (3 * (k2 + x1)) + k4 * (x2)/2 + k4 * x3;\n",
  "dx5/dt = -k3 * x5;\n",
  "y1 = x1;\n",
  "y2 = x2 + x3;\n",
  "y3 = x4"]],

  	"LV" = ["Lotka-Volterra System",[
  	"dx1/dt = a*x1 - b*x1*x2;\n",
  	"dx2/dt = -c*x2 + d*x1*x2;\n",
  	"y = x1;\n"]],
	"OralGlucose" = ["Example (with initial conditions assumed being unknown) from Section III of 'DAISY: an Efficient Tool to Test Global Identifiability. Some Case Studies'\nby G. Bellu, M.P. Saccomani",
	[
  "dG/dt = -(p1 + X) * G + p1 * Gb + v * R;\n",
  "dX/dt = -p2 * X + p3 * (u(t) - Ib);\n",
  "dR/dt = k;\n",
  "dIb/dt = 0;\n",
  "dGb/dt = 0;\n",
  "y1 = G;\n",
  "y2 = Ib;\n",
  "y3 = Gb;\n"]],

	"SEIR" = ["Taken from N. Tuncer, T. Le\n'Structural and practical identifiability analysis of outbreak models'\nhttps://doi.org/10.1016/j.mbs.2018.02.004\nEquation (2.2) with prevalence observations",
[
  "dS/dt = -b * S * In / N;\n",
  "dE/dt = b * S * In / N - nu * E;\n",
  "dIn/dt = nu * E - a * In;\n",
  "dN/dt = 0;\n",
  "y1 = In;\n",
  "y2 = N;\n"]],

	"SEIR2" = ["Taken from N. Tuncer, T. Le\n'Structural and practical identifiability analysis of outbreak models'\nhttps://doi.org/10.1016/j.mbs.2018.02.004\nEquation (2.2) with cumulative incidence observations",
	[
  "dS/dt = -b * S * In / N;\n",
  "dE/dt = b * S * In / N - nu * E;\n",
  "dIn/dt = nu * E - a * In;\n",
  "dN/dt = 0;\n",
  "dCu/dt = nu * E;\n",
  "y1 = Cu;\n",
  "y2 = N"]],

	"SIR_R0" = ["SIR R0",[
  "dS/dt = -b * In * S;\n",
  "dIn/dt = b * In * S - g * In;\n",
  "dR/dt = g * In;\n",
  "daux/dt = 0;\n",
  "y1 = In;\n",
  "y2 = b / g + aux;"]],

	"SIRSForced" = ["Taken from Capistran M., Moreles M., Lara B.\n'Parameter Estimation of Some Epidemic Models.\n The Case of Recurrent Epidemics Caused by Respiratory Syncytial Virus'\ndoi.org/10.1007/s11538-009-9429-3\nEquations (7)-(11)",
[
  "ds/dt = mu - mu * s - b0 * (1 + b1 * x1) * i * s + g * r;\n",
  "di/dt = b0 * (1 + b1 * x1) * i * s - (nu + mu) * i;\n",
  "dr/dt = nu * i - (mu + g) * r;\n",
  "dx1/dt = -M * x2;\n",
  "dx2/dt = M * x1;\n",
  "y1 = i;\n",
  "y2 = r;\n"]],

	"SlowFast" = ["Taken from Vajda S., Rabitz H.\n'Identifiability and Distinguishability of First-Order Reaction Systems', p. 701\nWe added an extra output x_C",
	[
  "dxA/dt = -k1 * xA;\n",
  "dxB/dt = k1 * xA - k2 * xB;\n",
  "dxC/dt = k2 * xB;\n",
  "deA/dt = 0;\n",
  "deC/dt = 0;\n",
  "y1 = eA * xA + eB * xB + eC * xC;\n",
  "y2 = xC;\n",
  "y3 = eA;\n",
  "y4 = eC"]],

	"Treatment" = ["Taken from N. Tuncer, T. Le\nStructural and practical identifiability analysis of outbreak models'\nhttps://doi.org/10.1016/j.mbs.2018.02.004\nEquation (2.3) with observed treatment",
	[
  "dS/dt = -b * S * In / N - d * b * S * Tr / N;\n",
  "dIn/dt = b * S * In / N + d * b * S * Tr / N - (a + g) * In;\n",
  "dTr/dt = g * In - nu * Tr;\n",
  "dN/dt = 0;\n",
  "y1 = Tr;\n",
  "y2 = N"]],

	"Tumor" = ["Example (with initial conditions assumed being unknown) from Section 3 of\n'Examples of testing global identifiability of biological and biomedical models with the DAISY software'\nby M.P. Saccomani, S. Audoly, G. Bellu, L. D'Angio",
[ "dx1/dt = -(k3 + k7) * x1 + k4 * x2;\n",
  "dx2/dt = k3 * x1 - (k4 + a * k5 + b * d1 * k5) * x2 + k6 * x3 + k6 * x4 + k5 * x2 * x3 + k5 * x2 * x4;\n",
  "dx3/dt = a * k5 * x2 - k6 * x3 - k5 * x2 * x3;\n",
  "dx4/dt = b * d1 * k5 * x2 - k6 * x4 - k5 * x2 * x4;\n",
  "dx5/dt = k7 * x1;\n",
  "da/dt = 0;\n",
  "db/dt = 0;\n",
  "dd1/dt = 0;\n",
  "y1 = x5;\n",
  "y2 = a;\n",
  "y3 = b;\n",
  "y4 = d1;\n"]]
]):

# Setup

timed_SIAN:=proc(sigma, params_to_assess, p, output_targets_sian, count_solutions, char)
	local output, data, start, finish:
	start:= time():
  #IdentifiabilityODE := proc(system_ODEs, params_to_assess, output_targets_sian, {char:=0, p := 0.99, count_solutions:=true, infolevel := 1, method := 2, num_nodes := 6})
	output := IdentifiabilityODE(sigma, params_to_assess, output_targets_sian, char=char, p=p, count_solutions=count_solutions, infolevel=1, method=2, num_nodes=1):
	finish:= time():
	DocumentTools:-SetProperty(output_targets_sian[runningtime], value, convert(finish-start, string), 'refresh'): # time
	return  output:
end proc:

timed_Multi:=proc(model, simplified_generators, no_bound, simplify_bound, max_perms, output_targets_multi)
	local start, output, finish, data, bound, generators:
	start:=time():
	bound, generators := op(MultiExperimentIdentifiableFunctions(model, simplified_generators, no_bound, simplify_bound, max_perms, output_targets_multi)):
	finish:=time():
	DocumentTools:-SetProperty(output_targets_multi[runningtime], value, convert(finish-start, string), 'refresh'):
	return [bound, finish-start, generators]: #table([output=bound, runtime=finish-start]):
end proc:

timed_Single:=proc(model, output_targets_single)
	local start, output, finish, data:
	start:=time():
	output := SingleExperimentIdentifiableFunctions(model, output_targets_single):
	finish:=time():
	DocumentTools:-SetProperty(output_targets_single[runningtime], value, convert(finish-start, string), 'refresh'):
	return output:
end proc:

with(StringTools):

sigmaParser := proc(sigma)
	local states, state_eqs, outputs, output_eqs;
	if SearchText("diff", sigma)>0 then
		sigma := [map(x->parse(x), Split(sigma, ";"))]:
	else
		LogExpression(sprintf("%q \n", Split(sigma, ";")), "LogAreaSIAN"):
		sigma := Split(sigma, ";"):

		states := map(x->Trim(RegSubs("d([a-zA-Z0-9]+)/dt(.*)" = "\\1", x)), select(x->SearchText(x, "/dt")>0, sigma)):
		state_eqs := select(x->Has(x, "/dt"), sigma):

		outputs :=  map(x->Trim(Split(x, "=")[1]), select(x->not Has(x, "/dt"), sigma)):
		output_eqs := select(x->not SearchText(x, "/dt")>0, sigma):

		state_eqs := map(x->convert(subs({seq(parse(states[i])=parse(cat(states[i],"(t)")), i=1..nops(states))}, parse(x)), string),  state_eqs):
		state_eqs := map(x->parse(RegSubs("d([a-zA-Z0-9]+)/dt(.*)" = "diff(\\1(t), t)\\2", x)), state_eqs):

		output_eqs := map(x->parse(Subs({seq(outputs[i]=cat(outputs[i],"(t)"), i=1..nops(outputs))}, x)), output_eqs):
		output_eqs := map(x->subs({seq(parse(states[i])=parse(cat(states[i],"(t)")), i=1..nops(states))}, x),  output_eqs):

		sigma := [op(state_eqs), op(output_eqs)]:
	end if:
	return sigma;
end proc:

DocumentTools:-SetProperty("RunningTimeSingle", value, ""):
DocumentTools:-SetProperty("RunningTimeMulti", value, ""):
DocumentTools:-SetProperty("RunningTimeSIAN", value, ""):

DocumentTools:-SetProperty("run_system", enabled, true):
DocumentTools:-SetProperty("Meter_sian", visible, true):
DocumentTools:-SetProperty("Meter_sian", value, 0):
DocumentTools:-SetProperty("sigma", enabled, true):
DocumentTools:-SetProperty("p", value, "0.99"):
DocumentTools:-SetProperty("params", value, ""):
DocumentTools:-SetProperty("replicas", value, "1"):

DocumentTools:-SetProperty("p", enabled, true):
DocumentTools:-SetProperty("params", enabled, true):
DocumentTools:-SetProperty("replicas", enabled, true):
DocumentTools:-SetProperty("LogAreaSE", value, ""):
DocumentTools:-SetProperty("LogAreaSIAN", value, ""):
DocumentTools:-SetProperty("LogAreaME", value, ""):

DocumentTools:-SetProperty("printSolutions", enabled, true):
DocumentTools:-SetProperty("printSolutions", value, true):

DocumentTools:-SetProperty("GlobalParams1", expression, NULL):
DocumentTools:-SetProperty("LocalParams1", expression, NULL):
DocumentTools:-SetProperty("NoIDParams1", expression, NULL):

DocumentTools:-SetProperty("Parameters", value, ""):

DocumentTools:-SetProperty("Bound", expression, NULL):
DocumentTools:-SetProperty("MultiFunctions", expression, NULL):
DocumentTools:-SetProperty("SingleFunctions", expression, NULL):

DocumentTools:-SetProperty("use_char", enabled, true):
DocumentTools:-SetProperty("use_char", value, false):
DocumentTools:-SetProperty("p_label", enabled, true):

DocumentTools:-SetProperty("ComputeId", enabled, true):
DocumentTools:-SetProperty("ComputeId", value, false):
DocumentTools:-SetProperty("bypass", enabled, false):

DocumentTools:-SetProperty("bypass", value, true):
DocumentTools:-SetProperty("SimplifiedGen", enabled, false):
DocumentTools:-SetProperty("SimplifiedGen", value, true):
DocumentTools:-SetProperty("SkipSingle", enabled, false):
DocumentTools:-SetProperty("SkipSingle", value, false):
DocumentTools:-SetProperty("Refine", enabled, false):
DocumentTools:-SetProperty("NoBound", enabled, false):
DocumentTools:-SetProperty("NoBound", value, false):
DocumentTools:-SetProperty("UsingUpTo", enabled, false):
DocumentTools:-SetProperty("MaxPermutations", enabled, false):
DocumentTools:-SetProperty("Permutations", enabled, false):
DocumentTools:-SetProperty("RunSIAN", enabled, true):
DocumentTools:-SetProperty("RunSIAN", value, true):
DocumentTools:-SetProperty("being_refined", caption, "");
DocumentTools:-SetProperty("sigma", value, "dx1/dt = a*x1 + x2*b + u(t);\ndx2/dt = x2*c + x1;\ny=x2"):
DocumentTools:-SetProperty("example_box", value, "Custom"):
DocumentTools:-SetProperty(reference, value, ""):
DocumentTools:-SetProperty("LocalLabel" , caption, "Locally Identifiable Paramters");
DocumentTools:-SetProperty(TxtOutput, visible, false);
DocumentTools:-SetProperty(TxtOutput, value, "");
DocumentTools:-SetProperty(SaveOutputLabel, visible, false);

readyToSave:=false:
counter:=0:
exname:="Custom"

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

$define gga_c_pbe_params
$include "gga_c_pbe.mpl"

pbeloc_b0 := 0.0375:
pbeloc_a  := 0.08:

(* we redefine beta here *)
mbeta := (rs, t) -> pbeloc_b0 + pbeloc_a*t^2*(1 - exp(-rs^2)):