text:
# Apply a domain to an expression.
# Apply will perform certain 'simplifications' to make sure the
# domain application is well-formed. e.g.,
#   DOMAIN( DBound( [x,y], [xt,yt] ), DConstrain(..) )
#      is the same as
#   DOMAIN(        ..               , DInto(x,xt,DInto(y,yt,DConstrain(..))) )
# basically, when we see an 'unbound' variable in the 'RHS' , we should bind
# it with the default 'DInto'.
export Apply := module ()
  uses Domain_Type, Hakaru, Utilities;
  # Make sure these refer to the global names, or Records don't work
  global context, weights, f_into, f_body, f_sum, f_constrain;

  # Validates input and calls `do_apply' which implements the main logic.  This
  # function returns a function which takes as input the program corresponding
  # to the 'body' (i.e. what is left after extracting the domain bounds and
  # shape) and produces a program equivalent to the original (i.e. what was
  # present before extracting domain bounds and shape), but simplified. In
  # equational language:
  #   let (b,e1) = Extract:-Bound(e0)
  #       (s,e2) = Extract:-Shape(e1)
  #       e3     = Apply(DOMAIN(s,b))(e2)
  #   in  e0 == e3
  # is an invariant of Domain.
  #
  # `Apply' also takes an additional argument, which is a record which should
  # contain handlers which determine how Domain constructors are applied to the
  # resulting expression, and a table of weights which appear outermost to each
  # domain bound variable. Some of these may have defaults (see *** below).
  #
  # The weights table is a mapping from variables appearing in the domain bounds
  # to weights which appear immediately outside the application of that domain
  # bound. It is passed to certain callbacks which can chose where to place said
  # weights.
  #
  # Each handler loosely corresponds to a Domain constructor, except for
  # `f_body'.
  # f_into : DomBoundKind -> LO x -> DomBoundVar -> DomBoundRange -> KB
  #        -> WeightsTable -> LO x
  #   The 1st, 3rd, and 4th arguments come directly from a DomInto;
  #   the 2nd argument is the recursive application of Apply;
  #   the 5th and 6th arguments come from the context
  # f_sum  : Sequence (LO x) -> LO x
  #   Sequence being a Maple sequence (i.e. `f_sum' is an n-ary function)
  # f_constrain : (forall x . LO x -> LO x) -> LO x -> LO x
  #   The 1st argument places a constraint on a LO, the 2nd arguement is the
  #   expression to be constrained.
  # f_body : LO x -> KB -> LO x
  #   This handler is called on the 'body', along with the context in which that
  #   'body' expression exists. This may be called multiple times, if for
  #   example the domain contains `DSum' or `DSplit'.
  export ModuleApply := proc(dom :: HDomain_mb, ctx0::record, $)
    local vs, sh, dbnd, ctx, defctx;
    if dom :: DomNoSol then
      error "cannot apply %1", dom;
    end if;
    dbnd, sh := op(dom);
    # *** Default handlers
    defctx :=
    Record('context'=Domain:-Bound:-contextOf(dbnd)
           ,'f_into'=do_mk
           ,'f_body'=do_body
           ,'f_sum'=`+`
           ,'f_constrain'=((f,x)->f(x)));
    ctx := merge_record_default(defctx,ctx0);
    if not(ctx::'apply_context') or
    {'context'} subset {exports(ctx0)} then
      error "invalid input: %1", ctx0;
    end if;
    vs := Domain:-Bound:-withVarsIxs(dbnd);
    (e->do_apply({}, e, vs, sh, ctx));
  end proc;

  local ModuleLoad := proc($)
    BindingTools:-load_types_from_table(
      # context := the context as a KB
      # weights := the weights outer to each integral
      # f_*     := handlers for subexpressions
      table([apply_context=
             '`record`(`context`,
                       `weights`,
                       `f_into`, `f_body`, `f_sum`, `f_constrain`)']));
  end proc;

  # Helpers for applying only a shape or a bounds.  Note that in the shape case,
  # it may not contain `DInto's as these are only valid in the presence of a
  # non-empty bounds. This is not checked.
  export Shape := proc(dsh :: DomShape, e, $) ModuleApply( DOMAIN( [], dsh ), e ) end proc;
  export Bound := proc(dbn :: DomBound, e, $) ModuleApply( DOMAIN( dbn, DConstrain() ), e ) end proc;

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

      shv, cond_outer :=
        select_cond_outer(shv, done_ union (indets(shv,DomBoundVar) intersect vars));

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

  # A default implementation for `f_into'
  export do_mk := proc(mk, e, vn, ty_, _kb)
    Domain:-ExtBound[mk]:-DoMk(e, vn, ty_);
  end proc;

  # A default implementation for `f_body'
  export do_body := proc(e, _kb) e end proc;

  # Like an `Indicator' but performs some early simplification which covers the
  # common case of `DConstrain()' (and is expressed in terms of Partition).
  local do_constrain := proc(cond0::{list,set,DomConstrain},$)
    local mk_cond, cond := cond0;
    if nops(cond)=0 then
      mk_cond := x->x;
    else
      cond := bool_And(op(cond));
      mk_cond := x->PARTITION([Piece(cond,x), Piece(Not(cond),0)]);
    end if;
    mk_cond;
  end proc;

  # A very cute way of pulling constraints out over other constructors
  # which relies on the fact that constraints are required to be innermost;
  # i.e. DomShape is a sum of products; i.e. this could also be called
  # 'factoring'. e.g.:
  # DSum( DConstrain(A, C), DConstrain(A, B) )
  # -> DSum( DConstrain(C), DConstrain(B) ), {A}
  # DSum( DConstrain(A), DConstrain(B) )
  # -> DSum( DConstrain(A), DConstrain(B) ), { }
  # DSum( DInto(x, .., DConstrain(A, B)), DInto(x, .., DConstrain(A, C))),
  # and `x' is not in `vars'
  # -> DSum( DInto(x, .., DConstrain(B)), DInto(x, .., DConstrain(C) ) ), {A}
  local select_cond_outer := proc(sh::DomShape, vars0::set({name,list(name)}), $)
    local csd, cs0, cs, ots, sh1, ins, vars := vars0;
    if sh=DConstrain() then return sh, {}; end if;
    vars := map(v->`if`(v::list,v,[v])[],vars); # handles(?) array variables

    # the constraints which appear in the shape and don't mention dependant variables
    csd := [op(indets(sh, And(DomConstrain,Not(satisfies(x->has(x,vars))))))];
    cs0 := map(x->{op(x)},csd);
    if cs0=[] then return sh, {}; end if; # if there are no such constraints

    # Split each constraint into a pair whose first component are
    # constraints common to all constraints in the shape
    cs := map(x->[selectremove(z->andmap(c->z in c,cs0),x)], cs0);

    # outers (common constraints) and inners (rest)
    ots := op([1,1],cs); ins := map(curry(op,2), cs);

    # substitute back into the shape the reduced constraints
    sh1 := subs(zip(`=`, csd, map(DConstrain@op,ins)), sh);
    userinfo(3, Domain, printf("select_cond_outer(%a, %a) = %a, %a\n", sh, vars, sh1, ots));
    sh1, ots;
  end proc;

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

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

    rn := `if`(evalb(vn=vn_), [], [vn=vn_]);
    upd_field(ctx, 'context', _->kb1), rn;
  end proc;

  # Performs multiple integrals, corresponding to those variables given in
  # `todo'. If `todo' is a set, then the order of integrals is one which matches
  # the order in the original expression. If `todo' is a list, this function
  # assumes that integrating in that order is valid (i.e. doesn't cause scope
  # extrusion). `kont' takes an expression and a KB corresponding to the
  # innermost context in which the input `e' will be used.
  #   do_mks(e, kont, [v0,v1..vn], dbnd, ctx{context=outerKB}) =
  #     Int(v0, RangeOf(v0,dbnd),
  #     Int(v1, RangeOf(v1,dnbd),
  #     ...
  #     Int(vn, RangeOf(vn,dnbd),
  #       kont(e, innerKB)) ... ))
  #   where innerKB is outerKB additionally augmented with
  #          the knowledge that
  #            `{v0 in RangeOf(v0,dbnd), v1 in RangeOf(v1,dbnd), .., vn in ..}'
  local do_mks := proc(e, kont, todo::({list,set})(DomBoundVar), dbnd :: DomBound, ctx, $)
    local vs, v_td, i, vt, kb0, v_mk, _, ctx1, rn, r := e;

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

    if nops(vs)=0 then
      return kont(r,ctx:-context);
    end if;

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

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

    r := do_mks(r, kont, vs, dbnd, ctx1);
    ctx:-f_into(v_mk, r, v_td, vt, ctx:-context, ctx:-weights[v_td]);
  end proc;
end module;


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

tm_lambda := 0.6866:
tm_beta   := 79.873:

(* below Equation (6) *)
tm_p  := x -> (X2S*x)^2:
tm_y  := x -> (2*tm_lambda - 1)^2 * tm_p(x):

(* Equation (7) *)
tm_f0 := x -> (1 + 10*(70*tm_y(x)/27) + tm_beta*tm_y(x)^2)^(1/10):

(* after Equation (9) *)
tm_R  := (x, t) -> 1 + 595*(2*tm_lambda - 1)^2 * tm_p(x)/54 \
   - (t - 3*(tm_lambda^2 - tm_lambda + 1/2)*(t - K_FACTOR_C - x^2/72))/K_FACTOR_C:

tm_fx_DME := (x, t) -> 1/tm_f0(x)^2 + 7*tm_R(x, t)/(9*tm_f0(x)^4):

tm_alpha := (x, t) -> (t - x^2/8)/K_FACTOR_C:

(* after Equation (11) *)
tm_qtilde := (x, t) -> 9/20*(tm_alpha(x, t) - 1) + 2*tm_p(x)/3:

(* Ratio tW/t; we have to make sure it's 1 at maximum *)
tm_tratio := (x, t) -> m_min(1.0, x^2/(8*t)):

tm_fx_SC := (x, t) -> (1 + 10*( \
       + (MU_GE + 50*tm_p(x)/729)*tm_p(x) + 146*tm_qtilde(x, t)^2/2025 \
       - 73*tm_qtilde(x,t)/405*(3/5*tm_tratio(x,t))*(1 - tm_tratio(x,t)))
       )^(1/10):

(* Equation 10 and below *)
tm_w := (x,t)-> (tm_tratio(x,t)^2 + 3*tm_tratio(x,t)^3)/(1 + tm_tratio(x,t)^3)^2:

tm_f := (x, u, t) -> tm_w(x,t)*tm_fx_DME(x, t) + (1 - tm_w(x,t))*tm_fx_SC(x, t):

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

ext: mpl


text:
input := FileTools:-Text:-ReadFile("AoC-2021-18-input.txt" ):

find5nested := proc(input::string,$)
local c, i, count;
    count := 0;
    for i from 1 to length(input) do
        if input[i] = "[" then count += 1; if count=5 then return i; end if;
        elif input[i] = "]" then count -= 1;
        end if;
    end do;
    return 0;
end proc:

PRINT := s->printf(cat(s,"\n"),_rest):

sfreduce := proc(inp::string,$)
local n, pin, regnums, input:=inp, oldinput;

    while oldinput <> input do
        oldinput := input;
    #check
        pin := parse(input);
        if not nops(pin) = 2 then
        error "bad input", input;
        end if;

        PRINT("reducing %s", input);

        n := find5nested(input);
        if n > 0 then
            input := ( sfexplode(input, n) );
            next;
        end if;

        regnums := indets(pin, integer);
        if max(regnums) >= 10 then
            input := ( sfsplit(input) );
            next;
        end if;
    end do;

    return input;

end proc:

sfexplode := proc(input::string, pos::posint:=find5nested(input), $)
local mid, toreplace, toreplacestr;
local lh := input[1..pos-1];
local rh := input[pos..-1]; # starting at the fifth [
local m,n := StringTools:-FirstFromLeft("]", rh);
local toexplode := rh[1..n];
    rh := rh[n+1..-1];
    toexplode := parse(toexplode); # get a list with ints
# righthand element
    n := StringTools:-FirstFromRight(",", lh); # any number will be followed by a ,

    PRINT("exploding %a at position %d", toexplode, pos);

    if n <> 0 then
        if lh[n-1] = "]" then
            while lh[n-1] = "]" do n -= 1 end do;
            mid := lh[n..-1];
            lh := lh[1..n-1]; # number still in lh
            n := max(StringTools:-FirstFromRight(",", lh),StringTools:-FirstFromRight("[", lh));
        else
            mid := lh[n..-1]; # number still in lh
            lh := lh[1..n-1];
            n := StringTools:-FirstFromRight("[", lh);
        end if;
        toreplace := lh[n+1..-1];
        lh := lh[1..n];
        toreplacestr := parse(toreplace) + toexplode[1];
        lh := cat(lh, toreplacestr, mid);
        end if;

    if toreplace = "" then
      error "bug"
    end if;

# lefthand element
    n := StringTools:-FirstFromLeft(",", rh); # any number will be proceded by a ,
    if n <> 0 then
        if rh[n+1] = "[" then
            while rh[n+1] = "[" do n += 1 end do;
            mid := rh[1..n];
            rh := rh[n+1..-1]; # number still in rh
            n := min(subs(0=NULL,[StringTools:-FirstFromLeft(",", rh), StringTools:-FirstFromLeft("]", rh)]));
        else
        mid := rh[1..n];
        rh := rh[n+1..-1]; # number still in rh
        n := StringTools:-FirstFromLeft("]", rh);
    end if;
        toreplace := rh[1..n-1];
        rh := rh[n..-1];
        toreplacestr := parse(toreplace) + toexplode[2];
        rh := cat(mid, toreplacestr, rh);
    end if;

    if toreplace = "" then
      error "bug"
    end if;

    return cat(lh,0,rh);

end proc:

sfsplit := proc(input::string, $)
local i,n,regnum,dig,diglocs;
    # find leftmost - this is slow
    # a better would be to just sweep to find a two digit number
    dig := map(d->cat("",d), [select(d->d>9,indets(parse(input), 'posint'))[]]);
    diglocs := map(d->StringTools:-Search(d, input), dig);
    i := min[index](diglocs);
    regnum := parse(dig[i]);
    n := diglocs[i];

    PRINT("splitting %d at position %d", regnum, n);
    if n > 0 then
        return cat(input[1..n-1], "[", floor(regnum/2), ",",
                   ceil(regnum/2) ,"]", input[n+length(rs)..-1]);
    else
        error "%1 is not in %2", regnum, input;
    end if;
end proc:

sfnorm := proc(input::{list,integer,string},$)
    if type(input, string) then return thisproc(parse(input)); end if;
    if type(input, integer) then return input; end if;
    if nops(input) <> 2 then error "bad sf number", input; end if;
    return thisproc(input[1])*3 + thisproc(input[2])*2;
end proc:

sfadd := proc(s1, s2)
   return cat("[",s1,",",s2,"]");
end proc:

# unit tests
sfexplode( "[[[[[9,8],1],2],3],4]")="[[[[0,9],2],3],4]";
sfexplode( "[7,[6,[5,[4,[3,2]]]]]") = "[7,[6,[5,[7,0]]]]";
sfexplode("[[6,[5,[4,[3,2]]]],1]") = "[[6,[5,[7,0]]],3]";
sfexplode("[[3,[2,[1,[7,3]]]],[6,[5,[4,[3,2]]]]]") = "[[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]";
sfexplode("[[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]")="[[3,[2,[8,0]]],[9,[5,[7,0]]]]";
sfadd("[[[[4,3],4],4],[7,[[8,4],9]]]", "[1,1]");
sfexplode(%);
sfexplode(%);
sfsplit(%);
sfsplit(%);
sfexplode(%)= "[[[[0,7],4],[[7,8],[6,0]]],[8,1]]";
sfreduce("[[[[[4,3],4],4],[7,[[8,4],9]]],[1,1]]") = "[[[[0,7],4],[[7,8],[6,0]]],[8,1]]";
sfsplit( "[[[[5,11],[13,0]],[[15,14],[14,0]]],[[2,[11,10]],[[0,8],[8,0]]]]") = "[[[[5,[5,6]],[13,0]],[[15,14],[14,0]]],[[2,[11,10]],[[0,8],[8,0]]]]";

PRINT := 'PRINT':
homework := StringTools:-Split(input,"\n"): nops(homework);

sofar := homework[1];
for i from 2 to nops(homework) do
    sofar := sfreduce(sfadd(sofar, homework[i]));
end do:
answer1 := sfnorm(sofar);
answer2 := max([seq(seq(ifelse(i<>j,sfnorm( sfreduce(sfadd(homework[i],homework[j])) ),NULL), j=1..nops(homework)),i=1..nops(homework))]);


text:
program numberlist;  { sample15 }
var e, n: integer;
begin
	writeln('Number list');
	writeln('    n  ', '   2n  ', ' n**2  ', ' n**3  ', ' 2**n  ');
	n := 0;
	e := 1;
	while n < 15 do begin
		writeln(n:6, ' ', 2*n:6, ' ', n*n:6, ' ', n*n*n:6, ' ', e:6);
		e := e * 2;
		n := n + 1;
	end
end.


text:
#Symetry operators#

#I need Kw, the superN2symetry operator

Kzij:= proc(poly, i, j, forceVars:=1)
    local p_poly;
    p_poly:= sort(subs({z||i = z||j, z||j = z||i}, poly));
    return p_poly;
end proc:

Kscrypt_ij:= proc(poly, ijvec)
    local KsijPoly, i,j;
    i:= ijvec[1];
    j:= ijvec[2];
	KsijPoly:= sort(subs({z||i = z||j, z||j = z||i, theta||i = theta||j, theta||j= theta||i, phi||i=phi||j, phi||j = phi||i}, poly));
    return KsijPoly;
end proc:

Ksigma_j:= proc(poly, someSet)
    local sigma_j, sigma_1, non_fixed_points, permVect, KsigmaPoly, RcycleIndex, theSubSeq, Ksigma_jPoly, table_index;
    sigma_j:= someSet;
    sigma_1:= [seq(i, i=1..nops(someSet))];
    non_fixed_points:= map(x -> if sigma_j[x]=x then NULL; else x; end if, [seq(i, i=1..nops(someSet))]);
    permVect:= [seq( [  z||(sigma_1[k]) = z||(sigma_j[k]), phi||(sigma_1[k]) = phi||(sigma_j[k]), theta||(sigma_1[k]) = theta||(sigma_j[k])  ], k in non_fixed_points )];
    permVect:= Flatten(permVect);
    KsigmaPoly:= subs(permVect, poly);
    return KsigmaPoly;
end proc:
Kw:= proc(poly, forceVars:=0)
    # alias K_omega
	#Somme sur toutes les permutations de S_N
#Ajout pour obtenir l<action diagonale seulement
    local NbVars, Vars, N, KwPoly, permutations;
	Vars:= giveVars(poly);
	#print(Vars);
	#print(poly);
	NbVars:= max(nops(Vars[1]), nops(Vars[2]), nops(Vars[3])); #Number of variables, (particles in CMS)
	if forceVars <> 0 then N:= forceVars; else N:=NbVars; end if;
	if N < NbVars then print(`Error, not enough variables`); N:= NbVars; end if;
	permutations:= permute(N);
    	#ici
	KwPoly:= map(onePerm->sort(expand(Ksigma_j(poly, onePerm))), permutations);
	#keeping only distinct permutations
	#KwPoly:= map(sort,KwPoly);
	KwPoly:= MakeUnique(KwPoly);
	KwPoly:= add(k, k in KwPoly);
	return KwPoly;
end proc:
kappa_ij:= proc(poly, i,j)
	local kappaij;
	kappaij:= sort(subs({phi||i = phi||j, phi||j=phi||i, theta||i=theta||j, theta||j=theta||i}, poly));
	return kappaij;
end proc:
MakeSubListFromList:= proc( identity_list, permutation_list)
	return [ seq( z||(identity_list[k]) = z||(permutation_list[k]), k=1..nops(identity_list))];
end proc:
AntiSymetrize:= proc(poly, variables_indices::list)
	local Apoly, perms_of_poly, subs_list, omega_sign, omega_perm, identity_perm, nbVars;
	if nops(variables_indices) < 2 then return poly; end if;
	nbVars:= nops(variables_indices);
	identity_perm:= variables_indices;
	omega_perm:= permute(variables_indices);
	omega_sign:= map(x-> sign_of_permutation(x), omega_perm);
	subs_list := map(x-> MakeSubListFromList(identity_perm, x), omega_perm);
	perms_of_poly:= map(x-> subs(x, poly), subs_list);
	Apoly:= add( omega_sign[k]*perms_of_poly[k], k=1..nops(omega_perm));
	return Apoly;
end proc:

Var_Symmetrize:= proc(poly, variables_indices::list)
	local Spoly, perms_of_poly, subs_list, omega_sign, omega_perm, identity_perm, nbVars;
	if nops(variables_indices) < 2 then return poly; end if;
	nbVars:= nops(variables_indices);
	identity_perm:= variables_indices;
	omega_perm:= permute(variables_indices);
	subs_list := map(x-> MakeSubListFromList(identity_perm, x), omega_perm);
	perms_of_poly:= map(x-> subs(x, poly), subs_list);
	Spoly:= add(perms_of_poly[k], k=1..nops(omega_perm));
	return Spoly;
end proc:
# Dunkl-Cherednik
Dunkl_i_n:= proc(i, n, poly, Nin:=0, parameter:= beta)
    local something, Vars, N, term1, term2, term3, term4, Dpoly, beta;
    option remember;
	#The Dunkl-like-Cherednik-... operator, still pretty clear code
	# i: the Dunkl index, n: the exponent, Nin: number of variables, parameter : beta, 1/alpha ...
	if n = 0 then return poly; end if;
	if n < 0 then return 0; end if;
	Vars:= giveVars(poly);
	N:= max(nops(Vars[1]), nops(Vars[2]), nops(Vars[3])); #Number of variables, (particles in CMS)
    	if Nin > N then N:= Nin; end if;
	if parameter <> beta then beta := parameter end if;
	term1:= z||i * diff(poly, z||i);
	term2:= beta*add(z||i/(z||i - z||j)*(poly- Kzij(poly,i,j)), j in seq(k, k=(i-1)..1, -1));
	term3:= beta*add(z||j/(z||i - z||j)*(poly- Kzij(poly,i,j)), j in seq(k, k=(i+1)..N, 1));
	term4:= -1*beta*(i-1)*poly;
	Dpoly:= factor(term1 + term2+ term3 + term4);

	if n=1 then return Dpoly; end if;
	if n>1 then return Dunkl_i_n(i,n-1, expand(sort(Dpoly)), Nin, beta); end if;
end proc:

sign_of_permutation:= proc(L::list)
	local n, N, i, j;

	n:=nops(L);
	if nops(convert(L,set))=n then
		N:=0;
		for i to n-1 do
			for j from i+1 to n do
				if L[i]>L[j] then N:=N+1; fi;
			end do;
		end do;
	else print(`Invalid input`);
	end if;

	return (-1)^N;
end proc:

#Operations regarding the representation of S_N only with the generators s_i
get_list_perm_rep:= proc(sigma, N)
	return map(x->PermApply(sigma,x), [seq(k, k=1..N)]);
end proc:

inversion_number:= proc(k, a_list)
	local inversion_number, chopped_list, inversions;
	if k = nops(a_list) then return 0; end if;
	if nops(a_list) = 0 then return 0; end if;
	chopped_list:= a_list[(k+1)..-1];
	inversions:= map(x-> if x < a_list[k] then 1; else 0; end if, chopped_list);
	return add(l, l in inversions);
end proc:

Total_inversions:= proc(a_list)
	local all_inversions;
	all_inversions:=map(x->inversion_number(x,a_list), [seq(k, k=1..(nops(a_list)-1))]);
	return add(k, k in all_inversions);
end proc:

gen_coxeter_rep_of_SN:=proc(N)
	local all_perms, adj_rep, generators, longest_chain, N_terms, permutations, number_of_inversions, perms_with_numb_of_inv, wanted_perms, choose_from, choices, numb_of_inv, new_elements, stop_it, a_choice, serie_of_perm, effective_perm, permuted, where_is_it; remember;
	#Trivial terms [0] is the identity and is not included in the generators;
	generators:= [ seq([k], k=1..(N-1))];
	adj_rep:= [[0],op(generators)];
	#This is the chain that reverses completly [1...N],
	#The idea is that you carry the 1 up to the end, then carry the 2 up to the 1, and so on and so forth
	#We see that this as exactly N(N-1)/2 terms, just like the maximal number of inversions.
	longest_chain:= [ seq(seq( j, j=1..i), i=(N-1)..1,-1)];
	adj_rep:=[op(adj_rep), longest_chain];
	N_terms:= [seq(k, k=1..N)];
	permutations:= permute(N);
	number_of_inversions:= map(Total_inversions, permutations);
	#We start the folling loop at 2 since we have the identity and the generators which have
	#total number of inversion of 0 and 1
	#We also skip the last length in the number of inversion since we have longest_chain
	for numb_of_inv from 2 to (N*(N-1)/2-1) do
		perms_with_numb_of_inv := ListTools:-SearchAll(numb_of_inv, number_of_inversions);
		wanted_perms := map(x-> op(x,permutations), [perms_with_numb_of_inv]);
		choose_from:= [ seq( op(generators), k=1..ceil(numb_of_inv/2)) ];
		choices:= permute(choose_from, numb_of_inv);
		choices:= map(void_redundant_list, choices);
		new_elements:=[];
		stop_it:= nops(wanted_perms);
		for a_choice in choices do
			serie_of_perm:= list_of_index_to_perm(a_choice);
			effective_perm:= prod_of_perm(serie_of_perm);
			permuted:= get_list_perm_rep(effective_perm, N);
			where_is_it:=ListTools:-Search(permuted, wanted_perms);
			if where_is_it <> 0 then
				new_elements:= [op(new_elements), a_choice];
				wanted_perms:= subsop(where_is_it=NULL, wanted_perms);
			end if;
			if nops(new_ele) = stop_it then break; end if;
		end do:
		adj_rep:=[op(adj_rep), op(new_elements)];
	end do:
		adj_rep:= map(x-> Flatten(x), adj_rep);
	return adj_rep;
end proc:

restricted_coxeter_rep:= proc(n1,n2,N)
	local SNrep, unallowed, valid_terms;
	SNrep:= gen_coxeter_rep_of_SN(N);
	#We remove the identity
	#SNrep:= SNrep[2..-1];
	unallowed:= { seq(k, k=1..(n1-1)), seq(k, k=(n2)..N)};
	valid_terms:=map(x-> if {op(x)} intersect unallowed = {} then x; else NULL; end if, SNrep);
	return valid_terms;
end proc:

SNcheckCoxeter:= proc(mbarbar, mbar, mubar, N)
	local S1, S2, S3, S4, Stot, doneS, Scheck;
	S1:= restricted_coxeter_rep(1,mbarbar,N);
	S2:= restricted_coxeter_rep(1+mbarbar,mbarbar+mbar,N);
	S3:= restricted_coxeter_rep(1+mbarbar+mbar,mbarbar+mbar+mubar,N);
	S4:= restricted_coxeter_rep(1+mbarbar+mbar+mubar,N,N);
	Stot:= {op(gen_coxeter_rep_of_SN(N))};
	doneS:= {op(S1), op(S2), op(S3), op(S4)};
	Scheck:= Stot minus doneS;
	Scheck:= [[0], op(Scheck)];
	return Scheck;
end proc:

void_redundant_list:= proc(list)
	local k;
	for k from 1 to nops(list)-1 do
		if list[k] = list[k+1] then return NULL;
		elif k+3 < nops(list) then
			if list[k..k+1] = list[k+2..k+3] then return NULL; end if;
		end if
	end do:
	return list;
end proc:

index_coxeter_gen_to_perm:= proc(a_gen)
	local the_val;
	the_val:= op(a_gen);
	return Perm([[the_val, the_val+1]]);
end proc:

list_of_index_to_perm:= proc(a_list)
	return map(index_coxeter_gen_to_perm, a_list);
end proc:

prod_of_perm:= proc(list_of_perm)
	local k, out;
	out:= Perm([]);
	for k from 1 to nops(list_of_perm) do
		out:=PermProduct(out,list_of_perm[k]);
	end do:
end proc:

#Operators related to nsMacdo
#T_sigma operation
opTi:= proc(i, expr_z, N)
	local expr, out, out1, out2, out2b;
	if i = N then return 0; end if;
	if i = 0 then return expr_z; end if;
	expr:= expr_z;
	out1:= t*expr;
	out2:= ((t*z||i - z||(i+1))/(z||i-z||(i+1)));
	out2b:= Kzij(expr, i, i+1) - expr;
	out2:= out2*out2b;
	out:= out1 + out2;
#	out:= t*expr + ((t*z||i - z||(i+1))/(z||i-z||(i+1)))*(Kzij(expr,i, i+1) - expr);
	#out:= simplify(factor(out));
	return out;
end proc:

inv_opTi:= proc(i, expr_z, N)
	local expr, out;
	return simplify(factor((t^(-1) - 1)*expr_z + t^(-1)*opTi(i, expr_z, N)));
end proc:

opT_sigma:= proc(sigma, expr, N)
	local out,k;
	out:= expr;
	for k from 1 to nops(sigma) do
		out:= factor(opTi(sigma[-k], out, N));
	end do:
	return out;
end proc:

tSymmetrize:= proc(expr_z, N, n1:=0, n2:=0)
	local coxeterSN, out, sigma;
	if n1 = 0 then
		coxeterSN:= gen_coxeter_rep_of_SN(N);
	else
		coxeterSN:= restricted_coxeter_rep(n1,n2,N);
	end if;
	out:= 0;
	for sigma in coxeterSN do
		out:= out + opT_sigma(sigma, expr_z,N);
	end do:
	out:= factor(out);
	return out;
end proc:
tVandermonde:= proc(the_range, remove_t:=0)
	local N, out;
	N:= the_range[2]-the_range[1];
	out:= t^(-N*(N-1)/2)*mul(mul( t*z||i - z||j, i=the_range[1]..(j-1)), j=the_range[1]..the_range[2]);
	if remove_t <> 0 then out:= subs(t=1, out); end if;
	return out;
end proc:
tantiSymmetrize:= proc(expr_z, n1:=0, n2:=0)
	local coxeterSN, out, sigma;
	out:= AntiSymetrize(expr_z, [seq(k,k=n1..n2)]);
	out:= tVandermonde([n1,n2])/tVandermonde([n1,n2],1)*out;
	out:= factor(out);
	return out;
end proc:
bkp_tantiSymmetrize:= proc(expr_z, N, n1:=0, n2:=0)
	local coxeterSN, out, sigma, permsign;
	if n1 = 0 then
		coxeterSN:= gen_coxeter_rep_of_SN(N);
	else
		coxeterSN:= restricted_coxeter_rep(n1,n2,N);
	end if;
	out:= 0;
	for sigma in coxeterSN do
		permsign:= nops(sigma);
		out:= out + (-1*t)^permsign*opT_sigma(sigma, expr_z,N);
	end do:
	out:= factor(out);
	return out;
end proc:
apply_sigma_si:= proc(expr, sigma)
	local out, k,i;
	if sigma = [0] then return expr; end if;
	out:= expr;
	for k from 1 to nops(sigma) do
		i := sigma[-k];
		out:= Kscrypt_ij(out, [i, i+1])
	end do:
	return out;
end proc:

var_omega:= proc(expr, N)
	local out, k;
	out:= opTaui(1,expr);
	for k from 1 to (N-1) do
		out:= Kzij(out,k,k+1);
	end do:
	return out;
end proc:

opTaui:= proc(i,expr)
	return factor(subs(z||i = q*z||i, expr));
end proc:

opYi:= proc(i, expr, N)
	local out, k;
	out:=expr;
	if i > 1 then
		for k from (i-1) to 1 by -1 do
			out:= inv_opTi(k, out, N);
		end do:
	end if;
	out:= var_omega(out, N);
	if i <> N then
		for k from (N-1) to i by -1 do
			out:= opTi(k,out, N);
		end do:
	end if;
	out:=factor(t^(-N+i)*out);
	return out;
end proc:

eigen_val_opYi:= proc(i, eta)
	local qexp,leta1,leta2,leta, out;
	qexp:= eta[i];
	if i >1 then
		leta1:= map(x-> if x>= eta[i] then 1; else 0; end if, eta[1..(i-1)]);
	else
		leta1:=0;
	end if;
	if i<nops(eta) then
		leta2:= map(x-> if x>eta[i] then 1; else 0; end if, eta[(i+1)..-1]);
	else
		leta2:= 0;
	end if;
	leta:= add(k, k in leta1) + add(k, k in leta2);
	out:= q^(qexp)*t^(-leta);
	return out;
end proc:


text:
# Acceleration Calculation for the Robot based on MDH frames
# Introduction
# Berechnung der Beschleunigung 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
# acceleration -> Berechnung der Beschleunigung aller Segmente
# linkframe -> Berechnung der Geschwindigkeit im Körper-KS (KSi)

# Sources
# [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":
printf("Generiere Beschleunigung für %s (Herleitung im Körper-KS)\n", 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(5)=-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
    if NJ <= 5 then # bei einfachen Systemen vereinfachungen vornehmen
      use_simplify := 1
    else # bei längeren kinematischen Ketten dauert es zu lange / bringt zu wenig
      use_simplify := 0:
    end if:
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(5): # fünfter Eintrag ist für Beschleunigung
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:
# Ergebnisse der Geschwindigkeit laden
read sprintf("../codeexport/%s/tmp/velocity_linkframe_floatb_%s_maple.m", robot_name, base_method_name):
omega_i_i:=omega_i_i:
rD_i_i:=rD_i_i:
# 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:
dD :=dD:
# Zeitableitungen der Geschwindigkeit laden.
# Die Berechnung soll nur an einer Stelle erfolgen. Siehe robot_tree_acceleration_mdh_angles.mw.
read sprintf("../codeexport/%s/tmp/acceleration_mdh_angles_maple.m", robot_name):
thetaDD :=thetaDD:
dDD := dDD:

# Calculate Acceleration
# First assume fixed base model with base velocity and acceleration set to zero
# Anfangsgeschwindigkeit definieren für floating base model
# Acceleration of Frames
#
# Beschleunigung der Basis der Koordinatensysteme der Körper (ausgedrückt im Körper-KS)
rDD_i_i := Matrix(3, NL):
# Winkelbeschleunigung von Körper i ausgedrückt im Körper-KS
omegaD_i_i := Matrix(3, NL):
# Anfangsgeschwindigkeit der Basis:
# twist: Basis-Geschwindigkeit bzgl Welt-KS ausgedrückt im Basis-KS
# eulangrpy: 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:
  rDD_i_i(1..3,1) := diff~(V_base_t(1..3,1),t):
  omegaD_i_i(1..3,1) := diff~(V_base_t(4..6,1),t):
end:
if base_method_name = "eulxyz" then:
  rDD_i_i(1..3,1) := Transpose(Trf_c(1..3, 1..3, 1)) . diff~(V_base_t(1..3,1),t):
  omegaD_i_i(1..3,1) := Transpose(Trf_c(1..3, 1..3, 1)) . T_basevel . diff~(V_base_t(4..6,1),t):
end:
# Erhöhe die Beschleunigung 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
  # [Ortmaier](7.14) (S.99) & (7.19) (S.99)
  if sigma(i) = 0 then # Drehgelenk
    omegaD_i_i(1 .. 3, i+1) := Matrix(Multiply(R_i_j,( Matrix(3,1,omegaD_i_i(1 .. 3, j))))) + Matrix( thetaDD(i,1)*<0;0;1>) + Matrix(CrossProduct(R_i_j.Matrix(3,1,omega_i_i(1 .. 3, j)),thetaD(i,1)*<0;0;1>)) :
  else: # Schubgelenk
    omegaD_i_i(1 .. 3, i+1) := Multiply(R_i_j, Matrix(3,1,omegaD_i_i(1 .. 3, j))):
  end if:
  # Terme vereinfachen (Teil 1)
  if use_simplify>=1 then
    tmp_t11:=time():
    tmp_l11 := length(omegaD_i_i(1 .. 3, i+1)):
    omegaD_i_i(1 .. 3, i+1) := simplify2(omegaD_i_i(1 .. 3, i+1)):
    tmp_l21 := length(omegaD_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
  # [Ortmaier](7.17) (S.99) & (7.22) (S.100)
  if sigma(i) = 0 then # Drehgelenk
    rDD_i_i(1 .. 3, i+1) := Multiply( R_i_j, ( rDD_i_i(1 .. 3, j) + omegaD_i_i(1 .. 3, j) &x r_j_j_i + CrossProduct(omega_i_i(1..3,j),(omega_i_i(1..3,j) &x r_j_j_i )))) :
  else: # Schubgelenk
    rDD_i_i(1 .. 3, i+1) := Matrix( Multiply( R_i_j, ( rDD_i_i(1 .. 3, j)  + CrossProduct(omega_i_i(1..3,j), (omega_i_i(1..3,j) &x r_j_j_i )) + CrossProduct(omegaD_i_i(1..3,j), r_j_j_i)))) + Matrix( 2* dD(i,1) * CrossProduct(R_i_j.omega_i_i(1..3,j),<0;0;1>)) + Matrix(dDD(i,1)*<0;0;1>)  :
  end if:
  printf("Beschleunigung für Körperkoordinatensystem %d aufgestellt (Herleitung im Körper-KS). %s\n", i, FormatTime("%Y-%m-%d %H:%M:%S")): #0=Basis
  # Terme vereinfachen (Teil 2)
  if use_simplify>=1 then
    tmp_t12:=time():
    tmp_l12 := length(rDD_i_i(1 .. 3, i+1)):
    rDD_i_i(1 .. 3, i+1) := simplify2(rDD_i_i(1 .. 3, i+1)):
    tmp_l22 := length(rDD_i_i(1 .. 3, i+1)):
    tmp_t22:=time():
    printf("%s: Terme für Beschleunigungen 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:

# Acceleration of Center of Mass

rDD_i_Si := Matrix(3, NL):

for i to NL do
  rDD_i_Si(1 .. 3, i) := Matrix(rDD_i_i(1 .. 3, i))+Matrix(CrossProduct(omegaD_i_i(1 .. 3, i), r_i_i_Si(1 .. 3, i)))+Matrix(CrossProduct(omega_i_i(1 .. 3, i), CrossProduct(omega_i_i(1 .. 3, i), r_i_i_Si(1 .. 3, i)))):
  printf("Beschleunigung für Körperschwerpunkt %d aufgestellt. %s\n", i-1, FormatTime("%Y-%m-%d %H:%M:%S")):
  # Terme vereinfachen
  if use_simplify>=1 then
    tmp_t1:=time(): tmp_l1 := length(rDD_i_Si(1 .. 3, i)):
    rDD_i_Si(1 .. 3, i) := simplify2(rDD_i_Si(1 .. 3, i)):
    tmp_t2:=time(): tmp_l2 := length(rDD_i_Si(1 .. 3, i)):
    printf("%s: Terme für Schwerpunkts-Beschleunigungen vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
end do:
# Export
# Maple Export
save omegaD_i_i, rDD_i_i,rDD_i_Si , sprintf("../codeexport/%s/tmp/acceleration_linkframe_floatb_%s_maple.m", robot_name, base_method_name):
printf("Maple-Ausdrücke exportiert. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):

# Matlab Export

if codegen_act then
  MatlabExport(convert_t_s(omegaD_i_i), sprintf("../codeexport/%s/tmp/acceleration_omegaDii_floatb_%s_linkframe_matlab.m", robot_name, base_method_name), codegen_opt):
  MatlabExport(convert_t_s(rDD_i_i), sprintf("../codeexport/%s/tmp/acceleration_rDDii_floatb_%s_linkframe_matlab.m", robot_name, base_method_name), codegen_opt):
  MatlabExport(convert_t_s(rDD_i_Si), sprintf("../codeexport/%s/tmp/acceleration_rDDiSi_floatb_%s_linkframe_matlab.m", robot_name, base_method_name), codegen_opt):
  printf("Beschleunigung in Matlab exportiert. %s\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_vxc *)
(* prefix:
  gga_x_lb_params *params;

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

lb_f0 := (rs, z, x) -> -my_piecewise3(x < 300,
         params_a_beta*x^2/(1 + 3*params_a_beta*x*arcsinh(params_a_gamma*x)),
         x/(3*log(2*params_a_gamma*x))):

lb_f := (rs, z, x) -> (params_a_alpha*(4/3)*LDA_X_FACTOR + lb_f0(rs, z, x))*n_spin(rs, z)^(1/3):

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


text:
# This makes the formal group law for integral Morava K-theory
# We truncate powers d+1 and higher

make_fgl_Morava := proc(p,n,d)
 local T,ld,x,y,u,i;

 T := table();
 T["p"] := p;
 T["n"] := n;
 T["d"] := d;
 T["height"] := n;
 T["degree"] := d;
 T["modulus"] := 0;

 ld := 0;
 while p^ld <= d/p do ld := ld+1; od;
 T["log_degree"] := ld;

 T["log"] := unapply(add(x^(p^(n*i))/p^i,i=0..ld),x);
 T["p_series"] := unapply(x^(p^n),x);

 Order := d+2;

 T["exp"] :=
  unapply(
   convert(
    solve(x = series(T["log"](y) + sin(y)^(d+1),y=0,d+1),y),
    polynom,x
   ),
   x
  );

 T["p_series"] :=
  unapply(expand(convert(series(T["exp"](p*T["log"](x)),x=0,d+1),polynom,x)),x);

 T["sum"] :=
  unapply(
   subs(u=1,expand(convert(series(T["exp"](T["log"](u*x)+T["log"](u*y)),u=0,d+1),polynom,u))),
   x,y
  ):

 T["sum0"] := unapply(mods(T["sum"](x,y),p),x,y);
 T["p_series0"] := unapply(mods(T["p_series"](x),p),x);

 return eval(T);
end:


text:
# This was originally the ggplot style from matplotlib
# from http://www.huyng.com/posts/sane-color-scheme-for-matplotlib/

# More options here
# https://matplotlib.org/users/customizing.html

patch.linewidth: 0.5
patch.facecolor: 348ABD  # blue
patch.edgecolor: EEEEEE
patch.antialiased: True
mathtext.default : regular

font.size: 6.0
font.family : sans-serif
font.sans-serif : DejaVu Sans #Arial also had a unicode problem
font.style : normal
#font.family: Arial # specifiying this because Bitstream doesnt show unicode

# Font Sizes
# xx-small, x-small, small, medium, large, x-large, xx-large, smaller, larger.
axes.facecolor: white #E5E5E5
#axes.edgecolor: white
axes.edgecolor: black
axes.spines.right : False
axes.spines.top : False
axes.linewidth: 0.5
axes.grid: False
axes.titlesize: large
axes.labelsize: medium
axes.labelcolor: 555555
axes.axisbelow: True       # grid/ticks are below elements (e.g., lines, text)

# matlab color cycle
axes.prop_cycle: cycler('color', ['0071BC', 'D85218', 'ECB01F', '7D2E8D', '76AB2F', '4CBDED', 'A1132E'])

#['E24A33', '348ABD', '988ED5', '777777', 'FBC15E', '8EBA42', 'FFB5B8'])
                   # E24A33 : red
                   # 348ABD : blue
                   # 988ED5 : purple
                   # 777777 : gray
                   # FBC15E : yellow
                   # 8EBA42 : green
                   # FFB5B8 : pink

xtick.color: 555555
xtick.major.width: 0.5
#xtick.direction: in
xtick.direction: out

ytick.color: 555555
ytick.major.width: 0.5
#ytick.direction: in
ytick.direction: out

#grid.color: white
grid.color: lightgray
grid.linestyle: :    # solid line

figure.facecolor: FFFFFF #white
figure.edgecolor: 0.50

figure.subplot.wspace  : 0.2

#legend.labelspacing
#legend.borderpad
legend.fontsize : small

max_stars_repo_name: npmurphy/CRNSynthesisFigures

max_issues_repo_name: npmurphy/CRNSynthesisFigures

max_forks_repo_name: npmurphy/CRNSynthesisFigures


text:
# Number of subspaces of dimension d in F^n, if |F|=q

grassmann_count := (q::posint) -> proc(n,d)
 local i;
 return mul(q^(n-i)-1,i=0..d-1)/mul(q^(d-i)-1,i=0..d-1);
end:

# Number of splittings of F^(n+m) as a sum of subspaces of
# dimensions n and m, if |F| = q.

splitting_count := (q::posint) -> (n,m) ->
 grassman_count(q)(n+m,n) * q^(n*m);

# tanabe_N(p,n,r)(k) is the number of irreducible F-linear representations
# of Z_p^n that have dimension p^k over F.

tanabe_N := (p,n,r) -> (k) -> `if`(k = 0,p^(n*r),(p^n-1) * p^((n-1)*k + n*(r-1)));

# The socle of K^0(BGL_{p^k}(F)) is generated by
# euler(V)^{N-1}, where N = tanabe_N_bar(n,v,k)

tanabe_N_bar := (p,n,r) -> proc(k)
 local j;
 p^(n*r) + (p^n-1)*p^(n*(r-1))*sum(p^((n-1)*j),j=1..k);
end:

tanabe_N_star := (p,n,r) -> proc(k)
 local j;
 p^(n*r-1) * (p-1) * sum(p^((n-1)*j),j=0..k-1);
end:

tanabe_N_lim := (p,n,r) -> (k) -> p^(n*r+n*k-k) * (1-p^(-n))/(1-p^(1-n));

# The natural numbers can be partitioned as the disjoint union of
# the sets tanabe_NN(p,n,r)(k), of size tanabe_N(p,n,r)(k).

tanabe_NN := (p,n,r) -> proc(k::nonnegint)
 local i;
 if k = 0 then
  return [seq(i,i=0..tanabe_N_bar(p,n,r)(0)-1)];
 else
  return [seq(i,i=tanabe_N_bar(p,n,r)(k-1)..tanabe_N_bar(p,n,r)(k)-1)];
 fi;
end:

# tanabe_s(p,n,r)(m) is the unique k such that m lies in tanabe_NN(p,n,r)(k).
tanabe_s := (p,n,r) -> proc(m)
 local k;
 k := 0;
 while m >= tanabe_N_bar(p,n,r)(k) do
  k := k+1;
 od;
 return k;
end:

# tanabe_X(p,n,r)(k) is the set of generators for E^0(BGL_*(F)) in
# degree p^k.

tanabe_X := (p,n,r) -> proc(k)
 local j;
 return [seq(c[p^k] ^ j, j = 0..tanabe_N(p,n,r)(k)-1)];
end:

# tanabe_prim_gen(p,n,r)(k) is a generator for the primitives
# in K^0(BGL_*(F)) in degree p^k
tanabe_prim_gen := (p,n,r) -> (k) -> c[p^k] ^ tanabe_N_bar(p,n,r)(k-1);

# tanabe_soc_gen(p,n,r)(k) is a generator for the socle of
# K^0(BGL_{p^k}(F))
tanabe_soc_gen  := (p,n,r) -> (k) -> c[p^k] ^ (tanabe_N_bar(p,n,r)(k) - 1);

# tanabe_m(p)(k) has p-adic valuation k and divides
# tanabe_m(p)(k+1).  Moreover, any positive integer n divides
# tanabe_m(p)(k) when k is sufficiently large.

tanabe_m := (p) -> proc(k)
  local m,i;
  m := 1;
  i := 1;
  while ithprime(i) <= p + k do
    m := m * ithprime(i) ^ k;
    i := i + 1;
  od:
  return m;
end:

check_tanabe_m := proc(p,k)
  _ASSERT(padic_val(tanabe_m(p)(k),p) = k,
          "v_p(m[k])");

  _ASSERT(type(tanabe_m(p)(k+1)/tanabe_m(p)(k),posint),
          "m[k] divides m[k+1]");
end:

# Given a polynomial f over the integers and a prime q, return
# the first irreducible factor of f mod q, in the standard
# Maple order.

first_factor := proc(f,q)
 local F;
 F := :-Factor(f) mod q;
 if type(F,`*`) then F := [op(F)] else F := [F]; fi:
 F := sort(map(g -> `if`(type(g,`^`),op(1,g),g),F));
 return F[1];
end:

# Return an irreducible polynomial f[k](x) of degree v0 * p^k over
# Z/q0.  This gives a field F(k) of degree p^k over F = F(0),
# where F(0) has order q = q0 ^ v0.  The polynomials are arranged
# so that the canonical primitive element x[k] in F(k) generates
# the group of units, and there is an inclusion F(k) -> F(k+1)
# sending x[k] to x[k+1] ^ m, where m = |F(k+1)^x|/|F(k)^x|.
tanabe_FF_poly := (q0,v0,p) -> proc(k::nonnegint)
 option remember;
 local q,m,f;

 q := q0 ^ v0;
 if k = 0 then
  f := NumberTheory[CyclotomicPolynomial](q0 ^ v0 - 1,x);
 else
  m := (q ^ (p ^ k) - 1)/(q ^ (p ^ (k - 1)) - 1);
  f := subs(x = x ^ m, tanabe_FF_poly(q0,v0,p)(k - 1));
 end:

 return first_factor(f,q0);
end:

# tanabe_F(q0,v0)(m) is the field denoted by F[m] in the paper,
# where q = q0 ^ v0 with q0 prime and v0 > 0.
# tanabe_FF(q0,v0)(p,k) is F(k) = F[p^k].

tanabe_F := (q0,v0) -> (m) -> GF(q0,v0 * m);
tanabe_FF := (q0,v0) -> (p,k) -> GF(q0,v0 * p ^ k,tanabe_FF_poly);

# Orders of various groups: G is the general linear group,
# N is the subgroup of monomial matrices, H is the block
# diagonal subgroup of p' index as in prop-non-p-power.

tanabe_G_order := (q) -> proc(d)
 local k;
 return mul(q ^ d - q ^ k,k = 0..d-1);
end:

tanabe_N_order := (q) -> (d) -> d! * (q - 1)^d;

tanabe_H_order := (q,p) -> proc(d)
 local dd,m,i;
 dd := digit_list(d,p);
 m := nops(dd);
 mul(tanabe_G_order(q)(p ^ i) ^ dd[i+1], i = 0..m-1);
end:

# The E^2 page of the AHSS converging to K_*(BV) has a trigrading
# with E^2_{ijk} = H_i(BV_k;K_j).  It has even generators
# b[i] and u, and odd generators e[i] with
#   |b[i]| = tanabe_b_trideg(p,n,r)(i) = [2i,0,1]
#   |e[i]| = tanabe_e_trideg(p,n,r)(i) = [2i+1,0,1]
#   |u|    = tanabe_u_trideg(p,n,r)    = [0,2,0].

tanabe_b_trideg  := (p,n,r) -> (i) -> [2*i,0,1];
tanabe_e_trideg  := (p,n,r) -> (i) -> [2*i+1,0,1];
tanabe_u_trideg  := (p,n,r) -> [0,2,0];

# The only nontrivial differentials have the form d_r, where
# r = tanabe_AHSS_page(p,n,r)(k) for some k >= 0
# [From now on we refer to this as the k'th differential,
#  and the page on which it occurs as the k'th page.]

tanabe_AHSS_page := (p,n,r) -> (k) -> 2*p^(n*(r+k)) - 1;

# Tridegree of the k'th differential
tanabe_AHSS_d_trideg := (p,n,r) -> (k) ->
 [1-2*p^(n*(r+k)),2*p^(n*(r+k))-2,0];

# The k'th page can be described in terms of generators
# bb[m,i] and ee[m,i]

tanabe_bb_defined := (p,n,r) -> (k) -> (m,i) ->
  m >= min(k,tanabe_s(p,n,r)(i));

tanabe_bb_indec := (p,n,r) -> (k) -> (m,i) ->
  m = min(k,tanabe_s(p,n,r)(i));

tanabe_ee_defined := (p,n,r) -> (k) -> (m,i) ->
  m = k and i >= tanabe_N_star(p,n,r)(k);

tanabe_bb_trideg := (p,n,r) -> (m,i) -> [2*p^m*i,0,p^m];
tanabe_ee_trideg := (p,n,r) -> (m,i) -> [2*p^m*i + 1,0,p^m];

tanabe_dbb0 := (p,n,r) -> (k) -> (m,i) ->
  ee[k,i - tanabe_N_bar(p,n,r)(k) + tanabe_N_star(p,n,r)(k)] *
    u ^ (p^(n * (k + r)) - 1);

tanabe_dbb := (p,n,r) -> (k) -> (m,i) ->
  `if`(m < k, 0, `if`(i < tanabe_N_bar(p,n,r)(k), 0, tanabe_dbb0(p,n,r)(k)(m,i)));

# This is the part of the inverse Poincare series contributed by
# irreducible representations of dimension at most p^k.

tanabe_inverse_poincare_series := (p,n,r) -> proc(k)
 local j;
 (1-t)^(p^(n*r)) * mul((1 - t^(p^j))^(p^(n*(r+j))*(1-1/p^n)),j=1..k);
end:

tanabe_poincare_series := (p,n,r) -> (k) ->
 convert(series(1/tanabe_inverse_poincare_series(p,n,r)(k),t=0,p^k+1),
         polynom,t);

tanabe_trideg := (p,n,r) -> proc(x)
 local d,d0,y;

 if type(x,indexed) then
  if op(0,x) = b then
   return tanabe_b_trideg(p,n,r)(op(x));
  elif op(0,x) = e then
   return tanabe_e_trideg(p,n,r)(op(x));
  elif op(0,x) = bb then
   return tanabe_bb_trideg(p,n,r)(op(x));
  elif op(0,x) = ee then
   return tanabe_ee_trideg(p,n,r)(op(x));
  else
   return [0,0,0];
  fi;
 elif x = u then
  return tanabe_u_trideg(p,n,r);
 elif type(x,`+`) then
  d := map(tanabe_trideg(p,n,r),{op(x)});
  if nops(d) = 1 then
   return op(d);
  else
   return FAIL;
  fi;
 elif type(x,`*`) then
  d := [0,0,0];
  for y in [op(x)] do
   d0 := tanabe_trideg(p,n,r)(y);
   if d0 = FAIL then return FAIL; fi;
   d := d +~ d0;
  od;
  return d;
 elif type(x,`^`) then
  return tanabe_trideg(p,n,r)(op(1,x)) *~ op(2,x);
 else
  return [0,0,0];
 fi;
end:

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

# Check v_p(q ^ j - 1) = r + v_p(j)
check_lem_vp := (p) -> proc(q::posint,j::posint)
 local r;
 r := padic_val(q - 1,p);
 _ASSERT(r > 0, "r > 0");

 _ASSERT(padic_val(q ^ j - 1,p) = r + padic_val(j,p),
         sprintf("v_{%d}(%d ^ %d - 1)",p,q,j));
end;

check_tanabe_numbers := proc()
 local n,r,k,p,L,N,N_bar,N_star,N_lim,P,err;

 assume(n::posint);
 assume(r::posint);
 assume(k::posint);
 assume(p::posint);

 N      := tanabe_N(p,n,r);
 N_bar  := tanabe_N_bar(p,n,r);
 N_star := tanabe_N_star(p,n,r);
 N_lim  := tanabe_N_lim(p,n,r);

 _ASSERT(
  simplify(N_bar(k) - N_bar(k-1) - N(k)) = 0,
  "N_bar(k) - N_bar(k-1) = N(k)"
 );

 _ASSERT(
  simplify(N_star(k) - (N_bar(k) - p^(n*k + n*r - k))) = 0,
  "N_star(k) = N_bar(k) - p^(n*k + n*r - k)"
 );

 _ASSERT(
  simplify(p^k * (N_bar(k) - N_star(k)) - p^(n * (k + r))) = 0,
  "p^k * (N_bar(k) - N_star(k)) = p^(n * (k + r))"
 );

 _ASSERT(
  simplify(p^(k + 1) * (N_bar(k) - N_star(k + 1)) - p^(n * (k + r))) = 0,
  "p^(k + 1) * (N_bar(k) - N_star(k + 1)) = p^(n * (k + r))"
 );

 _ASSERT(
  simplify(N(0) + sum(p^j * N(j),j=1..k) - p^(n*(r+k))) = 0,
  "sum(p^j * N(j),j=0..k) = p^(n*(r+k))"
 );

 _ASSERT(
  simplify((1 - p^(-k)) * N(0) +
           sum(expand((1 - p^(j - k)) * N(j)),j=1..k) - N_star(k)) = 0,
  "sum((1 - p^(j-k)) * N(j),j=0..k) = N_star(k)"
 );

 _ASSERT(
  simplify(p^(k + 1) * N_star(k + 1) - p^k * N_star(k) - (p - 1) * p^k * N_bar(k)) = 0,
  "p^(k + 1) * N_star(k + 1) = p^k * N_star(k) + (p - 1) * p^k * N_bar(k)"
 );

 _ASSERT(
  simplify(factor(N_bar(k) - (N_star(k)*(p^n-1)/(p-1) + p^(n*r)))) = 0,
  "N_bar(k) = N_star(k)*(p^n-1)/(p-1) + p^(n*r)"
 );

 _ASSERT(
  simplify(sum(p^j*(p-1)*N_bar(j),j=0..k-1) - p^k*N_bar(k) + p^(n*(r+k))) = 0,
  "sum(p^j*(p-1)*N_bar(j),j=0..k-1)"
 );

 _ASSERT(
  limit(simplify(N_bar(k)/N_lim(k)-1),k=infinity) = 0,
  "N_bar(k) is asymptotic to N_lim(k)"
 );
end:

check_factorial_congruences := proc()
 local p,i;
 p := ithprime(20);

 _ASSERT(
  {seq(modp((-1)^i * (i - 1)! * (p - i)!,p),i=1..p-1)} = {1},
  "(-1)^i * (i - 1)! * (p - i)! = 1 mod p"
 );

 _ASSERT(
  {seq(mods((-1)^i * binomial(p-1,i-1),p),i=1..p-1)} = {-1},
  "(-1)^i * binomial(p-1,i-1) = -1 mod p"
 );
end:

check_valuations := proc()
 local q0,v0,q,p,r,d,m;
 q0 := 19;
 v0 := 1;
 q  := q0 ^ v0;
 p  := 3;
 r  := padic_val(q - 1, p);
 d  := 20;

 _ASSERT(
  padic_val(tanabe_G_order(q)(d),p) - (d * r + (d - digit_sum(d,p))/(p - 1)) = 0,
  "valuation of |G(d)|"
 );

 m := tanabe_G_order(q)(d) / tanabe_N_order(q)(d):
 _ASSERT(type(m,posint) and padic_val(m,p) = 0,
  "N(d) has p' index in G(d)"
 );

 m := tanabe_G_order(q)(d) / tanabe_H_order(q,p)(d):
 _ASSERT(type(m,posint) and padic_val(m,p) = 0,
  "H(d) has p' index in G(d)"
 );

end:

check_diff_trideg := proc()
 local p,n,r,k,m,i,err,d0,d1,d2;

 d0 := tanabe_bb_trideg(p,n,r)(k,i);
 d1 := tanabe_AHSS_d_trideg(p,n,r)(k);
 d2 := tanabe_trideg(p,n,r)(tanabe_dbb0(p,n,r)(k)(m,i));

 err := simplify(d0 + d1 - d2, symbolic);

 _ASSERT(err = [0,0,0], "differential respects tridegrees");
end:


text:
(* type: work_lda *)

a :=  0.0311:
b := -0.048:
c :=  0.009:
d := -0.017:

f_ls := (rs, zeta) -> a*log(rs) + b + c*rs*log(rs) + d*rs:
f_ms := (rs, zeta) -> -0.06156 + 0.01898*log(rs):
f_hs := (rs, zeta) -> 0.438/rs + 1.325/rs^(3/2) - 1.47/rs^2 - 0.4/rs^(5/2):

f := (rs, zeta) -> convert(piecewise(
  rs < 0.7, f_ls(rs, zeta),
  rs < 10,  f_ms(rs, zeta),
  f_hs(rs, zeta)), 'Heaviside'):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
(*
 Copyright (C) 2017 M.A.L. Marques
               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: gga_exc *)
(* prefix:
  gga_c_lm_params *params;

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

$define lda_c_vbh_params
$include "lda_c_hl.mpl"

(* Constant after equation (2.25) in Hu-Langreth. It is equal to the
   numerical constant in equation (12) of the Langreth-Mehl paper,
   4.28e-3, once one accounts for the factor 2 from the conversion
   from Rydberg to Hartree *)
lm_J := Pi/(16*(3*Pi^2)^(4/3)):

(* Equation (2.23) in Hu-Langreth *)
lm_d := z -> sqrt(opz_pow_n(z, 5/3) + opz_pow_n(-z, 5/3))/sqrt(2):

(* F parameter, see after eqn (2.25) in Hu-Langreth; this yields
0.26181 which is rounds to 0.262 given by Langreth and Mehl. Hu and
Langreth say the external parameter f = 0.15 for comparison with LM,
but that f = 0.17 is preferable *)
lm_F := 2*sqrt(3)*params_a_lm_f / (2*(3/Pi)^(1/6)):

(* First term in eqn (2.25) in Hu-Langreth *)
lm_t1 := (z, xs0, xs1) ->
  -7/(9*2^(5/3)) * (xs0^2*opz_pow_n(z, 4/3) + xs1^2*opz_pow_n(-z, 4/3)):

(* Second term in eqn (2.25) in Hu-Langreth *)
lm_t2 := (rs, z, xt) ->
  2/lm_d(z) * exp(-lm_F*xt*n_total(rs)^(1/6)) * xt^2:

f := (rs, z, xt, xs0, xs1) ->
  + hl_f(rs, z)
  + lm_J*(lm_t1(z, xs0, xs1) + lm_t2(rs, z, xt))*n_total(rs)^(1/3):


text:
# MG10HL von Kawasaki
# Voraussetzung:
# fourbarpris muss mit der Toolbox berechnet worden sein (Arbeitsblatt "fourbarprisTE_kinematic_constraints.mw")
# fourbar1 muss mit der Toolbox berechnet worden sein (Arbeitsblatt "fourbar1TE_kinematic_constraints.mw")
#
# Offene Punkte / TODO:
# In [Bejaoui2018_S749] wurden die Systeme fourbar1 und fourbarpris etwas anders definiert. Die Änderung wurde hier noch nicht angepasst. Daher funktioniert das Modell aktuell nicht. Die kinematischen Ketten werden nicht geschlossen.
# Autor
# Abderahman Bejaoui
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-09
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Quellen
# [Bejaoui2018_S749]
#
# Initialisierung
# Import
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
# Konstante Winkel bereits hineinschreiben
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Variablen definieren für die Hilfswinkel
#
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
# Ähnliches Vorgehen wie in [1].
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:
# Einlesen der Ergebnisse der allgemeinen Viergelenkkette. TODO: Nehme die Ausgabe von
# fourbar1TE/tmp/kinematic_constraints_maple_inert.m und nutze kintmp_subsexp direkt von dort.
read sprintf("../codeexport/fourbar1TE/tmp/fourbar1TE_angles_trig_elim"):
# Liest Variable: "winkel"
;
read sprintf("../codeexport/fourbarprisTE/tmp/fourbarprisTE_angles_trig_elim"):
# Liest Variable: "winkel_fourbarpris"
;

# Viergelenkkette A-C-D-E
winkel1:=<<sin_rho910|sin(rho910_s)>;<cos_rho910|cos(rho910_s)>;<sin_rho19|sin(rho19_s)>;<cos_rho19|cos(rho19_s)>;<sin_rho3|sin(rho3_s)>;<cos_rho3|cos(rho3_s)>>:
for i from 1 to 6 do
  winkel1(i,2):=winkel(i,2):
  cos_phi_s := sin(qJ2s)*cos(phi23)+cos(qJ2s)*sin(phi23):
  sin_phi_s := cos(qJ2s)*cos(phi23)-sin(qJ2s)*sin(phi23):
  winkel1(i,2):=subs({sin(qJ1s)=sin_phi_s},winkel1(i,2)): # qJ1s in fourbar1TE ist der aktive Winkel
  winkel1(i,2):=subs({cos(qJ1s)=cos_phi_s},winkel1(i,2)):
  winkel1(i,2) := subs({l1=AC}, winkel1(i,2)):
  winkel1(i,2) := subs({l2=AE}, winkel1(i,2)):
  winkel1(i,2) := subs({l3=ED}, winkel1(i,2)):
  winkel1(i,2) := subs({l4=DC}, winkel1(i,2)):
end do:
winkel1(1..6,1):
sin_rho910_s:=winkel1(1,2): # für kintmp_subsexp
;
cos_rho910_s:=winkel1(2,2): # für kintmp_subsexp
;
sin_rho19_s:=winkel1(3,2): # für kintmp_subsexp
;
cos_rho19_s:=winkel1(4,2):  # für kintmp_subsexp
;
sin_rho3_s:=winkel1(5,2):
cos_rho3_s:=winkel1(6,2):
cos_rho23_s:=-cos_rho3_s: # für kintmp_subsexp
;
sin_rho23_s:=sin_rho3_s:  # für kintmp_subsexp
;
# Schubkurbel G-P-K
winkel2:=<<sin_eta1112|sin(eta1112_s)>;<cos_eta1112|cos(eta1112_s)>;<sin_eta3|sin(eta3_s)>;<cos_eta3|cos(eta3_s)>;<sin_eta45|sin(eta45_s)>;<cos_eta45|cos(eta45_s)>>:
for i from 1 to 6 do
	winkel2(i,2):=winkel_fourbarpris(i,2):
	winkel2(i,2):=subs({qJ1s=qJ6s},winkel2(i,2)):
	winkel2(i,2) := subs({GK=GK}, winkel2(i,2)):
     winkel2(i,2) := subs({GP=GP}, winkel2(i,2)):
     winkel2(i,2) := subs({HP=HP}, winkel2(i,2)):
end do:
winkel2
;
sin_eta1112_s:=winkel2(1,2): # für kintmp_subsexp
;
cos_eta1112_s:=winkel2(2,2): # für kintmp_subsexp
;
sin_eta3_s:=winkel2(3,2):
cos_eta3_s:=winkel2(4,2):
sin_eta45_s:=winkel2(5,2): # für kintmp_subsexp
;
cos_eta45_s:=winkel2(6,2): # für kintmp_subsexp
;
sin_eta411_s:=sin_eta3_s: # für kintmp_subsexp
;
cos_eta411_s:=-cos_eta3_s: # für kintmp_subsexp
;
cos_eta34_s:=-(cos_eta411_s*cos(phi3)-sin_eta411_s*sin(phi3)):  # für kintmp_subsexp
;
sin_eta34_s:=sin_eta411_s*cos(phi3)+cos_eta411_s*sin(phi3):  # für kintmp_subsexp
;
# Kintmp_subsexp
kintmp_subsexp(1,2):=sin_rho23_s:
kintmp_subsexp(2,2):=cos_rho23_s:
kintmp_subsexp(3,2):=sin_eta34_s:
kintmp_subsexp(4,2):=cos_eta34_s:
kintmp_subsexp(5,2):=sin_eta45_s:
kintmp_subsexp(6,2):=cos_eta45_s:
kintmp_subsexp(7,2):=sin_rho19_s:
kintmp_subsexp(8,2):=cos_rho19_s:
kintmp_subsexp(9,2):=sin_rho910_s:
kintmp_subsexp(10,2):=cos_rho910_s:
kintmp_subsexp(11,2):=sin_eta411_s:
kintmp_subsexp(12,2):=cos_eta411_s:
kintmp_subsexp(13,2):=sin_eta1112_s:
kintmp_subsexp(14,2):=cos_eta1112_s:
kintmp_subsexp(11..14,1);

# Kintmp_qs
kintmp_qs(1,1):=arctan(sin_rho23_s,cos_rho23_s):
kintmp_qs(2,1):=arctan(sin_eta34_s,cos_eta34_s):
kintmp_qs(3,1):=arctan(sin_eta45_s,cos_eta45_s):
kintmp_qs(4,1):=arctan(sin_rho19_s,cos_rho19_s):
kintmp_qs(5,1):=arctan(sin_rho910_s,cos_rho910_s):
kintmp_qs(6,1):=arctan(sin_eta411_s,cos_eta411_s):
kintmp_qs(7,1):=arctan(sin_eta1112_s,cos_eta1112_s):
# Export
kintmp_qt := convert_s_t(kintmp_qs):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
#printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d",robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert" ,robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
#printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
#kc_symbols :=Transpose(kc_symbols);
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);
#printf("Fertig. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
kc_symbols
;

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


text:
(*
 Copyright (C) 2019 Susi Lehtola

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

(* type: mgga_exc *)
(* prefix:
  mgga_x_rscan_params *params;

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

$include "mgga_x_scan.mpl"

(* Coefficients of the rSCAN switching function from SI, in reversed(!) order: 7, 6, ..., 0 *)
rscan_fx := [-0.023185843322, 0.234528941479, -0.887998041597, 1.451297044490, -0.663086601049, -0.4445555, -0.667, 1]:

np53 := (rs, z) -> n_spin(rs,z)^(5/3):

(* First regularization: tau^u -> tau^u + tau_r.
   tau_r gets scaled by 2^(2/3) due to the spin scaling. *)
rscan_alpha0 := (rs, z, x, t) -> (np53(rs,z)*m_max(t - x^2/8, 0))/(np53(rs,z)*K_FACTOR_C + params_a_taur/2):

(* Second regularization: alpha -> alpha^3/(alpha^2 + alpha_r) *)
rscan_alpha := (rs, z, x, t) -> rscan_alpha0(rs, z, x, t)^3/(rscan_alpha0(rs, z, x ,t)^2 + params_a_alphar):

(* f(alpha) replaced with a polynomial for alpha in [0, 2.5] *)
rscan_f_alpha_small := (a, ff) -> add(ff[8-i]*a^i, i=0..7):
rscan_f_alpha_large := a -> -params_a_d*exp(params_a_c2/(1 - a)):
rscan_f_alpha := (a, ff) -> my_piecewise3(a <= 2.5, rscan_f_alpha_small(m_min(a, 2.5), ff), rscan_f_alpha_large(m_max(a, 2.5))):
rscan_f   := (rs, z, x, u, t) -> (scan_h1x(scan_y(x, rscan_alpha(rs, z, x, t)))*(1 - rscan_f_alpha(rscan_alpha(rs, z, x, t), rscan_fx))
  + scan_h0x*rscan_f_alpha(rscan_alpha(rs, z, x, t), rscan_fx))*scan_gx(x):

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


text:
# MarSystemComposite
# Type = Series
# Name = extractNet

# MarControls = 8
# /Series/extractNet/mrs_bool/debug = 0
# /Series/extractNet/mrs_bool/mute = 0
# /Series/extractNet/mrs_natural/inObservations = 1
# /Series/extractNet/mrs_natural/inSamples = 512
# /Series/extractNet/mrs_natural/onObservations = 1
# /Series/extractNet/mrs_natural/onSamples = 1
# /Series/extractNet/mrs_real/israte = 22050
# /Series/extractNet/mrs_real/osrate = 43.0664

# Number of links = 1
# Synonyms of /Series/extractNet/mrs_bool/hasData =
# Number of synonyms = 1
# /Series/extractNet/NetworkTCPSource/src/mrs_bool/hasData

# nComponents = 6

# MarSystem
# Type = NetworkTCPSource
# Name = src

# MarControls = 12
# /NetworkTCPSource/src/mrs_bool/debug = 0
# /NetworkTCPSource/src/mrs_bool/mute = 0
# /NetworkTCPSource/src/mrs_bool/hasData = 1
# /NetworkTCPSource/src/mrs_natural/controlsPort = 5010
# /NetworkTCPSource/src/mrs_natural/dataPort = 5009
# /NetworkTCPSource/src/mrs_natural/inObservations = 1
# /NetworkTCPSource/src/mrs_natural/inSamples = 512
# /NetworkTCPSource/src/mrs_natural/onObservations = 1
# /NetworkTCPSource/src/mrs_natural/onSamples = 512
# /NetworkTCPSource/src/mrs_real/gain = 1
# /NetworkTCPSource/src/mrs_real/israte = 22050
# /NetworkTCPSource/src/mrs_real/osrate = 22050

# Number of links = 0

# MarSystem
# Type = Hamming
# Name = ham

# MarControls = 8
# /Hamming/ham/mrs_bool/debug = 0
# /Hamming/ham/mrs_bool/mute = 0
# /Hamming/ham/mrs_natural/inObservations = 1
# /Hamming/ham/mrs_natural/inSamples = 512
# /Hamming/ham/mrs_natural/onObservations = 1
# /Hamming/ham/mrs_natural/onSamples = 512
# /Hamming/ham/mrs_real/israte = 22050
# /Hamming/ham/mrs_real/osrate = 22050

# Number of links = 0

# MarSystem
# Type = Spectrum
# Name = spectrum

# MarControls = 8
# /Spectrum/spectrum/mrs_bool/debug = 0
# /Spectrum/spectrum/mrs_bool/mute = 0
# /Spectrum/spectrum/mrs_natural/inObservations = 1
# /Spectrum/spectrum/mrs_natural/inSamples = 512
# /Spectrum/spectrum/mrs_natural/onObservations = 512
# /Spectrum/spectrum/mrs_natural/onSamples = 1
# /Spectrum/spectrum/mrs_real/israte = 22050
# /Spectrum/spectrum/mrs_real/osrate = 43.0664

# Number of links = 0

# MarSystem
# Type = PowerSpectrum
# Name = psk

# MarControls = 9
# /PowerSpectrum/psk/mrs_bool/debug = 0
# /PowerSpectrum/psk/mrs_bool/mute = 0
# /PowerSpectrum/psk/mrs_natural/inObservations = 512
# /PowerSpectrum/psk/mrs_natural/inSamples = 1
# /PowerSpectrum/psk/mrs_natural/onObservations = 256
# /PowerSpectrum/psk/mrs_natural/onSamples = 1
# /PowerSpectrum/psk/mrs_real/israte = 43.0664
# /PowerSpectrum/psk/mrs_real/osrate = 43.0664
# /PowerSpectrum/psk/mrs_string/spectrumType = power

# Number of links = 0

# MarSystem
# Type = Centroid
# Name = cntrd

# MarControls = 8
# /Centroid/cntrd/mrs_bool/debug = 0
# /Centroid/cntrd/mrs_bool/mute = 0
# /Centroid/cntrd/mrs_natural/inObservations = 256
# /Centroid/cntrd/mrs_natural/inSamples = 1
# /Centroid/cntrd/mrs_natural/onObservations = 1
# /Centroid/cntrd/mrs_natural/onSamples = 1
# /Centroid/cntrd/mrs_real/israte = 43.0664
# /Centroid/cntrd/mrs_real/osrate = 43.0664

# Number of links = 0

# MarSystem
# Type = NetworkTCPSink
# Name = dest

# MarControls = 12
# /NetworkTCPSink/dest/mrs_bool/debug = 0
# /NetworkTCPSink/dest/mrs_bool/mute = 0
# /NetworkTCPSink/dest/mrs_natural/controlsPort = 5010
# /NetworkTCPSink/dest/mrs_natural/dataPort = 5009
# /NetworkTCPSink/dest/mrs_natural/inObservations = 1
# /NetworkTCPSink/dest/mrs_natural/inSamples = 1
# /NetworkTCPSink/dest/mrs_natural/onObservations = 1
# /NetworkTCPSink/dest/mrs_natural/onSamples = 1
# /NetworkTCPSink/dest/mrs_real/gain = 1
# /NetworkTCPSink/dest/mrs_real/israte = 43.0664
# /NetworkTCPSink/dest/mrs_real/osrate = 43.0664
# /NetworkTCPSink/dest/mrs_string/host = cyprus.cs.uvic.ca

# Number of links = 0

max_stars_repo_name: jaouahbi/VampPlugins

max_issues_repo_name: jaouahbi/VampPlugins

max_forks_repo_name: jaouahbi/VampPlugins


text:
var $test_var i32 public

func $test_extern_var ( var %i i32 ) i32 {
  return (dread i32 $test_var)
}

var $test_var i32 public = 0
 # 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:
repeat process ACV_lasers {
/tp @e[type=armor_stand,name=ACV_Cube] ~ ~ ~ ~ 0

/scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser] remove ACV_InCube
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser] remove ACV_EnteringCube
/scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser] remove ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ particle reddust ~ ~ ~ 0.1 0.1 0.1 0 10 force
/execute @e[type=armor_stand,name=ACV_Cube,tag=ACV_Redirecting] ~ ~2 ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,dy=0] add ACV_InCube
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=169,ry=-168,rxm=0,rx=0] ~0 ~ ~-0.5
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-167,ry=-146,rxm=0,rx=0] ~0.1913417161825448858642299920152 ~ ~-0.46193976625564337806409159469839
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-145,ry=-123,rxm=0,rx=0] ~0.35355339059327376220042218105242 ~ ~-0.35355339059327376220042218105242
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-122,ry=-101,rxm=0,rx=0] ~0.46193976625564337806409159469839 ~ ~-0.1913417161825448858642299920152
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-100,ry=-78,rxm=0,rx=0] ~0.5 ~ ~0
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-77,ry=-56,rxm=0,rx=0] ~0.46193976625564337806409159469839 ~ ~0.1913417161825448858642299920152
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-55,ry=-33,rxm=0,rx=0] ~0.35355339059327376220042218105242 ~ ~0.35355339059327376220042218105242
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-32,ry=-11,rxm=0,rx=0] ~0.1913417161825448858642299920152 ~ ~0.46193976625564337806409159469839
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=-10,ry=11,rxm=0,rx=0] ~0 ~ ~0.5
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=12,ry=33,rxm=0,rx=0] ~-0.1913417161825448858642299920152 ~ ~0.46193976625564337806409159469839
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=34,ry=56,rxm=0,rx=0] ~-0.35355339059327376220042218105242 ~ ~0.35355339059327376220042218105242
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=57,ry=78,rxm=0,rx=0] ~-0.46193976625564337806409159469839 ~ ~0.1913417161825448858642299920152
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=79,ry=101,rxm=0,rx=0] ~-0.5 ~ ~0
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=102,ry=123,rxm=0,rx=0] ~-0.46193976625564337806409159469839 ~ ~-0.1913417161825448858642299920152
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=124,ry=146,rxm=0,rx=0] ~-0.35355339059327376220042218105242 ~ ~-0.35355339059327376220042218105242
/tp @e[type=area_effect_cloud,name=ACV_Laser,rym=147,ry=168,rxm=0,rx=0] ~-0.1913417161825448858642299920152 ~ ~-0.46193976625564337806409159469839
/tp @e[type=area_effect_cloud,name=ACV_Laser,rxm=1,rx=90] ~ ~-0.5 ~
/tp @e[type=area_effect_cloud,name=ACV_Laser,rxm=-90,rx=-1] ~ ~0.5 ~
/execute @e[type=armor_stand,name=ACV_Cube,tag=!ACV_Redirecting] ~ ~1.5 ~ kill @e[type=area_effect_cloud,name=ACV_Laser,dy=0]
/execute @e[type=armor_stand,name=ACV_Cube,tag=ACV_Redirecting] ~ ~1.5 ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,tag=!ACV_InCube,dy=0] add ACV_EnteringCube
/execute @e[type=armor_stand,name=ACV_Cube,tag=ACV_Redirecting] ~ ~1.5 ~ tp @e[type=area_effect_cloud,name=ACV_Laser,tag=ACV_EnteringCube,dy=0] @e[type=armor_stand,name=ACV_Cube,r=2,c=1]
/execute @e[type=armor_stand,name=ACV_Cube,tag=ACV_Redirecting] ~ ~ ~ tp @e[type=area_effect_cloud,name=ACV_Laser,tag=ACV_EnteringCube,dy=0] ~ ~1.7 ~

