text:
`is_element/digraphs` := (V::set) -> proc(E)
 if not(`is_element/autorel`(V)(E)) then
  return false;
 fi;

 if not(`is_antisymmetric/autorel`(V)(E)) then
  return false;
 fi;

 return true;
end:

`is_leq/digraphs` := (V::set) -> (E0,E1) -> evalb(E0 minus E1 = {}):

`list_elements/digraphs` := proc(V::set)
 local n,EE,E,i,j,a,b;

 n := nops(V);
 EE := [{}];
 for i from 1 to n-1 do
  for j from i+1 to n do
   a := V[i];
   b := V[j];
   EE := [op(EE),
          seq(E union {[a,b]},E in EE),
          seq(E union {[b,a]},E in EE)];
  od;
 od;

 return EE;
end:

`count_elements/digraphs` := (V::set) -> 3^(nops(V)*(nops(V)-1)/2);

`random_element/digraphs` := (V::set) -> proc()
 local n,E,r,i,j,a,b,x;

 n := nops(V);
 E := NULL;
 r := rand(-1..1);
 for i from 1 to n-1 do
  for j from i+1 to n do
   a := V[i];
   b := V[j];
   x := r();
   if x = -1 then
    E := E,[b,a];
   elif x = 1 then
    E := E,[a,b];
   fi;
  od;
 od;

 return {E};
end:

`neighbour_table/digraphs` := (V) -> proc(E)
 local N,v,e;

 N := table();
 for v in V do N[v] := {}; od:
 for e in E do
  N[e[1]] := {op(N[e[1]]),e[2]};
  N[e[2]] := {op(N[e[2]]),e[1]};
 od:

 return eval(N);
end:

`in_neighbour_table/digraphs` := (V) -> proc(E)
 local N,v,e;

 N := table();
 for v in V do N[v] := {}; od:
 for e in E do
  N[e[2]] := {op(N[e[2]]),e[1]};
 od:

 return eval(N);
end:

`out_neighbour_table/digraphs` := (V) -> proc(E)
 local N,v,e;

 N := table();
 for v in V do N[v] := {}; od:
 for e in E do
  N[e[1]] := {op(N[e[1]]),e[2]};
 od:

 return eval(N);
end:

`component_relation/digraphs` := (V) -> (E) ->
 `transitive_closure/autorel`(V)(
  E union `op/autorel`(V)(E) union `id/autorel`(V)
 );

`components/digraphs` := (V) -> (E) ->
 `block_partition/equiv`(`component_relation/graphs`(V)(E));

`is_connected/digraphs` := (V) -> (E) ->
 evalb(`component_relation/digraphs`(V)(E) = `top/autorel`(V));

`skeleton/digraphs` := (V::set) -> proc(E)
 local n,V0,E0,ix,i;

 n := nops(V);
 V0 := {seq(i,i=1..n)};
 ix := table():
 for i from 1 to n do ix[V[i]] := i; od;
 E0 := map(e -> [ix[e[1]],ix[e[2]]],E);
 return [V0,E0];
end:

`is_element/dipaths` := (V) -> (E) -> proc(p)
 local i,n;

 if not type(p,list) then
  return false;
 fi;

 n := nops(p) - 1;

 if n < 0 then
  return false;
 fi;

 if {op(p)} minus V <> {} then
  return false;
 fi;

 for i from 1 to n do
  if not(member([p[i],p[i+1]],E)) then
   return false;
  fi;
 od;

 return true;
end:

`is_small_element/dipaths` := (V) -> (E) -> (n::nonnegint) -> proc(p)
 return evalb(`is_element/dipaths`(V)(E)(p) and (nops(p) = n+1));
end:

`random_small_element/dipaths` := (V) -> (E) -> (n::nonnegint) -> proc(num_tries := 10)
 local k,p,F;

 if n = 0 then
  if V = {} then
   return FAIL;
  else
   return [random_element_of(V)()];
  fi;
 elif n = 1 then
  if E = {} then
   return FAIL;
  else
   return random_element_of(E)();
  fi;
 fi;

 k := num_tries;

 while k > 0 do
  k := k-1;
  p := `random_small_element/dipaths`(V)(E)(n-1)();
  if p <> FAIL then
   F := select(e -> e[2] = p[1],E);
   if F <> {} then
    return [random_element_of(F)()[1],op(p)];
   fi;
  fi;
 od;

 return FAIL;
end:

`list_small_elements/dipaths` := (V) -> (E) -> proc(n::nonnegint)
 option remember;
 local P,Q,p,v,F,e;

 if n = 0 then
  return map(v -> [v],V);
 else
  P := `list_small_elements/dipaths`(V)(E)(n-1);
  Q := NULL;
  for p in P do
   v := p[1];
   F := select(e -> e[2] = v,E);
   Q := Q,seq([e[1],op(p)],e in F);
  od:
  return [Q];
 fi;
end:

`is_leq/dipaths` := NULL;
`count_elements/dipaths` := NULL:
`count_small_elements/dipaths` := NULL:

`length/dipaths` := (V) -> (E) -> (p) -> nops(p) - 1;

`is_trail/dipaths` := proc(p)
 local n,e,i;
 n := nops(p) - 1;
 e := {seq([p[i],p[i+1]],i=1..n)};
 return evalb(nops(e) = n);
end:

`is_cycle/dipaths` := proc(p)
 local n,i;
 n := nops(p) - 1;
 return evalb(p[1] = p[n+1] and nops({seq(p[i],i=1..n)}) = n);
end:

`is_small_element/ditrails` := (V) -> (E) -> (n::nonnegint) -> proc(p)
 if not `is_small_element/dipaths`(V)(E)(n)(p) then
  return false;
 fi;

 return `is_trail/dipaths`(p);
end:

`is_element/ditrails` := (V) -> (E) -> proc(p)
 if not `is_element/dipaths`(V)(E)(p) then
  return false;
 fi;

 return `is_trail/dipaths`(p);
end:

`list_small_elements/ditrails` := (V) -> (E) -> proc(n)
 option remember;
 local E0,E1,P,Q,p,v,F,i,e;

 if n = 0 then
  return map(v -> [v],V);
 else
  P := `list_small_elements/ditrails`(V)(E)(n-1);
  Q := NULL;
  for p in P do
   E0 := {seq([p[i],p[i+1]],i=1..n-1)};
   E1 := E minus E0;
   v := p[1];
   F := select(e -> e[2] = v,E1);
   Q := Q,seq([e[1],op(p)],e in F);
  od:
  return [Q];
 fi;
end:

`random_small_element/ditrails` := (V) -> (E) -> (n) -> proc(num_tries := 10)
 local i,k,p,E0,E1,v,F,e;

 k := num_tries;

 if n = 0 then
  if V = {} then
   return FAIL;
  else
   return [random_element_of(V)()];
  fi;
 elif n = 1 then
  if E = {} then
   return FAIL;
  else
   return random_element_of(E)();
  fi;
 fi;

 while k > 0 do
  p := `random_small_element/ditrails`(V)(E)(n-1)();
  k := k - 1;
  if p <> FAIL then
   E0 := {seq([p[i],p[i+1]],i=1..n-1)};
   E1 := E minus E0;
   v := p[1];
   F := select(e -> e[2] = v,E1);
   if F <> {} then
    e := random_element_of(F)();
    return [e[1],op(p)];
   fi;
  fi;
 od;

 return FAIL;
end:

`is_leq/ditrails` := NULL:
`count_elements/ditrails` := NULL:
`count_small_elements/ditrails` := NULL:

`roots/digraphs` := (V) -> proc(E)
 V minus map(e -> e[2],E);
end:

`root/digraphs` := (V) -> proc(E)
 local R;
 R := `roots/digraphs`(V)(E);
 if nops(R) = 1 then
  return R[1];
 else
  error "Digraph does not have a unique root";
 fi;
end:

`leaves/digraphs` := (V) -> proc(E)
 V minus map(e -> e[1],E);
end:

`is_forest/digraphs` := (V::set) -> proc(E)
 local NI,NO,R,v,U,U1;

 NI := `in_neighbour_table/digraphs`(V)(E);
 NO := `out_neighbour_table/digraphs`(V)(E);
 R := NULL;

 for v in V do
  if nops(NI[v]) > 1 then
   return false;
  elif NI[v] = {} then
   R := R,v;
  fi;
 od;

 R := {R};

 U := R;
 U1 := {op(U),seq(op(NO[v]),v in U)};
 while nops(U1) > nops(U) do
  U := U1;
  U1 := {op(U),seq(op(NO[v]),v in U)};
 od:

 if U <> V then return false; fi;

 return true;
end:

`is_tree/digraphs` := (V) -> (E) ->
 nops(`roots/digraphs`(V)(E)) = 1 and
  `is_forest/digraphs`(V)(E);

`level_table/forest_digraphs` := (V) -> proc(E)
 local X,N,L,k,v;

 X := `roots/digraphs`(V)(E);
 N := `out_neighbour_table/digraphs`(V)(E);
 L := table():
 k := 0;
 while X <> {} do
  for v in X do L[v] := k; od:
  X := map(v -> op(N[v]),X);
  k := k+1;
 od:

 return eval(L);
end:

`shadow_table/forest_digraphs` := (V) -> proc(E)
 local X,Y,Z,S,x,z,N,u;

 X := `leaves/digraphs`(V)(E);
 Y := V minus X;
 S := table():
 for x in X do S[x] := {x}; od:
 N := `out_neighbour_table/digraphs`(V)(E);
 while Y <> {} do
  Z := select(y -> ((N[y] minus X) = {}),Y);
  if Z = {} then
   error "Something wrong";
  fi;
  z := Z[1];
  S[z] := map(op,{seq(S[u],u in N[z])});
  X := X union {z};
  Y := Y minus {z};
 od:
 return eval(S);
end:

`shadow_set/forest_digraphs` := (V) -> proc(E)
 local S,v;
 S := `shadow_table/forest_digraphs`(V)(E);
 return {seq(S[v],v in V)};
end:

`is_element/stasheff_digraphs` := (V::set) -> proc(EL)
 local E,L,L0,r,i,TT,rTT;
 if not (type(EL,list) and nops(EL) = 2) then
  return false;
 fi;

 E,L := op(EL);

 if not(`is_tree/digraphs`(V)(E)) then
  return false;
 fi;

 L0 := `leaves/digraphs`(V)(E);
 if not({op(L)} = L0 and nops(L0) = nops(L)) then
  return false;
 fi;

 r := table();
 for i from 1 to nops(L) do r[L[i]] := i; od;

 TT := `shadow_set/forest_digraphs`(V)(E);
 rTT := map(T -> map(v -> r[v],T),TT);

 return `and`(op(map(`is_element/posint_intervals`,rTT)));
end:


text:
######################################################################
# W(N)(A) is Map(A,R^N)/Const(A,R^N).
#
# Elements are represented as tables indexed by A, whose values
# are lists of length N of real numbers.  The equivalence relation
# is built into the definition of the function `is_equal/W` rather
# than the representation of the elements.

`is_element/W` := (N::posint) -> (A::set) -> proc(x)
 `is_element/vector_functions`(N)(A)(x);
end;

######################################################################

`normalise/W` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,a;

 x0 := table();
 u := [0$N];

 for a in A do
  x0[a] := x[a];
  u := u +~ x[a];
 od;
 u := u /~ nops(A);
 for a in A do
  x0[a] := x0[a] -~ u;
 od;
 return(eval(x0));
end;

######################################################################

`bottom_normalise/W` := (N::posint) -> (A::set) -> proc(x)
 local x0,u,a,i;

 u := [seq(min(seq(x[a][i],a in A)),i=1..N)];
 x0 := table():
 for a in A do
  x0[a] := x[a] -~ u;
 od:
 return(eval(x0));
end;

######################################################################

`is_equal/W` := (N::posint) -> (A::set) -> proc(x,y)
 local x0,y0,a;

 x0 := `normalise/W`(N)(A)(x);
 y0 := `normalise/W`(N)(A)(y);

 return `is_equal/vector_functions`(N)(A)(x0,y0);
end;

######################################################################

`is_leq/W` := NULL;

######################################################################

`is_zero/W` := (N::posint) -> (A::set) -> proc(x)
 local x0,a;

 x0 := `normalise/W`(N)(A)(x);
 return `is_zero/vector_functions`(N)(A)(x0);
end;

######################################################################

`norm/W` := (N::posint) -> (A::set) -> proc(x)
 local x0,a;

 x0 := `normalise/W`(N)(A)(x);
 return `norm/vector_functions`(N)(A)(x0);
end;

######################################################################

`width/W` := (N::posint) -> (A::set) -> proc(x)
 local i,a;
 max(seq(max(seq(x[a][i],a in A)) -  min(seq(x[a][i],a in A)),i=1..N));
end;

######################################################################

`random_element/W` := (N::posint) -> (A::set) -> proc()
 local x,a,u;

 u := [0$N];
 x := table();
 for a in A do
  x[a] := `random_element/R`(N)();
  u := u +~ x[a];
 od:
 u := u /~ nops(A);
 for a in A do
  x[a] := x[a] -~ u;
 od;

 return eval(x);
end;

`list_elements/W` := NULL;
`count_elements/W` := NULL;

######################################################################
# Here C is assumed to be a subset of B, and rho is the usual map
# from W(N,B) to W(N,C).

`rho/W` := (B,C) -> proc(x)
 local y,c;

 y := table();
 for c in C do y[c] := x[c]; od;
 return eval(y);
end;

ext: mpl


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 *)

rge2_kappa := 0.8040:

rge2_den := s -> rge2_kappa + 1*MU_GE*s^2 + MU_GE^2*s^4/rge2_kappa:
rge2_f0  := s -> 1 + rge2_kappa * (1 - rge2_kappa/rge2_den(s)):
rge2_f   := x -> rge2_f0(X2S * x):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange(rge2_f, rs, zeta, xs0, xs1):


text:
   {  DakotaScen.mpl  }

TITLE
   DakotaScen;

INDEX
   product  := (Desk, Tble, Chair);
   resource := (Lumber, Finishing, Carpentry);

STOCHASTIC

TREE
   TWOSTAGE;

SCENARIOS
   scen := 1..3;

PROBABILITIES
   Prob[scen] := (0.3, 0.4, 0.3);

RANDOM DATA
   Demand[product,scen] := (50, 150, 250,
                            20, 110, 250,
                           200, 225, 500);

DATA
   Cost[resource] := (2, 4, 5.2);
   ProdReq[resource, product] := (8,   6,   1,
                                  4,   2,   1.5,
                                  2,   1.5, 0.5);
   Price[product] := (60, 40, 10);

   Budget := 8000;

DECISION VARIABLES
   Acquire[resource]

STAGE2 VARIABLES
   Produce[product];

MODEL
   Max Profit =  SUM(product: Price * Produce)
               - SUM(resource:  Cost * Acquire);

SUBJECT TO
   MaxBudget[resource]: Cost * Acquire <= Budget;

   MaxProduce[product]:  Produce <= Demand;

   MinAcquire[resource]: SUM(product: ProdReq * Produce) <= Acquire;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
