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

lp90_c0 := 0.80569:
lp90_d0 := 3.0124e-3:
lp90_k  := 4.0743e-3:

(* Equation (60) *)
lp90_f := (rs, z, xt, us0, us1) ->
  - (lp90_c0 + lp90_d0*t_vw(z, xt, us0, us1))/(rs/RS_FACTOR + lp90_k):

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


text:
# This file contains code related to the following situation.
# We have P0 in ICP_N(A) restricting to Q0 in ICP_N(B),
# where B = A \ {a}.  We also have Q1 in ICP_N(B) refining Q0.
# We let U denote the poset of elements M in ICP_N(A) such that
# the restriction of M to B is Q1, and M refines P0.
# The claim is that |U| is contractible.

`random_element/PQQ` := (N::posint) -> (A::set,a) -> proc()
 local B,P0,Q0,Q1;
 B := A minus {a};
 P0 := `random_element/ICP`(N)(A)():
 Q0 := `res/ICP`(N)(A,B)(P0):
 Q1 := NULL:
 while Q1 = NULL or not(`is_leq/ICP`(N)(B)(Q0,Q1)) do
  Q1 := `random_element/ICP`(N)(B)():
 od:
 return [P0,Q0,Q1];
end:

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

`analyse/PQQ` := (N::posint) -> (A::set,a) -> proc(PQQ)
 local B,P0,Q0,Q1,P0r,Q0r,Q1r,L1,n1,T,M,RR,SS,US,PN,PS,U,nU,U0,UR0,
  Q1B,AQ1B,RT,Bm,Bo,Bp,Bmr,Bpr,NB,NU,VV,b,i,j,k,l,beta,kk,ee;
 B := A minus {a};
 P0,Q0,Q1 := op(PQQ);
 T := table():
 T["N"] := N;
 T["A"] := A;
 T["B"] := B;
 T["a"] := a;
 T["P0"] := P0;
 T["Q0"] := Q0;
 T["Q1"] := Q1;
 T["P0_tot"] := `totalise/ICP/ord`(N)(A)(P0);
 T["Q0_tot"] := `totalise/ICP/ord`(N)(B)(Q0);
 T["Q1_tot"] := `totalise/ICP/ord`(N)(B)(Q1);
 P0r := table():
 Q0r := table():
 Q1r := table():
 for i from 1 to N do
  P0r[i] := `rank_table/preord`(A)(P0[i]);
  Q0r[i] := `rank_table/preord`(B)(Q0[i]);
  Q1r[i] := `rank_table/preord`(B)(Q1[i]);
 od:
 T["P0_rank_table"] := eval(P0r);
 T["Q0_rank_table"] := eval(Q0r);
 T["Q1_rank_table"] := eval(Q1r);
 RR := table():
 SS := table():
 US := table():
 for b in B do
  RR[b] := `list_elements/ICP`(N)({a,b}):
  SS[b] := map(`f/fibre_sphere/ICP`(N)(A)(a,b,Q1),RR[b]):
  US[b] := select(M -> `is_leq/ICP`(N)(A)(P0,M),SS[b]):
  PN[b] := evalb(nops(US[b]) > 0);
  PS[b] := evalb(modp(nops(US[b]),2) = 0);
  k := 1;
  while member([a,b],P0[k]) and member([b,a],P0[k]) do
   k := k+1;
  od:
  kk[b] := k;
  ee[b] := `if`(member([a,b],P0[k]),-1,+1);
 od:
 T["fibre_sphere"] := eval(SS);
 T["upper_fibre_part"] := eval(US);
 T["part_is_empty"] := eval(PE);
 T["part_is_sphere"] := eval(PS);
 T["P0_split_stage"] := eval(kk);
 T["P0_split_direction"] := eval(ee);
 T["fibre"] := {seq(op(SS[b]),b in B)};
 U := {seq(op(US[b]),b in B)};
 nU := nops(U);
 U0,UR0 := op(`skeleton/partord/from_comparator`(U)(`is_leq/ICP`(N)(A)));
 T["upper_fibre"] := U;
 T["U0"] := U0;
 T["upper_fibre_size"] := nU;
 T["upper_fibre_skeleton"] := UR0;
 T["upper_fibre_hasse_diagram"] :=
  `hasse_diagram/partord`(U0)(UR0);
 L1 := T["Q1_tot"];
 n1 := nops(L1);
 T["part_list"] := [seq(US[L1[i]],i=1..n1)];
 T["part_is_nonempty_list"] := [seq(PN[L1[i]],i=1..n1)];
 T["part_is_sphere_list"] := [seq(PS[L1[i]],i=1..n1)];
 T["left_nested_list"] :=
  [FAIL,seq(evalb({op(US[L1[i]])} minus {op(US[L1[i-1]])} = {}),i=2..n1)];
 T["right_nested_list"] :=
  [seq(evalb({op(US[L1[i]])} minus {op(US[L1[i+1]])} = {}),i=1..n1-1),FAIL];
 T["part_skeleton"] := [
  seq(select(i -> member(U[i],US[L1[j]]),U0),j=1..n1)
 ];
 Q1B := table():
 AQ1B := table():
 RT := table():
 Bm := table():
 Bo := table():
 Bp := table():
 Bmr := table():
 Bpr := table():
 NU  := table():
 for i from 1 to N do
  Q1B[i] := `block_partition/equiv`(B)(Q1[i] union `op/autorel`(B)(Q1[i]));
  AQ1B[i] := select(u -> member([u[1],a],P0[i]) or member([a,u[1]],P0[i]),Q1B[i]);
  for beta in AQ1B[i] do
   RT[i,beta] := `rank_table/preord`(beta)(Q1[i]);
   NB[i,beta] := max(0,seq(RT[i,beta][b]+1,b in beta));
   Bm[i,beta] := select(b -> not(member([a,b],P0[i])),beta);
   Bo[i,beta] := select(b -> member([a,b],P0[i]) and member([b,a],P0[i]),beta);
   Bp[i,beta] := select(b -> not(member([b,a],P0[i])),beta);
   Bmr[i,beta] := max(0,seq(RT[i,beta][b]+1,b in Bm[i,beta]));
   Bpr[i,beta] := min(NB[i,beta],seq(RT[i,beta][b],b in Bp[i,beta]));
   NU[i,beta] := Bpr[i,beta] - Bmr[i,beta] + 1;
   VV[i,beta] := {seq([i,beta,l-0.5],l=Bmr[i,beta]..Bpr[i,beta])};
  od:
 od:
 Q1B[N+1] := {seq({b},b in B)};
 AQ1B[N+1] := {};
 T["Q1_comparability_blocks"] := eval(Q1B):
 T["admissible_blocks"] := eval(AQ1B);
 T["admissible_tree"] := [seq(seq([i,beta],beta in AQ1B[i]),i=1..N)];
 T["Q1_rank_table"] := eval(RT);
 T["lower_set"]  := eval(Bm);
 T["middle_set"] := eval(Bo);
 T["upper_set"]  := eval(Bp);
 T["lower_rank"] := eval(Bmr);
 T["upper_rank"] := eval(Bpr);
 T["block_quotient_size"] := eval(NB);
 T["num_U"]      := eval(NU);
 T["V"] := {seq(seq(op(VV[i,beta]),beta in AQ1B[i]),i=1..N)};
 return eval(T):
end:

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

`f/PQQ/U/V` := (T) -> proc(M)
 local a,b,j,k,beta,beta0,gamma,rt,r;
 k := 1:
 beta := T["B"]:
 a := T["a"]:
 for j from 2 to T["N"] do
  beta0 := select(b -> member([a,b],M[j]) or member([b,a],M[j]),beta);
  if beta0 <> {} then
   beta := beta0;
   k := j;
  else
   break;
  fi;
 od:
 gamma := select(b -> not(member([a,b],M[k])),beta);
 rt := eval(T["Q1_rank_table"][k,beta]):
 r := max(-1,seq(rt[b],b in gamma)) + 0.5;
 return [k,beta,r];
end:

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

`g/PQQ/V/U` := (T) -> proc(x)
 local i,k,N,a,b,B,MT,QQ,MM,M,beta,r,b0,L,U,rt;
 k,beta,r := op(x);
 N := T["N"];
 a := T["a"];
 B := T["B"];
 MT := table():
 for i from 1 to N do
  QQ := T["Q1"][i];
  MM := QQ union {[a,a]};
  if i < k then
   b0 := beta[1];
   L := select(b -> member([b,b0],QQ),B);
   U := select(b -> member([b0,b],QQ),B);
   MM := {op(MM),seq([b,a],b in L),seq([a,b],b in U)};
  elif i = k then
   rt := T["Q1_rank_table"][k,beta];
   L := select(b -> rt[b] < r,beta);
   U := select(b -> rt[b] > r,beta);
   MM := {op(MM),seq([b,a],b in L),seq([a,b],b in U)};
  fi;
  MT[i] := MM;
 od:
 M := [seq(MT[i],i=1..N)];
 return M;
end:

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

`is_leq/PQQ/V` := (T) -> proc(x0,x1)
 local k0,k1,beta0,beta1,r0,r1,r2,b0;
 k0,beta0,r0 := op(x0);
 k1,beta1,r1 := op(x1);
 if k0 < k1 then return false; fi;
 if k0 = k1 then
  return evalb(beta0 = beta1 and r0 = r1);
 fi;
 b0 := beta0[1];
 if not(member(b0,beta1)) then return false; fi;
 r2 := T["Q1_rank_table"][k1,beta1][b0];
 if not(r1 = r2 - 0.5 or r1 = r2 + 0.5) then
  return false;
 fi;
 return true;
end:

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

`rho/PQQ` := (T) -> (j,e) -> proc(x)
 local k,beta,gamma,beta1,b0,QQ,rt,r1;
 k,beta,gamma := op(x);
 if j >= k then return FAIL; fi;
 b0 := beta[1];
 QQ := T["Q1"][j];
 beta1 := select(b -> member([b,b0],QQ) or member([b0,b],QQ),T["B"]);
 rt := T["Q1_rank_table"][j,beta1];
 r1 := rt[beta[1]] + 0.5 * e;
 return [j,beta1,r1];
end:

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

`describe/PQQ` := proc(T)
 cat(`describe/ICP`(T["N"])(T["A"])(T["P0"]),"\n",
     `describe/ICP`(T["N"])(T["B"])(T["Q0"]),"\n",
     `describe/ICP`(T["N"])(T["B"])(T["Q1"]),"\n");
end:

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

`check_UV/PQQ` := proc(T)
 local N,A,U,V,ff,gg,M,M0,M1,x,x0,x1,cU,cV;
 N := T["N"];
 A := T["A"];
 U := T["upper_fibre"]:
 V := T["V"]:
 ff := table():
 gg := table():
 for M in U do
  x := `f/PQQ/U/V`(T)(M);
  if not(member(x,V)) then return false; fi;
  M1 := `g/PQQ/V/U`(T)(x);
  if M <> M1 then return false; fi;
  ff[M] := x;
 od:
 for x in V do
  M := `g/PQQ/V/U`(T)(x);
  if not(member(M,U)) then return false; fi;
  x1 := `f/PQQ/U/V`(T)(M);
  if x <> x1 then return false; fi;
  gg[x] := M;
 od:
 for x0 in V do
  for x1 in V do
   M0 := gg[x0];
   M1 := gg[x1];
   cU := `is_leq/ICP`(N)(A)(M0,M1);
   cV := `is_leq/PQQ/V`(T)(x0,x1);
   if cU <> cV then return false; fi;
  od;
 od;
 return true;