/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ air -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ piston_extension -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ carpet -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ glass -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ stained_glass -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ glass_pane -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ stained_glass_pane -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir
/execute @e[type=area_effect_cloud,name=ACV_Laser] ~ ~ ~ detect ~ ~ ~ iron_bars -1 scoreboard players tag @e[type=area_effect_cloud,name=ACV_Laser,r=1,c=1] add ACV_InAir

/kill @e[type=area_effect_cloud,name=ACV_Laser,tag=!ACV_InAir]
}


text:
restart;
kernelopts(printbytes=false):
partitions := proc(c::list(integer),subtotal::integer) :: integer:
  local i,m,total,cards;
  cards := c;
  m := 0;
  # Hit
  for i from 1 to 10 do
    if (cards[i]>0) then
      total := subtotal+i;
      if (total < 21) then
        # Stand
        m := m+1;
	# Hit again
        cards[i] := cards[i]-1;
        m := m+partitions(cards, total);
        cards[i] := cards[i]+1;
      elif (total=21) then
        # Stand; hit again is an automatic bust
        m := m+1;
      end if;
    end if;
  end do;
  return m;
end proc:

deck := [4,4,4,4,4,4,4,4,4,16]:
d := 0:

for i from 1 to 10 do
  # Dealer showing
  p := 0:
  deck[i] := deck[i]-1:
  for j from 1 to 10 do
    deck[j] := deck[j]-1:
    p := p+partitions(deck,j):
    deck[j] := deck[j]+1:
  end do:
  deck[i] := deck[i]+1:
  printf("Dealer showing %d partitions = %d\n",i-1,p):
  d := d+p:
end do:
printf("Total partitions = %d\n",d):

max_stars_repo_name: octonion/recursion

max_issues_repo_name: octonion/recursion

max_forks_repo_name: octonion/recursion


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: mgga_exc *)
$include "mgga_xc_b98.mpl"

(* Fitted parameters, given below Eq. 50 *)
edmgga_mu := 0.21:
edmgga_a  := 0.704:

(* Eq. 40: F_x'(0) *)
edmgga_fp0   := 27*edmgga_mu/50:
(* Eq. 27: lim_(x -> -infinity) F_x(x) *)
edmgga_fminf := 1/3 * (4*Pi^2/3)^(1/3):
(* Eq. 44 *)
edmgga_b := 2*Pi/9 * sqrt(6/5):

(* Eq. 47 *)
edmgga_c1 := edmgga_fminf:
(* Eq. 48 *)
edmgga_c2 := 1 - edmgga_c1:
(* Eq. 49 *)
edmgga_c3 := edmgga_c2*sqrt(2*edmgga_a)/(2*edmgga_b):

(* Eq. 50. There's some problem in the equation, since the first line,
1/edmgga_c3^3 * (edmgga_c3^2 - edmgga_fp0/(2*edmgga_b^2)) * edmgga_c2,
does not agree with the second line. The first line evaluates to a
negative value, which makes the functional blow up. The second line
yields values that are close to the reported ones, but Tao might have
used a more accurate numerical value in his implementation since the
agreement is not perfect. Unfortunately there's no way to check this
as he has passed away. *)
edmgga_c4 := 1/edmgga_c3^3 * (edmgga_c3^2 - 0.09834*edmgga_mu):

(* Eq. 46 *)
(* Apparently, edmgga_x can sometimes be smaller than 0,
   as Qb can be very negative. The expansion for -infinity
   has an error of 1e-13 *)
edmgga_x_large := y -> -1/(2*y) + 1/(8*y^3):

edmgga_x := Qb -> my_piecewise3(
  Qb < -1e4/edmgga_a,
  edmgga_x_large(m_min(edmgga_a*Qb, -1e-4/edmgga_a)),
  edmgga_a*Qb + sqrt(1 + (edmgga_a*Qb)^2)
):

(* Eq. 45 *)
edmgga_f_x :=  x -> edmgga_c1 + edmgga_c2*x/(1 + edmgga_c3*sqrt(x)*arcsinh(edmgga_c4*(x - 1))):

(* the evalf is absolutely necessary to avoid exceptions *)
edmgga_f := (x, u, t) -> edmgga_f_x(edmgga_x(b98_q(x, u, t))):

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


text:
`is_element/prime_W` := (N::posint) -> (A::set) -> proc(x)
 `is_element/vector_functions`(N)(A)(x) and
 `sum/vector_functions`(N)(A)(x) = [0$N];
end;

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

`is_equal/prime_W` := (N::posint) -> (A::set) -> proc(x,y)
 return `is_equal/vector_functions`(N)(A)(x,y);
end;

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

`is_leq/prime_W` := NULL;

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

`is_zero/prime_W` := (N::posint) -> (A::set) -> proc(x)
 return `is_zero/vector_functions`(N)(A)(x);
end;

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

`norm/prime_W` := (N::posint) -> (A::set) -> proc(x)
 return `norm/vector_functions`(N)(A)(x);
end;

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

`normalise/prime_W` := (N::posint) -> (A::set) -> proc(x)
 local r,a;
 r := `norm/vector_functions`(N)(A)(x);
 if r = 0 then return FAIL; fi;
 return table([seq(a = simplify(x[a] /~ r),a in A)]);
end;

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

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

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

 return eval(x);
end;

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

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

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

 y := table();
 for c in C do y[c] := x[c]; od;
 z := `normalise/W`(N)(C)(y);
 return eval(z);
end;


text:
#
#     Test Suite for Reparam---integration and summation changes of variables
#
#Created by Carl 17Jun2016.

#common header for Hakaru-Maple test suites
kernelopts(assertlevel= 2): # be strict on all assertions while testing
kernelopts(opaquemodules= false): # allow testing of internal routines
if not NewSLO :: `module` then
  WARNING("loading NewSLO failed");
  `quit`(3)
end if:

with(Hakaru):
with(KB):
with(NewSLO):
########################

#common procedures for these tests.

TestReparam:= proc(
     e_in,
     e_out,   #Expected output, after postprocessing by `verify` and `&under`
     #`&under` is defined below.
     #Type `verification` is slightly modified (below) from its standard definition.
     ver::verification,
     #First infolevel is for initial test; second is for retest after failure.
     {infolevels::[nonnegative, nonnegative]:= [0, Float(infinity)]},
     {ctx::list:= []}   #Assumptions (not necessarily to be passed to `assuming`)
)
# For 'production' runs, we want things to be as quiet
# as possible.  We would like to do
# maple -q < ReparamT.mpl | grep -v "passed"
# and get 0 lines of output. If global `pretty`	is true, some extra lines are
# printed for readability. You may want to set this in a worksheet session.
local
     LOform:= toLO(e_in),  #preprocess input
     r
;
     if :-pretty then printf("\n") end if;
     :-infolevel[reparam]:= infolevels[1];
     if not CodeTools:-Test(reparam(LOform, :-ctx= ctx), e_out, ver, boolout, _rest) then
          #If test fails, do same test with diagnostic infolevel setting.
          if infolevels[2] > infolevels[1] then
               printf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"  );
               printf("\n*** Rerun at higher infolevel: ***\n");
               :-infolevel[reparam]:= infolevels[2];
               r:= CodeTools:-Test(reparam(LOform, :-ctx= ctx), e_out, ver, boolout, _rest);
               printf(  "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
               return r
          else
               return false
          end if
     else
          return
     end if
end proc:

#VerifyTools is an undocumented Maple library package similiar to TypeTools.
#Add a verification &under analogous to type &under.
VerifyTools:-AddVerification(`&under`= ((e1,e2,Ver,f)-> verify(f(e1,_rest), e2, Ver))):

#Need to overload Maple library `type/verification` so that expressions of the form 'v &under f'
#are considered verifications.
`type/verification/MapleLibrary`:= eval(`type/verification`):
`type/verification`:= overload([
     proc(Under::specfunc(`&under`))
     option overload;
          type(op(1,Under), `verification/MapleLibrary`)
     end proc,

     `type/verification/MapleLibrary`
]):

#Automatic simplifications: When an easy transformation can be used to
#transform a distribution that we don't support to one that we do, then a
#procedure in this section does it. These may be moved to NewSLO.

#Bernoulli
Bern:= proc(p)
     local i;
     Bind(Counting(0,1), i, Weight(idx([p,1-p],i), Ret(i)))
end proc:

#Exponential
Exponential:= mu-> GammaD(1,mu):

#Control printing of extra lines for readability.
if not assigned(pretty) then pretty:= false end if:

#
#   The tests
#

#Constant multiple in integral
TestReparam(
     Bind(Lebesgue(-infinity,infinity), x, Weight(2, Ret(2*x))),
     Lebesgue(-infinity,infinity),
     equal &under fromLO,
     label= "(t0a) Constant multiple with Lebesgue" #passing
):

#Logarithm with Uniform
TestReparam(
     Bind(Uniform(0,1), x, Ret(-log(x))),
     GammaD(1,1),
     equal &under fromLO,
     label= "(t0b) Logarithm with Uniform" #passing
):

#(t1) Symbolic affine transformation with Gaussian
TestReparam(
     Bind(Gaussian(mu,sigma), x, Ret((x-mu)/sigma)),
     Gaussian(0,1),
     equal &under fromLO,
     ctx= [sigma > 0],
     label= "(t1) Symbolic affine transformation with Gaussian" #passing
):

#(t2) ChiSquared to Gamma
#(inverse of relationship given in Relationships.hs)
# Special case of Gamma now recognized as being ChiSquared
TestReparam(
     ChiSquared(2*b),
     GammaD(b,2),
     equal &under (fromLO, _ctx= foldr(assert, empty, b > 0)),
     ctx= [b > 0],
     label= "(t2) ChiSquared to Gamma" #passing
):

#(t3) Symbolic constant multiple with Gamma
TestReparam(
     Bind(GammaD(alpha,beta), x, Ret(2*x/beta)),
     GammaD(alpha, 2),
     equal &under (fromLO, _ctx= foldr(assert, empty, alpha > 0, beta > 0)),
     ctx= [alpha > 0, beta > 0],
     label= "(t3) Symbolic constant multiple with Gamma" #passing
):

#(t4) Two-variable LFT with Gamma
#This fails due to multivariate changes of variable not being implemented yet.
TestReparam(
     Bind(GammaD(a,t), x1, Bind(GammaD(b,t), x2, Ret(x1/(x1+x2)))),
     BetaD(a,b),
     equal &under fromLO,
     infolevels= [2,2],
     label= "(t4) Two-variable LFT with Gamma (currently failing)"
):

#(t5) Logarithm with symbolic constant multiplier and Uniform
TestReparam(
     Bind(Uniform(0,1), x, Ret(-alpha*log(x))),
     GammaD(1,alpha),
     equal &under fromLO,
     ctx= [alpha >= 0],
     label= "(t5) Logarithm with symbolic constant multiplier and Uniform"
     #passing
):

#(t6) Poisson(mu) -> ? as mu -> infinity
#The concept of taking the limit at infinity of a distribution parameter isn't
#implemented yet.

#(t7) Quotient of standard normals to standard Cauchy
#Fails because it's multivariate.
#This should be the first multivariate integration to implement, due to
#simplicity; it might work as a simple multivariate substitution.
TestReparam(
     Bind(Gaussian(0,1), x1, Bind(Gaussian(0,1), x2, Ret(x1/x2))),
     Cauchy(0,1),
     equal &under fromLO,
     infolevels= [2,2],
     label=
          "(t7) Quotient of standard normals to standard Cauchy "
          "(currently failing)"
):

#(t8) Affine translation of quotient of standard normals to Cauchy
#This is a trivial variation of t7.
#Fails because it's multivariate.
TestReparam(
     Bind(Gaussian(0,1), x1, Bind(Gaussian(0,1), x2, Ret(a+alpha*x1/x2))),
     Cauchy(a,alpha),
     equal &under (fromLO, _ctx= foldr(assert, empty, alpha > 0)),
     ctx= [alpha > 0],
     infolevels= [2,2],
     label=
          "(t8) Affine translation of quotient of standard normals to Cauchy "
          "(currently failing)"
):

#(t9) Bernoulli to Binomial with n=1
#Fails because Binomial isn't implemented.
#Doesn't use reparam in a meaningful way.
TestReparam(
     Bern(p),
     Binomial(1,p),
     equal &under (fromLO, _ctx= foldr(assert, empty, 0 <= p, p <= 1)),
     ctx= [0 <= p, p <= 1],
     infolevels= [2,2],
     label= "(t9) Bernoulli to Binomial with n=1 (currently failing)"
):

#(t10) Beta(1,1) to Uniform(0,1)
#This one doesn't require reparam, but it passes using reparam.
TestReparam(
     #Note: No use of `Bind`. This passes by direct recognition by fromLO
     #without any meaningful use of reparam, although reparam is called.
     BetaD(1,1),
     Uniform(0,1),
     equal &under fromLO,
     label= "(t10) Beta(1,1) to Uniform(0,1)" #passing
):

#(t11) Laplace to difference of Gammas
#This one fails because Laplace isn't implemented. Also, there's nothing that
#reparam can do with this. It's up to fromLO to recognize the compound
#expression.
#Perhaps this should be an automatic simplification.
TestReparam(
     Laplace(a1,a2),
     Bind(GammaD(1,a1), x1, Bind(GammaD(1,a2), x2, Ret(x1-x2))),
     equal &under (fromLO, _ctx= foldr(assert, empty, a1 > 0, a2 > 0)),
     ctx= [a1 > 0, a2 > 0],
     infolevels= [2,2],
     label= "(t11) Laplace to difference of Gammas (currently failing)"
):

#(t12) Sum of iid Exponentials to Gamma
#This currently fails because it's multivariate.
TestReparam(
     Bind(Exponential(b), x1, Bind(Exponential(b), x2, Ret(x1+x2))),
     GammaD(2,b),
     equal &under (fromLO, _ctx= foldr(assert, empty, b > 0)),
     ctx= [b > 0],
     infolevels= [2,2],
     label= "(t12) Sum of iid Exponentials to Gamma (currently failing)"
):

#(t13) Weibull(1,b) to Exponential(1/b)
#This test is wrong in Relationships.hs!
#It should be Weibull(a,1) to Exponential(1/a).
#Fails because Weibull isn't implemented.
TestReparam(
     #Straightforward recognition; doesn't use reparam in a meaningful way.
     Weibull(a,1),
     Exponential(1/a),
     equal &under (fromLO, _ctx= foldr(assert, empty, a > 0)),
     ctx= [a > 0],
     infolevels= [2,2],
     label= "(t13) Weibull(a,1) to Exponential(1/a) (currently failing)"
):

#(t14) Standard normal over sqrt(ChiSquared) to StudentT
#Fails because it's multivariate.
TestReparam(
     Bind(Gaussian(0,1), x, Bind(ChiSquared(nu), y, Ret(x/sqrt(y/nu)))),
     StudentT(nu, 0, 1),
     equal &under (fromLO, _ctx= foldr(assert, empty, nu > 0)),
     ctx= [nu > 0],
     infolevels= [2,2],
     label= "(t14) Std normal over sqrt(Chi2) to StudentT (currently failing)"
):

#(t15) 1/InverseGamma(k,1/t) to GammaD(k,t)
TestReparam(
     Bind(InverseGammaD(k, 1/t), x, Ret(1/x)),
     GammaD(k,t),
     equal &under (fromLO, _ctx= foldr(assert, empty, k > 0, t > 0)),
     ctx= [k > 0, t > 0],
     label= "(t15) 1/InverseGammaD(k,1/t) to GammaD(k,t)" #passing
):

#(t16) Sum of std normals to normal(0,sqrt(2))
#Fails because it requires a multivariate change of vars.
TestReparam(
     Bind(Gaussian(0,1), x1, Bind(Gaussian(0,1), x2, Ret(x1+x2))),
     Gaussian(0,sqrt(2)),
     equal &under fromLO,
     infolevels= [2,2],
     label= "(t16) Sum of std normals to normal(0,sqrt(2)) (currently failing)"
):

#(t17) Sum of std normal and normal to normal
#Fails because it requires a multivariate change of vars.
TestReparam(
     Bind(Gaussian(0,1), x1, Bind(Gaussian(mu,sigma), x2, Ret(x1+x2))),
     Gaussian(mu, sqrt(1+sigma^2)),
     equal &under (fromLO, _ctx= foldr(assert, empty, sigma > 0)),
     ctx= [sigma > 0],
     infolevels= [2,2],
     label= "(t17) Sum of std normal and normal to normal (currently failing)"
):

#(t18) Sum of symbolic multiples of std normals to normal
#Fails because it requires a multivariate change of vars.
TestReparam(
     Bind(Gaussian(0,1), x1, Bind(Gaussian(0,1), x2, Ret(a1*x1 + a2*x2))),
     Gaussian(0, sqrt(a1^2 + a2^2)),
     equal &under fromLO,
     infolevels= [2,2],
     label=
          "(t18) Sum of symbolic multiples of std normals to normal "
          "(currently failing)"
):

#(t19) Sum of equiprobable binomials to binomial
#Fails for three reasons: it's multivariate, it's discrete, and Binomial isn't
#implemented.
TestReparam(
     Bind(Binomial(n1,p), x1, Bind(Binomial(n2,p), x2, Ret(x1+x2))),
     Binomial(n1+n2, p),
     equal &under (fromLO,
          _ctx= foldr(
               assert, empty,
               n1::integer, n1 > 0, n2::integer, n2 > 0, 0 <= p, p <= 1
          )
     ),
     ctx= [n1::integer, n1 > 0, n2::integer, n2 > 0, 0 <= p, p <= 1],
     infolevels= [2,2],
     label=
          "(t19) Sum of equiprobable binomials to binomial "
          "(currently failing)"
):

#(t20) Sum of n Bernoullis to Binomial
#Fails because Binomial isn't implemented.
TestReparam(
     Bind(Plate(n, i, Bern(p)), xs, Ret(sum(idx(xs,i),i=0..n-1))),
     Binomial(n,p),
     equal &under (fromLO,
          _ctx= foldr(assert, empty, n::integer, n > 0, 0 <= p, p <= 1)
     ),
     ctx= [n::integer, n > 0, 0 <= p, p <= 1],
     infolevels= [2,2],
     label= "(t20) Sum of n Bernoullis to Binomial (currently failing)"
):

#(t21) Sum of Poissons to Poisson
#Fails because reparam doesn't yet handle summations.
TestReparam(
     Bind(PoissonD(lambda1), x1, Bind(PoissonD(lambda2), x2, Ret(x1+x2))),
     PoissonD(lambda1 + lambda2),
     equal &under (fromLO,
          _ctx= foldr(assert, empty, lambda1 > 0, lambda2 > 0)
     ),
     ctx= [lambda1 > 0, lambda2 > 0],
     infolevels= [2,2],
     label= "(t21) Sum of Poissons to Poisson (currently failing)"
):

#(t22) Sum of Gammas to Gamma
#Fails because it's multivariate.
TestReparam(
     Bind(GammaD(a1,b), x1, Bind(GammaD(a2,b), x2, Ret(x1+x2))),
     GammaD(a1+a2, b),
     equal &under (fromLO, _ctx= foldr(assert, empty, a1 > 0, a2 > 0, b > 0)),
     ctx= [a1 > 0, a2 > 0, a3 > 0],
     infolevels= [2,2],
     label= "(t22) Sum of Gammas to Gamma (currently failing)"
):

#(t23) Sum of n Exponentials to Gamma
#Fails because it involves a Plate.
TestReparam(
     Bind(Plate(n, i, Exponential(b)), xs, Ret(sum(idx(xs,i), i= 0..n-1))),
     GammaD(n,b),
     equal &under (fromLO, _ctx= foldr(assert, empty, n::integer, n > 0, b > 0)),
     ctx= [n::integer, n > 0, b > 0],
     infolevels= [2,2],
     label= "(t23) Sum of n Exponentials to Gamma (currently failing)"
):

#(t24) Weibull "scaling" property.
#I can find no evidence that this is true; indeed, it seems trivial to disprove
#it.

#(t25) Product of lognormals is lognormal
#Fails because it's multivariate.
TestReparam(
     Bind(
          Gaussian(mu1,sigma1), x1,
          Bind(Gaussian(mu2,sigma2), x2, Ret(exp(x1)*exp(x2)))
     ),
     Bind(Gaussian(mu1+mu2, sqrt(sigma1^2 + sigma2^2)), x3, Ret(exp(x3))),
     equal &under (fromLO, _ctx= foldr(assert, empty, sigma1 > 0, sigma2 > 0)),
     ctx= [sigma1 > 0, sigma2 > 0],
     infolevels= [2,2],
     label= "(t25) Product of lognormals is lognormal (currently failing)"
):

#(t26) Cauchy of reciprocal to Cauchy
#Fails due to improper handling of the discontinuity in the substitution.
TestReparam(
     Bind(Cauchy(0, sigma), x, Ret(1/x)),
     Cauchy(0, 1/sigma),
     equal &under (fromLO, _ctx= foldr(assert, empty, sigma > 0)),
     ctx= [sigma > 0],
     infolevels= [2,2],
     label= "(t26) Cauchy of reciprocal to Cauchy (currently failing)"
):

#(t27) Symbolic constant multiple of Gamma to Gamma
TestReparam(
     Bind(GammaD(r,lambda), x, Ret(a*x)),
     GammaD(r, a*lambda),
     equal &under (fromLO,
          _ctx= foldr(assert, empty, r > 0, lambda > 0, a > 0)
     ),
     ctx= [r > 0, lambda > 0, a > 0],
     label= "(t27) Symbolic constant multiple of Gamma to Gamma" #passing
):

#(t28) Beta of 1-x to Beta
TestReparam(
     Bind(BetaD(a,b), x, Ret(1-x)),
     BetaD(b,a),
     equal &under (fromLO, _ctx= foldr(assert, empty, a > 0, b > 0)),
     ctx= [a > 0, b > 0],
     label= "(t28) Beta with 1-x to Beta" #passing
):

#(t29) Binomial of n-x to Binomial
#Fails because Binomial isn't implemented.
TestReparam(
     Bind(Binomial(n,p), x, Ret(n-x)),
     Binomial(n, 1-p),
     equal &under (fromLO,
          _ctx= foldr(assert, empty, n::integer, n > 0, 0 <= p, p <= 1)
     ),
     ctx= [n::integer, n > 0, 0 <= p, p <= 1],
     infolevels= [2,2],
     label= "(t29) Binomial of n-x to Binomial (currently failing)"
):

#(t30) Square of std normal to ChiSquared(1)
#Fails due to non-invertibility of substitution u= x^2.
TestReparam(
     Bind(Gaussian(0,1), x, Ret(x^2)),
     ChiSquared(1),
     equal &under fromLO,
     infolevels= [2,2],
     label= "(t30) Square of std normal to ChiSquared(1) (currently failing)"
):

#(t31) Sum of squares of n iid std normals to ChiSquared(n)
#Fails because it involves a Plate.
TestReparam(
     Bind(Plate(n, i, Gaussian(0,1)), xs, Ret(sum(idx(xs,i)^2, i= 0..n-1))),
     ChiSquared(n),
     equal &under (fromLO, _ctx= foldr(assert, empty, n::integer, n > 0)),
     ctx= [n::integer, n > 0],
     infolevels= [2,2],
     label=
          "(t31) Sum of squares of n iid std normals to ChiSquared(n) "
          "(currently failing)"
):

#(t32) Sum of ChiSquareds is ChiSquared
#Fails because it's multivariate.
TestReparam(
     Bind(ChiSquared(k1), x1, Bind(ChiSquared(k2), x2, Ret(x1+x2))),
     ChiSquared(k1+k2),
     equal &under (fromLO, _ctx= foldr(assert, empty, k1 > 0, k2 > 0)),
     ctx= [k1 > 0, k2 > 0],
     infolevels= [2,2],
     label= "(t32) Sum of ChiSquareds is ChiSquared (currently failing)"
):

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

a0 := 0.93222*RS_FACTOR:
kk := 9.47362e-3*RS_FACTOR:

f := (rs, zeta) -> -a0*(1 - kk*log(1 + rs/kk)/rs)/rs:


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

`is_element/P/itloc` := (N::set(nonnegint)) -> (AA) -> (T) -> proc(TT)
 `is_element/M1/itloc`(N)(AA)(TT) and `omega/itloc`(TT) = T;
end:

`is_leq/P/itloc` := (N::set(nonnegint)) -> (AA) -> (T) -> proc(TT0,TT1)
 `and`(seq(TT0[i] minus TT1[i] = {},i=1..nops(TT0)));
end:

`list_elements/P/itloc` := (N::set(nonnegint)) -> (AA) -> proc(T)
 local BB,L;

 BB := `cap/itloc`(AA,T);
 L := `list_elements/M1/itloc`(T)(BB);
 L := select(TT -> `omega/itloc`(TT) = T,L);
 return L;
end:

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

`is_element/P_star/itloc` := (N) -> (AA) -> (T) -> (iu) -> proc(TT)
 local i,u;

 i,u := op(iu);

 if not(`is_element/P/itloc`(N)(AA)(T)(TT)) then
  return false;
 fi;

 if max(u,op(TT[i])) <> u then
  return false;
 fi;

 return true;
end:

`Alpha/itloc` := (N) -> (AA) -> (T) -> (iu) -> proc(TT)
 local i,u,r,UU;

 i,u := op(iu);

 if max(op(TT[i])) < u then
  return TT;
 fi;

 r := nops(AA);
 UU := [
  seq(TT[k],k=1..i),
  TT[i+1] union {u},
  seq(TT[k],k=i+2..r)
 ];

 return UU;
end:

`Beta/itloc` := (N) -> (AA) -> (T) -> (iu) -> proc(TT)
 local i,u,r,UU;

 i,u := op(iu);

 if max(op(TT[i])) < u then
  return TT;
 fi;

 r := nops(AA);
 UU := [
  seq(TT[k],k=1..i-1),
  TT[i] minus {u},
  TT[i+1] union {u},
  seq(TT[k],k=i+2..r)
 ];

 return UU;
end:

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

`is_element/MM/itloc` := (N) -> (AA) -> (T) -> (i) -> proc(TT)
 local a,r,k,UU;

 if not(`is_element/P/itloc`(N)(AA)(T)(TT)) then
  return false;
 fi;

 a := `bot/threads/itloc`(`cap/itloc`(AA,T));
 r := nops(AA);
 UU := [
  select(t -> t <= a[1],T),
  seq(select(t -> a[k-1] <= t and t <= a[k],T),k=2..r)
 ];

 for k from 1 to i do
  if TT[k] <> UU[k] then return false; fi;
 od;

 return true;
end:

`is_element/NN/itloc` := (N) -> (AA) -> (T) -> (i) -> proc(TT)
 local a,r,k,UU;

 if not(`is_element/P/itloc`(N)(AA)(T)(TT)) then
  return false;
 fi;

 a := `bot/threads/itloc`(`cap/itloc`(AA,T));
 r := nops(AA);
 UU := [
  select(t -> t <= a[1],T),
  seq(select(t -> a[k-1] <= t and t <= a[k],T),k=2..r)
 ];

 for k from 1 to i do
  if TT[k] <> UU[k] then return false; fi;
 od;

 if i = 0 and not(member(min(op(T)),TT[1])) then
  return false;
 fi;

 if i > 0 and not(member(a[i],TT[i+1])) then
  return false;
 fi;

 return true;
end:

`is_element/LL/itloc` := (N) -> (AA) -> (T) -> (iu) -> proc(TT)
 local i,u,a,r,k,UU;

 i,u := op(iu);

 if not(`is_element/NN/itloc`(N)(AA)(T)(i-1)(TT)) then
  return false;
 fi;

 if max(op(TT[i])) > u then
  return false;
 fi;

 return true;
end:

`theta/itloc` := (N) -> (AA) -> (T) -> (i) -> proc(TT)
 local a,r;

 a := `bot/threads/itloc`(`cap/itloc`(AA,T));
 r := nops(AA);

 return [
  seq(TT[k],k=1..i),
  TT[i+1] union {a[i]},
  seq(TT[k],k=i+2..r)
 ];
end;