# This file sets up the 15-vertex triangulation of CP^2 introduced by
# Gaifullin in "A minimal triangulation of complex projective plane ..."
# (arxiv:0904.4222)

with(GroupTheory):

CP2_complex := proc()
 local T,V,F,p,a,b,d,s,S,om,ob,rt,
  id_S,o_S,inv_S,act_S,eq_S,eq_X,v,ix_v,i;

 T := table():

 # We take the vertices to be 1,...,15
 # In Gaifullin's presentation, the vertex set is
 # {transposition pairs} union {1,2,3,4} x {1,2,3}
 # Gaifullin's (i,j) is our 1 + 3 * (i-1) + (j-1)
 # The transposition pairs are our vertices 13-15.

 V := {seq(i,i=1..15)}:

 T["vertices"] := V;

 F := NULL:

 # This table encodes the three transposition pairs on {0,1,2,3}
 p := table():
 p[1,0] := 1; p[1,1] := 0; p[1,2] := 3; p[1,3] := 2;
 p[2,0] := 2; p[2,1] := 3; p[2,2] := 0; p[2,3] := 1;
 p[3,0] := 3; p[3,1] := 2; p[3,2] := 1; p[3,3] := 0;

 for a from 1 to 3 do
  for b[0] from 0 to 2 do
   for b[1] from 0 to 2 do
    for b[2] from 0 to 2 do
     for b[3] from 0 to 2 do
      if b[0] <> b[p[a,0]] and
	 b[1] <> b[p[a,1]] and
	 b[2] <> b[p[a,2]] and
	 b[3] <> b[p[a,3]] then
       F := F,{b[0]+1,b[1]+4,b[2]+7,b[3]+10,a+12};
      fi;
     od:
    od:
   od:
  od:
 od:
 F := {F}:

 T["max_simplices"] := F;
 T["all_simplices"] :=
   map(f -> op(combinat[powerset](f)),F) minus {{}}:

 T["simplices"] := table():
 for d from 0 to 4 do
  T["simplices"][d] := select(u -> nops(u) = d+1,T["simplices"]);
 od:

 s := table():
 T["s"] := eval(s);

 s[0] := Perm([]):
 s[1] := Perm([[ 1, 4],[ 2, 5],[ 3, 6],[14,15]]):
 s[2] := Perm([[ 4, 7],[ 5, 8],[ 6, 9],[13,14]]):
 s[3] := Perm([[ 7,10],[ 8,11],[ 9,12],[14,15]]):
 s[4] := Perm([[ 1, 2],[ 4, 5],[ 7, 8],[10,11]]):
 s[5] := Perm([[ 2, 3],[ 5, 6],[ 8, 9],[11,12]]):

 T["S4"] := Group(s[1],s[2],s[3]);
 T["S3"] := Group(s[4],s[5]);
 T["G"]  := Group(seq(s[i],i=1..5));

 S := table():
 T["S"] := eval(S);

 om := exp( 2*Pi*I/3);
 ob := exp(-2*Pi*I/3);
 rt := sqrt(3);

 S[0] := [ 1,[[ 1, 0, 0],[ 0, 1, 0],[ 0, 0, 1]]]:
 S[1] := [-1,[[ 0, 1, 0],[ 1, 0, 0],[ 0, 0, 1]]];
 S[2] := [-1,[[ 1, 0, 0],[ 0, 0, 1],[ 0, 1, 0]]];
 S[3] := [-1,[[ 0,-1, 0],[-1, 0, 0],[ 0, 0, 1]]];
 S[4] := [-1,[[ 1, 0, 0],[ 0, 1, 0],[ 0, 0, 1]]];
 S[5] := [-1,[[ 1, 0, 0],[ 0,om, 0],[ 0, 0,ob]]];

 id_S := S[0];

 o_S := proc(u,v)
  local m,n,mn,A,B,AB;
  m,A := op(u);
  n,B := op(v);
  mn := m*n;
  if n = -1 then A := conjugate(A); fi;
  AB := convert(Matrix(A) . Matrix(B),listlist);
  AB := simplify(expand(rationalize(AB)));
  return [mn,AB];
 end:

 inv_S := proc(u)
  local m,A;
  m,A := op(u);
  if m = -1 then A := conjugate(A); fi;
  A := convert(1/Matrix(A),listlist);
  A := simplify(expand(rationalize(A)));
  return [m,A];
 end:

 act_S := (u) -> proc(x)
  local y;
  y := Vector(x);
  y := simplify(convert(Matrix(u[2]) . y,list));
  if u[1] = -1 then y := map(conjugate,y); fi;
  return y;
 end:

 eq_S := proc(u,v)
  local w;
  if u[1] <> v[1] then return false; fi;
  w := o_S(u,inv_S(v));
  if w[2][1,1] = 0 then return false; fi;
  w := [w[1],map(simplify,map(expand,w[2] /~ w[2][1,1]))];
  return evalb(w = id_S);
 end:

 eq_X := proc(u,v)
  local i,j;
  `and`(seq(seq(evalb(simplify(u[i]*v[j]-u[j]*v[i])=0),j=i+1..3),i=1..2));
 end:

 v := table([
   1 = [  1,-om,-ob],
   2 = [  1,-ob,-om],
   3 = [  1, -1, -1],
   4 = [  1,-om, ob],
   5 = [  1,-ob, om],
   6 = [  1, -1,  1],
   7 = [  1, om,-ob],
   8 = [  1, ob,-om],
   9 = [  1,  1, -1],
  10 = [  1, om, ob],
  11 = [  1, ob, om],
  12 = [  1,  1,  1],
  13 = [  0,  0, rt],
  14 = [  0, rt,  0],
  15 = [ rt,  0,  0]
 ]):

 T["projective_embedding"] := eval(v);

 ix_v := proc(u)
  local i;
  for i from 1 to 15 do
   if eq_X(u,v[i]) then return i; fi;
  od;
  return FAIL;
 end:

 return eval(T):
end:


text:
func $singlefoo(var %a i32) i32 {
  return (dread i32 %a)
}

func $fact(var %n i32) i32 {
  var %m i32
  var %r1 i32
  var  %r2 i32
  var  %r3 i32
  callassigned &singlefoo(dread i32 %n) {
    dassign %m 0
  }
  # Not correct, but for parse and dump test.
  callassigned &singlefoo(dread i32 %n) {
    dassign %r1 1
    dassign %r2 2
    dassign %r3 3
  }
  return (dread i32 %n)
}
 # 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:
cube_complex := proc()
 local T,V,E,F,G0,G1,G,S,r,v,f,i,j,k,s,s1,t,g,eqs,sol;

 V := [seq(i,i=1..8)];

 T["vertices"] := V;

 T["short_edges"] := [
  [1,2],[1,3],[1,5],[2,4],[2,6],[3,4],[3,7],[4,8],[5,6],[5,7],[6,8],[7,8]
 ];

 T["long_edges"] := [[1,4],[1,6],[1,7],[2,8],[3,8],[4,8]];

 T["edges"] := sort([op(T["short_edges"]),op(T["long_edges"])]);

 T["squares"] := [
  [1,2,3,4],[1,2,5,6],[1,3,5,7],[2,4,6,8],[3,4,7,8],[5,6,7,8]
 ];

 T["faces"] := [
  [1,2,4],[1,3,4],[1,2,6],[1,5,6],[1,3,7],[1,5,7],
  [2,4,8],[2,6,8],[3,4,8],[3,7,8],[5,6,8],[5,7,8]
 ];

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

 T["vertex_index"]     := make_index(T["vertices"]):
 T["short_edge_index"] := make_index(T["short_edges"]):
 T["square_index"]     := make_index(T["squares"]):

 T["embedding_dim"] := 3;

 T["embedding"] := table([
  1 = [-1,-1,-1],
  2 = [ 1,-1,-1],
  3 = [-1, 1,-1],
  4 = [ 1, 1,-1],
  5 = [-1,-1, 1],
  6 = [ 1,-1, 1],
  7 = [-1, 1, 1],
  8 = [ 1, 1, 1]
 ]);

 v := eval(T["embedding"]);

 `plot/simplicial_complex`(T);

 S := NULL;
 for i from 1 to nops(T["squares"]) do
  f := T["squares"][i];
  f := [f[1],f[2],f[4],f[3]];
  S := S, polygon(map(j -> T["embedding"][j],f),colour=standard_colour(i));
 od:

 T["square_plot"] := display(S,scaling=constrained,axes=none);

 T["short_edge_centres"] :=
  map(e -> (v[e[1]] +~ v[e[2]]) /~ 2,T["short_edges"]);

 T["square_centres"] :=
  map(f -> (v[f[1]] +~ v[f[2]] +~ v[f[3]] +~ v[f[4]]) /~ 4, T["squares"]);

 f := (u) -> simplify(u /~ sqrt(add(u[i]^2,i=1..3)));

 T["poles"] := table([
  2 = map(f,T["short_edge_centres"]),
  3 = map(i -> f(v[i]),T["vertices"]),
  4 = map(f,T["square_centres"])
 ]);

 T["pole_plots"] := table([
  2 = map(u -> line(1.6 *~ u, -1.6 *~ u,color=green,thickness=4),T["poles"][2]),
  3 = map(u -> line(1.9 *~ u, -1.9 *~ u,color=red  ,thickness=4),T["poles"][3]),
  4 = map(u -> line(1.2 *~ u, -1.2 *~ u,color=blue ,thickness=4),T["poles"][4])
 ]);

 T["all_poles_plot"] := display(
  op(T["pole_plots"][2]),op(T["pole_plots"][3]),op(T["pole_plots"][4]),
  scaling=constrained,axes=none
 );

 G0 := [];
 G1 := [[2,4,1,3,6,8,5,7],[1,3,5,7,2,4,6,8]];
 G := [[1,2,3,4,5,6,7,8]];

 while G <> G0 do
  G0 := G;
  G := sort([op({seq(seq(`o/permutations`(8)(g,h),g in G1),h in G)})]);
 od:

 T["rotation_group"] := G;

 T["symmetry_group"] :=
   sort([op(G),op(map(s -> `o/permutations`(8)([8,7,6,5,4,3,2,1],s),G))]);

 r := table([1 = 1, 2 = 2, 3 = 3, 4 = 4, 5 = 4, 6 = 3, 7 = 2, 8 = 1]):

 T["to_axis_perm"]    := table():
 T["of_axis_perm"]    := table():
 T["symmetry_matrix"] := table():

 for s in T["rotation_group"] do
  s1 := `o/permutations`(8)([8,7,6,5,4,3,2,1],s);
  t := [r[s[1]],r[s[2]],r[s[3]],r[s[4]]];
  T["to_axis_perm"][s]  := t;
  T["to_axis_perm"][s1] := t;
  T["of_axis_perm"][t]  := s;

  g := Matrix(3,3,[seq(x[i],i=1..9)]);
  eqs := map(op,[seq(convert(g . Vector(v[i]) - Vector(v[s[i]]),list),i = 1..8)]);
  sol := solve(eqs);
  T["symmetry_matrix"][s ] := subs(sol,convert( g,listlist));
  T["symmetry_matrix"][s1] := subs(sol,convert(-g,listlist));
 od:

 return eval(T);
end:


text:
type $SS <struct {
        @real f32,
        @imag f32}>
var $sconst1 <struct {
        @f11 i32,
        @f22 <$SS>,
        @f33 i32,
        @f44 i32,
        @f55 f32,
        @f66 f32}> = [ 1=1007, 2= [1=11.11f, 2=22.22f], 5=-273.0f, 6=75.0f, 4=6023 ]
func $fact () i32 {
  return (dread i32 $sconst1 1)
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: openmaple/MapleCompiler

max_issues_repo_name: venshine/OpenArkCompiler

max_forks_repo_name: venshine/OpenArkCompiler


text:
rename_cohomology_gens := proc(T0::table,u::list)
 local u0,n0,r,T,i;

 u0 := [op(T0["cohomology_gens"])];
 n0 := nops(u0);

 if nops(u) <> n0 then
  error("Incorect number of cohomology generators");
 fi;

 r := [seq(u0[i] = u[i],i=1..n0)];

 T := copy(T0);
 T["cohomology_gens"] := subs(r,T0["cohomology_gens"]);
 T["cohomology_rels"] := subs(r,T0["cohomology_rels"]);
 T["cohomology_degrees"] := table([]):
 for i from 1 to n0 do
  T["cohomology_degrees"][u[i]] :=
   T0["cohomology_degrees"][u0[i]];
 od;
 T["cohomology_basis"] := subs(r,T0["cohomology_basis"]);
 for i from 0 to T["dimension"] do
  T["cohomology_graded_basis"][i] :=
   subs(r,T["cohomology_graded_basis"][i]);
 od:

 return eval(T);
end:

make_product_space := proc(T0,T1)
 local T,d0,d1,B0,B1,B,i0,i1,u0,u1,x,i;
 T := table():

 T["name"] := sprintf("%s x %s",T0["name"],T1["name"]);
 d0 := T0["dimension"];
 d1 := T1["dimension"];
 T["dimension"] := d0+d1;
 T["cohomology_gens"] := tdeg(op(T0["cohomology_gens"]),
                              op(T1["cohomology_gens"]));
 T["cohomology_rels"] := [op(T0["cohomology_rels"]),
                          op(T1["cohomology_rels"])];
 T["cohomology_degrees"] := table([]);
 for x in T0["cohomology_gens"] do
  T["cohomology_degrees"][x] := T0["cohomology_degrees"][x];
 od:
 for x in T1["cohomology_gens"] do
  T["cohomology_degrees"][x] := T1["cohomology_degrees"][x];
 od:

 T["cohomology_basis"] := [seq(seq(u0*u1,u0 in T0["cohomology_basis"]),
                                         u1 in T1["cohomology_basis"])];

 B0 := T0["cohomology_graded_basis"]:
 B1 := T1["cohomology_graded_basis"]:
 B := table():

 for i from 0 to d0+d1 do
  B[i] := NULL:

  for i0 from max(0,i-d1) to min(i,d0) do
   i1 := i - i0;
   B[i] := B[i],seq(seq(u0*u1,u0 in B0[i0]),u1 in B1[i1]);
  od:
 od;

 T["cohomology_graded_basis"] := eval(B):

 return eval(T);
end:


text:
# Berechne kinematische Zwangsbedingungen für den picker2Dm1 aus 2 Parallelogrammen und einer Fünfgelenkkette
# Einleitung
#
# Herleitung der Kinematischen Zwangsbedingungen für einen 2-FG hybriden Palettier-/Verpackungsroboter mit zwei Parallelogrammen und einer Fünfgelenkkette
# Roboter: Triowin (Hersteller), "TDR Series 2-Axis Parallel Robot"
# Link: http://en.triowin.com/tdr-series2-axisparallelrobot-15268685491267769.html
#
# Die kinematischen Zwangsbedingungen werden als Ersetzungsausdruck für die abhängigen Winkel aufgestellt.
# picker2Dm2TE -> modellierung der Zwangsbedingungen mit Ausdrücken für trigonometrische Elimination

# Voraussetzung
# Vier- und Fünf gelenkkette muss mit der Toolbox berechnet worden sein (Arbeitsblätter "fourbar1TE_kinematic_constraints.mw" und fivebar1TE_kinematic_constraints.mw)
#
# Autor
# Abderahman Bejaoui
#  Hiwi bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2019-01
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
# Import
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
# Konstante Winkel bereits hineinschreiben
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Variablen definieren für die Hilfswinkel
#
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:

# Sicherung der aktuellen Daten
backup_kintmp_qs := kintmp_qs:
backup_kintmp_qt := kintmp_qt:
backup_kintmp_subsexp := kintmp_subsexp:
# Aufruf der Viergelenkkette (als Parallelogramm)
read sprintf("../codeexport/fourbar2TE/tmp/kinematic_constraints_maple_inert.m"):
kintmp_subsexp_fourbar:= kintmp_subsexp:
kintmp_qs_fourbar := kintmp_qs:
kintmp_qt_fourbar := kintmp_qt:
# Aufrud der Fünfgelenkkette
read sprintf("../codeexport/fivebar1TE/tmp/kinematic_constraints_maple_inert.m"):
kintmp_subsexp_fivebar := kintmp_subsexp:
kintmp_qs_fivebar := kintmp_qs:
kintmp_qt_fivebar := kintmp_qt:
# Vorherige Werte für dieses System wieder herstellen
kintmp_qs := backup_kintmp_qs:
kintmp_qt := backup_kintmp_qt:
kintmp_subsexp := backup_kintmp_subsexp:
# Fünfgelenkkettte A-B-C-E-D
winkel1:=<<sin_eta3|sin(eta3_s)>;<cos_eta3|cos(eta3_s)>;<sin_eta12|sin(eta12_s)>;<cos_eta12|cos(eta12_s)>;<sin_eta24|sin(eta24_s)>;<cos_eta24|cos(eta24_s)>>:

for i from 1 to 6 do
  # Daten der allgemeinen Fünfgelenkkette laden
  winkel1(i,2):=kintmp_subsexp_fivebar(i,2):
  # Kinematikparameter der FGK an die des Roboters anpassen
  winkel1(i,2) := subs({BC=L1}, winkel1(i,2)):
  winkel1(i,2) := subs({AB=e}, winkel1(i,2)):
  winkel1(i,2) := subs({CD=L3}, winkel1(i,2)):
  winkel1(i,2) := subs({ED=L4}, winkel1(i,2)):
  winkel1(i,2) := subs({AE=L3}, winkel1(i,2)):
  # variable qJ bezogen auf FGK durch allgemeine Variable ersetzen
  winkel1(i,2) := subs({sin(qJ1s)= sin_phi_s},winkel1(i,2)):
  winkel1(i,2) := subs({cos(qJ1s)= cos_phi_s},winkel1(i,2)):
  winkel1(i,2) := subs({sin(qJ2s)= sin_psi_s},winkel1(i,2)):
  winkel1(i,2) := subs({cos(qJ2s)= cos_psi_s},winkel1(i,2)):
  # Allgemeine Variable durch qJ bezogen auf Roboter ersetzen
  winkel1(i,2) := subs({sin_phi_s = cos(qJ2s)},winkel1(i,2)):
  winkel1(i,2) := subs({cos_phi_s = -sin(qJ2s)},winkel1(i,2)):
  winkel1(i,2) := subs({sin_psi_s = sin(qJ1s)},winkel1(i,2)):
  winkel1(i,2) := subs({cos_psi_s = cos(qJ1s)},winkel1(i,2)):
end do:
# Umrechnung aus Zwangsbedingugen für Winkel aus MDH-Tabelle
sin_eta3_s:=winkel1(1,2):
cos_eta3_s:=winkel1(2,2):
sin_eta12_s:=winkel1(3,2):
cos_eta12_s:=winkel1(4,2):
sin_eta24_s:=winkel1(5,2):
cos_eta24_s:=winkel1(6,2):
sin_eta410_s:=sin_eta3_s:
cos_eta410_s:=-cos_eta3_s:

# Viergelenkkette B-H-L-C
winkel2:=<<sin_rho18|sin(rho18_s)>;<cos_rho18|cos(rho18_s)>;<sin_rho05|sin(rho05_s)>;<cos_rho05|cos(rho05_s)>;<sin_rho811|sin(rho811_s)>;<cos_rho811|cos(rho811_s)>>:
for i from 1 to 6 do
  # Daten des Parallelogramms laden
  winkel2(i,2) := kintmp_subsexp_fourbar(i,2):
  # Kinematikparameter des Parallelogramms an die des Roboters anpassen
  winkel2(i,2) := subs({l1=L5}, winkel2(i,2)):
  winkel2(i,2) := subs({l2=L1}, winkel2(i,2)):
  # variable qJ bezogen auf Parallelogramm durch allgemeine Variable ersetzen
  winkel2(i,2) := subs({sin(qJ1s)= sin_phi_s},winkel2(i,2)):
  winkel2(i,2) := subs({cos(qJ1s)= cos_phi_s},winkel2(i,2)):
  # Allgemeine Variable durch qJ bezogen auf Roboter ersetzen
  winkel2(i,2) := subs({cos_phi_s = -cos(phi05)*cos(qJ1s)-sin(phi05)*sin(qJ1s)},winkel2(i,2)):
  winkel2(i,2) := subs({sin_phi_s = -cos(phi05)*sin(qJ1s)+sin(phi05)*cos(qJ1s)},winkel2(i,2)):
end do:
sin_rho18_s:=winkel2(1,2):
cos_rho18_s:=winkel2(2,2):
sin_rho05_s:=winkel2(3,2):
cos_rho05_s:=winkel2(4,2):
sin_rho811_s:=winkel2(5,2):
cos_rho811_s:=winkel2(6,2):
# Viergelenkkette C-G-F-P
winkel3:=<<sin_xi23|sin(xi23_s)>;<cos_xi23|cos(xi23_s)>;<sin_xi612|sin(xi612_s)>;<cos_xi612|cos(xi612_s)>;<sin_xi39|sin(xi39_s)>;<cos_xi39|cos(xi39_s)>>:
cos_xi26_s:= - ( -cos(phi1)*(cos_rho18_s*cos_eta12_s+sin_rho18_s*sin_eta12_s)-sin(phi1)*(sin_eta12_s*cos_rho18_s-cos_eta12_s*sin_rho18_s)):
sin_xi26_s:= -cos(phi1)*(sin_eta12_s*cos_rho18_s-cos_eta12_s*sin_rho18_s)+sin(phi1)*(cos_eta12_s*cos_rho18_s+sin_eta12_s*sin_rho18_s):
for i from 1 to 6 do
  # Daten des Parallelogramms laden
  winkel3(i,2) := kintmp_subsexp_fourbar(i,2):
  # Kinematikparameter des Parallelogramms an die des Roboters anpassen
  winkel3(i,2) := subs({l1=L6}, winkel3(i,2)):
  winkel3(i,2) := subs({l2=L2}, winkel3(i,2)):
  # variable qJ bezogen auf Parallelogramm durch allgemeine Variable ersetzen
  winkel3(i,2) := subs({sin(qJ1s) = sin_phi_s},winkel3(i,2)):
  winkel3(i,2) := subs({cos(qJ1s) = cos_phi_s},winkel3(i,2)):
  # Allgemeine Variable durch qJ bezogen auf Roboter ersetzen
  winkel3(i,2) := subs({cos_phi_s = -cos(phi1)*(cos_rho18_s*cos_eta12_s+sin_rho18_s*sin_eta12_s)-sin(phi1)*(sin_eta12_s*cos_rho18_s-cos_eta12_s*sin_rho18_s)},winkel3(i,2)):
  winkel3(i,2) := subs({sin_phi_s = -cos(phi1)*(sin_eta12_s*cos_rho18_s-cos_eta12_s*sin_rho18_s)+sin(phi1)*(cos_eta12_s*cos_rho18_s+sin_eta12_s*sin_rho18_s)},winkel3(i,2)):
end do:
winkel3(1..6,1):
sin_xi23_s:=winkel3(1,2):
cos_xi23_s:=winkel3(2,2):
sin_xi612_s:=winkel3(3,2): # schon in MDH
;
cos_xi612_s:=winkel3(4,2): # schon in MDH
;
sin_xi39_s:=winkel3(5,2):
cos_xi39_s:=winkel3(6,2):
# Kintmp_subsexp
kintmp_subsexp(1,2) := sin_eta12_s:
kintmp_subsexp(2,2) := cos_eta12_s:
kintmp_subsexp(3,2) := sin_xi23_s:
kintmp_subsexp(4,2) := cos_xi23_s:
kintmp_subsexp(5,2) := sin_eta24_s:
kintmp_subsexp(6,2) := cos_eta24_s:
kintmp_subsexp(7,2) := sin_rho05_s:
kintmp_subsexp(8,2) := cos_rho05_s:
kintmp_subsexp(9,2) := sin_xi26_s:
kintmp_subsexp(10,2) := cos_xi26_s:
kintmp_subsexp(11,2) := sin_rho18_s:
kintmp_subsexp(12,2) := cos_rho18_s:
kintmp_subsexp(13,2) := sin_xi39_s:
kintmp_subsexp(14,2) := cos_xi39_s:
kintmp_subsexp(15,2) := sin_eta410_s:
kintmp_subsexp(16,2) := cos_eta410_s:
kintmp_subsexp(17,2) := sin_rho811_s:
kintmp_subsexp(18,2) := cos_rho811_s:
kintmp_subsexp(19,2) := sin_xi612_s:
kintmp_subsexp(20,2) := cos_xi612_s:
# Kintmp_qs
kintmp_qs(1,1):=%arctan(sin_eta12_s,cos_eta12_s):
kintmp_qs(2,1):=%arctan(sin_xi23_s,cos_xi23_s):
kintmp_qs(3,1):=%arctan(sin_eta24_s,cos_eta24_s):
kintmp_qs(4,1):=%arctan(sin_rho05_s,cos_rho05_s):
kintmp_qs(5,1):=%arctan(sin_xi26_s,cos_xi26_s):
kintmp_qs(6,1):=%arctan(sin_rho18_s,cos_rho18_s):
kintmp_qs(7,1):=%arctan(sin_xi39_s,cos_xi39_s):
kintmp_qs(8,1):=%arctan(sin_eta410_s,cos_eta410_s):
kintmp_qs(9,1):=%arctan(sin_rho811_s,cos_rho811_s):
kintmp_qs(10,1):=%arctan(sin_xi612_s,cos_xi612_s):
# Export
kintmp_qt := convert_s_t(kintmp_qs):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d",robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert" ,robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


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/.
*)

# Minimum and maximum functions that I can differentiate

m_min := (x1, x2) -> my_piecewise3(x1 > x2, x2, x1):
m_max := (x1, x2) -> my_piecewise3(x1 > x2, x1, x2):
m_abs := (x)      -> my_piecewise3(x > 0, x, -x):

# Teach maple to differentiate piecewise functions

`diff/my_piecewise3` :=
    proc(c, x1, x2, x) my_piecewise3(c, diff(x1, x), diff(x2, x)) end proc:

`diff/my_piecewise5` :=
    proc(c1, x1, c2, x2, x3, x) my_piecewise5(c1, diff(x1, x), c2, diff(x2, x), diff(x3,x)) end proc:

# This is the derivative of xc_E1_scaled = -exp(x)*Ei(-x) = exp(x)E1(x)

`diff/xc_E1_scaled` :=
  proc(y, x) (xc_E1_scaled(y) - 1/y) * diff(y, x) end proc:

# The derivative of xc_erfcx = exp(x^2)*erfc(x) = 2*x*xc_erfcx(x) - 2/sqrt(Pi)
`diff/xc_erfcx` :=
  proc(y, x) (2*y*xc_erfcx(y) - 2/sqrt(Pi)) * diff(y, x) end proc:

# a series of useful definitions

M_C         := 137.0359996287515: (* speed of light *)

X2S         := 1/(2*(6*Pi^2)^(1/3)):
X2S_2D      := 1/(2*(4*Pi)^(1/2)):

X_FACTOR_C    := 3/8*(3/Pi)^(1/3)*4^(2/3):
X_FACTOR_2D_C := 8/(3*sqrt(Pi)):
K_FACTOR_C    := 3/10*(6*Pi^2)^(2/3):

MU_GE       := 10/81:
MU_PBE      := 0.06672455060314922*(Pi^2)/3:
KAPPA_PBE   := 0.8040:

# generic conversion functions

$ifdef xc_dimensions_1d
DIMENSIONS   := 1:
RS_FACTOR    := 1/2:
$elif xc_dimensions_2d
DIMENSIONS   := 2:
RS_FACTOR    := 1/sqrt(Pi):
LDA_X_FACTOR := -X_FACTOR_2D_C:
$else
DIMENSIONS   := 3:
RS_FACTOR    := (3/(4*Pi))^(1/3):
LDA_X_FACTOR := -X_FACTOR_C:
$endif

r_ws       := n  -> RS_FACTOR/n^(1/DIMENSIONS):
n_total    := rs -> (RS_FACTOR/rs)^DIMENSIONS:

n_spin     := (rs, z) -> simplify((1 + z)*n_total(rs)/2):
sigma_spin := (rs, z, xs) -> simplify(xs^2*n_spin(rs, z)^(8/3)):
t_total    := (z, ts0, ts1) ->
  (ts0*((1 + z)/2)^(5/3) + ts1*((1 - z)/2)^(5/3)):
u_total    := (z, us0, us1) -> t_total(z, us0, us1):

# useful formulas that enter several functionals follow

# von Weizsäcker term
t_vw := (z, xt, us0, us1) -> (xt^2 - u_total(z, us0, us1))/8:

# Screening for extreme values of zeta
z_thr := z -> my_piecewise5(
  1 + z <= p_a_zeta_threshold, p_a_zeta_threshold - 1,
  1 - z <= p_a_zeta_threshold, 1 - p_a_zeta_threshold,
  z):

opz_pow_n := (z, n) -> my_piecewise3(1 + z <= p_a_zeta_threshold, (p_a_zeta_threshold)^n, (1 + z)^n):

# See Eq. (9) of Perdew1992_13244
f_zeta    := z -> (opz_pow_n(z,4/3) + opz_pow_n(-z,4/3) - 2)/(2^(4/3) - 2):
f_zeta_2d := z -> 1/2*(opz_pow_n(z,3/2) + opz_pow_n(-z,3/2)):

# used in several correlation functionals
mphi := z -> (opz_pow_n(z,2/3) + opz_pow_n(-z,2/3))/2:
tt   := (rs, z, xt) -> xt/(4*2^(1/3)*mphi(z)*sqrt(rs)):

# in the paper it is beta_a = 0.066725
beta_Hu_Langreth := rs -> 0.066724550603149220*(1 + 0.1*rs)/(1 + 0.1778*rs):

# Generate exchange and kinetic functionals from the expression for the
# enhancement factor
lda_x_spin := (rs, z) -> simplify(LDA_X_FACTOR*opz_pow_n(z,1 + 1/DIMENSIONS)*2^(-1-1/DIMENSIONS)*(RS_FACTOR/rs)):
lda_k_spin := (rs, z) -> simplify(K_FACTOR_C*opz_pow_n(z,5/3)*2^(-5/3)*(RS_FACTOR/rs)^2):

# Screen out small densities as well as the zero component from fully spin polarized densities due to terms of the form (1+z)^{-n}
screen_dens := (rs, z) -> (n_spin(rs,  z) <= p_a_dens_threshold):

# For most functionals zeta screening occurs inside the funcitonal,
#  but for B88 and B94 correlation need to screen out also outside
screen_dens_zeta := (rs, z) -> screen_dens(rs, z) or (1 + z <= p_a_zeta_threshold):

# non-separable GGA exchange
gga_exchange_nsp := (func, rs, z, xs0, xs1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, lda_x_spin(rs, z_thr( z))*func(rs, z_thr( z), xs0))
  + my_piecewise3(screen_dens(rs, -z), 0, lda_x_spin(rs, z_thr(-z))*func(rs, z_thr(-z), xs1)):
# GGA exchange
gga_exchange := (func, rs, z, xs0, xs1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, lda_x_spin(rs, z_thr( z))*func(xs0))
  + my_piecewise3(screen_dens(rs, -z), 0, lda_x_spin(rs, z_thr(-z))*func(xs1)):
# GGA kinetic energy
gga_kinetic := (func, rs, z, xs0, xs1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, lda_k_spin(rs, z_thr( z))*func(xs0))
  + my_piecewise3(screen_dens(rs, -z), 0, lda_k_spin(rs, z_thr(-z))*func(xs1)):

# non-separable meta-GGA exchange
mgga_exchange_nsp := (func, rs, z, xs0, xs1, u0, u1, t0, t1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, lda_x_spin(rs, z_thr( z))*func(rs, z_thr( z), xs0, u0, t0))
  + my_piecewise3(screen_dens(rs, -z), 0, lda_x_spin(rs, z_thr(-z))*func(rs, z_thr(-z), xs1, u1, t1)):
# meta-GGA exchange
mgga_exchange := (func, rs, z, xs0, xs1, u0, u1, t0, t1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, lda_x_spin(rs, z_thr( z))*func(xs0, u0, t0))
  + my_piecewise3(screen_dens(rs, -z), 0, lda_x_spin(rs, z_thr(-z))*func(xs1, u1, t1)):
# meta-GGA kinetic energy
mgga_kinetic := (func, rs, z, xs0, xs1, u0, u1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, lda_k_spin(rs, z_thr( z))*func(xs0, u0))
  + my_piecewise3(screen_dens(rs, -z), 0, lda_k_spin(rs, z_thr(-z))*func(xs1, u1)):

# This is the Stoll decomposition in our language
lda_stoll_par  := (lda_func, rs, z) ->
  my_piecewise3(screen_dens_zeta(rs,  z), 0, opz_pow_n(z, 1)/2 * lda_func(rs*2^(1/3)*opz_pow_n(z, -1/3), 1)):

lda_stoll_perp := (lda_func, rs, z) ->
  + lda_func(rs, z)
  - lda_stoll_par(lda_func, rs,  z)
  - lda_stoll_par(lda_func, rs, -z):

gga_stoll_par  := (gga_func, rs, z, xs, spin) ->
  my_piecewise3(screen_dens_zeta(rs, z), 0,
    gga_func(rs*2^(1/3)*opz_pow_n(z,-1/3), spin, xs, xs*(1 + spin)/2, xs*(1 - spin)/2) * opz_pow_n(z,1)/2):

# Curvature of the Fermi hole (without the current term)
Fermi_D := (xs, ts) -> 1 - xs^2/(8*ts):

# correction to Fermi_D similar to the one found in
#  JCP 127, 214103 (2007); doi: http://dx.doi.org/10.1063/1.2800011
Fermi_D_corrected := (xs, ts) -> (1 - xs^2/(8*ts)) * (1 - exp(-4*ts^2/params_a_Fermi_D_cnst^2)):

# Becke function used in several correlation functionals
b88_R_F := (f_x, rs, z, xs) -> 1/(2*X_FACTOR_C*n_spin(rs, z)^(1/3)*f_x(xs)):
b88_zss := (css, f_x, rs, z, xs) -> 2*css*b88_R_F(f_x, rs, z, xs):
b88_zab := (cab, f_x, rs, z, xs0, xs1) -> cab*(
  + my_piecewise3(screen_dens(rs,  z), 0, b88_R_F(f_x, rs, z_thr( z), xs0))
  + my_piecewise3(screen_dens(rs, -z), 0, b88_R_F(f_x, rs, z_thr(-z), xs1))
  ):

# The meta-GGA version
b94_R_F := (f_x, rs, z, xs, us, ts) -> 1/(2*X_FACTOR_C*n_spin(rs, z)^(1/3)*f_x(xs,us,ts)):
b94_zss := (css, f_x, rs, z, xs, us, ts) -> 2*css*b94_R_F(f_x, rs, z, xs, us, ts):
b94_zab := (cab, f_x, rs, z, xs0, xs1, us0, us1, ts0, ts1) -> cab*(
  + my_piecewise3(screen_dens(rs,  z), 0, b94_R_F(f_x, rs, z_thr( z), xs0, us0, ts0))
  + my_piecewise3(screen_dens(rs, -z), 0, b94_R_F(f_x, rs, z_thr(-z), xs1, us1, ts1))
  ):

# Power series often used in mggas
mgga_w := t -> (K_FACTOR_C - t)/(K_FACTOR_C + t):
mgga_series_w := (a, n, t) -> add(a[i]*mgga_w(t)^(i-1), i=1..n):

# Used in screened functionals
kF := (rs, z) -> (3*Pi^2)^(1/3) * opz_pow_n(z,1/3) * RS_FACTOR/rs:
nu := (rs, z) -> p_a_hyb_omega_0_/kF(rs, z):

(* Maple 2020 has a bug where series aren't computed to the order requested; this circumvents that *)
padding_order := 30:

(* Cap polynomial expansion to given order (throws out the padded terms, if any) *)
throw_out_large_n := proc(X,n) select(t -> abs(degree(t,{b}))<=n, X); end proc:

(* Function that makes f(a) smooth for a->infty *)
enforce_smooth_lr := proc(f, a, a_cutoff, expansion_order);
  (* Calculate large-a expansion *)
  f_large := a -> eval(throw_out_large_n(convert(series(f(b), b=infinity, expansion_order+padding_order), polynom), expansion_order), b=a):
  (* Return the series expansion for large a; also remove any numerical overflows from the original branch  *)
  my_piecewise3(a >= a_cutoff, f_large(m_max(a, a_cutoff)), f(m_min(a, a_cutoff))):
end proc:


text:
Digits := 120:

interface(quiet=true):

read "common-procedures.mpl":
read "triple-double.mpl":

#---------------------------------------------------------------------
# Procedure gal computes a number x and a boolean such that
#
# if the boolean is true
# - |x - xstart| is minimal
# - x = roundX(x) (x is a floating point number in format roundX)
# - |(roundY(f(x)) - f(x)) / f(x)| <= eps
#
# and the boolean is false if no such x could be found in iter iterations
#
# Limitations:
#
# - floating point numbers xp such that f(xp) = 0 are not considered
#
# Uses:
#
# - The successor function succX for floating point numbers roundX(t)

gal:=proc(xstart,roundX,f,roundY,epsi,succX,iter)
local xleft, xright, i, err, ima, d, x, eps, b;
d := Digits;
Digits := 120;
eps := abs(epsi);
xleft := roundX(xstart);
xright := roundX(succX(roundX(xstart)));
x := 0;
i := 0;
err := eps * 2;
while ((err > eps) and (i < iter)) do
	i := i + 1;
	if (type(i,odd)) then
		x := xright;
		xright := evalf(roundX(succX(roundX(xright))));
	else
		x := xleft;
		xleft := evalf(roundX(-succX(roundX(-xleft))));
	end if;
	ima := evalf(f(x));
	if (ima <> 0) then
		err := evalf(abs((ima - roundY(ima))/ima));
	end if;
end do;
b := false;
if (i < iter) then b := true; end if;
Digits := d;
(x,b);
end proc:

galDoubleToDouble:=proc(xstart,f,epsi,iter)
local g,b,xn;
(g,b) := gal(xstart,nearest,f,nearest,epsi,succDouble,iter);
if (b) then xn := g; else xn = xstart; end if;
xn;
end proc:

galDoubleToDoubleDouble:=proc(xstart,f,epsi,iter)
local g,b,xn;
(g,b) := gal(xstart,nearest,f,nearestDD,epsi,succDouble,iter);
if (b) then xn := g; else xn = xstart; end if;
xn;
end proc:

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


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: mgga_exc *)
(* replace: "mbrxc_x\(" -> "xc_mgga_x_mbrxc_get_x(" *)

$include "mgga_x_mbrxc_bg.mpl"

task_alpha := (x, t) -> (t/K_FACTOR_C) * m_max(1 - x^2/(8*t), 1e-10):

mggac_b1 := 3.712:
mggac_b2 := 2.0:
mggac_b4 := 0.1:
mggac_b3 := 2.595 + 0.5197*mggac_b4 + 0.559*mggac_b2:
mggac_b5 := -3*mggac_b3:

(* new definition of Q. The rest of the functional remains the same *)
(* We have Lambda = (32 Pi^2)^(2/3)/(6 Q) *)
mbrxc_Q := (x, t) ->
      + (32*Pi)^(2/3)/6
      * (1 + mggac_b4*task_alpha(x, t) + mggac_b5*task_alpha(x, t)^2)
      / (mggac_b1 + mggac_b2*task_alpha(x, t) + mggac_b3*task_alpha(x, t)^2):


text:
#@ Not autoload

with(LinearAlgebra):

p := 2;
chromatic_n_max := 4;  # NB degree(v[4]) = 30, degree(v[5]) = 62
chromatic_s_max := 20;

unprotect('v','m','t');
unassign('v','m','t');

v[0] := p;
m[0] := 1;
t[0] := 1;

protect('v','m','t');

# On BP_* we use a lexicographic ordering with v[1] >> v[2] >> v[3] >> ... >> 1.
# This seems to be the best order except that we really want to treat v[0]
# as a variable with v[0] >> v[1], so any monomial multiplied by p comes
# after all the bare monomials.  Then the terms p^t * monomial in any
# given degree have order type omega, and the terms after any given term
# span an invariant ideal.

BP_vars   := plex(seq(v[i],i=1..chromatic_n_max));

# On the other hand, it might be better to use some kind of order that treats
# the powers v[n]^(p^m) in a special way.  This might import more information
# from the chromatic spectral sequence into the ordering.  One option would be
# like
# .. v[1]^(p^2) >> v[2]^(p^2) >> .. >> v[1]^p >> v[2]^p >> .. v[1] >> v[2] >> ...

# On BP_*BP we use a lexicographic ordering with
# v[1] >> v[2] >> .... >> t[1] >> t[2] >> ... >> 1
# This privileges the ideal J = (t[1],...,t[k]); note that BP_*BP/J is
# the Hopf algebroid whose Ext calculates \pi_*(T(k)).

# Again, it might be better to use some kind of order that treats
# the powers t[n]^(p^m) in a special way.  This might import more information
# from the May spectral sequence into the ordering.

BPBP_vars := plex(seq(v[i],i=1..chromatic_n_max),seq(t[i],i=1..chromatic_n_max));

# The cobar complex involves variables t[n,s] corresponding to the copy of
# t[n] in the s'th tensor factor.  We use the ordering with
# v[1] >> v[2] >> ... >> t[1,1] >> t[2,1] >> ... >> t[1,2] >> t[2,2] >> ...
# so the variables from each tensor factor dominate those from the next tensor
# factor.  It is not clear whether this is the best choice.
BP_cobar_vars :=
 plex(op(BPBP_vars),seq(seq(t[i,j],i=1..chromatic_n_max),j=1..chromatic_s_max));

BP_cmp   := (a,b) -> TestOrder(a,b,BP_vars):
BPBP_cmp := (a,b) -> TestOrder(a,b,BPBP_vars):
BP_cobar_cmp := (a,b) -> TestOrder(a,b,BP_cobar_vars):

BP_degree := (u) -> degree(subs(BP_degree_rule,u),e);

# Basis for Z[v[i] : i >= n] in degree d
BP_basis := proc(d::integer,n::posint := 1)
 option remember;
 local m,r,i;
 m := 2*p^n-2;
 r := floor(d/m);
 if d < 0 or modp(d,2*p-2) <> 0 then
  return [];
 elif d = 0 then
  return [1];
 elif r = 0 then
  return [];
 else
  map(op,[seq(v[n]^i *~ BP_basis(d-m*i,n+1),i=0..r)]);
 fi;
end:

# Basis for Z[m[i] : i >= n] in degree d
HBP_basis := proc(d::integer,n::posint := 1)
 eval(subs(v = m,BP_basis(d,n)));
end:

# Basis for Z[t[i] : i >= n] in degree d
T_basis := proc(d::integer,n::posint := 1)
 eval(subs(v = t,BP_basis(d,n)));
end:

# Basis for Z[v[i] : i >= n] in degree <= d
BP_lower_basis := proc(d::integer,n::posint := 1)
 option remember;
 local m,r,i;
 m := 2*p^n-2;
 r := floor(d/m);
 if d < 0 or modp(d,2*p-2) <> 0 then
  return [];
 elif r = 0 then
  return [1];
 else
  map(op,[seq(v[n]^i *~ BP_lower_basis(d-m*i,n+1),i=0..r)]);
 fi;
end:

# Basis for Z[m[i] : i >= n] in degree d
HBP_lower_basis := proc(d::integer,n::posint := 1)
 eval(subs(v = m,BP_lower_basis(d,n)));
end:

# Basis for Z[t[i] : i >= n] in degree d
T_lower_basis := proc(d::integer,n::posint := 1)
 eval(subs(v = t,BP_lower_basis(d,n)));
end:

# Basis for BP_*BP in degree d
BPBP_basis := proc(d::integer)
 local B,m,u,v;
 B := NULL;
 for u in BP_lower_basis(d) do
  m := d - BP_degree(u);
  for v in T_basis(m) do
   B := B,(u*v);
  od;
 od;
 return [B];
end:

# Basis for the s-fold tensor power of Z[t_1,t_2,...] in degree d
# The copy of t[i] in the j'th tensor factor is represented by t[i,j]
T_power_basis := proc(s::nonnegint,d::integer)
 local B,m,R1,R2,u,v;
 if s = 0 then return `if`(d = 0,[1],[]); fi;

 R1 := {seq(t[n] = t[n,1],n=1..chromatic_n_max)};
 R2 := {seq(seq(t[n,i] = t[n,i+1],i=1..s-1),n=1..chromatic_n_max)};
 B := NULL;

 for u in subs(R1,T_lower_basis(d)) do
  m := d - BP_degree(u);
  for v in subs(R2,T_power_basis(s-1,m)) do
   B := B,(u*v);
  od;
 od;
 return [B];
end:

# Basis for the s-fold tensor power of the augmentation ideal in
# Z[t_1,t_2,...] in degree d
T_reduced_power_basis := proc(s::nonnegint,d::integer)
 local B,m,R1,R2,u,v;
 if s = 0 then return `if`(d = 0,[1],[]); fi;

 R1 := {seq(t[n] = t[n,1],n=1..chromatic_n_max)};
 R2 := {seq(seq(t[n,i] = t[n,i+1],i=1..s-1),n=1..chromatic_n_max)};
 B := NULL;

 for u in subs(R1,T_lower_basis(d)) do
  if u <> 1 then
   m := d - BP_degree(u);
   for v in subs(R2,T_reduced_power_basis(s-1,m)) do
    B := B,(u*v);
   od;
  fi;
 od;
 return [B];
end:

BP_cobar_basis := (s,d) ->
 [seq(seq(seq(a*b,a in BP_basis(i)),b in T_reduced_power_basis(s,d-i)),i=0..d)];

# Hazewinkel generators in terms of log coefficients
vm := proc(n::nonnegint)
 option remember;
 if n = 0 then
  return p;
 else
  return expand(p*m[n] - add(m[k]*vm(n-k)^(p^k),k=1..n-1));
 fi;
end:

# Log coefficients in terms of Hazewinkel generators
# This is inefficient; should use Ravenel's formulae instead
mv := proc(n::nonnegint)
 local err;
 option remember;

 if n = 0 then return 1; fi;

 err := vm(n) - v[n];
 err := expand(subs({seq(m[i] = mv(i),i=1..n-1)},err));

 return rhs(solve(err=0,{m[n]})[1]);
end:

# Right unit map on the log coefficients
eta_m := (k) -> add(m[i] * t[k-i]^(p^i),i=0..k);

# Right unit map on the Hazewinkel generators
eta_v := proc(k::posint,s::nonnegint)
 option remember;
 local u,R;

 if nargs = 1 then
  u := expand(subs({seq(m[i] = eta_m(i),i=1..k)},vm(k)));
  u := expand(subs({seq(m[i] = mv(i),i=1..k)},u));
  return u;
 else
  if s = 0 then
   return v[k];
  elif s = 1 then
   u := expand(subs({seq(m[i] = eta_m(i),i=1..k)},vm(k)));
   u := expand(subs({seq(m[i] = mv(i),i=1..k)},u));
   u := subs({seq(t[i] = t[i,1],i=1..chromatic_n_max)},u);
   return u;
  else
   u := eta_v(k,s-1);
   R := { seq(seq(t[i,j] = t[i,j+1],i=1..chromatic_n_max),j=1..s-1),
          seq(v[i] = eta_v(i,1),i=1..chromatic_n_max) };
   u := subs(R,u);
   return u;
  fi;
 fi;
end:

# Hopf algebroid coproduct on the generators t[n]
psi_t := proc(n)
 local a,b;
 option remember;
 if n = 0 then return 1; fi;

 a := add(add(mv(i)*t[j,1]^(p^i)*t[n-i-j,2]^(p^(i+j)),j=0..n-i),i=0..n);
 b := add(mv(i)*psi_t(n-i)^(p^i),i=1..n);
 a := subs({t[0,1]=1,t[0,2]=1},a);
 b := subs({t[0,1]=1,t[0,2]=1},b);
 return expand(a - b);
end:

d_BP_cobar_rule := proc(s,i)
 local R0,R1,R2,R3;
 if i = 0 then
  R0 := {seq(t[j]=t[j,1],j=1..chromatic_n_max)};
  R1 := {seq(t[j]=t[j,2],j=1..chromatic_n_max)};
  R2 := {seq(seq(t[j,k]=t[j,k+1],k=1..s),j=1..chromatic_n_max)};
  return
   {seq(v[j] = expand(subs(R0,eta_v(j))),j=1..chromatic_n_max),op(R1),op(R2)};
 else
  R0 := {seq(t[j,1]=t[j,i],j=1..chromatic_n_max),
         seq(t[j,2]=t[j,i+1],j=1..chromatic_n_max),
	 seq(v[j] = eta_v(j,i-1),j=1..chromatic_n_max)};
  R1 := {seq(seq(t[j,k]=t[j,k+1],k=i+1..s),j=1..chromatic_n_max)};
  R2 := {seq(t[j,i] = subs(R0,psi_t(j)),j=1..chromatic_n_max)};
  if i = 1 then
   R3 := {seq(t[j] = expand(subs(R0,psi_t(j))),j=1..chromatic_n_max)};
  else
   R3 := {seq(t[j] = t[j,1],j=1..chromatic_n_max)};
  fi;
  return {op(R1),op(R2),op(R3)};
 fi;
end:

d_BP_cobar := (s) -> (u) ->
 expand(add((-1)^i * subs(d_BP_cobar_rule(s,i),u),i=0..s+1));

d_BP_cobar_matrix := proc(s,d)
 local B1,B2,cf;
 B1 := BP_cobar_basis(s,d);
 B2 := BP_cobar_basis(s+1,d);
 cf := proc(u)
  local sol;
  sol := solve({coeffs(u - add(c[i]*B2[i],i=1..nops(B2)),indets(B2))});
  subs(sol,[seq(c[i],i=1..nops(B2))]);
 end;
 Transpose(Matrix(map(cf,map(d_BP_cobar(s),B1))));
end:

BP_cobar_cycles := proc(s,d)
 local M;
 M := d_BP_cobar_matrix(s,d);
 map(u -> Transpose(Vector(BP_cobar_basis(s,d))) . u,NullSpace(M));
end:

mu_BP_cobar := (s1,s2) -> proc(a,b)
 local i,b0;
 b0 := b;
 for i from 0 to s1-1 do
  b0 := expand(subs(d_BP_cobar_rule(s2+i,0),b0));
 od:
 return expand(a * b0);
end:

analyse_BP := proc(p_,n_)
 global p,chromatic_n_max;

 if nargs > 0 then p := p_; fi;
 if nargs > 1 then chromatic_n_max := n_; fi;

 vm(chromatic_n_max);
 mv(chromatic_n_max);
 eta_v(chromatic_n_max);
 psi_t(chromatic_n_max);
 NULL;
end:

analyse_BP_cobar := proc(s::nonnegint,d::integer)
 local i,R,R0,B1,B2,n1,n2,M,L0,P0,x,x1,L,P,Q,nx,LB2,y,HB,HE,T,U,L1,K,Ki;
 global BP_cobar_data;

 if d + s >= 2*(p^(chromatic_n_max + 1) - 1) then
  error("chromatic_n_max is too small");
 fi;

 if s > chromatic_s_max then
  error("chromatic_s_max is too small");
 fi;

 R := table():

 if s = 0 then
  if d = 0 then
   R["chain_basis"]     := [1];
   R["chain_rank"]      := 1;
   R["cycle_basis"]     := [1];
   R["cycle_rank"]      := 1;
   R["boundary_basis"]  := [];
   R["boundary_rank"]   := 0;
   R["pivot_data"]      := [[1,infinity]];
   R["non_cycle_basis"] := [];
   R["homology_basis"]  := [1];
   R["homology_exponents"] := [infinity];
  else
   R["chain_basis"]     := BP_cobar_basis(0,d);
   R["chain_rank"]      := nops(R["chain_basis"]);
   R["cycle_basis"]     := [];
   R["cycle_rank"]      := 0;
   R["boundary_basis"]  := [];
   R["boundary_rank"]   := 0;
   R["pivot_data"]      := [];
   R["non_cycle_basis"] := R["chain_basis"];
   R["homology_basis"]  := [];
   R["homology_exponents"] := [];
  fi;

  BP_cobar_data[s,d] := eval(R);
  return eval(R);
 fi;

 if type(BP_cobar_data[s-1,d+1],table) then
  R0 := BP_cobar_data[s-1,d+1];
  B1 := R0["non_cycle_basis"];
 else
  B1 := BP_cobar_basis(s-1,s+d):
 fi;

 B2 := BP_cobar_basis(s,s+d):
 n1 := nops(B1);
 n2 := nops(B2);
 R["chain_basis"] := B2;
 R["chain_rank"]  := n2;

 if n2 = 0 then
  R["cycle_basis"]     := [];
  R["cycle_rank"]      := 0;
  R["boundary_basis"]  := [];
  R["boundary_rank"]   := 0;
  R["pivot_data"]      := [];
  R["non_cycle_basis"] := [];
  R["homology_basis"]  := [];
  R["homology_exponents"] := [];
  BP_cobar_data[s,d] := eval(R);
  return eval(R);
 fi;

 M := Transpose(Matrix(map(coeff_list,map(d_BP_cobar(s-1),B1),B2)));
 L0,P0,x := op(Zpl_reduce(Transpose(M),p)):
 L := Transpose(L0):
 P := Transpose(P0):
 nx := nops(x):
 R["cycle_rank"] := nx;

 LB2 := convert(L0.Vector(B2),list):
 R["boundary_basis"] := [seq(LB2[i],i=1..nx)];
 R["cycle_basis"]    := [seq(LB2[i]/p^x[i][2],i=1..nx)];
 R["pivot_data"]     := x;

 y := sort([op({seq(i,i=1..n2)} minus {seq(x[i][1],i=1..nx)})]);
 R["non_cycle_basis"] := [seq(B2[i],i in y)];

 HB := [];
 HE := [];

 for i from 1 to nx do
  if x[i][2] > 0 then
   HB := [op(HB),R["cycle_basis"][i]];
   HE := [op(HE),p^x[i][2]];
  fi;
 od;

 R["homology_basis"] := HB;
 R["homology_exponents"] := HE;

 T := <IdentityMatrix(nx)|Matrix(nx,n1-nx)>;
 U := Matrix(nx,n2):
 for i from 1 to nx do U[i,x[i][1]] := 1; od:
 L1 := U.L.Transpose(T);
 Q := P.Transpose(T).(1/L1).U;
 K := <SubMatrix(L,1..n2,1..nx)|Matrix(n2,n2-nx)>;
 y := sort([op({seq(i,i=1..n2)} minus {seq(x[i][1],i=1..nx)})]);
 for i from 1 to n2 - nx do
  K[y[i],nx+i] := 1;
 od:
 Ki := 1/K;

 BP_cobar_data[s,d] := eval(R);
 return eval(R);
end:

BP_set_p := proc(p_)
 global p,v,BP_degree_rule,BP_basis,BP_lower_basis,vm,mv,eta_v,psi_t;

 p := p_;
 unprotect('v');
 v[0] := p;
 protect('v');

 BP_degree_rule := {
  seq(v[n] = e^(2*(p^n-1)) * v[n],n=1..chromatic_n_max),
  seq(t[n] = e^(2*(p^n-1)) * t[n],n=1..chromatic_n_max),
  seq(seq(t[n,i] = e^(2*(p^n-1)) * t[n,i],n=1..chromatic_n_max),i=1..chromatic_s_max)
 }:

 forget(BP_basis);
 forget(BP_lower_basis);
 forget(vm);
 forget(mv);
 forget(eta_v);
 forget(psi_t);
 NULL;
end:

save_BP_data := proc()
 local file;
 file := sprintf("%s/BP_%d.m",data_dir,p);
 save(p,chromatic_n_max,chromatic_s_max,vm,mv,eta_v,psi_t,BP_basis,BP_lower_basis,file);
end:

load_BP_data := proc(p)
 local file;
 file := sprintf("%s/BP_%d.m",data_dir,p);
 load(file);
end:

save_BP_cobar_data := proc()
 local file;
 file := sprintf("%s/BP_cobar_data_%d.m",data_dir,p);
 save(BP_cobar_data,file);
end:

load_BP_cobar_data := proc(p)
 local file;
 file := sprintf("%s/BP_cobar_data_%d.m",data_dir,p);
 load(file);
end:


text:
# https://matplotlib.org/users/customizing.html

#figure.titlesize : small        # size of the figure title (Figure.suptitle())
figure.dpi       : 120          # figure dots per inch

lines.linewidth : 2
lines.markersize : 10

legend.fontsize : 14
legend.handlelength : 1

axes.titlesize : 17
axes.labelsize : 17
axes.grid : False   # display grid or not --> Doesn't work. Just type "g"

xtick.direction : in
ytick.direction : in
xtick.labelsize : 14 # font size of tick labels
ytick.labelsize : 14
xtick.major.size     : 6      # major tick size in points
xtick.minor.size     : 4      # minor tick size in points
ytick.major.size     : 6
ytick.minor.size     : 4
xtick.minor.visible  : True   # turn minor ticks on by default
ytick.minor.visible  : True

# has to be a .ttf file in /Library/Fonts
# font.family      : sans-serif
# font.sans-serif : Helvetica Neue
#font.sans-serif : 'Comic Sans'
#font.style : italic
# font.serif : Optima
# font.family : Helvetica Neue Light
# font.style : light
#font.family : Optima

max_stars_repo_name: sweigart/pygama

max_issues_repo_name: sweigart/pygama

max_forks_repo_name: sweigart/pygama


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 *)
(* prefix:
  mgga_x_ms_params *params;

  assert(pt->params != NULL);
  params = (mgga_x_ms_params * ) (pt->params);
*)

fa := a -> (1 - a^2)^3 / (1 + a^3 + params_a_b*a^6):
f0 := (p, c) -> 1 + params_a_kappa*(1 - params_a_kappa/(params_a_kappa + MU_GE*p + c)):

f := (rs, x, t, u) -> f0(X2S^2*x^2, 0) + \
  fa((t - x^2/8)/K_FACTOR_C)*(f0(X2S^2*x^2, params_a_c) - f0(X2S^2*x^2, 0)):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
program parseformalparameters;
procedure func(a,b : integer; c,d : char);
begin
    a := b;
    c := d;
end;
var x1,x2 : integer;
var y1,y2 : char;
begin
    call func(x1,x2,y1,y2);
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


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 xc_dimensions_2d

_2d_b86_mgc_beta  := 0.003317:
_2d_b86_mgc_gamma := 0.008323:

_2d_b86_mgc_f := x -> 1 + _2d_b86_mgc_beta/X_FACTOR_2D_C*x^2/(1 + _2d_b86_mgc_gamma*x^2)^(3/4):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange(_2d_b86_mgc_f, rs, zeta, xs0, xs1):


text:
(*
 Copyright (C) 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_c_chachiyo_params *params;

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

(* Functional is based on Chachiyo correlation with modified spin scaling *)
$include "lda_c_chachiyo_mod.mpl"

(* Reduced gradient parameter *)
cha_t := (rs, xt) -> (Pi/3)^(1/6) / 4 * n_total(rs)^(1/6) * xt:

(* The full functional that agrees with the given reference values is *)
f_chachiyo_gga := (rs, z, xt, xs0, xs1) -> f_chachiyo(rs, z) * (1 + cha_t(rs, xt)^2)^(params_a_h / f_chachiyo(rs, z)):

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


text:
with(Threads:-Task):
printf("root exec\n");
Start(print,1,2,3);
printf("single task\n");
Start(null,Task=[print,1,2,3]);
printf("multi tasks\n");
Start(null,Task=[print,1],Task=[print,2],Task=[print,3]);
print("multi tasks with same function\n");
Start(null,Tasks=[print,[1],[2],[3]]);
print("complex parameters\n");
Start(passed,1,Task=[`+`,1,2,3],Tasks=[`*`,[2,3,3],[4,5,6]],x,y,z);
printf("执行完毕\n");


text:
with(Threads):
id:=Create(int(sin(x)^x,x),res);
res;
Wait(id);
res;


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 *)
(* prefix:
  mgga_x_scan_params *params;

  assert(pt->params != NULL);
  params = (mgga_x_scan_params * )(pt->params);
*)

p     := x -> X2S^2*x^2:
alpha := (x, t) -> (t - x^2/8)/K_FACTOR_C:

f_alpha := a -> piecewise(a <= 1, exp(-params_a_c1*a/(1 - a)), -params_a_d*exp(params_a_c2/(1 - a))):

h1x := x -> 1 + params_a_k1*(1 - params_a_k1/(params_a_k1 + x)):

b2 := sqrt(5913/405000):
b1 := (511/13500)/(2*b2):
b3 := 1/2:
b4 := MU_GE^2/params_a_k1 - 1606/18225 - b1^2:
y  := (x, a) -> MU_GE*p(x) + b4*p(x)^2*exp(-b4*p(x)/MU_GE) + (b1*p(x) + b2*(1 - a)*exp(-b3*(1 - a)^2))^2:

a1 := 4.9479:
gx := x -> 1 - exp(-a1/sqrt(X2S*x)):

h0x := 1.174:
f   := (rs, x, t, u) -> (h1x(y(x, alpha(x, t)))*(1 - f_alpha(alpha(x, t))) + h0x*f_alpha(alpha(x, t)))*gx(x):


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

`eta/height_functions`:= (A::set) -> `if`(nops(A)=1,table(A=0),FAIL);

######################################################################

`gamma/height_functions`:= (A::set,B::set) -> (p) -> proc(j,h)
 local k,P,phi,b0,b1,m,U,pU,b;

 k := table();
 P := `list_elements/nonempty_subsets`(A);

 phi := table();
 for b0 in B do
  m := 1;
  for b1 in B do
   if b1 <> b0 then m := min(m,j[{b0,b1}]); fi;
  od;
  phi[b0] := m;
 od;

 for U in P do
  pU := map(u -> p[u],U);
  if nops(pU) > 1 then
   k[U] := j[pU];
  else
   b := op(pU);
   k[U] := phi[b] * h[b][U];
  fi;
 od;

 return eval(k);
end;

######################################################################

`gamma_min/height_functions`:= (A::set,B::set) -> (p) -> proc(j,h)
 local k,P,phi,b0,b1,m,U,pU,b;

 k := table();
 P := `list_elements/nonempty_subsets`(A);

 phi := table();
 for b0 in B do
  m := 1;
  for b1 in B do
   if b1 <> b0 then m := min(m,j[{b0,b1}]); fi;
  od;
  phi[b0] := m;
 od;

 for U in P do
  pU := map(u -> p[u],U);
  if nops(pU) > 1 then
   k[U] := j[pU];
  else
   b := op(pU);
   k[U] := min(phi[b] , h[b][U]);
  fi;
 od;

 return eval(k);
end;


text:
`is_element/shuffles` := (n::nonnegint,m::nonnegint) -> proc(s)
 local i;

 if not(`is_element/permutations`(n+m)(s)) then
  return false;
 fi;

 for i from 1 to n+m-1 do
  if i <> n and s[i+1] < s[i] then
   return false;
  fi;
 od:

 return true;
end:

`is_equal/shuffles` := (n,m) -> (s,t) -> evalb(s = t);

`is_leq/shuffles` := NULL;

`from_subset/shuffles` := (n::nonnegint,m::nonnegint) -> proc(A)
 option remember;
 local X,B,i;

 X := {seq(i,i=1..n+m)};
 if A minus X <> {} then return FAIL; fi;
 if nops(A) <> n then return FAIL; fi;

 B := X minus A;
 return [op(sort(A)),op(sort(B))];
end:

`to_subset/shuffles` := (n::nonnegint,m::nonnegint) -> proc(s)
 local i;
 {seq(s[i],i=1..n)};
end:

`to_grid_path/shuffles` := (n::nonnegint,m::nonnegint) -> proc(s)
 local t,si,i;

 t := table():
 t[0] := [0,0];
 si := `inv/permutations`(n+m)(s);
 for i from 1 to n+m do
  t[i] := `if`(si[i] <= n,[1,0],[0,1]) +~ t[i-1];
 od:

 return eval(t);
end:

`from_grid_path/shuffles` := (n::nonnegint,m::nonnegint) -> proc(t)
 local A,i;

 A := select(i -> t[i][1] > t[i-1][1],{seq(i,i=1..n+m)});
 return `from_subset/shuffles`(n,m)(A);
end:

`random_element/shuffles` := (n::nonnegint,m::nonnegint) -> proc()
 local u,A,i;
 u := combinat[randperm](n+m);
 A := {seq(u[i],i=1..n)};
 return `from_subset/shuffles`(n,m)(A);
end:

`list_elements/shuffles` := proc(n::nonnegint,m::nonnegint)
 local AA;
 AA := map(A -> {op(A)},combinat[choose](n+m,n));
 return map(`from_subset/shuffles`(n,m),AA);
end:

`list_elements/inverse_shuffles` := proc(n::nonnegint,m::nonnegint)
 option remember;

 map(`inv/permutations`(n+m),`list_elements/shuffles`(n,m));
end:

`count_elements/shuffles` := (n::nonnegint,m::nonnegint) -> binomial(n+m,n);

`sgn/shuffles` := (n::nonnegint,m::nonnegint) -> proc(s)
 local i,j;
 return signum(mul(mul(s[j]-s[i],j=i+1..n+m),i=1..n+m-1));
end:


text:
$Amt:10.0
$maxAmt:120.0
$multipler:1.002
$additional:0.03
$interval:15.0

main()inline{

	once{
		RemoveImages(Self());
		SetImage(Self(), "main", "CustomEmitter");
		SetTextSize(0.6);
		SetDigitalis(CurrentCoords(), 1);
		SetUnitAttribute(Self(), CONST_NULLIFIERDAMAGES, False);
	}

	if(GetDigitalis(CurrentCoords())!=1) Amt*=0.85;

	SetDigitalisGrowth(CurrentCoords(), True);
	SetDigitalis(CurrentCoords(), 1);

	if(Amt<5)Destroy(Self(), 2);

	SetPopupText("Amt: " $ round(Amt, 2) $ endl() $ "Interval: " $ round(interval/30.0, 2));
	SetImageScale(Self(), "main", dup(Amt/maxAmt*1.7+0.3));
	AddCreeper(CurrentCoords(), Amt);

	if(Amt>=50)once enable = True;

	if(enable){
		[x, y] = RandCoordsInRange(CurrentCoords(), 10);
		noEnemy=True;
		do(GetEnemyUnitsInRange(x, y, 8)){
			pop(); // remove uid from stack
			noEnemy=False;
		}
		if(GetCreeper(x, y)>1 && noEnemy){
			Amt/=1.75;
			uid = CreateUnit("CRPLCORE", x, y);
			AddScriptToUnit(uid, "emitter.crpl");
			if(Amt<maxAmt/2)enable = False;
		}
	}

	if(Amt < maxAmt) Amt = Amt*multipler + additional;
	else once Amt = maxAmt;

	Delay(interval);
}


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 *)

$include "lda_c_rc04.mpl"

msigma := 1.43:
malpha := 2.30:

Bs := s -> 1/(1 + msigma*s^malpha):

f_tcs := (rs, z, xt) -> f_rc04(rs, z)*Bs(X2S*2^(1/3)*xt):

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


text:
#var $const1 <[13] i8> [ 104, 101, 108, 108, 111, 32, 119, 111, 114, 108, 100, 10, 0 ]
var $fconst1 <[18] f64> = [ 1007, 707, -273, 75, 0113, 0x4b, 75u, 75l, 75ul, 75lu, 3.1425926, 6.02e23, 6.02e+23, 1.6e-19, 3.0, 3.14159, 6.02e23, .233 ]
func $printf (var %p1 <* i8>)void
func $main ( ) i32 {
   call &printf (addrof a32 $fconst1)
   call &printf (conststr ptr "\"hello maple\"", conststr a32 "\nworld\n")
   # 你好: E4 BD A0 E5 A5 BD  中国: E4 B8 AD E5 9B BD
   call &printf (conststr16 ptr "\xE4\xBD\xA0\xE5\xA5\xBD", conststr16 a32 "\xE4\xB8\xAD\xE5\x9B\xBD")
   return (constval i32 0) }

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
with(plots):
with(plottools):

great_arc := proc(u,v,hue)
 local w,r,s,t,d,i;
 w := [u[2]*v[3]-u[3]*v[2],u[3]*v[1]-u[1]*v[3],u[1]*v[2]-u[2]*v[1]];
 r := evalf(sqrt(u[1]^2+u[2]^2+u[3]^2));
 d := u[1]*v[1]+u[2]*v[2]+u[3]*v[3];
 w := [v[1] - d * u[1]/r^2,v[2] - d * u[2]/r^2,v[3] - d * u[3]/r^2];
 s := evalf(sqrt(w[1]^2+w[2]^2+w[3]^2));
 w := [r*w[1]/s,r*w[2]/s,r*w[3]/s];
 t := evalf(arccos(d/r^2))/10;
 return(CURVES(
        [seq([cos(i*t)*u[1]+sin(i*t)*w[1],
              cos(i*t)*u[2]+sin(i*t)*w[2],
              cos(i*t)*u[3]+sin(i*t)*w[3]],i=0..10)],
        COLOR(HUE,hue)));
end:


text:
 func $foo (
  var %i i32
#var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     bnot i32(dread i32 %i))}

 func $foo1 (
  var %i i32, var %j i32, var %k i32,
  var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     bnot i32(dread i32 %i))}

 func $foo2 (
  var %i i32, var %j i32, var %k i32
  ) i32 {
   return (
     bnot i32(constval i32 0x111111111))}

 func $foo3 (
  var %i i64, var %j i64, var %k i32
  ) i64 {
   return (
     bnot i64(dread i64 %i))}

 func $foo5 (
  var %i i64, var %j i64, var %k i32
  ) i64 {
   return (
     bnot i64(constval i64 0x11111))}

# todo float bnot
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# Apply a domain to an expression.
# Apply will perform certain 'simplifications' to make sure the
# domain application is well-formed. e.g.,
#   DOMAIN( DBound( [x,y], [xt,yt] ), DConstrain(..) )
#      is the same as
#   DOMAIN(        ..               , DInto(x,xt,DInto(y,yt,DConstrain(..))) )
# basically, when we see an 'unbound' variable in the 'RHS' , we should bind
# it with the default 'DInto'.
export Apply := module ()
       export ModuleApply :=
         proc(dom :: HDomain_mb, kb0 :: t_kb
             ,f_into := "default"
             ,f_body := "default", $)
           local vs, sh, ctx;
           if dom :: DomNoSol then
               error "cannot apply %1", dom;
           end if;
           vs, sh := op(dom);
           vs := Domain:-Bound:-withVarsIxs(vs);
           ctx := [ kb0,
              `if`(f_into="default",`do_mk`,f_into),
              `if`(f_body="default",`do_body`,f_body) ];
           (e->do_apply({}, e, vs, sh, ctx));
       end proc;

       export Shape := proc(dsh :: DomShape, e, $) ModuleApply( DOMAIN( [], dsh ), e ) end proc;
       export Bound := proc(dbn :: DomBound, e, $) ModuleApply( DOMAIN( dbn, DConstrain() ), e ) end proc;

       # main loop of the apply function
       # done_  := variables already integrated
       # e      := expr
       # vs     := domain bounds
       # sh     := domain shape
       # ctx[0] := KB (context for f_into/f_body)
       # f_into := how to make a DomInto
       # f_body := how to make the expression body
       local do_apply := proc(done__, e, vs, sh_, ctx, $)
           local sh := sh_, done_ := done__
               , r, cond, cond_outer, vn, vt, shv, vars, deps, ctx1, rn;
           # If the solution is a constraint, and the constraint is true,
           # then just produce the expression. If it isn't necessarily true
           # (i.e. trivial) then produce a Partition with the constraint as a
           # guard.
           if sh :: DomConstrain then
               vars := Domain:-Bound:-varsOf(vs,"set") minus done_;
               cond, cond_outer := select_cond_outer(sh, vars);
               # if there are still integrals which have not been applied, apply
               # them now
               do_constrain(cond_outer)(do_mks(e, (r,kb1) -> do_constrain(cond)(op(3,ctx)(r, kb1)), vars, vs, ctx))
           # if the solution is a sum of solutions, produce the algebraic sum
           # of each summand of the solution applied to the expression.
           elif sh :: DomSum then
               `+`(seq(do_apply(done_, e, vs, s, ctx), s=sh))
           # if the solution is a split solution, just make `do_apply' over
           # the values of the Partition (the subsolutions)
           elif sh :: DomSplit then
               Partition:-Amap(
                 [ (c,_)->c
                 , (p,c)->do_apply(done_, e, vs, p, applyop(kb->KB:-assert(c,kb),1,ctx))
                 , c->c ], op(1, sh) );
           # performs the 'make' on the expression after recursively
           # applying the solution
           elif sh :: DomInto then
               # deconstruction
               vn, vt, shv := op(sh);
               # extract bounds which have not been applied upon which this
               # bound depends. Those are applied after this one, so are not
               # in 'scope' in the recursive call
               vars := Domain:-Bound:-varsOf(vs,"set");
               deps := (indets(vt, DomBoundVar) intersect vars) minus done_ ;
               deps := `union`(deps, {vn});
               done_ := `union`(done_, deps) ;

               shv, cond_outer := select_cond_outer(shv, done_);

               # build this integral, and the other this one depended on, then
               # recursively apply
               do_constrain(cond_outer)(
                 do_mks(e, (r,kb1) -> do_apply(done_, r, vs, shv, subsop(1=kb1,ctx)),
                      deps, Domain:-Bound:-upd(vs, vn, vt), ctx));
           else
               error "don't know how to apply %1", sh
           end if;
       end proc;

       export do_mk := proc(mk, e, vn, ty_, _kb)
           Domain:-ExtBound[mk]:-DoMk(e, vn, ty_);
       end proc;

       export do_body := proc(e, _kb) e end proc;

       local do_constrain := proc(cond0::{list,set,DomConstrain},$)
         local mk_cond, cond := cond0;
         if nops(cond)=0 then
           mk_cond := x->x;
         else
           cond := bool_And(op(cond));
           mk_cond := x->PARTITION([Piece(cond,x), Piece(Not(cond),0)]);
         end if;
         mk_cond;
       end proc;

       # A very cute way of pulling constraints out over other constructors
       # which relies on the fact that constraints are required to be innermost;
       # i.e. DomShape is a sum of products; i.e. this could also be called
       # 'factoring'.
       local select_cond_outer := proc(sh::DomShape, vars0::set({name,list(name)}), $)
         local csd, cs0, cs, ots, sh1, ins, vars := vars0;
         vars := map(v->`if`(v::list,v,[v])[],vars);
         csd := [op(indets(sh, And(DomConstrain,Not(satisfies(x->has(x,vars))))))];
         cs0 := map(x->{op(x)},csd);
         if cs0=[] then return sh, {}; end if;
         cs := map(x->[selectremove(z->andmap(c->z in c,cs0),x)], cs0);
         ots := op([1,1],cs); ins := map(curry(op,2), cs);
         sh1 := subs(zip(`=`, csd, map(DConstrain@op,ins)), sh);
         userinfo(3, Domain, printf("select_cond_outer(%a, %a) = %a, %a\n", sh, vars, sh1, ots));
         sh1, ots;
       end proc;

       # Move into an integral by augmenting the KB with the variables bound by
       # that integral.
       local into_mk := proc(dbnd::DomBound, vn, vt, ctx, $)
           local kb, v_i, v_k, vn_, kb1, rn;
           kb := op(1,ctx); if kb :: t_not_a_kb then return KB:-NotAKB() end if;

           v_i := Domain:-Bound:-varIx(dbnd, vn);
           v_k := op([1,v_i,3], dbnd);
           vn_, kb1 := Domain:-ExtBound[v_k]:-MakeKB(vn, Domain:-ExtBound[v_k]:-SplitRange(vt), kb);

           rn := `if`(evalb(vn=vn_), [], [vn=vn_]);
           subsop(1=kb1, ctx), rn;
       end proc;

       local do_mks := proc(e, kont, todo::({list,set})(DomBoundVar), dbnd :: DomBound, ctx, $)
           local vs, v_td, i, vt, kb0, v_mk, _, ctx1, rn, r := e;

           if todo :: set then
             vs := select(v->v in todo,Domain:-Bound:-varsOf(dbnd));
             vs := ListTools[Reverse](vs);
           else
             vs := todo;
           end if;

           if nops(vs)=0 then
             return kont(r,op(1,ctx));
           end if;

           v_td := op(1,vs);
           vs   := subsop(1=NULL,vs);
           i    := Domain:-Bound:-varIx(dbnd, v_td);
           _,vt,v_mk := op(op([1,i], dbnd));

           ctx1, rn := into_mk(dbnd, v_td, vt, ctx);
           if rn <> [] then
             r := subs(rn,r);
           end if;

           r := do_mks(r, kont, vs, dbnd, ctx1);
           op(2,ctx)(v_mk, r, v_td, vt, op(1,ctx));
       end proc;
end module;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
# This module forms part of NewSLO and is `$include`d there

RoundTrip := proc(e, t::t_type)
  local result;
  interface(screenwidth=9999, prettyprint=0, warnlevel=0,
    showassumed=0,quiet=true);
  kernelopts(assertlevel=0);
  result := eval(ToInert(Simplify(_passed)), _Inert_ATTRIBUTE=NULL);
  lprint(result)
end proc;

Simplify := proc(e, t::t_type, {ctx :: list := []}, $)
  subsindets(SimplifyKB(value(e), t, build_kb(ctx, "Simplify")),
             And({t_sum, t_product}, anyfunc(anything, anything=range)),
             e -> subsop(0 = `if`(e::t_sum, SumIE, ProductIE),
                         applyop(`+`, [2,2,2], e, 1)))
end proc;

SimplifyKB_ := proc(e, t::t_type, kb::t_kb, $)
  local patterns, x, kb1, ex;
  if t :: HMeasure(anything) then
    %fromLO(%improve(%toLO(e), _ctx=kb), _ctx=kb);
  elif t :: HFunction(anything, anything) then
    patterns := htype_patterns(op(1,t));
    if patterns :: Branches(Branch(PVar(name),anything)) then
      # Eta-expand the function type
      x := `if`(e::lam(name,anything,anything), op(1,e),
                op([1,1,1],patterns));
      x, kb1 := genType(x, op(1,t), kb, e);
      ex := app(e,x);
      lam(x, op(1,t), SimplifyKB_(ex, op(2,t), kb1))
    else
      # Eta-expand the function type and the sum-of-product argument-type
      x := `if`(e::lam(name,anything,anything), op(1,e), d);
      if depends([e,t,kb], x) then x := gensym(x) end if;
      ex := app(e,x);
      lam(x, op(1,t), 'case'(x,
        map(proc(branch)
            local eSubst, pSubst, p1, binds, y, kb1, i, pSubst1;
              eSubst, pSubst := pattern_match([x,e], x, op(1,branch));
              p1 := subs(pSubst, op(1,branch));
              binds := [pattern_binds(p1)];
              kb1 := kb;
              pSubst1 := table();
              for i from 1 to nops(binds) do
                y, kb1 := genType(op(i,binds), op([2,i],branch), kb1);
                pSubst1[op(i,binds)] := y;
              end do;
              pSubst1 := op(op(pSubst1));
              Branch(subs(pSubst1, p1),
                     SimplifyKB_(eval(eval(ex,eSubst),pSubst1), op(2,t), kb1))
            end proc,
            patterns)))
    end if
  else
    %simplify_assuming(e, kb)
  end if
end proc;

SimplifyKB := proc(e, t::t_type, kb::t_kb, $)
    eval(SimplifyKB_(args), [%fromLO=fromLO,%improve=improve,%toLO=toLO,%simplify_assuming=simplify_assuming]);
end proc;

# Testing

TestSimplify := proc(m, t, n::algebraic:=m, {verify:=simplify})
  local s, r;
  # How to pass keyword argument {ctx::list:=[]} on to Simplify?
  s, r := selectremove(type, [_rest], 'identical(ctx)=anything');
  CodeTools[Test](Simplify(m,t,op(s)), n, measure(verify), op(r))
end proc;

TestHakaru := proc(m, n::{set(algebraic),algebraic}:=m,
                   {simp:=improve, verify:=simplify, ctx::list:=[]})
  local kb := build_kb(ctx, "TestHakaru"), ver := measure(verify);
  ver := `if`(n::set, 'member'(ver), ver);

  CodeTools[Test](fromLO(simp(toLO(m), _ctx=kb), _ctx=kb), n,
    ver, _rest)
end proc;

TestDisint := module()
    export ModuleApply := proc(
      M::{t_Hakaru, list(anything)}, #args to disint, or just 1st arg
      n::set({t_Hakaru, identical(NULL)}), #desired return
      ctx::t_kb_atoms:= [], #context: assumptions, "knowledge"
      TLim::{positive, identical(-1)}:= 80 #timelimit
    )
      local disint_args, disint_var, expected := n;
      if M :: list then
        disint_args := [op(M),ctx];
      else
        disint_var := gensym('t');
        disint_args := [M,disint_var,ctx];
        expected := subs(:-`t`=disint_var,expected);
      end if;
      try
        do_test(disint_args, copy(expected), TLim, _rest);
      catch "time expired":
        error "Time expired while running: disint(%1)", disint_args;
      end try;
    end proc;

    # This is necessary because CodeTools seems to forgot the value
    # of our local variables (but only the `expected result' arg)
    # unless that arg is precisely a formal parameter, and we `copy' the
    # input to this function.
    local do_test := proc(disint_args, expected, tlim)
      timelimit( tlim, CodeTools[Test]
                 ( {disint(op(disint_args))}
                 , expected
                 , '`subset`(measure(simplify))'
                 , _rest));
    end proc;
end module;

TestEfficient := proc(e, t::t_type, t_kb := KB:-empty, {label::string := "Test" })
  local done_, result, todo;
  todo := SimplifyKB_(args[1..3]);
  done_ := eval(todo, [%fromLO=fromLO,%improve=improve,%toLO=toLO,%simplify_assuming=simplify_assuming]);

  _Env_TestTools_Try_printf := false;
  result := TestTools[Try](label, Efficient(done_),true);
  if result = NULL then
    printf("%s passed.\n", label);
  else
    error sprintf("%s FAILED.\n"
                  "The result of\n\t%a\n\tis not efficient.\n"
                 , label, todo );
  end if;
end proc;

# Test roughly for "efficient" Hakaru measure terms,
# i.e., those we want simplification to produce.
Efficient := proc(mm, $)
  local m, n;
  m := mm;
  if has(m, 'undefined') then
    return false;
  end if;
  while m :: '{lam(name, anything, anything),
               Context(anything, anything),
               case(anything, Branches(Branch(anything, anything))),
               And(specfunc(piecewise)
                  ,{anyfunc(anything, anything, Msum())
                   ,anyfunc(anything, anything, anything, Msum())})}' do
    m := op(`if`(op(0,m)='lam',3,`if`(op(0,m)='case',[2,1,2],2)),m);
  end do;
  if m :: 'Weight(anything, anything)' then m := op(2,m) end if;
  if has(m, '{infinity, Lebesgue, int, Int, Beta, GAMMA}') then
    return false;
  end if;
  for n in `if`(m :: 'specfunc(Msum)', m, [m]) do
    if n :: 'Weight(anything, anything)' then n := op(2,n) end if;
    if has(subsindets(n, 'specfunc(Weight(anything, anything), Msum)',
                      s -> `if`(Testzero(`+`(op(map2(op, 1, s))) - 1),
                                map2(op, 2, s), s)),
           '{Msum, Weight}') then
      return false;
    end if;
  end do;
  return true;
end proc;

