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

`is_element/SEM` := (N::posint) -> (A::set) -> proc(eta)
 local AA,a,b,c,ok;
 global reason;

 if not(type(eta,table)) then
  reason := [convert(procname,string),"eta is not a table",eta];
  return false;
 fi;

 AA := [seq(seq([a,b],b in A),a in A)];
 if {indices(eta)} <> {op(AA)} then
  reason := [convert(procname,string),"eta is not indexed by A^2",eta,AA];
  return false;
 fi;

 for a in A do
  for b in A do
   if eta[a,b] = 0 and a <> b then
    reason := [convert(procname,string),"eta does not separate a and b",a,b,eta];
    return false;
   fi;

   if eta[a,b] + eta[b,a] <> 0 then
    reason := [convert(procname,string),"eta is not antisymmetric at [a,b]",a,b,eta];
    return false;
   fi;

   for c in A do
    ok := `is_leq/E`(N)(0,eta[a,c]) or
          `is_leq/E`(N)(eta[a,b],eta[a,c]) or
          `is_leq/E`(N)(eta[b,c],eta[a,c]);
    if not(ok) then
     reason := [convert(procname,string),"min(0,eta[a,b],eta[a,c]) > eta[a,c]",a,b,c,eta];
     return false;
    fi;

    ok := `is_leq/E`(N)(eta[a,c],0) or
          `is_leq/E`(N)(eta[a,c],eta[a,b]) or
          `is_leq/E`(N)(eta[a,c],eta[b,c]);
    if not(ok) then
     reason := [convert(procname,string),"eta[a,c] > max(0,eta[a,b],eta[a,c])",a,b,c,eta];
     return false;
    fi;
   od;
  od;
 od;

 return true;
end;

`is_equal/SEM` := (N::posint) -> (A::set) -> proc(eta,zeta)
 local a,b;
 global reason;

 for a in A do
  for b in A do
   if eta[a,b] <> zeta[a,b] then
    reason := [convert(procname,string),"eta[a,b] <> zeta[a,b]",a,b,eta[a,b],zeta[a,b]];
    return false;
   fi;
  od;
 od;

 return true;
end:

`is_leq/SEM` := (N::posint) -> (A::set) -> proc(eta,zeta)
 local a,b;

 for a in A do
  for b in A do
   if not(`preceq/E`(N)(A)(eta[a,b],zeta[a,b])) then
    return false;
   fi;
  od;
 od;

 return true;
end:

`build/SEM` := (N::posint) -> (A::set) -> proc(s::list,u::list)
 local n,L,i,j,k,m,eta,a,b;

 n := nops(A);
 if nops(s) <> n or nops(u) <> n-1 then return FAIL; fi;
 if {op(s)} <> A then return FAIL; fi;

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

 eta := table;

 for a in A do
  eta[a,a] := 0;
  for b in A do
   i := L[a];
   j := L[b];
   if i < j then
    m := min(seq(u[k],k=i..j-1));
    eta[a,b] :=  epsilon^m;
    eta[b,a] := -epsilon^m;
   fi;
  od;
 od;

 return eval(eta);
end:

`random_element/SEM` := (N::posint) -> (A::set) -> proc()
 local s,m,u,i;

 s := `random_element/ord`(A)();
 m := nops(A);
 u := [seq(rand(0..N-1)(),i=1..m-1)];
 return `build/SEM`(N)(A)(s,u);
end:

`list_elements/SEM` := (N::posint) -> proc(A::set)
 option remember;
 local S,U,s,u,m,i,j;

 S := `list_elements/ord`(A);
 U := [[]];
 m := nops(A);
 for i from 1 to m-1 do
  U := map(u -> seq([op(u),j],j=0..N-1),U);
 od:

 return [seq(seq(`build/SEM`(N)(A)(s,u),u in U),s in S)];
end:

`count_elements/SEM` := (N::posint) -> (A::set) ->
  nops(A)! * N^(nops(A) - 1);

`mu/F/SEM` := (N::posint) -> (A::set) -> proc(x)
 local eta,a,b,u,v,i;

 eta := table();
 for a in A do
  for b in A do
   u := x[b] -~ x[a];
   v := add(u[i] * epsilon^(i-1),i = 1..N);
   eta[a,b] := `lambda/P/E`(N)(v);
  od;
 od;

 return eval(eta);
end:

`sigma/SEM/F` := (N::posint) -> (A::set) -> proc(eta)
 local n,m,a,b,A1,x,u,i;

 n := table():

 for a in A do
  for b in A do
   n[a,b] := `is_leq/E`(N)(0,eta[a,b]);
  od;
 od;

 A1 := sort([op(A)],(a,b) -> n[a,b]);

 x := table():
 u := 0;
 m := nops(A1);
 x[A1[1]] := [0$N];

 for i from 2 to m do
  u := u + eta[A1[i-1],A1[i]];
  x[A1[i]] := [seq(coeff(u,epsilon,i),i=0..N-1)];
 od;

 return(eval(x));
end:


text:
flavor 1
srclang 3
id 65535
numfuncs 5
import "RCCycleTest.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &LRCCycleTest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "RCCycleTest.jar"}
srcfileinfo {
  1 "RCCycleTest.java"}
javaclass $LRCCycleTest_3B <$LRCCycleTest_3B> public
javaclass $LTest__A_3B <$LTest__A_3B>
javaclass $LTest__B_3B <$LTest__B_3B>
func &LRCCycleTest_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LRCCycleTest_3B>>) void
func &LRCCycleTest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg3_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void
func &LTest__A_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LTest__A_3B>>) void
func &LTest__A_3B_7CgenCycle_7C_28_29V virtual (var %_this <* <$LTest__A_3B>>) void
func &LTest__B_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LTest__B_3B>>) void
var $_C_STR_d4b6a4788952b3ae6554bbc29994ed26 fstatic <[5] u64> readonly = [0, 0x2900000000, 0x73616c63ae9a5a99, 0x6320736920422073, 0x64657463656c6c6f]
var $_C_STR_a360f9874cf7de12fce9dffe2c4c052a fstatic <[6] u64> readonly = [0, 0x3100000000, 0x73616c634a44d58c, 0x6e20736920422073, 0x656c6c6f6320746f, 0x64657463]
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
func &LRCCycleTest_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$LRCCycleTest_3B>>) void {
  funcid 48153
  var %Reg1_R43696 <* <$LRCCycleTest_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43696 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43696))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &LRCCycleTest_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg3_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48154
  var %Reg0_R43699 <* <$LTest__A_3B>>
  var %Reg2_R43699 <* <$LTest__A_3B>>
  var %Reg0_R43698 <* <$LTest__B_3B>>
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>>
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161337 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_161338 <* <$Ljava_2Flang_2FString_3B>>

  intrinsiccallwithtype <$LRCCycleTest_3B> JAVA_CLINIT_CHECK ()
  #INSTIDX : 0||0000:  new
  intrinsiccallwithtype <$LTest__A_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R43699 0 (gcmalloc ref <$LTest__A_3B>)
  #INSTIDX : 3||0003:  dup
  #INSTIDX : 4||0004:  invokespecial
  callassigned &LTest__A_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg0_R43699) {}
  #INSTIDX : 7||0007:  astore_1
  dassign %Reg2_R43699 0 (dread ref %Reg0_R43699)
  #INSTIDX : 8||0008:  aload_1
  #INSTIDX : 9||0009:  invokevirtual
  virtualcallassigned &LTest__A_3B_7CgenCycle_7C_28_29V (dread ref %Reg2_R43699) {}
  #INSTIDX : 12||000c:  aload_1
  #INSTIDX : 13||000d:  getfield
  dassign %Reg0_R43698 0 (iread ref <* <$LTest__A_3B>> 4 (dread ref %Reg2_R43699))
  #INSTIDX : 16||0010:  ifnonnull
  brtrue @label0 (ne i32 ref (dread ref %Reg0_R43698, constval ref 0))
  #INSTIDX : 19||0013:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 22||0016:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_d4b6a4788952b3ae6554bbc29994ed26) { dassign %L_STR_161337 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_161337)
  #INSTIDX : 24||0018:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 27||001b:  goto
  goto @label1
@label0   #INSTIDX : 30||001e:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 33||0021:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_a360f9874cf7de12fce9dffe2c4c052a) { dassign %L_STR_161338 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_161338)
  #INSTIDX : 35||0023:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
@label1   #INSTIDX : 38||0026:  return
  return ()
}
func &LTest__A_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LTest__A_3B>>) void {
  funcid 48155
  var %Reg1_R43699 <* <$LTest__A_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43699 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43699))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &LTest__A_3B_7CgenCycle_7C_28_29V virtual (var %_this <* <$LTest__A_3B>>) void {
  funcid 48156
  var %Reg3_R43699 <* <$LTest__A_3B>>
  var %Reg0_R43698 <* <$LTest__B_3B>>
  var %Reg2_R43698 <* <$LTest__B_3B>>

  dassign %Reg3_R43699 0 (dread ref %_this)
  #INSTIDX : 0||0000:  new
  intrinsiccallwithtype <$LTest__B_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R43698 0 (gcmalloc ref <$LTest__B_3B>)
  #INSTIDX : 3||0003:  dup
  #INSTIDX : 4||0004:  invokespecial
  callassigned &LTest__B_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg0_R43698) {}
  #INSTIDX : 7||0007:  astore_1
  dassign %Reg2_R43698 0 (dread ref %Reg0_R43698)
  #INSTIDX : 8||0008:  aload_0
  #INSTIDX : 9||0009:  aload_1
  #INSTIDX : 10||000a:  putfield
  iassign <* <$LTest__A_3B>> 4 (dread ref %Reg3_R43699, dread ref %Reg2_R43698)
  #INSTIDX : 13||000d:  aload_1
  #INSTIDX : 14||000e:  aload_0
  #INSTIDX : 15||000f:  putfield
  iassign <* <$LTest__B_3B>> 4 (dread ref %Reg2_R43698, dread ref %Reg3_R43699)
  #INSTIDX : 18||0012:  return
  return ()
}
func &LTest__B_3B_7C_3Cinit_3E_7C_28_29V constructor (var %_this <* <$LTest__B_3B>>) void {
  funcid 48157
  var %Reg1_R43698 <* <$LTest__B_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43698 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43698))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}

ext: mpl

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[v, x, z, w, y], [y2, y1] ], derivations = [t]):
sys := [
diff(v(t), t) - (k*y(t) - v(t)*u),
diff(x(t), t) - (lm - x(t)*d - x(t)*v(t)*beta),
diff(z(t), t) - (c*w(t)*q*y(t) - h*z(t)),
diff(w(t), t) - (-b*w(t) + c*w(t)*x(t)*y(t) - c*w(t)*q*y(t)),
diff(y(t), t) - (x(t)*v(t)*beta - a*y(t)),
y2(t) - (z(t)),
y1(t) - (w(t))
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));


text:
#@ Not autoload

E := proc(u)
 local n,i;
 if type(u,`+`) then
  return(map(E,u));
 else
  return('E'(u));
 fi;
end:

H := proc(u)
 local n,i,J;
 if type(u,`+`) then
  return(map(H,u));
 else
  return('H'(u));
 fi;
end:

P := proc(u)
 local n,i;
 if type(u,`+`) then
  return(map(P,u));
 else
  return('P'(u));
 fi;
end:

o := proc(a,b)
 if a=0 or b=0 then
  return(0);
 elif type(b,`+`) then
  return(map2(o,a,b));
 elif type(a,indexed) and op(0,a)=x and op(2,a)=[] then
  return(b);
 elif type(b,indexed) and op(0,b)=x and op(2,b)=[] then
  return(a);
 else
  return('o'(args));
 fi;
end:

is_o := proc(u)
 type(u,specfunc(anything,o));
end:

SpherePi := proc(n,m)
 local stem,is_stable,G;
 stem := n-m;
 is_stable := evalb(m>stem+1);

 if stem<0 or m<=0 or (n>1 and m=1) then
  return(Group());
 elif stem=0 then
  return(Group([],[x[n,[]]]));
 elif is_stable then
  G := StableSpherePi(stem);
  return(Group(map(a -> x[m,op(2,a)],op(1,G))));
 else
  return('SpherePi'(n,m));
 fi;
end:

Source := proc(a)
 if type(a,indexed) and op(0,a) = x then
  return(op(1,a) + `+`(op(op(2,a))));
 elif is_o(a) and nops(a) > 0 then
  return Source(op(nops(a),a))
 else
  return('Source'(args));
 fi;
end:

Target := proc(a)
 if type(a,indexed) and op(0,a) = x then
  return(op(1,a));
 elif is_o(a) and nops(a) > 0 then
  return Target(op(1,a))
 else
  return('Target'(args));
 fi;
end:

Stem := proc(a) Source(a) - Target(a); end:

IsAdmissible := proc(a)
 local m,J,r,i;
 if type(a,indexed) and op(0,a) = x then
  m := op(1,a);
  J := op(2,a);
  r := nops(J);
  if r>0 and J[1] >= m then
   return(false);
  fi;
  for i from 1 to r-1 do
   if J[i+1] > 2*J[i] then
    return(false);
   fi;
  od;
  return(true);
 else
  return('IsAdmissible'(args));
 fi;
end:

# 1-stem

SpherePi( 3, 2) := Group([],[x[2,[1]]]):
SpherePi( 4, 3) := Group([x[3,[1]]]):

StableSpherePi( 1) := Group([x[infinity,[1]]]):

# 2-stem

SpherePi( 4, 2) := Group([x[2,[1,1]]]):
SpherePi( 5, 3) := Group([x[3,[1,1]]]):

StableSpherePi( 2) := Group([x[infinity,[1,1]]]):

# 3-stem

SpherePi( 5, 2) := Group([x[2,[1,1,1]]]):
SpherePi( 6, 3) := Group([x[3,[2,1]],
                          x[3,[1,1,1]]]):
SpherePi( 7, 4) := Group([x[4,[2,1]],
                          x[4,[1,1,1]]],
                         [x[4,[3]]]):

StableSpherePi( 3) := Group([x[infinity,[3]],
		             x[infinity,[2,1]],
		             x[infinity,[1,1,1]]]):

# 4-stem

SpherePi( 6, 2) := Group([x[2,[1,1,1,1]],
                          x[2,[1,2,1]]]):
SpherePi( 7, 3) := Group([x[3,[2,1,1]]]):
SpherePi( 8, 4) := Group([x[4,[3,1]],
                          x[4,[2,1,1]]]):
SpherePi( 9, 5) := Group([x[5,[3,1]]]):

StableSpherePi( 4) := Group([]):

# 5-stem

SpherePi( 7, 2) := Group([x[2,[1,2,1,1]]]):
SpherePi( 8, 3) := Group([x[3,[2,1,1,1]]]):
SpherePi( 9, 4) := Group([x[4,[3,1,1]],
                          x[4,[2,1,1,1]]]):
SpherePi(10, 5) := Group([x[5,[3,1,1]]]):
SpherePi(11, 6) := Group([],[x[6,[5,0]]]):

StableSpherePi( 5) := Group([]):

# 6-stem

SpherePi( 8, 2) := Group([x[2,[1,2,1,1,1]]]):
SpherePi( 9, 3) := Group([]):
SpherePi(10, 4) := Group([x[4,[3,3]],
                          x[4,[3,2,1]],
                          x[4,[3,1,1,1]]]):
SpherePi(11, 5) := Group([x[5,[3,3]]]):
SpherePi(12, 6) := Group([x[6,[3,3]]]):
SpherePi(13, 7) := Group([x[7,[3,3]]]):

StableSpherePi( 6) := Group([x[infinity,[3,3]]]):

# 7-stem

SpherePi( 9, 2) := Group([]):
SpherePi(10, 3) := Group([]):
SpherePi(11, 4) := Group([]):
SpherePi(12, 5) := Group([x[5,[4,1,1,1]]]):
SpherePi(13, 6) := Group([x[6,[5,1,1]],
                          x[6,[4,1,1,1]]]):
SpherePi(14, 7) := Group([x[7,[6,1]],
                          x[7,[5,1,1]],
                          x[7,[4,1,1,1]]]):
SpherePi(15, 8) := Group([x[8,[6,1]]],[x[8,[7]]]):

StableSpherePi( 7) := Group([x[infinity,[7]],
		             x[infinity,[6,1]],
                             x[infinity,[5,1,1]],
                             x[infinity,[4,1,1,1]]]):

# 8-stem

SpherePi(10, 2) := Group([]):
SpherePi(11, 3) := Group([x[3,[2,3,3]]]):
SpherePi(12, 4) := Group([x[4,[2,3,3]]]):
SpherePi(13, 5) := Group([x[5,[2,3,3]]]):
SpherePi(14, 6) := Group([x[6,[5,3]],
                          x[6,[5,2,1]],
                          x[6,[5,1,1,1]],
                          x[6,[2,3,3]]]):
SpherePi(15, 7) := Group([x[7,[6,1,1]],
                          x[7,[5,3]],
                          x[7,[2,3,3]]]):
SpherePi(16, 8) := Group([x[8,[7,1]],
                          x[8,[6,1,1]],
                          x[8,[5,3]],
                          x[8,[2,3,3]]]):
SpherePi(17, 9) := Group([x[9,[7,1]],
                          x[9,[5,3]],
                          x[9,[2,3,3]]]):

StableSpherePi( 8) := Group([x[infinity,[5,3]],
                             x[infinity,[2,3,3]]]):

# 9-stem

SpherePi(11, 2) := Group([x[2,[1,2,3,3]]]):
SpherePi(12, 3) := Group([x[3,[2,4,1,1,1]],
                          x[3,[1,2,3,3]]]):
SpherePi(13, 4) := Group([x[4,[3,3,3]],
                          x[4,[2,4,1,1,1]],
                          x[4,[1,2,3,3]]]):
SpherePi(14, 5) := Group([x[5,[3,3,3]],
                          x[5,[2,4,1,1,1]],
                          x[5,[1,2,3,3]]]):
SpherePi(15, 6) := Group([x[6,[3,3,3]],
                          x[6,[2,4,1,1,1]],
                          x[6,[1,2,3,3]]]):
SpherePi(16, 7) := Group([x[7,[6,1,1,1]],
                          x[7,[3,3,3]],
                          x[7,[2,4,1,1,1]],
                          x[7,[1,2,3,3]]]):
SpherePi(17, 8) := Group([x[8,[7,1,1]],
                          x[8,[6,1,1,1]],
                          x[8,[3,3,3]],
                          x[8,[2,4,1,1,1]],
                          x[8,[1,2,3,3]]]):
SpherePi(18, 9) := Group([x[9,[7,1,1]],
                          x[9,[3,3,3]],
                          x[9,[2,4,1,1,1]],
                          x[9,[1,2,3,3]]]):
SpherePi(19,10) := Group([x[10,[3,3,3]],
                          x[10,[2,4,1,1,1]],
                          x[10,[1,2,3,3]]],
                         [x[10,[9,0]]]):

StableSpherePi( 9) := Group([x[infinity,[3,3,3]],
                             x[infinity,[2,4,1,1,1]],
                             x[infinity,[1,2,3,3]]]):

# 10-stem

SpherePi(12, 2) := Group([x[2,[1,2,4,1,1,1]],
                          x[2,[1,1,2,3,3]]]):
SpherePi(13, 3) := Group([x[3,[2,2,3,3]],
                          x[3,[1,2,4,1,1,1]],
                          x[3,[1,1,2,3,3]]]):
SpherePi(14, 4) := Group([x[4,[3,6,1]],
                          x[4,[3,5,1,1]],
                          x[4,[3,4,1,1,1]],
                          x[4,[2,2,3,3]],
                          x[4,[1,2,4,1,1,1]],
                          x[4,[1,1,2,3,3]]]):
SpherePi(15, 5) := Group([x[5,[4,3,3]],
                          x[5,[2,2,3,3]],
                          x[5,[1,2,4,1,1,1]],
                          x[5,[1,1,2,3,3]]]):
SpherePi(16, 6) := Group([x[6,[4,3,3]],
                          x[6,[2,2,3,3]],
                          x[6,[1,2,4,1,1,1]],
                          x[6,[1,1,2,3,3]]]):
SpherePi(17, 7) := Group([x[7,[4,3,3]],
                          x[7,[2,2,3,3]],
                          x[7,[1,2,4,1,1,1]],
                          x[7,[1,1,2,3,3]]]):
SpherePi(18, 8) := Group([x[8,[7,3]],
                          x[8,[7,2,1]],
                          x[8,[7,1,1,1]],
                          x[8,[4,3,3]],
                          x[8,[2,2,3,3]],
                          x[8,[1,2,4,1,1,1]],
                          x[8,[1,1,2,3,3]]]):
SpherePi(19, 9) := Group([x[9,[7,3]],
                          x[9,[7,2,1]],
                          x[9,[7,1,1,1]],
                          x[9,[1,2,4,1,1,1]]]):
SpherePi(20,10) := Group([x[10,[7,3]],
                          x[10,[7,2,1]],
                          x[10,[1,2,4,1,1,1]]]):
SpherePi(21,11) := Group([x[11,[7,3]],
                          x[11,[1,2,4,1,1,1]]]):

StableSpherePi(10) := Group([x[infinity,[1,2,4,1,1,1]]]):

# 11-stem

SpherePi(13, 2) := Group([x[2,[1,2,2,3,3]],
                          x[2,[1,1,1,2,3,3]],
                          x[2,[1,1,2,4,1,1,1]]]):
SpherePi(14, 3) := Group([x[3,[2,3,3,3]],
                          x[3,[2,2,4,1,1,1]],
                          x[3,[2,1,2,3,3]],
                          x[3,[1,1,2,4,1,1,1]]]):
SpherePi(15, 4) := Group([x[4,[3,6,1,1]],
                          x[4,[3,5,3]],
                          x[4,[3,2,3,3]],
                          x[4,[2,3,3,3]],
                          x[4,[2,2,4,1,1,1]],
                          x[4,[2,1,2,3,3]],
                          x[4,[1,1,2,4,1,1,1]]]):
SpherePi(16, 5) := Group([x[5,[4,4,1,1,1]],
                          x[5,[3,5,3]],
                          x[5,[3,2,3,3]],
                          x[5,[2,2,4,1,1,1]],
                          x[5,[1,1,2,4,1,1,1]]]):
SpherePi(17, 6) := Group([x[6,[5,3,3]],
                          x[6,[4,4,1,1,1]],
                          x[6,[2,2,4,1,1,1]],
                          x[6,[1,1,2,4,1,1,1]]]):
SpherePi(18, 7) := Group([x[7,[5,3,3]],
                          x[7,[4,4,1,1,1]],
                          x[7,[2,2,4,1,1,1]],
                          x[7,[1,1,2,4,1,1,1]]]):
SpherePi(19, 8) := Group([x[8,[5,3,3]],
                          x[8,[4,4,1,1,1]],
                          x[8,[2,2,4,1,1,1]],
                          x[8,[1,1,2,4,1,1,1]]]):
SpherePi(20, 9) := Group([x[9,[5,3,3]],
                          x[9,[4,4,1,1,1]],
                          x[9,[2,2,4,1,1,1]],
                          x[9,[1,1,2,4,1,1,1]]]):
SpherePi(21,10) := Group([x[10,[4,4,1,1,1]],
                          x[10,[2,2,4,1,1,1]],
                          x[10,[1,1,2,4,1,1,1]]]):
SpherePi(22,11) := Group([x[11,[4,4,1,1,1]],
                          x[11,[2,2,4,1,1,1]],
                          x[11,[1,1,2,4,1,1,1]]]):
SpherePi(23,12) := Group([x[12,[4,4,1,1,1]],
                          x[12,[2,2,4,1,1,1]],
                          x[12,[1,1,2,4,1,1,1]]]):

StableSpherePi(11) := Group([x[infinity,[4,4,1,1,1]],
                             x[infinity,[2,2,4,1,1,1]],
                             x[infinity,[1,1,2,4,1,1,1]]]):

# 12-stem

StableSpherePi(12) := Group([]):

# 13-stem

StableSpherePi(13) := Group([]):

# 14-stem

StableSpherePi(14) := Group([x[infinity,[7,7]],
                             x[infinity,[6,2,3,3]]]):

N := 11:
degrees_above := (k) -> [seq(k+i,i=0..N),infinity];

for n in degrees_above(1) do
 TodaName(x[n,[]]) := iota[n]:
od:

for n in degrees_above(2) do
 TodaName(x[n,[1]]) := eta[n]:
 TodaName(x[n,[1,1]]) := o(eta[n],eta[n+1]):
 TodaName(x[n,[1,2,3,3]]) := o(eta[n],epsilon[n+1]):
 TodaName(x[n,[1,2,4,1,1,1]]) := o(eta[n],mu[n+1]):
od:

for n in degrees_above(3) do
 TodaName(x[n,[2,3,3]]) := epsilon[n]:
 TodaName(x[n,[2,4,1,1,1]]) := mu[n]:
od:

for n in degrees_above(4) do
 TodaName(x[n,[3]]) := nu[n]:
 TodaName(x[n,[3,3,3]]) := o(nu[n],nu[n+3],nu[n+6]):
od:

for n in degrees_above(5) do
 TodaName(x[n,[4,4,1,1,1]]) := zeta[n]:
 TodaName(x[n,[2,2,4,1,1,1]]) := 2 * zeta[n]:
 TodaName(x[n,[1,1,2,4,1,1,1]]) := 4 * zeta[n]:
 TodaName(x[n,[3,3]]) := o(nu[n],nu[n+3]):
 TodaName(x[n,[2,1]]) :=  2 * nu[n]:
 TodaName(x[n,[1,1,1]]) :=  4 * nu[n]:
od:

for n in degrees_above(6) do
 TodaName(x[n,[5,3]]) := nubar[n]:
od:

for n in degrees_above(8) do
 TodaName(x[n,[7]]) := sigma[n]:
od:

for n in degrees_above(9) do
 TodaName(x[n,[6,1]]) := 2*sigma[n]:
 TodaName(x[n,[5,1,1]]) := 4*sigma[n]:
 TodaName(x[n,[4,1,1,1]]) := 8*sigma[n]:
od:

TodaName(x[2,[1,1,1,1]]) := - o(eta[2],nuprime):
TodaName(x[8,[7]]) := sigma[8]:
TodaName(x[2,[1,2,2,3,3]]) := 3 * o(eta[2],epsilonprime):
TodaName(x[2,[1,2,2,4,1,1,1]]) := 3 * o(eta[2],muprime):
TodaName(x[2,[1,1,1]]) := o(eta[2],eta[3],eta[4]):
TodaName(x[3,[1,1,1]]) := 2 * nuprime:
TodaName(x[3,[2,1]]) := nuprime:
TodaName(x[4,[1,1,1]]) := 2 * E(nuprime):
TodaName(x[4,[2,1]]) := E(nuprime):
TodaName(x[2,[1,2,1]]) := 2 * o(eta[2],nuprime):
TodaName(x[3,[2,1,1]]) := o(nuprime,eta[6]):
TodaName(x[4,[2,1,1]]) := o(E(nuprime),eta[7]):
TodaName(x[4,[3,1]]) := o(nu[4],eta[7]):
TodaName(x[5,[3,1]]) := o(nu[5],eta[8]):
TodaName(x[2,[1,2,1,1]]) := o(eta[2],nuprime,eta[6]):
TodaName(x[3,[2,1,1,1]]) := o(nuprime,eta[6],eta[7]):
TodaName(x[4,[3,1,1]]) := o(nu[4],eta[7],eta[8]):
TodaName(x[4,[2,1,1,1]]) := o(E(nuprime),eta[7],eta[8]):
TodaName(x[5,[3,1,1]]) := o(nu[5],eta[8],eta[9]):
TodaName(x[6,[5,0]]) := w[6]:
TodaName(x[2,[1,2,1,1,1]]) := o(eta[2],nuprime,eta[6],eta[7]):
TodaName(x[4,[3,3]]) := o(nu[4],nu[7]):
TodaName(x[4,[3,2,1]]) := 2*o(nu[4],nu[7]):
TodaName(x[4,[3,1,1,1]]) := 4*o(nu[4],nu[7]):
TodaName(x[5,[4,1,1,1]]) := sigmathird:
TodaName(x[6,[5,1,1]]) := sigmasecond:
TodaName(x[6,[4,1,1,1]]) := 2*sigmasecond:
TodaName(x[7,[6,1]]) := sigmaprime:
TodaName(x[7,[5,1,1]]) := 2*sigmaprime:
TodaName(x[7,[4,1,1,1]]) := 4*sigmaprime:
TodaName(x[8,[6,1]]) := E(sigmaprime):
TodaName(x[8,[5,1,1]]) := 2*E(sigmaprime):
TodaName(x[8,[4,1,1,1]]) := 4*E(sigmaprime):
TodaName(x[6,[5,2,1]]) := 2*nubar[6]:
TodaName(x[6,[5,1,1,1]]) := 4*nubar[6]:
TodaName(x[7,[6,1,1]]) := o(sigmaprime,eta[14]):
TodaName(x[8,[6,1,1]]) := o(E(sigmaprime),eta[15]):
TodaName(x[8,[7,1]]) := o(sigma[8],eta[15]):
TodaName(x[9,[7,1]]) := o(sigma[9],eta[16]):
TodaName(x[7,[6,1,1,1]]) := o(sigmaprime,eta[14],eta[15]):
TodaName(x[8,[6,1,1,1]]) := o(E(sigmaprime),eta[15],eta[16]):
TodaName(x[8,[7,1,1]]) := o(sigma[8],eta[15],eta[16]):
TodaName(x[9,[7,1,1]]) := o(sigma[9],eta[16],eta[17]):
TodaName(x[2,[1,1,2,3,3]]) := o(eta[2],eta[3],epsilon[4]):
TodaName(x[3,[2,2,3,3]]) := epsilonprime:
TodaName(x[3,[1,1,2,3,3]]) := 2*epsilonprime:
TodaName(x[4,[2,2,3,3]]) := E(epsilonprime):
TodaName(x[4,[1,1,2,3,3]]) := 2*E(epsilonprime):
TodaName(x[4,[3,6,1]]) := o(nu[4],sigmaprime):
TodaName(x[4,[3,5,1,1]]) := 2*o(nu[4],sigmaprime):
TodaName(x[4,[3,4,1,1,1]]) := 4*o(nu[4],sigmaprime):
TodaName(x[5,[4,3,3]]) := o(nu[5],sigma[8]):
TodaName(x[5,[2,2,3,3]]) := 2*o(nu[5],sigma[8]):
TodaName(x[5,[1,1,2,3,3]]) := 4*o(nu[5],sigma[8]):
TodaName(x[6,[4,3,3]]) := o(nu[6],sigma[9]):
TodaName(x[6,[2,2,3,3]]) := 2*o(nu[6],sigma[9]):
TodaName(x[6,[1,1,2,3,3]]) := 4*o(nu[6],sigma[9]):
TodaName(x[7,[4,3,3]]) := o(nu[7],sigma[10]):
TodaName(x[7,[2,2,3,3]]) := 2*o(nu[7],sigma[10]):
TodaName(x[7,[1,1,2,3,3]]) := 4*o(nu[7],sigma[10]):
TodaName(x[8,[4,3,3]]) := o(nu[8],sigma[11]):
TodaName(x[8,[2,2,3,3]]) := 2*o(nu[8],sigma[11]):
TodaName(x[8,[1,1,2,3,3]]) := 4*o(nu[8],sigma[11]):
TodaName(x[8,[7,3]]) := o(sigma[8],nu[15]):
TodaName(x[8,[7,2,1]]) := 2*o(sigma[8],nu[15]):
TodaName(x[8,[7,1,1,1]]) := 4*o(sigma[8],nu[15]):
TodaName(x[9,[7,3]]) := o(sigma[9],nu[16]):
TodaName(x[9,[7,2,1]]) := 2*o(sigma[9],nu[16]):
TodaName(x[9,[7,1,1,1]]) := 4*o(sigma[9],nu[16]):
TodaName(x[10,[9,0]]) := w[10]:
TodaName(x[10,[7,3]]) := o(sigma[10],nu[17]):
TodaName(x[10,[7,2,1]]) := 2*o(sigma[10],nu[17]):
TodaName(x[11,[7,3]]) := o(sigma[11],nu[18]):
TodaName(x[2,[1,1,1,2,3,3]]) := 2*o(eta[2],epsilonprime):
TodaName(x[2,[1,1,2,4,1,1,1]]) := o(eta[2],eta[3],mu[4]):
TodaName(x[3,[2,2,4,1,1,1]]) := muprime:
TodaName(x[3,[1,1,2,4,1,1,1]]) := 2*muprime:
TodaName(x[3,[2,3,3,3]]) := o(epsilon[3],nu[11]):
TodaName(x[3,[2,1,2,3,3]]) := o(nuprime,epsilon[6]):
TodaName(x[4,[2,2,4,1,1,1]]) := E(muprime):
TodaName(x[4,[1,1,2,4,1,1,1]]) := 2*E(muprime):
TodaName(x[4,[2,3,3,3]]) := o(epsilon[4],nu[12]):
TodaName(x[4,[2,1,2,3,3]]) := o(E(nuprime),epsilon[7]):
TodaName(x[4,[3,6,1,1]]) := o(nu[4],sigmaprime,eta[14]):
TodaName(x[4,[3,5,3]]) := o(nu[4],nubar[7]):
TodaName(x[4,[3,2,3,3]]) := o(nu[4],epsilon[7]):
TodaName(x[5,[3,5,3]]) := o(nu[5],nubar[8]):
TodaName(x[5,[3,2,3,3]]) := o(nu[5],epsilon[8]):
TodaName(x[6,[5,3,3]]) := o(nubar[6],nu[14]):
TodaName(x[6,[3,2,3,3]]) := 2*o(nubar[6],nu[14]):
TodaName(x[7,[5,3,3]]) := o(nubar[7],nu[15]):
TodaName(x[8,[5,3,3]]) := o(nubar[8],nu[16]):
TodaName(x[9,[5,3,3]]) := o(nubar[9],nu[17]):

generators  :=  proc(G)
 map(op,[op(G)]):
end:

for i from 0 to N do
 for j from 1 to i+2 do
  G  :=  SpherePi(i+j,j):
  for a in generators(G) do
   if not(Source(a)=i+j and
          Target(a)=j and
          IsAdmissible(a)) then
    print([i+j,j,a]):
   fi:
  od:
 od:
od:

all_gens  :=  [seq(seq(op(generators(SpherePi(i,j))),i=1..j+N),j=1..N+2)]:

for a in all_gens do
 GN(TodaName(a)) := a:
od:

for a in all_gens do
 m := Source(a):
 n := Target(a):
 b := two(a):
 if type(b,specfunc(anything,two)) then
  two(a) := 0:
 fi:
 b := E(a):
 if type(b,specfunc(anything,E)) then
  E(a) := x[op(1,a)+1,op(2,a)]:
 fi:
 b := H(a):
 if type(b,specfunc(anything,H)) then
  if (2*m-1>n or n=m) then
   H(a) := 0:
  else
   J := op(2,a):
   if J[1] > n-1 then
    H(a) := 0:
   else
    H(a) := x[2*n-1,J[2..-1]]:
   fi:
  fi:
 fi:

 o(a,x[m,[]]) := a;
 o(x[n,[]],a) := a;
od:

for n in degrees_above(5) do
 two(x[n,[3]]) := x[n,[2,1]]:
 two(x[n,[2,1]]) := x[n,[1,1,1]]:
 two(x[n,[4,4,1,1,1]]) := x[n,[2,2,4,1,1,1]]:
 two(x[n,[2,2,4,1,1,1]]) := x[n,[1,1,2,4,1,1,1]]:
od:

for n in degrees_above(6) do
 two(x[n,[5,1,1]])       := x[n,[4,1,1,1]]:
od:

for n in degrees_above(9) do
 two(x[n,[7]]) := x[n,[6,1]]:
 two(x[n,[6,1]]) := x[n,[5,1,1]]:
 two(x[n,[5,1,1]]) := x[n,[4,1,1,1]]:
od:

two(x[2,[1]])           := x[2,[1,0]]:
two(x[2,[1,1,1,1]])     := x[2,[1,2,1]]:
two(x[2,[1,1,1,1]])     := x[2,[1,2,1]]:
two(x[2,[1,2,2,3,3]])   := x[2,[1,1,1,2,3,3]]:
two(x[3,[2,1]])         := x[3,[1,1,1]]:
two(x[4,[3]])           := x[4,[3,0]]:
two(x[4,[3,2,1]])       := x[4,[3,1,1,1]]:
two(x[4,[3,3]])         := x[4,[3,2,1]]:
two(x[4,[3,6,1]])       := x[4,[3,5,1,1]]:
two(x[5,[2,2,3,3]])     := x[5,[1,1,2,3,3]]:
two(x[5,[3,5,1,1]])     := x[5,[3,4,1,1,1]]:
two(x[5,[4,3,3]])       := x[5,[2,2,3,3]]:
two(x[6,[2,2,3,3]])     := x[6,[1,1,2,3,3]]:
two(x[6,[4,3,3]])       := x[6,[2,2,3,3]]:
two(x[6,[5,1,1]])       := x[6,[4,1,1,1]]:
two(x[6,[5,2,1]])       := x[6,[5,1,1,1]]:
two(x[6,[5,3,3]])       := x[6,[3,2,3,3]]:
two(x[6,[5,3]])         := x[6,[5,2,1]]:
two(x[7,[2,2,3,3]])     := x[7,[1,1,2,3,3]]:
two(x[7,[4,3,3]])       := x[7,[2,2,3,3]]:
two(x[7,[5,1,1,1]])     := x[7,[4,1,1,1,1]]:
two(x[7,[6,1]])         := x[7,[5,1,1]]:
two(x[8,[2,2,3,3]])     := x[8,[1,1,2,3,3]]:
two(x[8,[4,3,3]])       := x[8,[2,2,3,3]]:
two(x[8,[5,1,1,1]])     := x[8,[4,1,1,1,1]]:
two(x[8,[6,1]])         := x[8,[5,1,1]]:
two(x[8,[7]])           := x[8,[7,0]]:
two(x[8,[7,2,1]])       := x[8,[7,1,1,1]]:
two(x[8,[7,3]])         := x[8,[7,2,1]]:
two(x[9,[7,2,1]])       := x[9,[7,1,1,1]]:
two(x[9,[7,3]])         := x[9,[7,2,1]]:
two(x[10,[7,3]])        := x[10,[7,2,1]]:

for n in degrees_above(2) do
 o(x[n,[1]],x[n+1,[1]]) := x[n,[1,1]]:
 o(x[n,[1]],x[n+1,[1,1]]) := x[n,[1,1,1]]:
 o(x[n,[1,1]],x[n+2,[1]]) := x[n,[1,1,1]]:
 o(x[n,[1]],x[n+1,[2,3,3]]) := x[n,[1,2,3,3]]:
 o(x[n,[1]],x[n+1,[2,4,1,1,1]]) := x[n,[1,2,4,1,1,1]]:
od:

for n in degrees_above(4) do
 o(x[n,[3]],x[n+3,[3]]) := x[n,[3,3]]:
 o(x[n,[3]],x[n+3,[3,3]]) := x[n,[3,3,3]]:
 o(x[n,[3,3]],x[n+6,[3]]) := x[n,[3,3,3]]:
od:

# Predictable ternary+ composites
o(x[2,[1]],    x[3,[2,1]],   x[6,[1]],   x[7,[1]]) := x[2,[1,2,1,1,1]]:
o(x[3,[2,1]],  x[6,[1]],     x[7,[1]])             := x[3,[2,1,1,1]]:
o(x[4,[2,1]],  x[7,[1]],     x[8,[1]])             := x[4,[2,1,1,1]]:
o(x[4,[3]],    x[7,[1]],     x[8,[1]])             := x[4,[3,1,1]]:
o(x[4,[3]],    x[7,[6,1]],   x[14,[1]])            := x[4,[3,6,1,1]]:
o(x[4,[3]],    x[7,[3]],     x[10,[3]])            := x[4,[3,3,3]]:
o(x[5,[3]],    x[8,[1]],     x[9,[1]])             := x[5,[3,1,1]]:
o(x[7,[6,1]],  x[14,[1]],    x[15,[1]])            := x[7,[6,1,1,1]]:
o(x[8,[6,1]],  x[15,[1]],    x[16,[1]])            := x[8,[6,1,1,1]]:
o(x[8,[7]],    x[15,[1]],    x[16,[1]])            := x[8,[7,1,1]]:
o(x[9,[7]],    x[16,[1]],    x[17,[1]])            := x[9,[7,1,1]]:

# Predictable composites
o(x[2,[1]],    x[3,[1,2,3,3]])     := x[2,[1,1,2,3,3]]:
o(x[2,[1]],    x[3,[1,1,2,3,3]])   := x[2,[1,1,1,2,3,3]]:
o(x[2,[1]],    x[3,[1,2,4,1,1,1]]) := x[2,[1,1,2,4,1,1,1]]:
o(x[2,[1]],    x[3,[2,1]])         := x[2,[1,2,1]]:
o(x[2,[1]],    x[3,[2,1,1]])       := x[2,[1,2,1,1]]:
o(x[2,[1]],    x[3,[2,1,1,1]])     := x[2,[1,2,1,1,1]]:
o(x[2,[1]],    x[3,[1,1,1]])       := x[2,[1,1,1,1]]:
o(x[2,[1,1]],  x[4,[2,1]])         := 0:
o(x[2,[1,1]],  x[4,[1,1]])         := x[2,[1,1,1,1]]: #?
o(x[2,[1,1,1]],x[5,[1]])           := x[2,[1,1,1,1]]:
o(x[2,[1]],    x[3,[2,1]])         := x[2,[1,2,1]]:
o(x[2,[1,2,1]],x[6,1])             := x[2,[1,2,1,1]]:
o(x[3,[1]],    x[4,[2,1]])         := 0:
o(x[3,[1]],    x[4,[2,1,1]])       := 0:
o(x[3,[1]],    x[4,[2,4,1,1,1]])   := x[3,[1,2,4,1,1,1]]:
o(x[3,[2,1]],  x[6,[1]])           := x[3,[2,1,1]]:
o(x[3,[2,1]],  x[6,[1,1]])         := x[3,[2,1,1,1]]:
o(x[3,[2,1,1]],x[7,[1]])           := x[3,[2,1,1,1]]:
o(x[3,[2,1]],  x[6,[2,3,3]])       := x[3,[2,1,2,3,3]]:
o(x[3,[2,3,3]],x[11,[3]])          := x[3,[2,3,3,3]]:
o(x[4,[2,1]],  x[7,[1]])           := x[4,[2,1,1]]:
o(x[4,[2,1]],  x[7,[2,3,3]])       := x[4,[2,1,2,3,3]]:
o(x[4,[2,3,3]],x[12,[3]])          := x[4,[2,3,3,3]]:
o(x[4,[3]],    x[7,[1]])           := x[4,[3,1]]:
o(x[4,[3]],    x[7,[2,3,3]])       := x[4,[3,2,3,3]]:
o(x[4,[3]],    x[7,[3]])           := x[4,[3,3]]:
o(x[4,[3]],    x[7,[5,3]])         := x[4,[3,5,3]]:
o(x[4,[3]],    x[7,[6,1]])         := x[4,[3,6,1]]:
o(x[5,[3]],    x[8,[1]])           := x[5,[3,1]]:
o(x[5,[3]],    x[8,[2,3,3]])       := x[5,[3,2,3,3]]:
o(x[5,[3]],    x[8,[5,3]])         := x[5,[3,5,3]]:
o(x[6,[5,3]],  x[14,[3]])          := x[6,[5,3,3]]:
o(x[7,[5,3]],  x[15,[3]])          := x[7,[5,3,3]]:
o(x[7,[6,1]],  x[14,[1]])          := x[7,[6,1,1]]:
o(x[8,[5,3]],  x[16,[3]])          := x[8,[5,3,3]]:
o(x[8,[6,1]],  x[15,[1]])          := x[8,[6,1,1]]:
o(x[8,[7]],    x[15,[1]])          := x[8,[7,1]]:
o(x[8,[7]],    x[15,[3]])          := x[8,[7,3]]:
o(x[9,[5,3]],  x[17,[3]])          := x[9,[5,3,3]]:
o(x[9,[7]],    x[16,[1]])          := x[9,[7,1]]:
o(x[9,[7]],    x[16,[3]])          := x[9,[7,3]]:
o(x[10,[7]],   x[17,[3]])          := x[10,[7,3]]:
o(x[11,[7]],   x[18,[3]])          := x[11,[7,3]]:
o(x[2,[1,1]],  x[4,[1]])           := x[2,[1,1,1]]:
o(x[2,[1,1]],  x[4,[2,3,3]])       := x[2,[1,1,2,3,3]]:

o(x[2, [1, 1]], x[4, [2, 4, 1, 1, 1]]) := x[2, [1, 1, 2, 4, 1, 1, 1]]:
o(x[2, [1, 1]], x[4, [2, 1, 1]]) := 0:
o(x[2, [1, 2, 1, 1]], x[7, [1]]) := x[2, [1, 2, 1, 1, 1]]:

# Unpredictable composites
o(x[2,[1]],    x[3,[2,1]])   := x[2,[1,1,1,1]] + x[2,[1,2,1]]:
o(x[2,[1]],    x[3,[2,2,3,3]]) := x[2,[1,2,2,3,3]]+x[2,[1,1,1,2,3,3]]:
o(x[5,[3]],    x[8,[7]])     := x[5,[4,3,3]]:
o(x[6,[3]],    x[9,[7]])     := x[6,[4,3,3]]:
o(x[7,[3]],    x[10,[7]])    := x[7,[4,3,3]]:
o(x[8,[3]],    x[11,[7]])    := x[8,[4,3,3]]:

# Suspensions
E(x[2,[1,1,1,1]]) := 0:
E(x[2,[1,2,1]])   := 0:

# (************************************************************************)
# (* Suspensions                                                          *)

# TodaName(x[6,[5,0]])               := w[6]:
# TodaName(x[6,[5,1,1]])             := sigmasecond:
# TodaName(x[7,[6,1]])               := sigmaprime:
# TodaName(x[8,[6,1]])               := E(sigmaprime):
# TodaName(x[10,[9,0]])              := w[10]:

# TodaName(x[2,[1,1,1,2,3,3]])       := 2*o(x[2,[1]],x[3,[2,2,3,3]]):
# TodaName(x[2,[1,2,1]])             := 2*o(x[2,[1]],x[3,[2,1]]):
# TodaName(x[3,[1,1,1]])             := 2*x[3,[2,1]]:
# TodaName(x[3,[1,1,2,3,3]])         := 2*x[3,[2,2,3,3]]:
# TodaName(x[3,[1,1,2,4,1,1,1]])     := 2*x[3,[2,2,4,1,1,1]]:
# TodaName(x[4,[1,1,1]])             := 2*x[4,[2,1]]:
# TodaName(x[4,[1,1,2,3,3]])         := 2*x[4,[2,2,3,3]]:
# TodaName(x[4,[1,1,2,4,1,1,1]])     := 2*x[4,[2,2,4,1,1,1]]:
# TodaName(x[4,[3,2,1]])             := 2*o(x[4,[3]],x[7,[3]]):
# TodaName(x[4,[3,5,1,1]])           := 2*o(x[4,[3]],x[7,[6,1]]):
# TodaName(x[5,[2,2,3,3]])           := 2*o(x[5,[3]],x[8,[7]]):
# TodaName(x[6,[2,2,3,3]])           := 2*o(x[6,[3]],x[9,[7]]):
# TodaName(x[6,[3,2,3,3]])           := 2*o([6,[5,3]],x[14,[3]]):
# TodaName(x[6,[4,1,1,1]])           := 2*x[6,[5,1,1]]:
# TodaName(x[6,[5,2,1]])             := 2*[6,[5,3]]:
# TodaName(x[7,[2,2,3,3]])           := 2*o(x[7,[3]],x[10,[7]]):
# TodaName(x[7,[5,1,1]])             := 2*x[7,[6,1]]:
# TodaName(x[8,[2,2,3,3]])           := 2*o(x[8,[3]],x[11,[7]]):
# TodaName(x[8,[5,1,1]])             := 2*x[8,[6,1]]:
# TodaName(x[8,[7,2,1]])             := 2*o(x[8,[7]],x[15,[3]]):
# TodaName(x[9,[7,2,1]])             := 2*o(x[9,[7]],x[16,[3]]):
# TodaName(x[10,[7,2,1]])            := 2*o(x[10,[7]],x[17,[3]]):

# TodaName(x[4,[3,1,1,1]])           := 4*o(x[4,[3]],x[7,[3]]):
# TodaName(x[4,[3,4,1,1,1]])         := 4*o(x[4,[3]],x[7,[6,1]]):
# TodaName(x[5,[1,1,2,3,3]])         := 4*o(x[5,[3]],x[8,[7]]):
# TodaName(x[6,[1,1,2,3,3]])         := 4*o(x[6,[3]],x[9,[7]]):
# TodaName(x[6,[5,1,1,1]])           := 4*[6,[5,3]]:
# TodaName(x[7,[1,1,2,3,3]])         := 4*o(x[7,[3]],x[10,[7]]):
# TodaName(x[7,[4,1,1,1]])           := 4*x[7,[6,1]]:
# TodaName(x[8,[1,1,2,3,3]])         := 4*o(x[8,[3]],x[11,[7]]):
# TodaName(x[8,[4,1,1,1]])           := 4*x[8,[6,1]]:
# TodaName(x[8,[7,1,1,1]])           := 4*o(x[8,[7]],x[15,[3]]):
# TodaName(x[9,[7,1,1,1]])           := 4*o(x[9,[7]],x[16,[3]]):

# TodaName(x[2,[1,1,1,1]])           := - o(x[2,[1]],x[3,[2,1]]):
# TodaName(x[2,[1,2,2,3,3]])         := 3 * o(x[2,[1]],x[3,[2,2,3,3]]):
# TodaName(x[2,[1,2,2,4,1,1,1]])     := 3 * o(x[2,[1]],x[3,[2,2,4,1,1,1]]):