end:

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

`check_rho/PQQ` := proc(T)
 global reason;
 local VV,i,j,k,x,y,z,C,y1,y2,y3,y4,y5,y6;
 VV := table():
 for k from 1 to T["N"] do
  VV[k] := select(x -> x[1] = k,T["V"]):
 od:
 for k from 2 to T["N"] do
  for j from 1 to k-1 do
   for x in VV[k] do
    y := `rho/PQQ`(T)(j,-1)(x);
    z := `rho/PQQ`(T)(j,+1)(x);
    C := select(w -> `is_leq/PQQ/V`(T)(x,w),VV[j]);
    if (y = z) or (C <> {y,z}) then
     return false;
    fi;
   od;
  od;
 od;
 for k from 3 to T["N"] do
  for j from 2 to k-1 do
   for i from 1 to j-1 do
    for x in VV[k] do
     y1 := `rho/PQQ`(T)(i,-1)(x);
     y2 := `rho/PQQ`(T)(i,+1)(x);
     y3 := `rho/PQQ`(T)(i,-1)(`rho/PQQ`(T)(j,-1)(x));
     y4 := `rho/PQQ`(T)(i,-1)(`rho/PQQ`(T)(j,+1)(x));
     y5 := `rho/PQQ`(T)(i,+1)(`rho/PQQ`(T)(j,-1)(x));
     y6 := `rho/PQQ`(T)(i,+1)(`rho/PQQ`(T)(j,+1)(x));
     if not(y3 = y1 and y4 = y1 and y5 = y2 and y6 = y2) then
      reason := [i,j,k,x,y1,y2,y3,y4,y5,y6];
      return false;
     fi;
    od;
   od;
  od;
 od;

 return true;
end:

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

`check_UL/PQQ` := proc(T)
 global reason;
 local N,A,B,a,b,Q,kk,k,L,U,S,i,e;
 N := T["N"];
 a := T["a"];
 A := T["A"];
 B := T["B"];
 Q := T["Q1"];
 kk := T["P0_split_stage"];
 for b in B do
  k := kk[b];
  L := {seq(seq(`f_alt/fibre_sphere/ICP`(N)(A)(a,b,Q)([i,e]),e in {-1,1}),i=1..k-1)};
  U := {seq(seq(`f_alt/fibre_sphere/ICP`(N)(A)(a,b,Q)([i,e]),e in {-1,1}),i=k+1..N)};
  S := {op(T["upper_fibre_part"][b])};
  if not((U intersect S = {}) and (L minus S = {})) then
   return false;
  fi;
 od:
 return true;
end:

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

`check_is_sphere/PQQ` := proc(T)
 global reason;
 local N,A,B,a,b,P0,Q0,Q1,kk,k,ee,e,C,D,E,F,u;
 N := T["N"];
 a := T["a"];
 A := T["A"];
 B := T["B"];
 P0 := T["P0"];
 Q0 := T["Q0"];
 Q1 := T["Q1"];
 kk := eval(T["P0_split_stage"]);
 ee := eval(T["P0_split_direction"]);
 for b in B do
  k := kk[b];
  e := ee[b];
  if e = -1 then
   C := select(x -> member([x,b],Q0[k]) and member([b,x],Q0[k]),B);
   D := select(x -> member([x,b],Q1[k]) and
                    not(member([b,x],Q1[k])),C);
   E := select(x -> member([a,x],P0[k]) and
                    member([x,b],P0[k]) and
                    not(member([x,a],P0[k])) and
                    not(member([b,x],P0[k])),B);
   F := select(x -> member([x,b],Q1[k]) or member([b,x],Q1[k]),E);
  else
   C := select(x -> member([x,b],Q0[k]) and member([b,x],Q0[k]),B);
   D := select(x -> member([b,x],Q1[k]) and
                    not(member([x,b],Q1[k])),C);
   E := select(x -> member([x,a],P0[k]) and
                    member([b,x],P0[k]) and
                    not(member([a,x],P0[k])) and
                    not(member([x,b],P0[k])),B);
   F := select(x -> member([x,b],Q1[k]) or member([b,x],Q1[k]),E);
  fi;
  u := not(evalb(D = {} and F = {}));
  if u <> T["part_is_sphere"][b] then
   reason := [b,k,e,C,D,E,F]:
   return false;
  fi;
 od:
 return true;
end:

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

`check_sphere_in_ball/PQQ` := proc(T)
 global reason;
 local N,A,B,a,b,P0,Q0,Q1,US,kk,k,ee,e,C,C1,b1;
 N := T["N"];
 a := T["a"];
 A := T["A"];
 B := T["B"];
 P0 := T["P0"];
 Q0 := T["Q0"];
 Q1 := T["Q1"];
 kk := eval(T["P0_split_stage"]);
 ee := eval(T["P0_split_direction"]);
 US := eval(T["upper_fibre_part"]);
 for b in B do
  if T["part_is_sphere"][b] then
   k := kk[b];
   e := ee[b];
   if e = -1 then
    C := select(x -> member([x,b],Q1[k]) and not(member([x,a],P0[k])),B);
    if C = {} then
     reason := [b,k,e,C];
     return false;
    fi;
    # Choose b1 to be Q1[k]-minimal in C
    C1 := C;
    while C1 <> {} do
     b1 := C1[1];
     C1 := select(x -> not(member([b1,x],Q1[k])),C1);
    od:
    if T["part_is_sphere"][b1] or
       ({op(US[b])} minus {op(US[b1])} <> {}) or
       member([b,b1],Q1[k]) then
     reason := [b,k,e,C,b1];
     return false;
    fi;
   else
    C := select(x -> member([b,x],Q1[k]) and not(member([a,x],P0[k])),B);
    if C = {} then
     reason := [b,k,e,C];
     return false;
    fi;
    # Choose b1 to be Q1[k]-maximal in C
    C1 := C;
    while C1 <> {} do
     b1 := C1[1];
     C1 := select(x -> not(member([x,b1],Q1[k])),C1);
    od:
    if T["part_is_sphere"][b1] or
       ({op(US[b])} minus {op(US[b1])} <> {}) or
       member([b1,b],Q1[k]) then
     reason := [b,k,e,C,b1];
     return false;
    fi;
   fi;
  fi;
 od:
 return true;
end:

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

`check_nonempty_interval/PQQ` := proc(T)
 local L,J,i,j,j0,j1;
 L := T["part_is_nonempty_list"];
 if not(`or`(op(L))) then return false; fi;
 J := select(i -> L[i],[seq(i,i=1..nops(L))]);
 j0 := min(op(J));
 j1 := max(op(J));
 return `and`(seq(L[j],j=j0..j1));
end:

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

`check_adjacent_nonempty/PQQ` := proc(T)
 local L,P,i;
 L := T["part_is_nonempty_list"];
 P := T["part_list"];
 for i from 1 to nops(L)-1 do
  if L[i] and L[i+1] and ({op(P[i])} intersect {op(P[i+1])} = {}) then
   return false;
  fi;
 od;
 return true;
end:

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

`check_incomparable_min/PQQ` := proc(T)
 local n,LN,RN,P,i,b,P1;
 n := nops(T["B"]);
 LN := T["left_nested_list"];
 RN := T["right_nested_list"];
 P := T["part_list"];
 for i from 1 to n-1 do
  if not(LN[i+1]) and not(RN[i]) then
   b := T["Q1_tot"][i];
   P1 := `m/fibre_sphere/ICP`(T["N"])(T["A"])(a,b,T["Q1"]):
   if not (member(P1,{op(P[i])}) and member(P1,{op(P[i+1])})) then
    return false;
   fi;
  fi;
 od;
 return true;
end:

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

`check_intersections_nested/PQQ` := proc(T)
 global reason;
 local n,P,PS,i,j,k;
 n := nops(T["B"]);
 P := T["part_list"];
 for i from 1 to n do PS[i] := {op(P[i])}; od:
 for i from 1 to n-2 do
  for j from i+1 to n-1 do
   for k from j+1 to n do
    if (PS[i] intersect PS[k]) minus PS[j] <> {} then
     reason := [i,j,k];
     return false;
    fi;
   od;
  od;
 od;
 return true;
end:

ext: mpl


text:
######################################################################
# triples(A) is the set of all distinct triples (a,b,c) in A^3

`is_element/triples` := (A::set) -> proc(abc)
 type(abc,list) and nops(abc) = 3 and
 member(abc[1],A) and member(abc[2],A) and member(abc[3],A) and
 abc[1] <> abc[2] and abc[1] <> abc[3] and abc[2] <> abc[3];
end;

`is_equal/triples` := (A::set) -> (abc,def) -> evalb(abc = def):

`is_leq/triples` := NULL;

`random_element/triples` := (A::set) -> proc()
 local n,r,i,j,k,abc;
 n := nops(A);

 if n < 3 then return FAIL; fi;

 abc := [];

 while nops({op(abc)}) < 3 do
  i := rand(1..n)();
  j := rand(1..n)();
  k := rand(1..n)();
  abc := [A[i],A[j],A[k]];
 od;

 return abc;
end;

`list_elements/triples` := proc(A::set)
 local n,L,i,j,k;

 n := nops(A);
 L := NULL;

 for i from 1 to n do
  for j from 1 to n do
   for k from 1 to n do
    if i <> j and i <> k and j <> k then
     L := L,[A[i],A[j],A[k]];
    fi;
   od;
  od;
 od;

 return [L];
end:

`count_elements/triples` := (A::set) ->
 nops(A) * (nops(A) - 1) * (nops(A) - 2);


text:
(*
    交换子表计算：
    + 计算交换子的表达式形式。
    + 将交换子用生成元进行线性表出。
    + 计算伴随变换矩阵。
    + 生成偏微分方程组，并求解不变量。
*)
$ifndef _BASIC_
$define _BASIC_

$include "Logout.mpl"
$include "InvSimplify.mpl"

$ifndef _HEADERS_
PDETools:-declare(quiet):
macro(Pa=`\x26\x50\x61\x72\x74\x69\x61\x6C\x44\x3B`);# 偏导符号作为函数名
$endif

SymsHolder:=module()
    option object;
    export  default_syms:={x,y,z,t,u,v,w},      # 默认符号集
            syms:=default_syms,                 # 当前符号集
            pnames:={a,c,d,epsilon,Delta,phi};  # 受保护的名字
end module:

# 修改微分算子的符号集合
setSymbols:=proc(s::set(name):=SymsHolder:-default_syms)
    description "设置函数的变量名集合";
    if ( s intersect SymsHolder:-pnames <> {}) then
        error "变量不能包含%1",SymsHolder:-pnames;
    end if;
    SymsHolder:-syms:=s;
end proc;

# 获取变量名集合
getSymbols:=proc()
    description "获取变量名集合";
    SymsHolder:-syms;
end proc;

# 自定义微分算子操作，作用到函数f上
d:=proc()
    description "用于生成微分算子表达式";
    if not {_passed} subset SymsHolder:-syms then
        error "表达式只能包含以下变量 %1, 可以使用 setSymbols 设置变量集合，但是不能包含 %2",SymsHolder:-syms,SymsHolder:-pnames;
    end if;
    return diff(Pa(SymsHolder:-syms[]),_passed);
end proc;

# 自定义交换子计算符
`&c`:=proc(a,b)
    description "计算两个生成微分算子的交换子";
    expand(eval(subs(Pa(SymsHolder:-syms[])=b,a)-subs(Pa(SymsHolder:-syms[])=a,b)));
end proc:

(*
* 将表达式分解为非线性项并提取系数
* 输入：
*      f 表达式
* 输出：
*     T 非线性项->系数 的映射表
*)
getKd:=proc(f)
    local T:=table(),kd;

    (*
    * 将表达式分解为非线性基并提取系数的递归子函数
    * 输入：
    *     f 表达式
    *     T 保存结果的表
    * 输出：
    *     T T会被修改
    *)
    kd:=proc(f,T)
    local i,p,v,x;
    if type(f,`+`) then
        for i from 1 to nops(f) do
            thisproc(op(i,f),T);
        end do;
    elif type(f,`*`) then
        p:=1;
        v:=1;
        for i from 1 to nops(f) do
            x:=op(i,f);
            if type(x,extended_numeric) then
                p:=p*x;
            else
                v:=v*x;
            end if;
        end do;
        T[v]:=p;
    else
        T[f]:=1;
    end if;
    return;
    end proc:

    kd(f,T);
    return eval(T);
end proc:

(*
* 获取表达式关于给定非线性项集的系数向量
* 输入：
*     f    表达式
*     s    非线性项集
* 输出：
*     v    系数向量
*)
getPmVec:=proc(f,s)
    local n,v,i,tb;
    tb:=getKd(f);
    n:=numelems(s);
    v:=Vector(n);
    for i from 1 to n do
        if assigned(tb[s[i]]) then
            v[i]:=tb[s[i]];
        end if;
    end do;
    return eval(v);
end proc:

(*
* 求解表达式关于基的线性表出
* 输入：
*     f    表达式
*     A    基关于非线性项集的系数矩阵
*     stbs    非线性项集
*     sbs    基的符号表示
* 输出：
*     r    表达式关于基的线性表出，求解失败返回原表达式
*)
ans2v:=proc(f,A,stbs)
    local r;
    try
        r:=LinearAlgebra[LinearSolve](A,getPmVec(f,stbs));
    catch:
        error "所给生成元不能构成一组基";
    end try;
end proc:

(*
* 计算所有结果
* 输入：一组生成元
* 输出：
*     AD    伴随变换矩阵的数组
*     ADA    A[1]*A[2]*...*A[n]
*     dts    不变量数组
*)
getTransMatAndPDE:=proc(vv::list)
    local tbs,stbs,vvv,M,n,sbs,i,j,A,tmpv,MK,AD,ADA,ADT,BA,pPhi,eq,AList,dts,eqs;

    if not (indets(vv,name) subset SymsHolder:-syms) then
        error "表达式只能包含以下变量 %1, 可以使用 setSymbols 设置变量集合，但是不能包含 %2",
            SymsHolder:-syms,SymsHolder:-pnames;
    end if;

    vvv:=expand(vv):
    n:=numelems(vvv):
    sbs:=Vector[row](1..n,i->v[i]):# 生成元的表示符号
    flogf[2]("Input:");
    flog[2](seq(sbs[i]=vv[i],i=1..n));

    # 计算交换子矩阵，这里得到的是关于f的结果，需要进一步用基表示
    M:=Matrix(1..n, 1..n, (i, j)->vvv[i] &c vvv[j]):# 交换子的表达式形式
    MK:=Matrix(1..n,1..n);# 交换子关于生成元的系数向量矩阵

    # 将原交换子表用基表出
    tbs:=getKd~(vvv):# 非线性项及其系数映射表
    stbs:={map(indices,tbs,nolist)[]}:# 非线性项集
    # 生成元关于非线性项集的系数矩阵
    A:=Matrix(1..numelems(stbs),1..numelems(tbs),
    (i,j)->`if`(assigned(tbs[j][stbs[i]]),tbs[j][stbs[i]],0)):
    # 计算每个交换子关于生成元的系数
    # 这里将交换子的表达式转化成了生成元的线性表出
    for i from 1 to n do
        for j from 1 to n do
            if (M(i,j)<>0) then
                if (i<=j) then
                    tmpv:=ans2v(M(i,j),A,stbs);
                    M(i,j):=sbs.tmpv;
                    MK(i,j):=convert(tmpv,list);
                else
                    M(i,j):=-M(j,i);
                    MK(i,j):=-MK(j,i);
                end if;
            else
                MK(i,j):=convert(Vector[row](1..n),list);
            end if;
        end do;
    end do;
    flogf[2]("Commutator table:");
    flog[2](M);

    # 伴随矩阵
    AD:=Array(1..n);
    ADA:=LinearAlgebra[IdentityMatrix](n);
    ADT:=Matrix(1..n,1..n);
    flogf[2]("Adjoint transformation matrixes :");
    for i from 1 to n do
        AD[i]:=LinearAlgebra[MatrixExponential](Matrix(convert(MK[i],list)),-epsilon[i]);
        ADA:=ADA.AD[i];
        ADT(i,1..n):=subs~(epsilon[i]=epsilon,(AD[i].sbs^%T)^%T);
        flogf[2]("Adjoint transformation matrix of %a",sbs[i]);
        flog[2](AD[i]);
    end do;
    flogf[2]("General adjoint transformation matrix");
    flog[2](ADA);

    flogf[2]("Adjoint representation table:");
    flog[2](ADT);

    # 计算不变量
    BA:=Matrix(1..n,1..n,(i,j)->b[i]*a[j]);
    pPhi:=add(BA*~M);
    eqs:=getPDE(pPhi,sbs);
    return AD,ADA,eqs;
end proc:

# 计算不变量
# 已化简并排序
getInvariants:=proc(eqs)
    local res;
    res:=pdsolve(eqs);
    res:=res[];
    res:=[op(op(2,res))];
    res:=sortByComplexity(res);# 按照复杂度升序输出
    flogf[0]("解得的不变量");
    map(x->flog[0]('Delta'[x]=res[x]),[seq(i,i=1..numelems(res))]);
    res:=simplifyInvariants(res);# 不变量化简
    res:=sortByComplexity(res);# 按照复杂度升序输出
    flogf[0]("化简后的不变量");
    map(x->flog[0]('Delta'[x]=res[x]),[seq(i,i=1..numelems(res))]);
    return res;
end proc;

(*
* 生成不变量的偏微分方程组
* 输入：
*     p        p=[w,v],w=sum(b[j]*v[j]),v=sum(a[i]*v[i])
*     sbs      生成元符号集 v[1],...,v[n]
* 输出：
*     偏微分方程组
*)
getPDE:=proc(p,sbs)
    local n:=numelems(sbs),i,eq,eqs;
    uses phi=phi(seq(a[i],i=1..n));
    eq:=add(coeff(p,sbs[i])*diff(phi,a[i]),i=1..n);
    eqs:={seq(coeff(eq,b[i]),i=1..n)} minus {0};
end proc:

$endif


text:
`is_element/nat_partitions` := (total::nonnegint,num_parts::nonnegint) -> proc(u)
 return type(u,list(nonnegint)) and nops(u) = num_parts and `+`(op(u)) = total;
end:

`is_leq/nat_partitions` := NULL;

`list_elements/nat_partitions` := proc(total,num_parts)
 local L,M,i,j,m,u;

 if num_parts = 0 then
  if total = 0 then
   return [[]];
  else
   return [];
  fi;
 fi;

 if num_parts = 1 then return [[total]]; fi;

 L := [seq([i],i=0..total)];
 for j from 2 to num_parts-1 do
  M := NULL;
  for u in L do
   m := total - `+`(op(u));
   M := M,seq([op(u),j],j=0..m);
  od:
  L := [M];
 od:
 L := map(u -> [op(u),total - `+`(op(u))],L);

 return L;
end:

`count_elements/nat_partitions` := (total::nonnegint,num_parts::nonnegint) ->
 binomial(total+num_parts-1,total);

`random_element/nat_partitions` := (total::nonnegint,num_parts::nonnegint) -> proc()
 local X,i;
 if total = 0 and num_parts = 0 then
  return [];
 fi;

 X := random_subset_of({seq(i,i=1..total+num_parts-1)},num_parts-1);
 X := [0,op(sort(X)),total+num_parts];
 return [seq(X[i+1]-X[i]-1,i=1..num_parts)];
end:

`is_element/pos_partitions` := (total::nonnegint,num_parts::nonnegint) -> proc(u)
 return type(u,list(posint)) and nops(u) = num_parts and `+`(op(u)) = total;
end:

`is_leq/pos_partitions` := NULL;

`random_element/pos_partitions` := (total::nonnegint,num_parts::nonnegint) -> proc()
 local p;
 if num_parts > total then return FAIL; fi;
 p := `random_element/nat_partitions`(total - num_parts,num_parts)();
 return p +~ 1;
end:

`list_elements/pos_partitions` := proc(total,num_parts)
 map(u -> u +~ 1,`list_elements/nat_partitions`(total - num_parts,num_parts));
end:

`count_elements/pos_partitions` := (total::nonnegint,num_parts::nonnegint) ->
 binomial(total-1,num_parts-1);


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

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

 global reason;

 procname_ := "is_element/ACP";

 if not (type(Q,list) and nops(Q) = N) then
  reason := [procname_,"Q is not a list of length N",Q,N];
  return false;
 fi;

 for i from 1 to N do
  if not `is_element/preord`(A)(Q[i]) then
   reason := [procname_,"Q[i] is not a preorder on A",Q[i],A,reason];
   return false;
  fi;
 od;

 if not `is_total/preord`(A)(Q[1]) then
  reason := [procname_,"Q[1] is not total",Q[1]];
  return false;
 fi;

 for i from 1 to N-1 do
  U := Q[i] intersect `op/preord`(A)(Q[i]);
  V := Q[i+1] union `op/preord`(A)(Q[i+1]);
  if U <> V then
   reason := [procname_,"Q[i] and Q[i+1] do not match",Q[i],Q[i+1]];
   return false;
  fi;
 od;

 return true;
end:

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

 procname_ := "is_equal/ACP";

 if Q1 <> Q2 then
  reason := [procname_,"Q1 <> Q2",Q1,Q2];
  return false;
 fi;

 return true;
end:

`is_leq/ACP` := (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/ACP` := (N::posint) -> proc(A::set)
 local W,X,Y,Z,Q,R,S,pi,B;

 if nops(A) = 0 then
  return FAIL;
 fi;

 if N = 1 then
  X := `list_elements/total_preord`(A);
  return map(Q -> [Q],X);
 fi;

 X := `list_elements/ACP`(N-1)(A);
 Y := NULL;

 for Q in X do
  pi := `block_partition/preord`(A)(Q[N-1]);
  Z := [{}];
  for B in pi do
   W := `list_elements/total_preord`(B);
   Z := [seq(seq(R union S,S in W),R in Z)];
  od:
  Y := Y,seq([op(Q),R],R in Z);
 od;
 return [Y];
end:

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

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

 if nops(A) = 0 then
  return FAIL;
 fi;

 R := `random_element/total_preord`(A)();
 Q := [R];
 for i from 2 to N do
  pi := `block_partition/preord`(A)(R);
  R := {};
  for B in pi do
   S := `random_element/total_preord`(B)();
   R := R union S;
  od;
  Q := [op(Q),R];
 od;

 return Q;
end:

`rank_vector/ACP` := (N::posint) -> (A::set) -> proc(Q)
 local i;
 [seq(`rank/preord`(N)(A)(Q[i]),i=1..N)];
end:

`rank/ACP` := (N::posint) -> (A::set) -> proc(Q)
 local i;
 add(`rank/preord`(N)(A)(Q[i]),i=1..N);
end:

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

`is_constant/ACP` := (N::posint) -> (A::set) -> proc(Q)
 return evalb(nops(Q[N]) = nops(A)^2);
end;

`is_separated/ACP` := (N::posint) -> (A::set) -> proc(Q)
 return evalb(Q[N] minus `id/autorel`(A) = {});
end;

`bot/ACP` := (N::posint) -> (A::set) -> [`top/autorel`(A) $ N];

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

`mu/W/ACP` := (N::posint) -> (A::set) -> proc(x)
 local QT,ET,AA,i,a,b;

 QT := table():
 ET := table():
 AA := {seq(seq([a,b],b in A),a in A)};
 ET[0] := AA;

 for i from 1 to N do
  QT[i] := select(ab -> is(x[ab[1]][i] <= x[ab[2]][i]),ET[i-1]);
  ET[i] := QT[i] intersect `op/autorel`(A)(QT[i]);
 od:

 return [seq(QT[i],i=1..N)];
end:

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

 u := [seq(`rank_table/preord`(A)(Q[i]),i=1..N)];
 x := table():
 for a in A do
  x[a] := [seq(u[i][a],i=1..N)];
 od:

 return eval(x);
end;

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

`is_element/realisation/ACP` := (N::posint) -> (A::set) ->
 `is_element/realisation/generic`("ACP",
  eval(`is_element/ACP`(N)(A)),
  eval(`is_leq/ACP`(N)(A))
 );

`is_equal/realisation/ACP` := (N::posint) -> (A::set) ->
 `is_equal/realisation/generic`("ACP",
  eval(`is_equal/ACP`(N)(A))
 );

`sigma/realisation/ACP/W` := (N::posint) -> (A::set) -> proc(u)
 local x,y,m,a,Q,S;

 x := table();
 for a in A do
  x[a] := [0$N];
 od:

 S := map(op,[indices(u)]);
 for Q in S do
  y := `sigma/ACP/W`(N)(A)(Q);
  m := u[Q];
  for a in A do
   x[a] := x[a] +~ (m *~ y[a]);
  od:
 od:

 return eval(x);
end;

`mu_aux/ACP` := (N::posint) -> (A::set) -> proc(x)
 local Q,t,i,ab,a,b,s,y,u;

 Q := `mu/W/ACP`(N)(A)(x);

 if nops(Q[N]) = nops(A)^2 then
  return(table());
 fi;

 t := 1;
 for i from 1 to N do
  for ab in Q[i] minus `op/autorel`(A)(Q[i]) do
   a,b := op(ab);
   t := min(t,x[b][i] - x[a][i]);
  od;
 od;

 s := `sigma/ACP/W`(N)(A)(Q);
 y := table();
 for a in A do
  y[a] := x[a] -~ (t *~ s[a]);
 od;

 u := `mu_aux/ACP`(N)(A)(y);
 if `is_element/RR`(u[Q]) then
  u[Q] := u[Q] + t;
 else
  u[Q] := t;
 fi;

 return eval(u);
end;

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

`rank_vector/ACP` := (N) -> (A) -> proc(Q)
 local i;
 return [seq(`rank/preord`(A)(Q[i])-1,i=1..N)];
end;

`rank/ACP` := (N) -> (A) -> proc(Q)
 return `+`(op(`rank_vector/ACP`(N)(A)(Q)));
end;

`cofunctor/ACP` := (N::posint) -> (A::set,B::set) -> (f) -> proc(Q)
 map(`cofunctor/autorel`(A,B)(f),Q);
end;

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

`is_element/split_spots/ACP` := (N::posint) -> (A::set) -> (Q) -> proc(iCD)
 global reason;
 local procname_,i,B,C,D,B1,E,C1;

 procname_ := "is_element/split_spots/ACP";

 if not(type(iCD,list) and nops(iCD) = 3) then
  reason := [procname_,"iCD is not a list of length 3",iCD];
  return false;
 fi;

 i,C,D := op(iCD);
 if not(type(i,posint) and i <= N) then
  reason := [procname_,"i is not in [1,N]",i,N];
  return false;
 fi;

 if not(type(C,set) and nops(C) > 0 and C minus A = {}) then
  reason := [procname_,"C is not a nonempty subset of A",C,A];
  return false;
 fi;

 if not(type(D,set) and nops(D) > 0 and D minus A = {}) then
  reason := [procname_,"D is not a nonempty subset of A",D,A];
  return false;
 fi;

 if C intersect D <> {} then
  reason := [procname_,"C and D are not disjoint",C,D];
  return false;
 fi;

 B := C union D;

 B1 := select(a -> member([a,B[1]],Q[i]) and member([B[1],a],Q[i]),A);
 if B1 <> B then
  reason := [procname_,"C u D is not an equivalence class for Q[i]",C,D,i,Q[i]];
  return false;
 fi;

 if i < N then
  E := select(ab -> member([ab[2],ab[1]],Q[i+1]),Q[i+1]);
  C1 := map(ab -> ab[1],select(ab -> member(ab[1],C),E));
  if C minus C1 <> {} then
   reason := [procname_,"C is not a union of equivalence classes for Q[i+1]",C,D,i,Q[i+1]];
   return false;
  fi;
 fi;

 return true;
end:

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

`list_elements/split_spots/ACP` := (N::posint) -> (A::set) -> proc(Q)
 local E,L,i,B,m,UU,U,CC,C,a;

 E := [seq(`block_partition/preord`(A)(Q[i]),i=1..N),{seq({a},a in A)}];
 L := NULL;
 for i from 1 to N do
  for B in E[i] do
   UU := select(U -> U minus B = {},E[i+1]);
   m := nops(UU);
   UU := select(U -> nops(U) < m,`list_elements/nonempty_subsets`(UU));
   CC := map(U -> map(op,U),UU);
   L := L,seq([i,C,B minus C],C in CC);
  od:
 od:
 return [L];
end:

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

`split/ACP` := (N::posint) -> (A::set) -> (Q) -> proc(iCD)
 local i,j,C,D,CD,DC,c,d;

 i,C,D := op(iCD);
 CD := {seq(seq([c,d],d in D),c in C)};
 DC := {seq(seq([d,c],d in D),c in C)};
 return([
  seq(Q[j],j=1..i-1),
  Q[i] minus DC,
  seq(Q[j] minus DC minus CD,j=i+1..N)
 ]);
end:

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

`is_element/glue_spots/ACP` := (N::posint) -> (A::set) -> (Q) -> proc(iCDS)
 global reason;
 local procname_,i,C,D,S,Qi,Qis,Ei,C1,D1,E,nC,nD;

 procname_ := "is_element/glue_spots/ACP";

 if not(type(iCDS,list) and nops(iCDS) = 4) then
  reason := [procname_,"iCDS is not a list of length 4",iCDS];
  return false;
 fi;

 i,C,D,S := op(iCDS);
 if not(type(i,posint) and i <= N) then
  reason := [procname_,"i is not in [1,N]",i,N];
  return false;
 fi;

 if not(type(C,set) and nops(C) > 0 and C minus A = {}) then
  reason := [procname_,"C is not a nonempty subset of A",C,A];
  return false;
 fi;

 if not(type(D,set) and nops(D) > 0 and D minus A = {}) then
  reason := [procname_,"D is not a nonempty subset of A",D,A];
  return false;
 fi;

 if C intersect D <> {} then
  reason := [procname_,"C and D are not disjoint",C,D];
  return false;
 fi;

 Qi := Q[i];
 Ei,Qis := selectremove(ab -> member([ab[2],ab[1]],Qi),Qi);

 if not(member([C[1],D[1]],Qis)) then
  reason := [procname_,"not(C < D)",C,D,Qi];
  return false;
 fi;

 C1 := select(a -> member([a,C[1]],Ei),A);
 if C1 <> C then
  reason := [procname_,"C is not an equivalence class",C,C1,Qi];
  return false;
 fi;

 D1 := select(a -> member([a,D[1]],Ei),A);
 if D1 <> D then
  reason := [procname_,"D is not an equivalence class",D,D1,Qi];
  return false;
 fi;

 E := select(a -> member([C[1],a],Qis) and member([a,D[1]],Qis),A);
 if E <> {} then
  reason := [procname_,"C and D are not adjacent",C,D,E,Qi];
  return false;
 fi;

 if i = N then
  if S = {} then
   return true;
  else
   reason := [procname_,"i=N but S <> {}",i,N,S];
   return false;
  fi;
 fi;

 nC := nops(`block_partition/preord`(C)(`res/preord`(A,C)(Q[i+1])));
 nD := nops(`block_partition/preord`(D)(`res/preord`(A,D)(Q[i+1])));

 if not (type(S,set(posint)) and
         nops(S) = nC and
         max(op(S)) <= nC+nD) then
  reason := [procname_,"S is not a subset of size nC in [1,nC+nD]",S,C,D,Q[i+1]];
  return false;
 fi;

 return true;
end:

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

`list_elements/glue_spots/ACP` := (N::posint) -> (A::set) -> proc(Q)
 local L,i,j,Qi,Qis,E,E0,E1,B,U,C,a,b,r0,nn,SS,S;

 L := NULL;
 for i from 1 to N do
  Qi := Q[i];
  Qis := remove(ab -> member([ab[2],ab[1]],Qi),Qi);
  E := `block_partition/preord`(A)(Qi);
  E0 := NULL;
  E1 := table():
  for B in E do
   E0 := E0,B[1];
   E1[B[1]] := B;
  od:
  E0 := {E0}:
  U := table():
  C := table():
  for a in E0 do
   U[a] := select(b -> member([a,b],Qis),E0);
  od:
  for a in E0 do
   C[a] := U[a] minus {seq(op(U[b]),b in U[a])};
  od:
  if i < N then
   r0 := `rank_table/preord`(A)(Q[i+1]);
   for a in E0 do
    nn[a] := 1 + max(seq(r0[b],b in E1[a]));
   od;
  fi;

  for a in E0 do
   for b in C[a] do
    if i < N then
     SS := combinat[choose]({seq(j,j=1..nn[a]+nn[b])},nn[a]);
    else
     SS := {{}};
    fi;
    L := L,seq([i,E1[a],E1[b],S],S in SS);
   od;
  od;
 od:
 return [L];
end:

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

`glue/ACP` := (N::posint) -> (A::set) -> (Q) -> proc(iCDS)
 local i,C,D,S,CD,DC,c,d,P,j,k,rc,rd,nc,nd,ic,id,n,Cn,Dn,BB,x,y;

 i,C,D,S := op(iCDS);
 CD := {seq(seq([c,d],d in D),c in C)};
 DC := {seq(seq([d,c],d in D),c in C)};
 P := table():
 for j from 1 to i-1 do
  P[j] := Q[j];
 od:
 P[i] := Q[i] union DC;
 if i < N then
  rc := `rank_table/preord`(C)(Q[i+1]);
  rd := `rank_table/preord`(D)(Q[i+1]);
  nc := max(seq(rc[c],c in C)) + 1;
  nd := max(seq(rd[d],d in D)) + 1;
  for n from 1 to nc do Cn[n] := select(c -> rc[c] = n-1,C); od;
  for n from 1 to nd do Dn[n] := select(d -> rd[d] = n-1,D); od;
  ic := 0;
  id := 0;
  BB := NULL;
  for j from 1 to nc + nd do
   if member(j,S) then
    ic := ic+1;
    BB := BB,Cn[ic];
   else
    id := id+1;
    BB := BB,Dn[id];
   fi;
  od;
  BB := [BB];
  P[i+1] := Q[i+1] union {seq(seq(seq(seq([x,y],y in BB[k]),x in BB[j]),k=j..nc+nd),j=1..nc+nd)};
 fi;
 for j from i+2 to N do
  P[j] := Q[j];
 od:
 return [seq(P[j],j=1..N)];
end:

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

`describe/ACP` := (N::posint) -> (A::set) -> proc(Q)
 local s,E,i,j,k,r,s0,s1,s2,B,C,m;

 s := "";
 E := {A};
 for i from 1 to N do
  s0 := "";
  for B in E do
   if nops(B) > 1 then
    r := `rank_table/preord`(B)(Q[i] intersect `top/autorel`(B));
    m := max(map(b -> r[b],B));
    s1 := "";
    for j from 0 to m do
     C := select(b -> r[b]=j,B);
     s2 := "";
     for k from 1 to nops(C) do
      s2 := cat(s2,`if`(k>1,"=",""),sprintf("%A",C[k]));
     od:
     s1 := cat(s1,`if`(j>0,"<",""),s2);
    od;
    s0 := cat(s0,`if`(s0="","",","),s1);
   fi;
  od:
  s := cat(s,`if`(s="","",";"),s0);
  E := `block_partition/preord`(A)(Q[i]);
 od;
 s := cat("(",s,")");
end:


text:
[plugin]
pluginname = RebootMadDevice
description = Reboot MAD Devices if they are not responding
author = GhostTalker
url = https://github.com/GhostTalker/RebootMadDevice/tree/V2
version = 2.4.0

max_stars_repo_name: Debaucus/RebootMadDevice

max_issues_repo_name: Debaucus/RebootMadDevice

max_forks_repo_name: Debaucus/RebootMadDevice


text:
`is_element/ICP_alt` := (N::posint) -> (A::set) -> proc(Ru)
 local n,R,u;
 global reason;

 n := nops(A);
 if n = 0 then
  reason := [convert(procname,string),"A is empty"];
  return false;
 fi;

 if not(type(Ru,list) and nops(Ru) = 2) then
  reason := [convert(procname,string),"Ru is not a list of length two"];
  return false;
 fi;

 R,u := op(Ru);

 if not `is_element/ord`(A)(R) then
  reason := [convert(procname,string),"R is not an order on A",reason];
  return false;
 fi;

 if not(type(u,list(posint)) and
        nops(u) = n-1 and
        max(op(u)) <= N) then
  reason := [convert(procname,string),"u is not in [1,N]^(n-1)"];
  return false;
 fi;

 return true;
end:

`is_equal/ICP_alt` := (N::posint) -> (A::set) -> proc(Ru1,Ru2)
 global reason;

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

 return true;
end:

`list_elements/ICP_alt` := (N::posint) -> proc(A::set)
 local RR,uu,n,R,u,i;

 RR := `list_elements/ord`(A);
 n := nops(A);
 uu := [[]];
 for i from 1 to n-1 do
  uu := [seq(seq([op(u),i],i=1..N),u in uu)];
 od:

 return [seq(seq([R,u],u in uu),R in RR)];
end:

`random_element/ICP_alt` := (N::posint) -> (A::set) -> proc()
 local R,u,i;

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

`count_elements/ICP_alt` := (N::posint) -> proc(A::set)
 local n;

 n := nops(A);
 return n! * N^(n-1);
end:


text:
type $RR <class {
  @e1 i32,
  @e2 f32,
  @classname "foo",
  @e3 f64}>
type $SS <class <$RR> { @accessflag 1234,
                        &method1(agg)agg,
                        &method2(void)void}>
func $foo (
	var %x <$SS>) i32 {
  dassign %x 2 ( constval i32 32 )
  return ( dread i32 %x 2 ) }
 # 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:
type $I1 <interface {@g0 i32}>
type $I2 <interfaceincomplete {@g1 i32}>
type $RR <interface <$I1> <$I2> {@e1 i32, &method1(agg)agg, &method2(void)void}>
javainterface $UU <$RR> public static
type $CC <interface {@h0 i32}>
type $TT <class <$CC> {@f1 i32, @f2 i8, @f3 i8, $I1, $I2, $RR}>
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# This function assumes that u is a symmetric polynomial in
# x[1],...,x[n] and returns a polynomial in c[1],...,c[n]
# that agrees with u after replacing the variables c[i] by
# the elementary symmetric functions in x[1],...,x[n].
# If the optional argument p_ is supplied, then the calculation
# is done mod p_.

newton_rewrite := proc(u,n,x,c,p_)
 local p,f,i,v,w,a,m,k,cx,cxp,vars;

 p := `if`(nargs > 4,p_,0);

 f := expand(mul(t+x[i],i=1..n));
 for i from 1 to n do
  cx[i] := coeff(f,t,n-i);
 od:

 cxp := proc(i,j)
  option remember;
  local u;

  if j = 0 then
   return 1;
  else
   u := expand(cx[i]*cxp(i,j-1));
   if p > 0 then u := mods(u,p); fi;
   return u;
  fi;
 end:

 v := u;
 w := 0;
 vars := plex(seq(x[i],i=1..n));
 while v <> 0 do
  a,m := LeadingTerm(v,vars);
  k := [seq(degree(m,x[i]),i=1..n),0];
  w := w + a * mul(c[i]^(k[i]-k[i+1]),i=1..n);
  v := expand(v - a * mul(cxp(i,k[i]-k[i+1]),i=1..n));
  if nargs > 4 then v := modp(v,p_); fi;
 od:

 return w;
end:

orbit_sum := proc(m,n,x)
 local S,R,v,C,c,r,i;
 S := combinat[permute](n);
 R := map(s -> [seq(x[i]=x[s[i]],i=1..n)],S);
 v := add(subs(r,m),r in R);
 C := {coeffs(v,{seq(x[i],i=1..n)})};
 c := igcd(op(C));
 v := expand(v/c);
 return v;
end:


text:
# # Draw weak and strong Bruhat diagram
# # by Shizuo Kaji

# it requires "coxeter and weyl" and "poset" packages by J.Stembridge

# http://www.math.lsa.umich.edu/~jrs/maple.html

with(posets);
with(coxeter); with(weyl);
with(WeylOps):

# Draw weak (left or right) bruhat diagram
#
# Example:
# H,L := wbruhat(F4, {}, 'L')
# plot_poset(H, labels = L, stretch = 1)

wbruhat:=proc(R,roots::set:={}, lr::symbol:='L', upto::integer :=0)
             local H,i,j,w,newweyl,oldweyl,p,pos,wlabel,udim;
  udim := `if`(upto>0,upto,nops(pos_roots(R)));
  wlabel:=[[]]; oldweyl:=[[]];
  H:=NULL;
  while nops(wlabel[-1])<udim and nops(oldweyl)>0 do
    newweyl:=[];
  	for i  from 1 to nops(oldweyl) do
	    for j in {$1..rank(R)} minus roots do
	    	w:=`if`(lr='L',reduce([j,op(oldweyl[i])],R),reduce([op(oldweyl[i]),j],R));
			if nops(w)=nops(oldweyl[i])+1 then
				if not member(w, newweyl, 'pos') then
					H:=H, [nops(wlabel)-nops(oldweyl)+i, nops(wlabel)+nops(newweyl)+1];
					newweyl:=[op(newweyl), w];
				else
					H:=H, [nops(wlabel)-nops(oldweyl)+i, nops(wlabel)+pos];
				end if;
			end if;
	    end do;
	end do;
    wlabel:=[op(wlabel), op(newweyl)]; oldweyl:=newweyl;
  end do;
  return {H},table([seq(i = convert(wlabel[i], 'string'), i = 1 .. nops(wlabel))]):
end proc:

# Draw strong bruhat diagram
#
# Example:
# H,L,E := wbruhat(F4, {}, 3)
# plot_poset(H, labels = L, stretch = 1)
# E contains edge labeling

sbruhat := proc (R, roots::set:={}, upto::integer:=0)
local H, w, r, v, i, j, beta, pos, newvects, newreps, vects,
    len, reps, oldvects, oldreps, edge, udim,B,P;

    B := base(R); r := nops(B); P:=pos_roots(R);
    udim := `if`(upto>0,upto,nops(P));
	newvects := [add(weights(R)[i], i={$1..r} minus roots)];
    newreps := []; reps := [[]]; vects := newvects;

    for len from 1 to udim do;
		oldvects := newvects; oldreps := [newreps];
		newvects := []; newreps := NULL;
		for j to nops(oldvects) do
			for i to r do
				if 0 < iprod(B[i], oldvects[j]) then
					v := reflect(B[i], oldvects[j]);
					if not member(v, newvects) then
						newvects := [op(newvects), v];
						newreps := newreps, reduce([i, op(oldreps[j])], R)
					end if
				end if
			end do
		end do;
		if newreps=NULL then break; end if;
        reps := [op(reps), newreps];
        vects := [op(vects), op(newvects)];
    end do;

    H:= NULL;
    for j from 1 to nops(vects) do
        for beta in P do
            if iprod(beta,vects[j])>0 then
                v:=reflect(beta,vects[j]);
                if member(v,vects,'pos') then
                    edge[j,pos]:=root_coords(beta,R);
                    H:= H, [j,pos];
                end if;
            end if;
        end do;
    end do;
    return {H},table([seq(i=[i,reps[i]],i=1..nops(reps))]),edge;
end proc;

# check if v is minimal w.r.t. a reflection subgroup;
has_descent:=proc(w::list,refsub_gen::set,R) local ref;
  for ref in refsub_gen do
    if nops(reduce([op(w),op(ref)],R))<nops(w) then return true; end if;
 end do;
 return false;
end proc:

# reduced word w.r.t. a reflection subgroup
ref_reduce:=proc(w::list,refsub_gen::set,R) local u,v,neww,is_min;
    neww:=reduce(w,R);
    is_min:=false;
    while is_min=false do;
        is_min:=true;
        for u in refsub_gen do
            v:=reduce([op(neww),op(u)],R);
            if nops(v)<nops(neww) then
                neww:=v;
                is_min:=false;
            end if;
        end do;
    end do;
    return neww;
end proc:

# strong bruhat graph w.r.t. a reflection subgroup geenrated by
# (non-simple) roots
refbruhat:=proc(R,refsub::set:={},upto::integer:=0) local H,i,j,r,u,v,w,len,
    reps,newreps,oldreps,vects,newvects,oldvects,refsub_gen,udim,pos,ref,P,B,edge;
    B:=base(R); r:=rank(R);
    P:=pos_roots(R);
    udim := `if`(upto>0,upto,nops(P));

    refsub_gen:={};
    for ref in refsub do;
        vec2fc(reflect(ref,interior_pt(R)),B,'w');
        refsub_gen:={op(refsub_gen),w};
    end do;

    # enumerate reduced words
	newvects := {interior_pt(R)};
    reps := [[]]; newreps := reps; vects := newvects;

    for len from 1 to udim do;
		oldvects := newvects; oldreps := newreps;
		newvects := []; newreps := [];
		for j to nops(oldvects) do
			for i to r do
				if 0 < iprod(B[i], oldvects[j]) then
                    v := reflect(B[i], oldvects[j]);
                    w := reduce([i, op(oldreps[j])], R);
					if not member(v, newvects) then
                        if not has_descent(w,refsub_gen,R) then
                            newreps := [op(newreps), w];
                            newvects := [op(newvects), v];
                        end if;
					end if
				end if
			end do
		end do;
		if newreps={} then break; end if;
        reps := [op(reps), op(newreps)];
        vects := [op(vects), op(newvects)];
    end do;

    # edges
    H:=NULL;
    for j from 1 to nops(vects) do;
        for ref in P do
            if iprod(ref,vects[j])<0 then next; end if;
            vec2fc(reflect(ref,vects[j]),R,'w');
            w:=ref_reduce(w,refsub_gen,R):
            member(w,reps,'pos');
            if j<pos then
                if assigned(edge[j,pos]) then
                    edge[j,pos]:=[op(edge[j,pos]),root_coords(ref,R)];
                else
                    edge[j,pos]:=[root_coords(ref,R)];
                    H:=H, [j,pos];
                end if
            end if;
        end do;
    end do;
    return {H},table([seq(i = [i,reps[i]], i = 1 .. nops(reps))]),edge:
end proc:

max_stars_repo_name: shizuo-kaji/Maple-flag-cohomology

max_issues_repo_name: shizuo-kaji/Maple-flag-cohomology

max_forks_repo_name: shizuo-kaji/Maple-flag-cohomology


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

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

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

(* These come from pbe correlation *)
params_a_gamma := (1 - log(2))/Pi^2:
mgamma := params_a_gamma:

(* r2scan values *)
params_a_dp2 := 0.361:

(* Equation (S6) *)
r2scan_alpha := (z, xt, ts0, ts1) -> (t_total(z, ts0, ts1) - xt^2/8) / (K_FACTOR_C*t_total(z, 1, 1) + params_a_eta*xt^2/8):

(* Equation (S26) *)
r2scan_f_alpha_neg := a -> exp(-params_a_c1*a/(1 - a)):
r2scan_f_alpha := (a, ff) -> my_piecewise5(
  a <= 0, r2scan_f_alpha_neg(m_min(a, 0)),
  a <= 2.5, rscan_f_alpha_small(m_min(a, 2.5), ff),
  rscan_f_alpha_large(m_max(a, 2.5))):

(* Equation (S28) *)
r2scan_d := z -> (opz_pow_n(z,5/3) + opz_pow_n(-z,5/3))/2:

(* Equation (S33): beta(rs), this is the same as in gga_c_regtpss *)
beta_a := 0.066724550603149220:
beta_b := 0.1:
beta_c := 0.1778:
mbeta := (rs) -> beta_a*(1 + beta_b*rs)/(1 + beta_c*rs):

(* Equation (S30) *)
w1 := (rs, z) -> exp(-f_pw(rs, z)/(mgamma*mphi(z)^3)) - 1:

(* Equation (S27); note that the paper indexes starting from zero *)
r2scan_dfc2 := ff -> add(i*ff[8-i], i=1..7):

(* According to James Furness, this is LSDA0 - see also Equation (S25) *)
r2scan_elsda0 := (rs, z) -> scan_eclda0(rs)*scan_Gc(z):
(* while LSDA1 is just Perdew-Wang *)
r2scan_elsda1 := (rs, z) -> f_pw(rs, z):

(* Derivatives wrt rs *)
r2scan_delsda0 := (rs, z) -> eval(diff(r2scan_elsda0(x1, x2), x1), [x1=rs, x2=z]):
r2scan_delsda1 := (rs, z) -> eval(diff(r2scan_elsda1(x1, x2), x1), [x1=rs, x2=z]):

(* Equation (S34) *)
r2scan_dy := (rs, z, s) -> r2scan_dfc2(rscan_fc)/(27 * mgamma * r2scan_d(z) * mphi(z)^3 * w1(rs, z)) * (
    + 20*rs*(r2scan_delsda0(rs, z) - r2scan_delsda1(rs, z))
    - 45*params_a_eta*(r2scan_elsda0(rs, z) - r2scan_elsda1(rs, z))
  ) * s^2*exp(-s^4/params_a_dp2^4):

(* Equation (S32) *)
r2scan_y := (rs, z, t) -> mbeta(rs)*t^2/(mgamma*w1(rs, z)):

(* Equation (S31) *)
r2scan_g := (rs, z, s, t) -> 1/(1 + 4*(r2scan_y(rs, z, t) - r2scan_dy(rs, z, s)))^(1/4):

(* Equation (S29) *)
fH := (rs, z, s, t) -> mgamma*mphi(z)^3*log(1 + w1(rs, z)*(1 - r2scan_g(rs, z, s, t))):

(* Now we can build ec1 from (S24) *)
r2scan_ec1 := (rs, z, s, t) -> f_pw(rs, z) + fH(rs, z, s, t):

(* Equation (S35)-(S41) are same as SCAN *)
r2scan_ec0 := (rs, z, s) -> scan_e0(rs, z, s):

(* and the functional itself *)
r2scan_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  r2scan_ec1(rs, z, X2S*2^(1/3)*xt, tt(rs, z, xt)) + r2scan_f_alpha(r2scan_alpha(z, xt, ts0, ts1), rscan_fc)*(
    + r2scan_ec0(rs, z, X2S*2^(1/3)*xt) - r2scan_ec1(rs, z, X2S*2^(1/3)*xt, tt(rs, z, xt))):

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


text:
type $ClassA <class {@e1 i32, @e2 f32, @e3 f64}>
type $ClassB <class {@e1 i32}>

func $foo() void {
  var %Reg1 <* <$ClassB>>
  dassign %Reg1 0 (intrinsicopwithtype ptr <* <$ClassB>> DEX_CONST_CLASS ())
  intrinsiccallwithtype <* <$ClassA>> DEX_CHECK_CAST (dread ptr %Reg1 0)
}

max_stars_repo_name: MapleSystem/OpenArkCompiler

max_issues_repo_name: MapleSystem/OpenArkCompiler

max_forks_repo_name: MapleSystem/OpenArkCompiler


text:
redundant_DIntos := module()
  export ModuleApply := proc(vs :: DomBound, sh :: DomShape, $)
    # This 'simplification' removes redundant information, but it is
    # entirely pointless as the result should be the same anyways. This
    # is mainly here as an assertion that Apply properly
    # re-applies integrals when the domain shape does not explicitly
    # state them.
    subsindets( sh, DomInto
              , proc (x, $)
                  local x_vn, x_t0, x_rest, x_t, x_mk;
                  x_vn, x_t0, x_rest := op(x);
                  x_t, x_mk := Domain:-Bound:-get(vs, x_vn);
                  if x_t = x_t0 then
                      x_rest
                  else
                      x
                  end if;
                end proc );
  end proc;

  export SimplName  := "Obviously redundant 'DInto's";
  export SimplOrder := 2;
end module;

empty_DIntos := module()
  uses Domain;

  export ModuleApply := proc(vs,sh,$)
    subsindets(sh, satisfies(is_empty), _->DSum());
  end proc;

  local is_empty := proc(e,$)
    evalb(e = DSum()) or (evalb(op(0,e) = DInto) and is_empty(op(3,e)));
  end proc;

  export SimplName  := "Empty DIntos";
  export SimplOrder := 15;
end module;

max_stars_repo_name: zaxtax/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 *)

thakkar_f0 := x -> 1 + 0.0055*x^2/(1 + 0.0253*x*arcsinh(x)):
thakkar_f1 := x -> -0.072*x/(1 + 2*4^(1/3)*x):

thakkar_f := x -> thakkar_f0(x) + thakkar_f1(x):

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


text:
with(LinearAlgebra):
with(Groebner):

# Exterior powers of matrices (implicitly using lexicographic ordering
# of the standard bases of exterior powers of R^n).

exterior_power := proc(k::nonnegint,A)
 local m,n,S,T,s,t;
 m,n := Dimension(A);
 S := combinat[choose](m,k);
 T := combinat[choose](n,k);
 Matrix([seq([seq(Determinant(SubMatrix(A,s,t)),t in T)],s in S)]);
end:

# If S is a subset of {1,..,n} then we have an associated shuffle
# permutation of {1,..,n}.  The next function computes the signature.

shuffle_sgn := proc(n,S)
 local k,T,ST;
 k := nops(S);
 T := sort([op({seq(i,i=1..n)} minus S)]);
 ST := [op(S),op(T)];
 return signum(mul(mul(ST[j]-ST[i],j=i+1..n),i=1..n));
end:

# Given a matrix A we get an induced map of Koszul algebras.  Each
# Koszul algebra has a natural Frobenius form, with using which we
# get a kind of dual map in the opposite direction, which is computed
# by the function below.  For a square matrix A, the composite of the
# k'th exterior power and the k'th koszul adjoint (in either order)
# is multiplication by det(A).

koszul_adjoint := proc(k::nonnegint,A)
 local m,n,M,N,S,T,Sc,Tc,Se,Te,s,t,i,j;
 m,n := Dimension(A);
 M := {seq(i,i=1..m)};
 N := {seq(j,j=1..n)};
 S := combinat[choose](m,k);
 T := combinat[choose](n,k);
 Sc := map(s -> sort([op(M minus {op(s)})]),S);
 Tc := map(t -> sort([op(N minus {op(t)})]),T);
 Se := [seq([op(S[i]),op(Sc[i])],i=1..nops(S))];
 Se := map(u -> signum(mul(mul(u[j]-u[i],j=i+1..m),i=1..m)),Se);
 Te := [seq([op(T[i]),op(Tc[i])],i=1..nops(T))];
 Te := map(u -> signum(mul(mul(u[j]-u[i],j=i+1..n),i=1..n)),Te);

 M := Matrix(nops(T),nops(S));
 for i from 1 to nops(S) do
  for j from 1 to nops(T) do
   M[j,i] := Determinant(SubMatrix(A,Sc[i],Tc[j])) * Se[i] * Te[j];
  od:
 od:

 return M;
end:

# A list u of length n in a ring R gives a differential on the
# Koszul algebra of R^n.  The function below gives the matrix of
# this differential with respect to the standard bases of the
# k'th and (k-1)'th exterior powers of R^n.

koszul_d := proc(k::posint,u)
 local n,S,T,s,t,i,j,MT,M;
 n := nops(convert(u,list));
 S := combinat[choose](n,k);
 T := combinat[choose](n,k-1);
 MT := table():
 for s in S do
  for t in T do
   MT[t,s] := 0;
  od:
  for i from 1 to k do
   t := [seq(s[j],j=1..i-1),seq(s[j],j=i+1..k)];
   MT[t,s] := (-1)^(i-1) * u[s[i]];
  od;
 od:
 M := Matrix([seq([seq(MT[T[i],S[j]],j=1..nops(S))],i=1..nops(T))]);
 return M;
end:

# koszul_p can be set equal to a prime number if we want to work in that
# characteristic.

koszul_p := infinity;

koszul_pmod := (x) -> `if`(koszul_p=infinity,x,mods(x,koszul_p)):

# The function below sets up a number of global variables.  The
# original v can be supplied as a list.  Then v_vec will be the corresponding
# row vector.  Also u_vec will be the vector of variables x[i], and w_vec
# will be the vector of powers x[i]^h, where h is minimal such that x[i]^h
# lies in I = (v[1],...,v[n])  for all i.  Next, A and B will be matrices with
# u_vec = v_vec.A and v_vec = w_vec.B, which witnesses the inclusions
# (x[1]^(h),...,x[n]^(h)) <= I <= (x[1],...,x[n]).  These matrices induce
# morphisms of Koszul complexes alpha : K(u) ->K(v) and beta: K(v)->K(w).  It
# is easy to understand K(u) and K(w) in particular,  they only have  homology
# in degree zero.  All three complexes can be regarded as commutative
# Frobenius algebras in the category of chain complexes, so there are adjoint
# morphisms alpha^! : K(v) -> K(u)  and beta^! : K(w) -> K(v).  The idea is
# to use these to show that K(v) also has homology concentrated in degree zero.

koszul_setup := proc(v)
 global n,u_vec,v_vec,w_vec,x_vec,x_vars,Iv,Mv,Nv,
        soc_u,soc_v,soc_w,A,B,p;
 local a,b,c,i,h,v0,v1,x0,x1,char;
 n := nops(v);
 v_vec := Transpose(Vector(v));
 x_vec := Transpose(<seq(x[i],i=1..n)>);
 w_vec := x_vec;
 soc_w := 1;
 x_vars := tdeg(seq(x[i],i=1..n));

 a := table();
 v0 := v_vec;
 v1 := v_vec;
 for i from n to 1 by -1 do
  v0 := v1;
  v1 := subs(x[i]=0,v1);
  a[i] := koszul_pmod((v0 - v1)/x[i]);
 od;

 B := Matrix([seq(convert(a[i],list),i=1..n)]);
 soc_v := koszul_pmod(Determinant(B));
 char := `if`(koszul_p=infinity,NULL,characteristic = koszul_p);

 Iv,Mv := Basis(convert(v,list),x_vars,output=extended,char);
 Mv := Transpose(Matrix(Mv));
 Nv := (z) -> koszul_pmod(NormalForm(z,Iv,x_vars,char));

 h := 1;
 x0 := map(Nv,convert(x_vec,list));
 x1 := x0;
 while h < 100 and x1 <> [0$n] do
  h := h+1;
  x1 := map(Nv,x0 *~ x1);
 od;

 if h >= 100 then
  error("Variables are not nilpotent mod given ideal");
 fi;

 u_vec := Transpose(Vector([seq(x[i]^h,i=1..n)]));
 soc_u := mul(x[i]^(h-1),i=1..n);

 a := NULL;
 b := NULL;
 for i from 1 to n do
  unassign('b');
  c := NormalForm(x[i]^h,Iv,x_vars,b,char);
  a := a,b;
 od;

 A := map(expand,Mv . Transpose(Matrix([a])));
 NULL;
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: work_gga_c *)

$define gga_c_pbe_params
$include "gga_c_pbe.mpl"

$include "lda_c_2d_amgb.mpl"

rs2D_factor := 1.704:
q2d_dd := 1e6:

q2d_rs2D := (rs, xt) -> rs2D_factor*rs*sqrt(X2S*xt)/RS_FACTOR:

fac := t -> t^4*(1 + t^2)/(q2d_dd + t^6):

f_q2d := (rs, z, xt, xs0, xs1) ->
  (1 - fac(tt(rs, z, xt)))*f_pbe(rs, z, xt, xs0, xs1) + fac(tt(rs, z, xt))*f_amgb(q2d_rs2D(rs, xt), z):

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

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
# Energy Calculation for Atlas Robot based on MDH frames
# Introduction
# Berechnung von potentieller 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)
# par2 -> Parametersatz 2 (erstes Moment anstelle von Schwerpunkt als Parameter)

# Anmerkung
# Nur potentielle Energie möglich
# Kinetische Energie stimmt nicht überein. Funktioniert nicht mit diesem Parametervektor und Geschwindigkeiten im Basis-KS
# Für die kinetische Energie sollte robot_tree_floatb_twist_rotmat_energy_linkframe_par2.mw genommen werden!
# 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
# 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(einzelne), 2=auch E_pot (komplett)
if not assigned(simplify_options) or simplify_options(6)=-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 := 2: # standardmäßig alle simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(6): # sechster Eintrag ist für Energie
end if:
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf_c := Trf_c:
read sprintf("../codeexport/%s/tmp/kinematics_com_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
r_W_i_Si := r_W_i_Si:
mr_W_i_Si := mr_W_i_Si:

# Potential Energy
U_b := Matrix(NL, 1):
U_grav := 0:
for i to NL do
  r_W_W_i := Matrix(3,1,Trf_c(1 .. 3, 4, i)):
  U_b[i] := -Multiply(Transpose(g_world), M[i, 1]*r_W_W_i+Matrix(3,1,mr_W_i_Si(1 .. 3, i))):
  printf("%s. Potentielle Energie aus Gravitation für Körper %d berechnet (im Welt-KS, mit Parametersatz 2).\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. 2; 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. 2) 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 nochmal in Gesamtsumme zusammenfassen (zusätzlich)
  if use_simplify>=2 then
    tmp_t1:=time(): tmp_l1 := length(U_grav):
    printf("%s. Beginne Vereinfachung für potentielle Energie (Param. 2; 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 (Param. 2) vereinfacht (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_par2_maple.m", robot_name, base_method_name):
# Kinetic Energy
# Using velocities in world frame and MX,MY,MZ as a parameter does not work.
# Only calculation of potential energy in this worksheet
# Matlab Export
printf("%s: Potentielle Energie (Param. 2) 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_par2_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:
for i from 1 to 6 do
  U_s_fixb := subs({V_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_par2_matlab.m", robot_name), codegen_opt):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

`is_element/partord` := (A::set) -> proc(R)
 global reason;

 if not `is_element/autorel`(A)(R) then
  reason := ["is_element/partord","R is not a relation on A",R,A,reason];
  return false;
 fi;

 if not(`is_reflexive/autorel`(A)(R)) then
  reason := ["is_element/partord","R is not reflexive",R];
  return false;
 fi;

 if not(`is_transitive/autorel`(A)(R)) then
  reason := ["is_element/partord","R is not transitive",R];
  return false;
 fi;

 if not(`is_separated/autorel`(A)(R)) then
  reason := ["is_element/partord","R is not separated",R];
  return false;
 fi;

 return true;
end;

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

`is_equal/partord`     := eval(`is_equal/autorel`);
`is_leq/partord`       := eval(`is_leq/autorel`);
`is_separated/partord` := eval(`is_separated/autorel`);

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

`op/partord` := eval(`op/autorel`);

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

`random_element/partord` := (A::set) -> proc()
 local M,B,E,RB,RA,DS,xy,x,y,d,m,a,e;
 if nops(A) <= 1 then return {seq([a,a],a in A)}; fi;

 M := `random_element/nonempty_subsets`(A)();
 B := A minus M;
 RB := `random_element/partord`(B)();
 DS := `list_elements/subsets`(B);
 for xy in RB do
  x,y := op(xy);
  if x <> y then
   DS := select(S -> (member(x,S) or not(member(y,S))),DS);
  fi;
 od;
 d := nops(DS);
 RA := RB union {seq([m,m],m in M)};
 for m in M do
  E := DS[rand(1..d)()];
  RA := RA union {seq([e,m],e in E)};
 od;

 return RA;
end:

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

`top_element/partord` := (A) -> proc(R)
 local n,a,U;
 n := nops(A);
 for a in A do
  U := select(e -> e[2] = a,R);
  if nops(U) = n then
   return a;
  fi;
 od:

 return FAIL;
end:

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

`bottom_element/partord` := (A) -> proc(R)
 local n,a,U;
 n := nops(A);
 for a in A do
  U := select(e -> e[1] = a,R);
  if nops(U) = n then
   return a;
  fi;
 od:

 return FAIL;
end:

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

`rank_table/partord` := (A::set) -> proc(R)
 local pi,pi0,a,b,r,k,B;

 r := table():
 for a in A do
  r[a] := nops(select(x -> member([x,a],R),A)) - 1;
 od;

 return eval(r);
end;

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

`hasse_diagram/partord` := (A::set) -> proc(R)
 local R1,R2;
 R1 := R minus `id/autorel`(A);
 R2 := `o/autorel`(A)(R1,R1);
 return R1 minus R2;
end:

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

`maximal_chains/partord` := (A::set) -> proc(R)
 local H,LC,UC,a,C,C0,extra,c,u;

 H := `hasse_diagram/partord`(A)(R);
 LC := table():
 UC := table():
 for a in A do
  LC[a] := map(e -> e[1],select(e -> e[2] = a,H));
  UC[a] := map(e -> e[2],select(e -> e[1] = a,H));
 od:

 C := map(a -> [a],select(a -> LC[a] = {},A));
 extra := true;
 while extra do
  extra := false;
  C0 := C;
  C := NULL;
  for c in C0 do
   u := UC[c[-1]];
   if u = {} then
    C := C,c;
   else
    C := C,seq([op(c),a],a in u);
    extra := true;
   fi;
  od:
  C := {C};
 od;

 return C;
end:

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

`partord_is_leq` := (A::set) -> (R) -> (a,b) -> member([a,b],R);

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

`partord_is_equiv` := (A::set) -> (R) -> (a,b) ->
  member([a,b],R) and member([b,a],R);

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

`partord_is_less` := (A::set) -> (R) -> (a,b) ->
  member([a,b],R) and not(member([b,a],R));

######################################################################
# Returns [A0,R0] where A0 is of the form {1,...,n} and
# (A0,R0) is isomorphic to (A,R).

`skeleton/partord` := (A::set) -> proc(R)
 local n,i,A0,Ai,R0;
 n := nops(A);
 A0 := {seq(i,i=1..n)};
 Ai := table():
 for i from 1 to n do Ai[A[i]] := i; od;
 R0 := map(u -> [Ai[u[1]],Ai[u[2]]],R);
 return [A0,R0];
end:

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

`partord/from_comparator` := (A::set) -> proc(C)
 local R,a,b;

 R := {seq(seq([a,b],b in A),a in A)};
 R := select(ab -> C(op(ab)),R);

 return [A,R];
end:

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

`skeleton/partord/from_comparator` := (A::{set,list}) -> proc(C)
 local n,i,j,A0,R0;

 n := nops(A);
 A0 := {seq(i,i=1..n)};
 R0 := {seq(seq([i,j],j = 1..n),i = 1..n)};
 R0 := select(ij -> C(A[ij[1]],A[ij[2]]),R0);

 return [A0,R0];
end:

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

`homology/partord` := (A::set) -> proc(R,p::nonnegint := 0)
 local U,C,Cn,Ci,dM,rdM,BN,H,c,c1,d,d_max,a,i,i1,j,k,t;

 U := table():
 for a in A do
  U[a] := select(b -> (b <> a and member([a,b],R)),A):
 od:
 C[0] := [seq(chain(a),a in A)]:
 d := 0;
 while C[d] <> [] do
  d := d+1;
  C[d] := [seq(seq(chain(op(c),a),a in U[op(d,c)]),c in C[d-1])]:
 od:
 d_max := d - 1;
 for d from 0 to d_max do
  Cn[d] := nops(C[d]):
  Ci[d] := table():
  for i from 1 to nops(C[d]) do
   Ci[d][C[d][i]] := i;
  od:
 od:

 dM := table():
 rdM := table():
 rdM[0] := 0;
 rdM[d_max+1] := 0;

 for d from 1 to d_max do
  dM[d] := Matrix(Cn[d-1],Cn[d]);
  for i from 1 to Cn[d] do
   c := [op(C[d][i])];
   for j from 0 to d do
    c1 := chain(seq(c[k+1],k=0..j-1),seq(c[k+1],k=j+1..d));
    i1 := Ci[d-1][c1];
    dM[d][i1,i] := (-1)^j;
   od:
  od:
  if p = 0 then
   rdM[d] := LinearAlgebra[Rank](dM[d]);
  else
   rdM[d] := LinearAlgebra[Modular][Rank](p,dM[d]);
  fi;
 od:

 BN := table():
 for d from 0 to d_max do
  BN[d] := Cn[d] - rdM[d] - rdM[d+1];
 od:

 H := table():
 H["characteristic"] := p;
 H["dim"] := d_max;
 H["chains"] := eval(C);
 H["chain_index"] := eval(Ci);
 H["differential_matrix"] := eval(dM);
 H["chain_number"] := eval(Cn);
 H["betti_number"] := eval(BN):
 H["poincare_series"] := unapply(add(BN[i]*t^i,i=0..d_max),t);
 H["euler_characteristic"] := add((-1)^d * BN[d],d = 0..d_max):

 return eval(H):
end:

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

`simplex/partord` := proc(n)
 local i,j;
 return [{seq(i,i=0..n)},{seq(seq([i,j],j=i..n),i=0..n)}];
end:

`crown/partord` := proc(n)
 local A,R,i;
 A := {seq(i,i=1..2*n)};
 R := {seq([i,i],i=1..2*n),
       seq([2*i-1,2*i],i=1..n),
       seq([2*i+1,2*i],i=1..n-1),
       [1,2*n]};
 return [A,R];
end:

`saw/partord` := proc(n)
 local A,R,i;
 A := {seq(i,i=1..n)};
 R := {seq([i,i],i=1..n),
       seq([2*i-1,2*i],i=1..floor(n/2)),
       seq([2*i+1,2*i],i=1..floor((n-1)/2))};
 return [A,R];
end:

`diamond/partord` := proc(n)
 local A,R,i,j;

 A := {seq(i,i=1..2*n)};
 R := {seq([i,i],i=1..2*n),
       seq(seq([2*i-1,j],j=2*i+1..2*n),i=1..n),
       seq(seq([2*i  ,j],j=2*i+1..2*n),i=1..n)};
 return [A,R];
end:


text:
type $RR <class {@e1 i32, @e2 f32, @e3 f64}>
type $SS <classincomplete <$RR> {@f1 i32, @f2 i8, @f3 i8}>
javaclass $TT <$SS> public static
func &foo (
	var %x <$SS>) i32 {
  dassign %x 2 ( constval i32 32 )
  return ( dread i32 %x 2 ) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
######################################################################
# CP(N) is the projective space CP^N

`is_element/CP` := (N::posint) -> proc(x,t := constant)
 if not (type(x,[t $ (N+1)])) then
  return false;
 fi;

 if simplify(x) = [0$(N+1)] then
  return false;
 fi;

 return true;
end:

`is_equal/CP` := (N::posint) -> proc(x,y)
 local i,j;
 for i from 1 to N do
  for j from i+1 to N do
   if simplify(x[i] * y[j] - x[j] * y[i]) <> 0 then
    return false;
   fi;
  od;
 od;

 return true;
end:

`is_leq/CP` := NULL;

`random_element/CP` := (N::posint) -> proc() local i; [seq(`random_element/CC`(),i=1..N+1)] end;

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


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_x_br89_params *params;

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

(* replace: "br89_x\(" -> "xc_mgga_x_br89_get_x(" *)

(* This is the derivative of f = x*exp(-2.0/3.0*x)/(x - 2) = y = 2*Pi^(2/3)/(3*Q) *)
br89_aux_dfdx := x -> -2/3 * exp(-2*x/3) * (x^2 - 2*x + 3) / (x - 2)^2:

`diff/br89_x` := proc(Q, g)
  -2/3 * Pi^(2/3) * diff(Q, g)/(Q^2 * br89_aux_dfdx(br89_x(Q)))
end proc:

br89_Q := (x, u, t) ->
  (u - 4*params_a_gamma*t + params_a_gamma*x^2/2)/6:

br89_min_Q := 5.0e-13:
br89_cQ := Q -> my_piecewise3(abs(Q) < br89_min_Q,
  my_piecewise3(Q > 0, br89_min_Q, -br89_min_Q), Q):

br89_v := x ->
  -2*Pi^(1/3)/X_FACTOR_C * exp(x/3)*(1 - exp(-x)*(1 + x/2))/x:

br89_mx := Q -> br89_x(Q):

br89_f := (x, u, t) ->
  - br89_v(br89_mx(br89_cQ(br89_Q(x, u, t))))/2 *
  (1 + params_a_at*mgga_series_w([0, 1, 0, -2, 0, 1], 6, t)):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  mgga_exchange(br89_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: work_mgga_c *)

par_n := 6:

gamma_x := 0.004:
par_x := [
    [  0.85,  0, 0],
    [  1.007, 0, 1],
    [  0.259, 1, 0],
    [  0,   0, 0],
    [  0,   0, 0],
    [  0,   0, 0]
]:

gamma_ss := 0.2:
par_ss := [
    [  0.443,  0, 0],
    [ -1.437,  0, 4],
    [ -4.535,  1, 0],
    [ -3.39,   2, 0],
    [  4.278,  4, 3],
    [  0,    0, 0]
]:

gamma_os := 0.006:
par_os := [
    [  1.000,  0, 0],
    [  1.358,  1, 0],
    [  2.924,  2, 0],
    [ -8.812,  2, 1],
    [ -1.39,   6, 0],
    [  9.142,  6, 1]
]:

$include "lda_x_erf.mpl"
$include "b97mv.mpl"

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


text:
var $g_i8_u <i8>
var $g_i8 <i8> = 8
var $g_u8_u <u8>
var $g_u8 <u8> = 8

var $g_i16_u <i16>
var $g_i16 <i16> = 16
var $g_u16_u <u16>
var $g_u16 <u16> = 16

var $g_i32_u <i32>
var $g_i32 <i32> = 32
var $g_u32_u <u32>
var $g_u32 <u32> = 32

var $g_i64_u <i64>
var $g_i64 <i64> = 0x123456789a
var $g_u64_u <u64>
var $g_u64 <u64> = 0x123456789a

var $g_f32_u <f32>
var $g_f32 <f32> = 1.2f

var $g_f64_u <f64>
var $g_f64 <f64> = 2.4

# i8/u8
func $add_i8_u () i8 {
   return (
     add i32(dread i32 $g_i8_u, constval i32 0))}

func $add_i8 () i8 {
   return (
     add i32(dread i32 $g_i8, constval i32 0))}

func $add_u8_u () u8 {
   return (
     add u32(dread u32 $g_u8_u, constval u32 0))}

func $add_u8 () u8 {
   return (
     add u32(dread u32 $g_u8, constval u32 0))}

# i16/u16
func $add_i16_u () i16 {
   return (
     add i32(dread i32 $g_i16_u, constval i32 0))}

func $add_i16 () i16 {
   return (
     add i32(dread i32 $g_i16, constval i32 0))}

func $add_u16_u () u16 {
   return (
     add u32(dread u32 $g_u16_u, constval u32 0))}

func $add_u16 () u16 {
   return (
     add u32(dread u32 $g_u16, constval u32 0))}

# i32/u32
func $add_i32_u () i32 {
   return (
     add i32(dread i32 $g_i32_u, constval i32 0))}

func $add_i32 () i32 {
   return (
     add i32(dread i32 $g_i32, constval i32 0))}

func $add_u32_u () u32 {
   return (
     add u32(dread u32 $g_u32_u, constval u32 0))}

func $add_u32 () u32 {
   return (
     add u32(dread u32 $g_u32, constval u32 0))}

# i64/u64
func $add_i64_u () i64 {
   return (
     add i64(dread i64 $g_i64_u, constval i64 0))}

func $add_i64 () i64 {
   return (
     add i64(dread i64 $g_i64, constval i64 0))}

func $add_u64_u () u64 {
   return (
     add u64(dread u64 $g_u64_u, constval u64 0))}

func $add_u64 () u64 {
   return (
     add u64(dread u64 $g_u64, constval u64 0))}

# # f32
# func $add_f32_u () f32 {
   # return (
     # add f32(dread f32 $g_f32_u, constval f32 0.0f))}

# func $add_f32 () f32 {
   # return (
     # add f32(dread f32 $g_f32, constval f32 0.0f))}

# # f64
# func $add_f64_u () f64 {
   # return (
     # add f64(dread f64 $g_f64_u, constval f64 0))}

# func $add_f64 () f64 {
   # return (
     # add f64(dread f64 $g_f64, constval f64 0))}
 # 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: gga_exc *)
(* prefix:
  gga_x_ak13_params *params;

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

ak13_f0 := s -> 1 + params_a_B1*s*log(1 + s) + params_a_B2*s*log(1 + log(1 + s)):
ak13_f  := x -> ak13_f0(X2S*x):

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


text:
`is_element/young_partitions` := (k,n) -> proc(lambda)
 local i;

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

 if nops(lambda) <> k then
  return false;
 fi;

 if k > 0 and lambda[1] > n-k then
  return false;
 fi;

 for i from 1 to k-1 do
  if lambda[i] < lambda[i+1] then
   return false;
  fi;
 od;

 return true;
end:

`young_partition/from_set` := (n,k) -> proc(P)
 local i;
 return [seq(n-k+i-P[i],i=1..k)];
end:

`list_elements/young_partitions` := proc(n,k)
 map(`young_partition/from_set`(n,k),combinat[choose](n,k));
end:

`random_element/young_partitions` := (n,k) -> proc()
 if k > n then return FAIL; fi;
 `young_partition/from_set`(n,k)(combinat[randcomb](n,k));
end:

`count_elements/young_partitions` := (n,k) -> binomial(n,k);

`weight/young_partitions` := (lambda) -> `+`(op(lambda)):

`length/young_partitions` := (lambda) -> nops(lambda):

`cell_dim/young_partitions` := (n,k) ->
 k * (n - k) - `weight/young_partitions`(lambda):

`matrix/young_partitions` := (n,k) -> (lambda) -> proc(a)
 local A,P,S,u,i,j;
 A := Matrix(k,n);
 P := [seq(n-k+i-lambda[i],i=1..k)];
 S := {op(P)};
 u := 1;
 for i from 1 to k do
  for j from 1 to P[i] - 1 do
   if not(member(j,S)) then
    A[i,j] := a[u];
    u := u+1;
   fi;
  od:
  A[i,P[i]] := 1;
 od:
 return A;
end:

`giambelli_matrix/young_partitions` := proc(lambda)
 local i,j,k,f;
 k := nops(lambda);
 f := (i) -> `if`(i < 0,0,`if`(i = 0,1,c[i]));
 Matrix([seq([seq(f(lambda[i]-i+j),j=1..k)],i=1..k)]);
end:

`schur_function/young_partitions` := (lambda) ->
 Determinant(`giambelli_matrix/young_partitions`(lambda)):