# Load a file in concrete Hakaru syntax (using the "momiji" command)
# and return its term (in which Sum and Product are inert) and type.
Concrete := proc(path::string, $)
  local cmd, res, dangerous_chars;
  cmd := FileTools:-AbsolutePath(path,
    (FileTools:-ParentDirectory@@2)(LibraryTools:-FindLibrary(Hakaru)));
  dangerous_chars := [ " ", "'", """", `if`(kernelopts(platform)="windows", [], ["\\"])[] ];
  if ormap((c->StringTools:-Has(cmd,c)), dangerous_chars) then
    error "Dangerous characters in path: %1", cmd;
  end if;
  cmd := cat("momiji ", cmd);
  res := ssystem(cmd);
  if res :: [0, string] then
    parse(cat("use Hakaru in ", op(2,res), " end use"));
  else
    error "ssystem %1: %2", cmd, res;
  end if;
end proc;


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 *)

$include "gga_x_ityh.mpl"

ityh_attenuation := a  -> attenuation_yukawa(a):


text:
# This checks some formulae in the LaTeX file Out_SL2.tex

check_aut_SL2 := proc()
 local i2,alpha,beta,beta0,gamma,mu,sigma,rho,
  m,n,g,h,a,b,c,d,t,u,v,rels,
  is_in_U,is_in_T,is_in_B,err;

 i2 := IdentityMatrix(2);
 alpha := (v) -> <<1|v>,<0|1>>;
 mu := (u) -> <<u|0>,<0|1/u>>;

 is_in_U := (g) -> g[2,1] = 0 and g[1,2] = 0 and g[1,1] * g[2,2] = 1;
 is_in_T := (g) -> g[2,1] = 0 and g[1,1] = 1 and g[2,2] = 1;
 is_in_B := (g) -> g[2,1] = 0 and g[1,1] * g[2,2] = 1;
 rho := (u) -> g[1,1];
 sigma := (a) -> <<0|1>,<1|a>>;
 beta := (a,b) -> alpha(b) . sigma(a) . (1/alpha(b)) . (1/sigma(a));
 beta0 := (b) -> <<1+b^2|b>,<b|1>>;
 gamma := (b,u) -> <<1/u|0>,<(1+b^2)/u|u>>;

 _ASSERT(
  Equal(map(modp,simplify(beta(a,b)),2),beta0(b)),
  "Formula for beta(b)"
 ):

 err := gamma(b,u).beta0(b).(1/gamma(b,u)) - sigma(b^2);
 err := map(modp,simplify(subs(b = u^2,err)),2);
 _ASSERT(
  Equal(err,<<0|0>,<0|0>>),
  "beta(b) is conjugate to sigma(b^2)"
 ):

 m := <<1/u|0>,<0|u>>;
 n := <<u|u>,<0|1/u>>;
 _ASSERT(
  Equal(m.n,alpha(1)),
  "m n = alpha(1)"
 );

 g := <<a|b>,<c|d>>;
 g := subs(solve(map(op,convert(g.alpha(1) - alpha(1).g,listlist)),{c,d}),g);
 _ASSERT(
  Equal(g,<<a|b>,<0|a>>),
  "Centraliser of alpha(1)"
 );

 g := <<u|v>,<0|1/u>>;
 h := alpha(v/(1/u-u));
 _ASSERT(
  Equal(simplify(h.mu(u).(1/h)),g),
  "D_u is conjugate to mu(u) in B"
 );

 g := <<a|b>,<c|d>>;
 rels := map(numer,map(op,convert(mu(u).g.mu(u).g - i2,listlist)));
 rels := map(coeffs,expand(rels),u);
 g := subs(solve(rels,{a,c,d}),g);
 _ASSERT(
  Equal(g,<<0|b>,<1/b|0>>),
  "(mu(u) . g)^2 = 1 only on U"
 );

 _ASSERT(
  Equal(mu(u) . alpha(t) . (1/mu(u)),alpha(u^2*t)),
  "Conjugate of alpha(t) by mu(u)"
 );
end:


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: mgga_exc *)
(* prefix:
  mgga_c_m05_params *params;

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

$define lda_c_pw_params
$define lda_c_pw_modified_params
$include "lda_c_pw.mpl"

$include "b97.mpl"

(* The parallel and perpendicular components of the energy *)
m05_comp := (rs, z, spin, xs, t) ->
  + lda_stoll_par(f_pw, rs,  z,  1)
  * b97_g(params_a_gamma_ss, params_a_css, xs)
  * Fermi_D_corrected(xs, t):

m05_fpar  := (rs, z, xs0, xs1, t0, t1) ->
  + m05_comp(rs,  z,  1, xs0, t0)
  + m05_comp(rs, -z, -1, xs1, t1):

m05_fperp := (rs, z, xs0, xs1, t0, t1) ->
  + lda_stoll_perp(f_pw, rs,  z)
  * b97_g(params_a_gamma_ab, params_a_cab, sqrt(xs0^2 + xs1^2)):

m05_f := (rs, z, xs0, xs1, t0, t1) ->
  + m05_fpar (rs, z, xs0, xs1, t0, t1)
  + m05_fperp(rs, z, xs0, xs1, t0, t1):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  m05_f(rs, z, xs0, xs1, t0, t1):


text:
type $ss1 <struct {@b <!T>,
                       @x <* f32>}>
type $ss2 <struct {@b1 <struct {@b <!SS>,
                                    @x <* f32>}>,
                       @pa <* <struct {@b f64,
                                       @x <* f32>}>>,
                       @x1 <[10] f64>,
                       @dx <!U>}>
type $ss3 <[] <!AAA>>
func &foo (var %x <!W>) i32 {
  var %y <!W>
  dassign %y (dread agg %x)
  return (dread i32 %y 1)
}
var $g <$ss1{!T=i32}>
var $h <$ss2{!U=i32,!SS=<$ss1{!T=f32}>}>
func &bar () i32 {
  var %i i32
  callinstant &foo<{!W=<$ss2{!U=f32,!SS=i32}>}>( dread agg $g )
  callinstantassigned &foo<{!W=<$ss1{!T=i32}>}>( dread agg $g ) { dassign %i 0 }
  virtualcallinstant &foo<{!W=<$ss2{!U=f32,!SS=i32}>}>( dread agg $g )
  virtualcallinstantassigned &foo<{!W=<$ss1{!T=i32}>}>( dread agg $g ) { dassign %i 0 }
  superclasscallinstant &foo<{!W=<$ss2{!U=f32,!SS=i32}>}>( dread agg $g )
  superclasscallinstantassigned &foo<{!W=<$ss1{!T=i32}>}>( dread agg $g ) { dassign %i 0 }
  interfacecallinstant &foo<{!W=<$ss2{!U=f32,!SS=i32}>}>( dread agg $g )
  interfacecallinstantassigned &foo<{!W=<$ss1{!T=i32}>}>( dread agg $g ) { dassign %i 0 }
  return ( dread i32 %i )
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


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_x *)

mbeta  := 0.003317:
mgamma := 0.008323:

f := x -> 1 + mbeta/X_FACTOR_C*x^2/(1 + mgamma*x^2)^(3/4):


text:
# Energy Calculation for the Robot based on MDH frames
# Introduction
# Berechnung von potentieller und kinetischer Energie für den Roboter.
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb -> floating base wird durch base twist (Geschwindigkeit der Basis) oder vollständige Orientierung (Euler-Winkel) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# energy -> Berechnung der Energie
# worldframe -> Berechnung der Positionen und Geschwindigkeiten im Welt-KS (KS W) anstatt im Basis-KS (KS 0)
# par1 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ)
# Authors
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
#
# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I (WS 2014/15)
# [Ott2008] Cartesian Impedance Control of Redundant and Flexible-Joint Robots
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_opt := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # Laden zur Einstellung der Term-Vereinfachung
;
# Einstellungen für Term-Vereinfachungen: 0=keine, 1=E_pot, 2=auch E_kin (einzelne), 3=auch E_kin (komplett)
if not assigned(simplify_options) or simplify_options(6)=-1 then
  # Die Vereinfachung ist ein Kompromiss aus Laufzeit des Programms und Güte der Ergebnisse
  # Bei zu komplizierten Termen stürzt Maple bei der Vereinfachung unvorhersehbar ab.
  # Die folgenden Regeln sind konservative Einstellungen, damit die Code-Generierung nicht zu lange dauert.
  if not kin_constraints_exist then # ohne Zwangsbedingungen. Annahme: Bringt dann sowieso nicht so viel
    use_simplify := 0:
  else # für Systeme mit Zwangsbedingungen
    if NJ <= 7 then
      use_simplify := 3: # standardmäßig bei einfachen Systemen alle simplify-Befehle anwenden
    elif NJ <= 11 then # das ist nur eine grobe Abschätzung
      use_simplify := 2: # Simplify reduzieren (Annahme: Optimierung Gesamt-Energie dauer länger)
    else
      use_simplify := 1: # simplify-Befehle weiter begrenzen (sonst hängt Maple sich auf)
    end if:
  end if:
else
  use_simplify := simplify_options(6): # sechster Eintrag ist für Energie
end if:

# Ergebnisse der Kinematik laden (Rotationsmatrizen, Schwerpunktskoordinaten, Geschwindigkeiten)
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf := Trf:
Trf_c := Trf_c:
read sprintf("../codeexport/%s/tmp/kinematics_com_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
r_W_W_Si := r_W_W_Si:
r_W_i_Si := r_W_i_Si:
read sprintf("../codeexport/%s/tmp/velocity_worldframe_floatbase_%s_par1_maple.m", robot_name, base_method_name):
omega_W_i := omega_W_i:
rD_W_i := rD_W_i:
rD_W_Si := rD_W_Si:

# Potential Energy
U_b := Matrix(NL, 1):
U_grav := 0:
for i to NL do
  U_b[i] := -M[i, 1]*Multiply(Transpose(g_world), r_W_W_Si(1 .. 3, i)):
  printf("%s. Potentielle Energie aus Gravitation für Körper %d berechnet (im Welt-KS, Parametersatz 1).\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), i-1): #0=Basis
  # Terme vereinfachen (nur einzelne Energien und nicht die Summe.
  # Annahme: Zusammenfassung verschiedener Körper nicht zielführend
  if use_simplify>=1 then
    tmp_t1:=time(): tmp_l1 := length(U_b[i,1]):
    printf("%s. Beginne Vereinfachung für potentielle Energie (Param. 1; Körper %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1):
    U_b[i,1] := simplify2(U_b[i,1]):
    tmp_t2:=time(): tmp_l2 := length(U_b[i,1]):
    printf("%s. Terme für potentielle Energie (Param. 1) vereinfacht (Körper %d). Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
  U_grav := U_grav+U_b[i, 1]:
  # Terme vereinfachen
  if use_simplify>=1 then
    tmp_t1:=time(): tmp_l1 := length(U_grav):
    printf("%s. Beginne Vereinfachung für potentielle Energie (Param. 1; Summe bis Körper %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1):
    U_grav := simplify2(U_grav):
    tmp_t2:=time(): tmp_l2 := length(U_grav):
    printf("%s. Terme für potentielle Energie vereinfacht (Param. 1; Summe bis Körper %d). Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
end do:
# Maple Export
save U_grav, sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_worldframe_par1_maple.m", robot_name, base_method_name):
# Kinetic Energy
# Use velocities in base frame and SX,SY,SZ as a parameter
# Berechne kinetische Energie aller durch Gelenkwinkel bewegter Körper.
T_b := Matrix(NL, 1):
T := 0:
for i to NL do
  # Trägheitstensor (3x3) um den Körperschwerpunkt in Körper-KS
  I_i_Si_Tensor := Matrix([[I_i_Si[1, i], I_i_Si[2, i], I_i_Si[3, i]], [I_i_Si[2, i], I_i_Si[4, i], I_i_Si[5, i]], [I_i_Si[3, i], I_i_Si[5, i], I_i_Si[6, i]]]):
  # Rotation vom Welt-KS zum jeweiligen Körper-KS (Die Basis ist dabei laufende Nummer 1)
  R_W_i := Matrix(Trf_c(1 .. 3, 1 .. 3, i)):
  # Trägheitstensor (3x3) um den Körperschwerpunkt in Basis-KS
  I_i_W := Multiply(R_W_i, Multiply(I_i_Si_Tensor, Transpose(R_W_i))):
  T_b[i] := (1/2)*M[i, 1]*Multiply(Transpose(rD_W_Si(1 .. 3, i)), rD_W_Si(1 .. 3, i))+(1/2)*Multiply(Transpose(omega_W_i(1 .. 3, i)), Multiply(I_i_W, omega_W_i(1 .. 3, i))):
  printf("%s. Kinetische Energie für Körper %d berechnet (im Welt-KS, Parametersatz 1)\n", FormatTime("%Y-%m-%d %H:%M:%S"), i-1):#0=Basis
  # Terme vereinfachen (nur einzelne Energien und nicht die Summe.
  # Annahme: Zusammenfassung verschiedener Körper nicht zielführend
  if use_simplify>=2 then
    tmp_t1:=time(): tmp_l1 := length(T_b[i,1]):
    printf("%s. Beginne Vereinfachung für kinetische Energie (Param. 1; Körper %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1):
    T_b[i,1] := simplify2(T_b[i,1]):
    tmp_t2:=time(): tmp_l2 := length(T_b[i,1]):
    printf("%s. Terme für kinetische Energie (Param. 1) vereinfacht (Körper %d). Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
  T := T+T_b[i, 1]:

  # Terme vereinfachen (optional, Standardmäßig nur kinetische Energie einzelner Körper optimieren)
  # Für komplizierte Systeme hängt sich Maple hier auf.
  if use_simplify>=3 then
    tmp_t1:=time(): tmp_l1 := length(T):
    printf("%s. Beginne Vereinfachung für kinetische Energie (Param. 1; Summe bis Körper %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1):
    T := simplify2(T):
    tmp_t2:=time(): tmp_l2 := length(T):
    printf("%s. Terme für kinetische Energie (Param. 1; Summe bis Körper %d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i-1, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
end do:
# Maple Export
save T, sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_worldframe_par1_maple.m", robot_name, base_method_name):
# Matlab Export
printf("%s: Energie (Param. 1) berechnet und gespeichert. Beginne Matlab-Export.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S")):
# Potential Energy
# Floating Base
U_s:=convert_t_s(U_grav):
if codegen_act then
  MatlabExport(U_s, sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_worldframe_par1_matlab.m", robot_name, base_method_name), codegen_opt):
end if:
# Fixed Base
U_s_fixb:=U_s:
for i from 1 to NQB do
  U_s_fixb := subs({X_base_s[i,1]=0},U_s_fixb):
end do:
if codegen_act then
  MatlabExport(U_s_fixb, sprintf("../codeexport/%s/tmp/energy_potential_fixb_worldframe_par1_matlab.m", robot_name), codegen_opt):
end if:
# Maple Export (Fixed-Base)
save U_s_fixb, sprintf("../codeexport/%s/tmp/energy_potential_fixb_worldframe_par1_maple.m", robot_name):
# Kinetic Energy
# Floating Base
# Eliminiere die Basis-Position (nicht: Orientierung). Die kinetische Energie kann nicht davon abhängen.
# Falls die Ausdrücke doch vorkommen, liegt es daran, dass Maple diese nicht wegoptimieren konnte.
T_s:=convert_t_s(T):
for i from 1 to 3 do
  T_s := subs({X_base_s[i,1]=0},T_s):
end do:
if codegen_act then
  MatlabExport(T_s, sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_worldframe_par1_matlab.m", robot_name, base_method_name), codegen_opt):
end if:
# Fixed Base
T_s_fixb:=T_s:
for i from 1 to NQB do
  T_s_fixb := subs({X_base_s[i,1]=0},T_s_fixb):
end do:
for i from 1 to 6 do
  T_s_fixb := subs({V_base_s[i,1]=0},T_s_fixb):
end do:
if codegen_act then
  MatlabExport(T_s_fixb, sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_worldframe_par1_matlab.m", robot_name), codegen_opt):
end if:
# Maple Export (Fixed-Base)
save T_s_fixb, sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_worldframe_par1_maple.m", robot_name):

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

`is_element/SCP` := (N::posint) -> (A::set) -> proc(Q)
 local i,U,V;

 global reason;

 if not `is_element/ACP`(N)(A)(Q) then
  reason := [convert(procname,string),"Q in ACP(N)(A)",reason];
  return false;
 fi;

 if nops(Q[N]) = nops(A)^2 then
  reason := [convert(procname,string),"Q[N] has only one block",Q[N]];
  return false;
 fi;

 return true;
end:

`is_equal/SCP` := (N::posint) -> (A::set) -> proc(Q1,Q2)
 global reason;

 if Q1 <> Q2 then
  reason := [convert(procname,string),"Q1 <> Q2",Q1,Q2];
  return false;
 fi;

 return true;
end:

`is_leq/SCP` := (N::posint) -> (A::set) -> proc(Q1,Q2)
 local i;

 for i from 1 to N do
  if Q2[i] minus Q1[i] <> {} then
   return false;
  fi;
 od;

 return true;
end:

`list_elements/SCP` := (N::posint) -> proc(A::set)
 local X,n;
 n := nops(A);
 X := `list_elements/ACP`(N)(A);
 X := select(Q -> nops(Q[N]) < n^2,X);
 return X;
end:

`count_elements/SCP` := (N::posint) -> proc(A::set)
 local d;
 return add(Stirling2(nops(A),d)*d!*N^(d-1),d=2..nops(A));
end:

`random_element/SCP` := (N::posint) -> (A::set) -> proc()
 local i,n,pi,Q,R,S,B,C,ok;

 if nops(A) <= 1 then
  return FAIL;
 fi;

 ok := false;

 while not ok do
  Q := `random_element/ACP`(N)(A)();
  if nops(Q[N]) < nops(A)^2 then
   ok := true;
  fi;
 od;

 return Q;
end:

`res/SCP` := (N::posint) -> (A::set,B::set) -> proc(Q)
 return map(`intersect`,Q,`top/autorel`(B));
end:

######################################################################

`gamma/SCP` := (N::posint) -> (A::set) -> proc(Q)
 local i;
 [seq(`op/autorel`(A)(Q[i]) minus Q[i],i=1..N)];
end:

######################################################################

`mu/SW/SCP` := (N::posint) -> (A::set) -> proc(x)
 return `mu/W/ACP`(N)(A)(x);
end:

`sigma/SCP/SW` := (N::posint) -> (A::set) -> proc(Q)
 return `bottom_normalise/SW`(N)(A)(`sigma/ACP/W`(N)(A)(Q));
end;

######################################################################

`describe/SCP` := eval(`describe/ACP`):


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: mgga_exc *)

zlp_c := 0.828432   *RS_FACTOR:
zlp_d := 2.15509e-2 *RS_FACTOR:
zlp_k := 2.047107e-3*RS_FACTOR:

zlp_f := (rs, z, xt, us0, us1) ->
  - (zlp_c + zlp_d*t_vw(z, xt, us0, us1))
  * (1 - zlp_k*log(1 + rs/zlp_k)/rs)/rs:

f := (rs, z, xt, xs0, xs1, us0, us1, ts0, ts1) ->
  zlp_f(rs, z, xt, us0, us1):


text:
fd := fopen("nonlinfit.txt", APPEND):
with(plots):
with(Statistics):
with(LinearAlgebra):

dta := ImportMatrix("/tmp/dta.csv"):
len := RowDimension(dta):
x := dta[1 .. len, 1]; y := dta[1 .. len, 2]:
#pointplot([x, y]);
N:=15;
tmp := NonlinearFit(c*exp(a*ln(binomial(N, v))-(b-a)*ln(factorial(v))), x, y, v, output = [parametervalues, residualmeansquare], initialvalues = [a = 1, b = 1, c = 0.1e-5]):
print(tmp);
fprintf(fd,"%f , %f, %f,%f\n", eval(a, tmp[1]), eval(b, tmp[1]), eval(c, tmp[1]), tmp[2]):
fclose(fd):

max_stars_repo_name: jkobject/CodonUsageBias

max_issues_repo_name: jkobject/CodonUsageBias

max_forks_repo_name: jkobject/CodonUsageBias


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[x5, x6, x4, x2, x1, x3], [y1, y2] ], derivations = [t]):
sys := [
diff(x5(t), t) - (k5*x6(t) + k4*x6(t) - k6*x5(t)*x3(t)),
diff(x6(t), t) - (-k5*x6(t) - k4*x6(t) + k6*x5(t)*x3(t)),
diff(x4(t), t) - (-k3*x4(t) - k2*x4(t) + k1*x1(t)*x2(t)),
diff(x2(t), t) - (k3*x4(t) + k2*x4(t) + k1*x1(t)*x2(t)),
diff(x1(t), t) - (k4*x6(t) + k2*x4(t) - k1*x1(t)*x2(t)),
diff(x3(t), t) - (k5*x6(t) + k3*x4(t) - k6*x5(t)*x3(t)),
y1(t) - (x3(t)),
y2(t) - (x2(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));