# TodaName(x[2,[1,1,2,3,3]])         := o(x[2,[1]],x[3,[1]],x[4,[2,3,3]]):
# TodaName(x[2,[1,1,2,4,1,1,1]])     := o(x[2,[1]],x[3,[1]],x[4,[2,4,1,1,1]]):
# TodaName(x[2,[1,2,1,1,1]])         := o(x[2,[1]],x[3,[2,1]],x[6,[1]],x[7,[1]]):
# TodaName(x[2,[1,2,1,1]])           := o(x[2,[1]],x[3,[2,1]],x[6,[1]]):
# TodaName(x[3,[2,1,1,1]])           := o(x[3,[2,1]],x[6,[1]],x[7,[1]]):
# TodaName(x[4,[2,1,1,1]])           := o(x[4,[2,1]],x[7,[1]],x[8,[1]]):
# TodaName(x[4,[3,1,1]])             := o(x[4,[3]],x[7,[1]],x[8,[1]]):
# TodaName(x[4,[3,6,1,1]])           := o(x[4,[3]],x[7,[6,1]],x[14,[1]]):
# TodaName(x[7,[6,1,1,1]])           := o(x[7,[6,1]],x[14,[1]],x[15,[1]]):
# TodaName(x[8,[6,1,1,1]])           := o(x[8,[6,1]],x[15,[1]],x[16,[1]]):
# TodaName(x[8,[7,1,1]])             := o(x[8,[7]],x[15,[1]],x[16,[1]]):
# TodaName(x[9,[7,1,1]])             := o(x[9,[7]],x[16,[1]],x[17,[1]]):
# TodaName(x[5,[3,1,1]])             := o(x[5,[3]],x[8,[1]],x[9,[1]]):


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

mu    := 0.0617:
kappa := 1.245:
alpha := 0.0483:

f0 := s -> 1 + mu*s^2*exp(-alpha*s^2) + kappa*(1 - exp(-1/2*alpha*s^2)):
f  := x -> f0(X2S*x):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
# Velocity Calculation for the Robot based on MDH frames
# Introduction
# Berechnung der Geschwindigkeit von Koordinatensystemen und Schwerpunkten
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb_twist -> floating base wird durch base twist (Geschwindigkeit der Basis) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# velocity -> Berechnung der Geschwindigkeit aller Segmente
# linkframe -> Berechnung der Geschwindigkeit im Körper-KS (KSi)
# 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":
printf("%s. Generiere Geschwindigkeit für %s (Herleitung im Körper-KS)\n", FormatTime("%Y-%m-%d %H:%M:%S"), robot_name):
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: # nur zum Abschätzen der Komplexität
;
# Term-Vereinfachungen einstellen
if not assigned(simplify_options) or simplify_options(4)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(4): # vierter Eintrag ist für Geschwindigkeit
end if:

# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf := Trf:
Trf_c := Trf_c:
# Zeitableitungen der MDH-Drehwinkel laden.
# Die Berechnung soll nur an einer Stelle erfolgen. Siehe robot_tree_velocity_mdh_angles.mw.
read sprintf("../codeexport/%s/tmp/velocity_mdh_angles_maple.m", robot_name):
thetaD := thetaD:

# Calculate Velocities
# First assume fixed base model with base velocity and acceleration set to zero
# Anfangsgeschwindigkeit definieren für floating base model
# Velocities of Frames
#
# Analog zu [Ortmaier2014], Gl. (7.6)
# Geschwindigkeit der Basis der Koordinatensysteme der Körper (ausgedrückt im Körper-KS)
rD_i_i := Matrix(3, NL):
# Winkelgeschwindigkeit von Körper i ausgedrückt im Körper-KS
omega_i_i := Matrix(3, NL):
# Anfangsgeschwindigkeit der Basis:
# twist: Basis-Geschwindigkeit bzgl Welt-KS ausgedrückt im Basis-KS
# eulxyz: V_base_t beinhaltet die Geschwindigkeit der Basis im Welt-KS, ausgedrückt im Welt-KS. Daher ist für eine Darstellung im Körper-KS noch die Rotation erforderlich.
if base_method_name = "twist" then:
  rD_i_i(1..3,1) := V_base_t(1..3,1):
  omega_i_i(1..3,1) := V_base_t(4..6,1):
end:
if base_method_name = "eulxyz" then:
  rD_i_i(1..3,1) := Transpose(Trf_c(1..3, 1..3, 1)) . V_base_t(1..3,1):
  omega_i_i(1..3,1) := Transpose(Trf_c(1..3, 1..3, 1)) . T_basevel . V_base_t(4..6,1):
end:
# Erhöhe die Geschwindigkeit jedes Körpers
# Betrachte dazu nur die durch Gelenke angetriebenen Körper, nicht die Basis
for i from 1 to NJ do # Gelenke durchgehen
  # Körper der von Gelenkwinkel i bewegt wird: Körperindex i+1
  # Vorgängerkörper bestimmen
  j := v(i) + 1:
  # Geschwindigkeit des Vorgängers; Trf_c(...,i+1) enthält z-Achse des Körperkoordinatensystems, das von qi bewegt wird.
  # [Ortmaier2014] (7.7) (S.115) (dort falsche Indizes für MDH), [KhalilDombre2002] (9.14)
  R_j_i := Trf(1..3,1..3,i): # Rotation vom Vorgänger-Körper (j) zu diesem Körper (i+1)
  R_i_j := Transpose(R_j_i):
  # [GautierKhalil1988], equ.7: omega_jj aus [GautierKhalil1988] entspricht omega_i_i(1 .. 3, i+1) hier
  if sigma(i) = 0 then # Drehgelenk
    omega_i_i(1 .. 3, i+1) := Multiply(R_i_j,Matrix(3,1,omega_i_i(1 .. 3, j))) + thetaD(i,1)*<0;0;1>:
  else: # Schubgelenk
    omega_i_i(1 .. 3, i+1) := Multiply(R_i_j,Matrix(3,1,omega_i_i(1 .. 3, j)))
  end if:
  # Terme vereinfachen (Teil 1)
  if use_simplify>=1 then
    tmp_t11:=time():
    tmp_l11 := length(omega_i_i(1 .. 3, i+1)):
    omega_i_i(1 .. 3, i+1) := simplify2(omega_i_i(1 .. 3, i+1)):
    tmp_l21 := length(omega_i_i(1 .. 3, i+1)):
    tmp_t21:=time():
  end if:
  # Vektor vom Ursprung des vorherigen Koordinatensystems zu diesem KS
  r_j_j_i := Trf(1 .. 3, 4, i):
  # [GautierKhalil1988], equ.8: v_jj aus [GautierKhalil1988] entspricht rD_i_i(1 .. 3, i+1) hier
  if sigma(i) = 0 then # Drehgelenk
    rD_i_i(1 .. 3, i+1) := Multiply( R_i_j, ( rD_i_i(1 .. 3, j) + CrossProduct(omega_i_i(1 .. 3, j), r_j_j_i) ) ):
  else: # Schubgelenk
    rD_i_i(1 .. 3, i+1) := Matrix(Multiply( R_i_j, ( rD_i_i(1 .. 3, j) + CrossProduct(omega_i_i(1 .. 3, j), r_j_j_i) ) ) ) + Matrix(dD(i,1)*<0;0;1>):
  end if:
  printf("%s. Geschwindigkeit für Körperkoordinatensystem %d aufgestellt (Herleitung im Körper-KS).\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), i-1): #0=Basis
  # Terme vereinfachen (Teil 2)
  if use_simplify>=1 then
    tmp_t12:=time():
    tmp_l12 := length(rD_i_i(1 .. 3, i+1)):
    rD_i_i(1 .. 3, i+1) := simplify2(rD_i_i(1 .. 3, i+1)):
    tmp_l22 := length(rD_i_i(1 .. 3, i+1)):
    tmp_t22:=time():
    printf("%s: Terme für Geschwindigkeiten vereinfacht. Länge: %d->%d / %d->%d. Rechenzeit %1.1fs und %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l11, tmp_l21, tmp_l12, tmp_l22, tmp_t21-tmp_t11, tmp_t22-tmp_t12):
  end if:
end do:
# Export
# Maple Export
save omega_i_i, rD_i_i, sprintf("../codeexport/%s/tmp/velocity_linkframe_floatb_%s_maple.m", robot_name, base_method_name):
printf("%s. Maple-Ausdrücke exportiert.\n", FormatTime("%Y-%m-%d %H:%M:%S")):

# Matlab Export
if codegen_act then
  MatlabExport(convert_t_s(omega_i_i), sprintf("../codeexport/%s/tmp/velocity_omegaii_floatb_%s_linkframe_matlab.m", robot_name, base_method_name), codegen_opt):
  MatlabExport(convert_t_s(rD_i_i), sprintf("../codeexport/%s/tmp/velocity_rDii_floatb_%s_linkframe_matlab.m", robot_name, base_method_name), codegen_opt):
  printf("%s. Geschwindigkeiten in Matlab exportiert.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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_c_tca.mpl"

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

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

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


text:
# Taken from
# N. Tuncer, T. Le
# "Structural and practical identifiability analysis of outbreak models"
# https://doi.org/10.1016/j.mbs.2018.02.004
# Equation (2.2) with cumulative incidence observations
read "../ComputeIdentifiableFunctions.mpl";

# Prevalence observations
model := [
  diff(S(t), t) = b * S(t) * In(t) / N,
  diff(E(t), t) = b * S(t) * In(t) / N - nu * E(t),
  diff(In(t), t) = nu * E(t) - a * In(t),
  y(t) = In(t),
  y2(t) = N
]:

me := MultiExperimentIdentifiableFunctions(model, simplified_generators=true):
printf("The bound for the number of experiments is %a, this means that the fields of single-experiment and multi-experiment identifiable functions coincide, and they are equal to: %a\n", me[1], me[3]):

# Cummulative incidence observations
# using change of variables Ninv := 1 / N to simplify the computation
model := [
  diff(S(t), t) = b * S(t) * In(t) * Ninv,
  diff(E(t), t) = b * S(t) * In(t) * Ninv - nu * E(t),
  diff(In(t), t) = nu * E(t) - a * In(t),
  diff(Cu(t), t) = nu * E(t),
  y(t) = Cu(t),
  y2(t) = Ninv
]:

me := MultiExperimentIdentifiableFunctions(model, simplified_generators=true, infolevel=1):
printf("The bound for the number of experiments is %a, this means that the field of multi-experiment identifiable functions is the same, %a\n", me[1], me[3]):

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


text:
p := 5;

# b(i) is the image in H_{2i}(QS^0) of the standard generator of H_{2i}(BC_p)
# This is zero unless i is divisible by p-1.  An expression like b(8,4,6)
# represents the circle product of b(8), b(4) and b(6).

b := proc()
 if map(mods,{0,args},p-1) <> {0} or min(0,args) < 0 then
  return(0);
 fi;
 'b'(op(sort([args])));
end:

# Signature of permutations of {1,...,n}, represented as lists of values.
# Sends non-permutations to zero.
signature := proc(s)
 local i,j,x;
 x := 1;
 for i from 1 to nops(s)-1 do
  for j from i+1 to nops(s) do
   if op(i,s) > op(j,s) then
    x := -x;
   elif op(i,s) = op(j,s) then
    return(0);
   fi;
  od;
 od;
 x;
end:

# a(i) is the image in H_{2i+1}(QS^0) of the standard generator of H_{2i+1}(BC_p)
# This is zero unless i+1 is divisible by p-1.  An expression like a(11,3,7)
# represents the circle product of a(11), a(3) and a(7).
a := proc()
 if map(modp,{-1,args},p-1) <> {p-2} or min(0,args) < 0 then
  return(0);
 fi;
 signature([args]) * 'a'(op(sort([args])));
end:

# The elements b(i) satisfy \sum b(i,j) s^i t^j=\sum b(i,j) s^i (s+t)^j.
# By expanding this out we see that b_rel(i,j) = 0 for all i and j,
# where b_rel(i, j) is as defined below.

b_rel := proc(i,j)
 add(mods(binomial(i+j-u,j),p) * b(u,i+j-u),u=0..i-1);
end:

b_rels := proc(d) [seq(b_rel(i,d-i),i=0..d)]; end:

# The elements a(i) satisfy \sum a(i,j) s^i t^j=\sum a(i,j) s^i (s+t)^j.
# By expanding this out we see that a_rel(i,j) = 0 for all i and j,
# where a_rel(i, j) is as defined below.

a_rel := proc(i,j)
 add(mods(binomial(i+j-u,j),p) * a(u,i+j-u),u=0..i-1);
end:

a_rels := proc(d) [seq(a_rel(i,d-i),i=0..d)]; end:

# Degree function for expressions in a's, b's, P's and Q's.
# Should be refactored to use apply_deg
# Not sure what the P's and Q's are
deg := proc(x)
 local d,n;
 if type(x,`+`) then
  d := map(deg,{op(x)});
  if nops(d) = 1 then
   return(op(1,d));
  else
   error("inhomogeneous sum");
  fi;
 elif type(x,`*`) then
  return(`+`(op(map(deg,[op(x)]))));
 elif type(x,function) and op(0,x) = o then
  return(`+`(op(map(deg,[op(x)]))));
 elif type(x,integer) then
  return(0);
 elif type(x,function) and op(0,x) = a then
  return(`+`(op(map(i->2*i+1,[op(x)]))));
 elif type(x,function) and op(0,x) = b then
  return(`+`(op(map(i->2*i,[op(x)]))));
 elif type(x,function) and op(0,x) = P then
  n := nops(x)/2;
  return( -2*(p-1)*add(op(2*i,x),i=1..n) - add(op(2*i-1,x),i=1..n));
 elif type(x,function) and op(0,x) = Q then
  n := nops(x)/2;
  return( 2*(p-1)*add(op(2*i,x),i=1..n) - add(op(2*i-1,x),i=1..n));
 fi;
end:

# Circle product.  Shoud be refactored to use apply_*
o := proc()
 local xx;
 xx := map(mods,[args],p);
 if member(0,xx) then
  return(0);
 else
  return('o'(op(xx)));
 fi;
end:

# Bockstein operation.  Recognises star and circle products, a's, b's and Q's
beta := proc(x)
 local xx,n,s,y,i,j;
 if type(x,`+`) then
  return(map(beta,x));
 elif type(x,`*`) then
  xx := [op(x)];
  n := nops(xx);
  s := 1;
  y := 0;
  for i from 1 to n do
   y := y + s * mul(xx[j],j=1..(i-1)) * beta(xx[i]) * mul(xx[j],j=i+1..n);
   s := s * (-1)^deg(xx[i]);
  od;
  return(mods(expand(y),p));
 elif type(x,function) and op(0,x) = o then
  xx := [op(x)];
  n := nops(xx);
  s := 1;
  y := 0;
  for i from 1 to n do
   y := y + s * o(seq(xx[j],j=1..(i-1)),beta(xx[i]),seq(xx[j],j=i+1..n));
   s := s * (-1)^deg(xx[i]);
  od;
  return(mods(expand(y),p));
 elif type(x,integer) then
  return(0);
 elif type(x,function) and op(0,x) = a then
  return(0);
 elif type(x,function) and op(0,x) = b then
  xx := [op(x)];
  n := nops(xx);
  y := 0;
  for i from 1 to n do
   if xx[i] > 0 then
    y := y + o(a(xx[i]-1),b(seq(xx[j],j=1..i-1),seq(xx[j],j=i+1..n)));
   fi;
  od;
  return(mods(expand(y),p));
 elif type(x,function) and op(0,x) = Q then
  if nops(x) = 0 or op(1,x) = 1 then return(0); fi;
  return(Q(1,op(2..-1,x)));
 else
  return('beta'(x));
 fi;
end:

# Dyer-Lashof operation.  Only knows linearity and the Cartan rule
QQ := proc(i,x)
 local y,z;
 if type(x,`+`) then
  return(map2(QQ,i,x));
 elif type(x,`*`) then
  y := op(1,x);
  z := `*`(op(2..-1,x));
  if type(y,integer) then
   return(mods(expand(y*QQ(i,z)),p));
  else
   return(mods(expand(add(QQ(j,y)*QQ(i-j,z),j=0..i)),p));
  fi;
 elif type(x,integer) then
  if i = 0 then
   return(x);
  else
   return(0);
  fi;
 else
  return('QQ'(i,x));
 fi;
end:

# Steenrod operation.  Knows the Cartan rule for star and circle products,
# and behaviour on a's, b's and Q's.
PP := proc(i,x)
 local y,z,j,s;
 if i<0 then
  return(0);
 elif i=0 then
  return(x);
 elif type(x,`+`) then
  return(map2(PP,i,x));
 elif type(x,`*`) then
  y := op(1,x);
  z := `*`(op(2..-1,x));
  if type(y,integer) then
   return(mods(expand(y*PP(i,z)),p));
  else
   return(mods(expand(add(PP(j,y)*PP(i-j,z),j=0..i)),p));
  fi;
 elif type(x,function) and op(0,x) = o then
  y := op(1,x);
  z := o(op(2..-1,x));
  if type(y,integer) then
   return(mods(expand(y * PP(i,z)),p));
  else
   return(mods(expand(add(o(PP(j,y),PP(i-j,z)),j=0..i)),p));
  fi;
 elif type(x,integer) then
  if i = 0 then
   return(x);
  else
   return(0);
  fi;
 elif type(x,function) and op(0,x) = b then
  j := op(1,x);
  if nops(x) = 1 then
   return(mods(binomial(j-(p-1)*i,i),p)*b(j-(p-1)*i));
  else
   z := b(op(2..-1,x));
   return(mods(expand(add(binomial(j-(p-1)*k,k)*o(b(j-(p-1)*k),PP(i-k,z)),k=0..j/(p-1))),p));
  fi;
 elif type(x,function) and op(0,x) = a then
  j := op(1,x);
  if nops(x) = 1 then
   return(mods(binomial(j-(p-1)*i,i),p)*a(j-(p-1)*i));
  else
   z := a(op(2..-1,x));
   return(mods(expand(add(binomial(j-(p-1)*k,k)*o(a(j-(p-1)*k),PP(i-k,z)),k=0..j/(p-1))),p));
  fi;
 elif type(x,function) and op(0,x) = Q then
  if nops(x) = 0 then
   if i=0 then return(x) else return(0); fi;
  fi;
  s := op(2,x);
  y := Q(op(3..-1,x));
  if op(1,x) = 0 then
   return(add((-1)^(i+k)*binomial((s-i)*(p-1),i-p*k)*oQ(Q(0,s-i+k),PP(k,y)),k=0..i/p));
  else
   return(zap(
    add((-1)^(i+k)*binomial((s-i)*(p-1)-1,i-p*k)*oQ(Q(1,s-i+k),PP(k,y)),k=0..i/p) +
    add((-1)^(i+k)*binomial((s-i)*(p-1)-1,i-p*k-1)*oQ(Q(0,s-i+k),PP(k,beta(y))),k=0..i/p)
   ));
  fi;
 else
  return('PP'(i,x));
 fi;
end:

# Length function for Steenrod or Dyer-Lashof words.
len := proc(x)
 if type(x,function) and (op(0,x) = P or op(0,x) = Q) then
  return(nops(x)/2);
 else
  return('len'(x));
 fi;
end:

# Excess for Dyer-Lashof words
excess := proc(x)
 local n;
 if type(x,function) and op(0,x) = Q then
  n := nops(x)/2;
  if n = 0 then return(0); fi;
  return(2*op(2,x) - op(1,x) - add(2*(p-1)*op(2*i,x)-op(2*i-1,x),i=2..n));
 else
  return('excess'(x));
 fi;
end:

# Admissibility criterion for Dyer-Lashof words
is_admissible := proc(x)
 local i,n;
 if type(x,function) and op(0,x) = Q then
  n := nops(x)/2;
  if n = 0 then
   return(true);
  fi;
  if op(2,x) < op(1,x) then
   return(false);
  fi;
  for i from 2 to n do
   if p*op(2*i,x) - op(2*i-1,x) < op(2*i-2,x) then
    return(false);
   fi;
  od;
  return(true);
 else
  return('is_admissible'(x));
 fi;
end:

Q_adem := proc()
 local i,j,u,ep,r,dl,s,v;
 if nargs = 0 then
  return(Q());
 fi;
 if args[2] < args[1] then
  # (we must have args[1] = 1, args[2] = 0)
  if nargs = 2 then
#   error("naked beta");
   return(Q(1,0));
  fi;
  if args[3] = 1 then
   return(0);
  else
   return(Q(1,args[4..-1]));
  fi;
 fi;
 for i from 2 to nargs/2 do
  u  := args[1..2*i-4];
  ep := args[2*i-3];
  r  := args[2*i-2];
  dl := args[2*i-1];
  s  := args[2*i  ];
  v  := args[2*i+1..-1];
  if r > p*s-dl then
   if dl = 0 then
    return(
     mods(add(
      (-1)^(r+j)*binomial((p-1)*(j-s)-1,p*j-r)*
       Q(u,ep,r+s-j,0,j,v),
      j=ceil(r/p)..(r-(p-1)*s-1)
     ),p)
    );
   else
    if ep = 0 then
     return(mods(add(
      (-1)^(r+j)*binomial((p-1)*(j-s),p*j-r)*
       Q(u,1,r+s-j,0,j,v),
      j=ceil(r/p)..(r-(p-1)*s)
     ),p) -
     mods(add(
      (-1)^(r+j)*binomial((p-1)*(j-s)-1,p*j-r-1)*
       Q(u,0,r+s-j,1,j,v),
       j = ceil((r+1)/p)..(r-(p-1)*s)
     ),p));
    else
     return(mods(add(
      (-1)^(r+j+1)*binomial((p-1)*(j-s)-1,p*j-r-1)*
       Q(u,1,r+s-j,1,j,v),
       j = ceil((r+1)/p)..(r-(p-1)*s)
     ),p));
    fi;
   fi;
  fi;
 od;
 return(Q(args));
end:

Q_exc0 := proc()
 if excess(Q(args)) < 0 then
  return(0);
 else
  return(Q(args));
 fi;
end:

Q_comp := proc(a,b)
 local aa,bb,n,m,i,d;
 aa := a;
 bb := b;
 if type(aa,`*`) then
  aa := select(type,a,specfunc(integer,Q));
 fi;
 if type(bb,`*`) then
  bb := select(type,b,specfunc(integer,Q));
 fi;
 if not(type([aa,bb],[specfunc(integer,Q)$2])) then
  error("invalid arguments");
 fi;
 n := nops(aa)/2;
 m := nops(bb)/2;
 for i from 1 to min(n,m) do
  d := 2*op(2*i,aa) - op(2*i-1,aa) - 2*op(2*i,bb) + op(2*i-1,bb);
  if d < 0 then
   return(true);
  elif d > 0 then
   return(false);
  fi;
 od;
 if n<m then return(true); fi;
 return(false);
end:

Q_terms := proc(x)
 if type(x,specfunc(integer,Q)) then
  return([x]);
 elif type(x,`*`) then
  return(select(type,[op(x)],specfunc(integer,Q)));
 elif type(x,`+`) then
  return(map(op,map(Q_terms,[op(x)])));
 else
  return([]);
 fi;
end:

Q_bot := proc(x)
 local t;
 t := sort(Q_terms(x),Q_comp);
 if t = [] then
  return(0);
 else
  return(t[1]);
 fi;
end:

Q_top := proc(x)
 local t;
 t := sort(Q_terms(x),Q_comp);
 if t = [] then
  return(0);
 else
  return(t[nops(t)]);
 fi;
end:

oQ := proc(a,b)
 local n,aa,bb;
 if type(a,`+`) then
  return(map(oQ,a,b));
 elif type(b,`+`) then
  return(map2(oQ,a,b));
 fi;
 if type(a,`*`) then
  n,aa := selectremove(type,a,integer);
  if n<>1 then
   return(expand(n*oQ(aa,b)));
  fi;
 fi;
 if type(b,`*`) then
  n,bb := selectremove(type,b,integer);
  if n<> 1 then
   return(expand(n*oQ(a,bb)));
  fi;
 fi;
 if type(a,integer) or type(b,integer) then
  return(expand(a*b));
 fi;
 if type(a,specfunc(integer,Q)) and type(b,specfunc(integer,Q)) then
  return(Q(op(a),op(b)));
 fi;
end:

zap := proc(x) mods(eval(subs(Q=Q_exc0,eval(subs(Q = Q_adem,x)))),p); end:


text:
INF = #7ff0000000000000
NAN = #7ff8000000000000
EPS = #3cb0000000000000

TRUE = 1
FALSE = 0

E = 2.71828182845904523536
PI = 3.141592653589793236462643
GAMMA = 0.577215664901632860606512 ~ Euler-Mascheroni constant
I = (0,1)

function string(x)
    return x*" "
end

function value(x)
    return x
end

REAL_TYPE = 1
COMPLEX_TYPE = 2
VECTOR_TYPE = 3
CVECTOR_TYPE = 4
MATRIX_TYPE = 5
CMATRIX_TYPE = 6
STRING_TYPE = 7
LIST_TYPE = 8

function isreal(x)
    return type(x) == REAL_TYPE
end

function iscomplex(x)
    return type(x) == COMPLEX_TYPE
end

function isvector(x)
    return type(x) == VECTOR_TYPE
end

function iscvector(x)
    return type(x) == CVECTOR_TYPE
end

function ismatrix(x)
    return type(x) == MATRIX_TYPE
end

function iscmatrix(x)
    return type(x) == CMATRIX_TYPE
end

function isstring(x)
    return type(x) == STRING_TYPE
end

function islist(x)
    return type(x) == LIST_TYPE
end

function rows(m)
    s = size(m)
    return s[0]
end

function cols(m)
    s = size(m)
    return s[1]
end

function left(s, n)
    if n <= size(s) then
        return s
    end
    return s + (n-size(s))*" "
end

function right(s, n)
    if n <= size(s) then
        return s
    end
    return (n-size(s))*" " + s
end

function center(s, n)
    if n <= size(s) then
        return s
    end
    d = n - size(s)
    l = round(d/2)
    r = d - l
    return l*" " + s + r*" "
end

function min(x, y)
    if x <= y then
        return x
    end
    return y
end

function max(x, y)
    if x >= y then
        return x
    end
    return y
end

function cis(x)
    return (cos(x), sin(x))
end

function pow(x, y)
    return x^y
end

function csc(x)
    return 1/sin(x)
end

function sec(x)
    return 1/cos(x)
end

function cot(x)
    return 1/tan(x)
end

function acsc(x)
    return asin(1/x)
end

function asec(x)
    return acos(1/x)
end

function acot(x)
    return atan(1/x)
end

function csch(x)
    return 1/sinh(x)
end

function sech(x)
    return 1/cosh(x)
end

function coth(x)
    return 1/tan(x)
end

function acsch(x)
    return asinh(1/x)
end

function asech(x)
    return acosh(1/x)
end

function acoth(x)
    return atanh(1/x)
end

function mag(x)
    m = 0
    for i in 0:length(x)
        m = m + x[i]*x[i]
    end
    return sqrt(m)
end

function cross(x, y)
    return [x[1]*y[2] - x[2]*y[1], x[2]*y[0] - x[0]*y[2], x[0]*y[1] - x[1]*y[0]]
end

function islower(s)
    c = s[0]
    if "a" <= c and c <= "z" then
        return TRUE
    end
    return FALSE
end

function isupper(s)
    c = s[0]
    if "A" <= c and c <= "Z" then
        return TRUE
    end
    return FALSE
end

function tolower(s)
    rv = ""
    for i in 0:size(s)-1
        if isupper(s[i]) then
	    rv = rv + char(ascii(s[i]) + 32)
	else
	    rv = rv + s[i]
	end
    end
    return rv
end

function toupper(s)
    rv = ""
    for i in 0:size(s)-1
        if islower(s[i]) then
	    rv = rv + char(ascii(s[i]) - 32)
	else
	    rv = rv + s[i]
	end
    end
    return rv
end

max_stars_repo_name: yugyfoog/MPL

max_issues_repo_name: yugyfoog/MPL

max_forks_repo_name: yugyfoog/MPL


text:
######################################################################
# W(N)(A) is Map(A,RR^N) modulo translation, except that we take
# W(N)(emptyset) to be empty.
#
# 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)
 local a;
 global reason;

 if A = {} then
  reason := [convert(procname,string),"W(N)(emptyset) = emptyset"];
  return false;
 fi;

 if not(type(x,table)) then
  reason := [convert(procname,string),"x is not a table",x];
  return false;
 fi;

 if map(op,{indices(x)}) <> A then
  reason := [convert(procname,string),"x is not indexed by A",x,A];
  return false;
 fi;

 for a in A do
  if not `is_element/R`(N)(x[a]) then
   reason := [convert(procname,string),"x[a] is not in R^N",a,x[a],N];
   return false;
  fi;
 od;

 return true;
end;

`is_equal/W` := (N::posint) -> (A::set) -> proc(x,y)
 local a0,u0,a;
 global reason;

 if A = {} then
  reason := [convert(procname,string),"W(N)(emptyset) = emptyset"];
  return false;
 fi;

 a0 := A[1];
 u0 := x[a0] -~ y[a0];

 for a in A do
  if not(`is_equal/R`(N)(x[a] -~ y[a],u0)) then
   reason := [convert(procname,string),"x[a] - y[a] <> x[a0] - y[a0]",a,a0];
   return false;
  fi;
 od;

 return true;
end;

`is_leq/W` := NULL;

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

 x := table();
 for a in A do
  x[a] := `random_element/R`(N)();
 od:

 return eval(x);
end;

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


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

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

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

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

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

XX := s -> 1/(1 + params_a_alpha*s^2):
ff := s -> XX(s) + params_a_gamma*(1 - XX(s)):

f := (rs, z, xt, xs0, xs1) -> f_pw(rs, z)*(
  + opz_pow_n( z,1)/2 * ff(X2S*xs0)
  + opz_pow_n(-z,1)/2 * ff(X2S*xs1)
):


text:
# Konvertierung Vektor zu symmetrischer Matrix
# Einleitung
# Generiere die Umwandlung von Vektor mit oberem rechten Teil einer Symmetrischen Matrix zur Matrix selbst für die Dimensionen des Roboters.
# Diese Umwandlung verursacht weniger Rechenoperationen als die Verwendung der numerischen Funktion vec2symmat in Matlab, da dort jedes Mal die Indizes beim Zusammenbauen der Matrix getestet werden müssen.
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-12
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# 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):
# Einstellungen für Code-Export: Optimierungsgrad (2=höchster) und Aktivierung jedes Terms.
read "../helper/proc_MatlabExport":
read "../helper/proc_vector2symmat":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
printf("Generiere Symmat2Vector-Funktionen für %s\n", robot_name):
codegen_opt := 0: # Soll nicht von Einstellung in robot_env überschrieben werden.
;
# Funktion symmat2vector für den Roboter definieren
# Erstelle eine Dummy-Variable (mv), die als temporäre Variable in Matlab dient (zum Zusammensetzen der Matrix).
clear('mv'):
M_NQ:= vec2symmat(mv, NQ):
M_NQJ:= vec2symmat(mv, NQJ):
M_6:= vec2symmat(mv, 6):
# Warnungen bei Code-Generierung unterdrücken. Die Meldung das der Ausdruck mv() in Matlab nicht bekannt ist, spielt keine Rolle, da diese Variable nach dem Einsetzen des Codes vorher definiert sein wird.
interface(warnlevel=0):
MatlabExport(M_NQ, sprintf("../codeexport/%s/tmp/vec2symmat_%d_matlab.m", robot_name, NQ), codegen_opt):
MatlabExport(M_NQJ, sprintf("../codeexport/%s/tmp/vec2symmat_%d_matlab.m", robot_name, NQJ), codegen_opt):
MatlabExport(M_6, sprintf("../codeexport/%s/tmp/vec2symmat_%d_matlab.m", robot_name, 6), codegen_opt):


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

sys := [
diff(s(t), t) = -b1*b0*s(t)*x1(t)*i(t) - b0*s(t)*i(t) - s(t)*mu + mu + g*r(t),
diff(x2(t), t) = M*x1(t),
diff(i(t), t) = -nu*i(t) + b1*b0*s(t)*x1(t)*i(t) + b0*s(t)*i(t) - mu*i(t),
diff(r(t), t) = nu*i(t) - mu*r(t) - g*r(t),
diff(x1(t), t) = -M*x2(t),
y1(t) = i(t),
y2(t) = r(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
# Euler found rational curves on the surface given by
# w+x+y+z = wxyz + a = 0.  The formula here is similar to Euler's
# but simpler, and is taken from slides of Elkies

`f/euler_quartic` := proc(a,x) local i; [add(x[i],i=1..4),mul(x[i],i=1..4)+a]; end:

`g/euler_quartic` := (a,s) -> [
 1/2*s^(-3)*(s^4-4*a)^( 2)*(s^4+12*a)^(-1),
 2*a*s^(-3)*(s^4-4*a)^(-1)*(s^4+12*a)^(-1)*(3*s^4+4*a)^2,
 1/2*s     *               (s^4+12*a)     *(3*s^4+4*a)^(-1),
  -2*s^5   *(s^4-4*a)^(-1)*(s^4+12*a)     *(3*s^4+4*a)^(-1)
];

check_euler_quartic := proc()
 local a,s;

 _ASSERT(
  simplify(`f/euler_quartic`(a,`g/euler_quartic`(a,s))) = [0,0],
  "Euler quartic relation"
 );
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: lda_exc *)
(* replace: "int1\(" -> "xc_integrate(func1, NULL, 1e-20, " *)
(* replace: "int2\(" -> "xc_integrate(func2, NULL, 1e-20, " *)
(* prefix:
  lda_x_1d_exponential_params *params;

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

$define xc_dimensions_1d

`diff/int1` := proc(g, x) diff(g, x) * x1d_inter(g)   end proc:
`diff/int2` := proc(g, x) diff(g, x) * x1d_inter(g)*g end proc:

x1d_inter := x -> xc_E1_scaled(x^2):

x1d_R := rs -> Pi*params_a_beta/(2*rs):
x1d_fs := (rs, z) ->
  -((1 + z)*int1((1 + z)*x1d_R(rs)) - int2((1 + z)*x1d_R(rs))/x1d_R(rs))
    / (4.0*Pi*params_a_beta):

x1d_f := (rs, z) ->
  + my_piecewise3(screen_dens_zeta(rs,  z), 0, x1d_fs(rs, z_thr( z)))
  + my_piecewise3(screen_dens_zeta(rs, -z), 0, x1d_fs(rs, z_thr(-z)))
:

f := (rs, z) -> x1d_f(rs, z):


text:
# An extended Stasheff tree is like a standard Stasheff tree, except that we
# are allowed vertices with only one child, which we call nodes.  By deleting
# the nodes we obtain a standard Stasheff tree TT.  For each T in TT we let
# m[T] be the number of nodes lying just below T.  The function m : TT -> N
# is arbitrary and determines the original tree.  We will write n for the
# number of leaves (which are never nodes) and k for the number of vertices
# that are neither leaves nor the root.  Note that the root may or may not
# be a node.
#
# We represent an extended Stasheff tree as the set of pairs [T,m[T]].

`is_element/extended_stasheff_trees` := (n::posint,k::nonnegint) -> proc(UU)
 global reason;
 local TT,T,U,m,kk;

 if not (type(UU,set)) then
  reason := ["is_element/extended_stasheff_trees","UU is not a set",UU];
  return false;
 fi;

 kk := 0;
 TT := NULL;

 for U in UU do
  if not(type(U,list) and nops(U) = 2) then
   reason := ["is_element/extended_stasheff_trees","U is not a list of length 2",U,UU];
   return false;
  fi;
  T,m := op(U);
  if not(type(m,nonnegint)) then
   reason := ["is_element/extended_stasheff_trees","m is not a natural number",m,U,UU];
   return false;
  fi;
  kk := kk + m + 1;
  TT := TT,T;
 od:

 TT := {TT};

 if not(`is_element/standard_stasheff_trees`(n)(TT)) then
  reason := ["is_element/extended_stasheff_trees","TT is not a standard Stasheff tree",TT,reason];
  return false;
 fi;

 kk := kk - n - 1;

 if kk <> k then
  reason := ["is_element/extended_stasheff_trees","wrong number of internal vertices",kk,k,UU];
  return false;
 fi;

 return true;
end:

`is_equal/extended_stasheff_trees` := (n::posint,k::nonnegint) -> proc(UU0,UU1)
 return evalb(UU0 = UU1);
end:

`is_leq/extended_stasheff_trees` := NULL:

`list_elements/extended_stasheff_trees` := proc(n::posint,k::nonnegint)
 local TTT,TT,T,P,p,i,j,UU,L,u;

 TTT := `list_elements/standard_stasheff_trees`(n);

 L := NULL;
 for TT in TTT do
  p := nops(TT);
  j := k - (p - n - 1);
  if j >= 0 then
   P := `list_elements/nat_partitions`(j,p);
   for u in P do
    UU := {seq([TT[i],u[i]],i=1..p)};
    L := L,UU;
   od:
  fi;
 od:
 return [L];
end:

`random_element/extended_stasheff_trees` := (n::posint,k::nonnegint) -> proc()
 local TT,m,P,UU,i;

 TT := `random_element/standard_stasheff_trees`(n)();
 while (nops(TT) > n + k + 1) do
  TT := `random_element/standard_stasheff_trees`(n)();
 od;
 m := nops(TT);

 P := `random_element/nat_partitions`(n+k+1-m,m)();
 UU := {seq([TT[i],P[i]],i=1..m)};
 return UU;
end:

`count_elements/extended_stasheff_trees` := (n::posint,k::nonnegint) ->
 binomial(n+k,n)*binomial(n+k,n-1)/(n+k);

`node_count/extended_stasheff_trees` := proc(UU)
 local U;
 return add(U[2],U in UU);
end:

# This removes the nodes from an extended Stasheff tree to give an
# ordinary Stasheff tree.

`strip/extended_stasheff_tree` := proc(UU)
 map(U -> U[1],UU);
end:

`narayana_path/extended_stasheff_tree` := proc(UU)
 local TT,n,CC,C,m,u,l,UU0;

 TT := `strip/extended_stasheff_tree`(UU);
 n := max(op(map(op,TT)));
 CC := `root_children/standard_stasheff_tree`(n)(TT);

 # Number of nodes under the root
 m := select(U -> nops(U[1])=n,UU)[1][2];

 u := 1$m;

 for C in CC do
  l := min(op(C));
  UU0 := select(U -> U[1] minus C = {},UU);
  UU0 := map(U -> [U[1] -~ (l-1),U[2]],UU0);
  u := u,1,op(`narayana_path/extended_stasheff_tree`(UU0)),-1;
 od:

 u := u,(-1)$m;
 u := [u];
 return u;
end:

# This takes a pair of extended Stasheff trees and adds an edge
# joining the roots, taking the root of the first one as the root
# of the combined tree.
`splice/extended_stasheff_tree` := (n0,k0) -> (n1,k1) -> proc(UU0,UU1)
 local UU2,R2,m,i;

 UU2 := map(U -> [U[1] +~ n0,U[2]],UU1);
 R2 := select(U -> nops(U[1]) = n1,UU2)[1];
 m := R2[2];
 if m > 0 then
  UU2 := UU2 minus {R2} union {[R2[1],m-1]};
 else
  UU2 := UU2 minus {R2};
 fi;
 return {op(UU0),op(UU2),[{seq(i,i=1..n0+n1)},0]};
end:

# This adds a new left-most leaf and connects it to the root.
# This is essentially the same as splice(ZZ,UU), where ZZ is
# the 0-corolla, which we have not allowed as an extended
# Stasheff tree.  However, it would be a bad idea to allow ZZ
# here, because sprouting adds a leaf, which is inconsistent with
# the usual effect of splicing on the leaf count.

`sprout/extended_stasheff_tree` := (n,k) -> proc(UU)
 local UU2,R2,m,i;

 UU2 := map(U -> [U[1] +~ 1,U[2]],UU);
 R2 := select(U -> nops(U[1]) = n,UU2)[1];
 m := R2[2];
 if m > 0 then
  UU2 := UU2 minus {R2} union {[R2[1],m-1]};
 else
  UU2 := UU2 minus {R2};
 fi;
 return {[{1},0],op(UU2),[{seq(i,i=1..n+1)},0]};
end:

# This adds an extra node at the bottom
# This is essentially the same as splice(UU,ZZ), where ZZ is
# the 0-corolla, which we have not allowed as an extended
# Stasheff tree.
`grow/extended_stasheff_tree` := (n::posint,k::nonnegint) -> proc(UU)
 local R;
 R := select(U -> nops(U[1]) = n,UU)[1];
 return (UU minus {R}) union {[R[1],R[2]+1]};
end:


text:
  banish := module()
    export ModuleApply :=
         proc(g, h :: name, kb :: t_kb_mb, levels :: extended_numeric,
                 x :: name, make, $)
    # banish(g, h, kb, levels, x, make), where the integrand h and the
    # integration variable x take scope over the integral g patently linear
    # in h, should be equivalent to make(kb, g), where the integration operator
    # make binds x in g, except
    #   - integration over x is performed innermost rather than outermost;
    #   - if levels < infinity then levels controls how deeply to banish g;
    #   - make is invoked with the KB in the first argument extended.
    local subintegral, w, y, kb1, lo, hi, m, loops, xx, less;
    if g = 0 or kb :: t_not_a_kb then
      0
    elif levels <= 0 then
      make(kb, g)
    elif not depends(g, x) then
      make(kb, 1) * g
    elif g :: `+` then
      map(banish, _passed)
    elif g :: `*` then
      (subintegral, w) := selectremove(depends, g, h);
      if subintegral :: `*` then error "Nonlinear integral %1", g end if;
      banish(subintegral, h, kb, levels, x, banish_weight(make, w));
    elif g :: 'And'('specfunc({Int,int,Sum,sum})',
                    'anyfunc'('anything','name'='range'('freeof'(h)))) then
      lo, hi      := op(op([2,2],g));
      m           := make;
      if depends(lo, x) then
        m  := banish_guard(eval(m), lo<y);
        lo := -infinity;
      end if;
      if depends(hi, x) then
        m  := banish_guard(eval(m), y<hi);
        hi := infinity;
      end if;
      y, kb1 := `if`(op(0,g) in '{Int,int}',
        genLebesgue(op([2,1],g), lo, hi, kb, make),
        genSummation(op([2,1],g), lo, hi, kb, make));
      subintegral := subs(op([2,1],g)=y, op(1,g));
      op(0,g)(banish(subintegral, h, kb1, levels-1, x, m), y=lo..hi)

    # if `g' is an 'n-' integral or sum of the form
    # {Ints/Sums}(f, nm, x..y, [nm0=x0..y0, .., nmk=xk..yk])
    #   where x..y is free of `h'
    elif g :: 'And'('specfunc({Ints,ints,Sums,sums})',
                    'anyfunc'('anything', 'name', 'range'('freeof'(h)),
                              'list(name=range)')) then
      lo, hi      := op(op(3,g));       # 'primary' (?) integration bounds
      loops       := op(4,g);           # rest of the integration bounds
      xx          := map(lhs, loops);   # [nm0 .. nmk]
      m           := make;
      less        := `if`(op(0,g) in '{Ints,ints}', `<`, `<=`);

      # banish the primary bounds if they depend on the primary integration
      # variable
      if depends(lo, x) then
        m  := banish_guard(m, forall(xx, less(lo, mk_idx(y,loops))));
        lo := -infinity;
      end if;
      if depends(hi, x) then
        m  := banish_guard(m, forall(xx, less(mk_idx(y,loops), hi)));
        hi := infinity;
      end if;

      # augment the KB with a new binding (y) corresponding to the primary
      # variable of integration.
      y, kb1 := genType(op(2,g),
                        mk_HArray(`if`(op(0,g) in '{Ints,ints}',
                                       HReal(open_bounds(lo..hi)),
                                       HInt(closed_bounds(lo..hi))),
                                  op(4,g)),
                        kb);

      # if there are any more integration bounds, then add the information
      # about the last (?) integration bound to the KB.
      if nops(op(4,g)) > 0 then
        kb1 := assert(size(y)=op([4,-1,2,2],g)-op([4,-1,2,1],g)+1, kb1);
          ASSERT(type(kb1,t_kb), "banish/{Ints,Sums}: invalid integration bounds found "
                 "in input.");
      end if;

      subintegral := subs(op(2,g)=y, op(1,g));
      op(0,g)(banish(subintegral, h, kb1, levels-1, x, m), y, lo..hi, op(4,g));
    elif g :: {t_pw, Partition} then
      `if`(g::t_pw,foldr_piecewise,Partition:-Foldr)(banish_pw_cons(x, h, levels), banish_pw_nil, g)(make, kb)
    # elif g ::  then
      # Partition:-Foldr(banish_pw_cons(x, h, levels), banish_pw_nil, g)(make, kb)
    elif g :: t_case then
      map_piecewiselike(banish, _passed)
    elif g :: 'toLO:-integrate(freeof(x), Integrand(name, anything), list)' then
      y           := gensym(op([2,1],g));
      subintegral := subs(op([2,1],g)=y, op([2,2],g));
      subsop(2=Integrand(y, banish(subintegral, h, kb, levels-1, x, make)), g)
    else
      make(kb, g)
    end if
  end proc;

  local banish_pw_nil := proc(make, kb, $) 0 end proc;

  local banish_pw_cons := proc(x, h, levels, $)
        proc(cond, th, el, $) proc(make, kb, $)
          if depends(cond, x) then
            banish(th, h, kb, levels-1, x, banish_guard(make, cond))
              + el(banish_guard(make, Not(cond)), kb)
          else
            piecewise(cond,
              banish(th, h, assert_mb(cond, kb), levels-1, x, make),
              el(make, assert_mb(Not(cond), kb)))
          end if
        end proc end proc;
  end proc;

  local banish_guard := proc(make, cond, $)
    if cond :: 'And(specfunc(Not), anyfunc(anything))' then
      # Work around simplify/piecewise bug:
      #   > simplify(piecewise(Not(i=0), 1, 0))
      #   a
      # (due to PiecewiseTools:-ImportImplementation:-UseSolve calling
      # solve(Not(i=0), {i}, 'AllSolutions', 'SolveOverReals'))
      proc(kb,g,$) make(kb, piecewise(op(1,cond),0,g)) end proc
    else
      proc(kb,g,$) make(kb, piecewise(cond,g,0)) end proc
    end if
  end proc;

  local banish_weight := proc(make, w, $)
    proc(kb,g,$) make(kb, w*g) end proc
  end proc;

end module;

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


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

airy_a5  := 133.983631:
airy_a6  :=   3.217063:
airy_a7  := 136.707378:
airy_a8  :=   3.223476:
airy_a9  :=   2.675484:
airy_a10 :=   3.473804:

airy_f1 := s -> (1 - airy_a5*s^airy_a6 + airy_a7*s^airy_a8)/(1 + airy_a9*s^airy_a10):

$include "gga_x_lag.mpl"

airy_f := x -> lag_f0(X2S*x) + airy_f1(X2S*x):

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


text:
(*
 Copyright (C) 2020 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: mgga_exc *)

(* This is the definition in the paper *)
th_f0 := (x, u, t) -> -27*Pi/(10*t) * (1 + 7*x^2/(108*t)):

(* Since we write this as an enhancement functional, we need to divide
   out the LDA prefactor. The paper also defines tau without one half *)
th_f := (x, u, t) -> -th_f0(x,u,2*t) / X_FACTOR_C:

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


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

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

dk_f := x ->
  add(params_a_aa[i]*x^(2*(i-1)), i=1..5) /
  add(params_a_bb[i]*x^(2*(i-1)), i=1..5):

f := (rs, zeta, xt, xs0, xs1) -> gga_kinetic(dk_f, rs, zeta, xs0, xs1):


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

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

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

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

  assert(params->N > 1);
*)

$define xc_dimensions_2d

prm_q := 3.9274: (* 2.258 *)

(* Equation (4) *)
beta := rs -> prm_q/(sqrt(Pi)*rs):

phi  := rs -> beta(rs)/(beta(rs) + sqrt(Pi)/2):

f0 := rs ->
  + sqrt(Pi)*beta(rs)*(phi(rs) - 1)/(2*sqrt(2 + params_a_c)) (* original version has (phi-1)^2 *)
  + phi(rs)*(phi(rs) - 1)/(2 + params_a_c)
  + sqrt(Pi)*phi(rs)*phi(rs)/(4*beta(rs)*(2 + params_a_c)^1.5)
  + sqrt(Pi)*beta(rs)*(phi(rs) - 1)/sqrt(1 + params_a_c)
  + phi(rs)/(1 + params_a_c):

f  := (rs, z) -> f0(rs)*Pi/(2*prm_q*prm_q):


text:
#struct aa{
#  double b;
#  float *x;
#};

#struct bb {
#  struct aa b1;
#  struct aa *pa;
#  double x1[10];
#  double dx;
#};

#typedef struct ss {
#  int f1;
#  char f2:3;
#  struct aa a3;
#  struct aa *pa1;
#  char f3:5;
#  struct bb *b3;
#  struct bb *b2;
#  struct aa xt[23];
#  struct aa(* funcc)(int, int, double);
#}SS;

type $unname1 <struct {@b f64,
                       @x <* f32>}>
type $unname2 <struct {@b1 <struct {@b f64,
                                    @x <* f32>}>,
                       @pa <* <struct {@b f64,
                                       @x <* f32>}>>,
                       @x1 <[10] f64>,
                       @dx f64}>
type $SS      <struct {@f1 i32,
                       @f2:3 i8,
                       @a3 <structincomplete {@b f64,
                                    @x <*f32>}>,
                       @pa1 <* <struct {@b f64,
                                       @x <* f32> } >>,
                       @f3:5 i32,
                       @b3 <* <struct {@b1 <struct {@b f64,
                                                    @x <* f32>}>,
                                       @pa <* <struct {@b f64,
                                                       @x <* f32>}>>,
                                       @x1 <[10] f64>,
                                       @dx f64}>>,
                       @b2 <* $unname2>,
                       @xt <[23] <struct {@b f64,
                                          @a <*f32>}>>,
                       @funcc <func(i32, i32, f64)<struct {@b f64,
                                                           @x <* f32>}>>}>

 # 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:
program sample11;
var n, sum, data : integer;
begin
	writeln('input the number of data');
	readln(n);
	sum := 0;
	while n > 0 do begin
		readln(data);
		sum := sum + data;
		n := n - 1
	end;
	writeln('Sum of data = ', sum)
end.

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing