text:
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(NewSLO):
with(Partition):

# covers primitive constructs
model1 :=
  Bind(Gaussian(0,1), x,
  Bind(Msum(Ret(0),Weight(1,Lebesgue(-infinity,infinity))), y,
  Ret(1/exp(x^2+y^2)))):

# simplifies to
model1s :=
  Bind(Gaussian(0,1), x,
  Msum(Ret(exp(-x^2)),
       Bind(Lebesgue(-infinity,infinity), y,
       Ret(exp(-(x^2+y^2)))))):

CodeTools[Test](value(integrate(model1,z->z)), (sqrt(Pi)+1)/sqrt(3), equal,
  label="primitive constructs + integrate + value");

TestHakaru(model1, model1s, label = "primitive constructs simplification");

# Unknown measures -- no changes
u1 := Bind(m, x, Ret(x^2)):
u2 := Bind(Gaussian(0,1), x, m(x)):

TestHakaru(u1, label = "binding unknown m");
TestHakaru(u2, u2, label = "sending to unknown m");

# hygiene for Bind
TestHakaru(Bind(Bind(Uniform(0,1),x,Ret(x^2)),y,Ret(x^3)),
           Ret(x^3),
           label = "hygiene for Bind");

# example with an elaborate simplifier to do reordering of
# integration, which in turn integrates out a variable
model3 := Bind(Gaussian(0,1),x,Gaussian(x,1)):

TestHakaru(model3, Gaussian(0,sqrt(2)),
 simp = (lo -> subsindets(lo, 'specfunc(Int)',
  i -> subsindets(IntegrationTools[CollapseNested](
                  IntegrationTools[Combine](i)), 'Int(anything,list)',
  i -> subsop(0=int, applyop(ListTools[Reverse], 2, i))))),
  label = "use simplifier to integrate out variable");

TestHakaru(Bind(GammaD(1,1),lambda,PoissonD(lambda)), NegativeBinomial(1,1/2),
           label = "integrate out GammaD with PoissonD likelihood to Geometric");
TestHakaru(Bind(GammaD(r,1/(1/p-1)),lambda,PoissonD(lambda)), NegativeBinomial(r,p),
           ctx = [0<p,p<1],
           label = "integrate out GammaD with PoissonD likelihood to NegativeBinomial");

# Kalman filter; note the parameter + assumption
module()
  local y, kalman;
  kalman :=
    Bind(Gaussian(0,1),x,Weight(NewSLO:-density[Gaussian](x,1)(y),Ret(x))):
  TestHakaru(kalman,
    Weight(exp(-y^2/4)/2/sqrt(Pi),Gaussian(y/2,1/sqrt(2))),
    label = "Kalman filter", ctx = [y::real]);
end module:

# piecewise
model4 :=
  Bind(Gaussian(0,1),x,
  Bind(piecewise(x<0,Ret(0),x>4,Ret(4),Ret(x)),y,
  Ret(y^2))):
model4s := { Bind(Gaussian(0,1),x,piecewise(x<0,Ret(0),4<x,Ret(16),
  0<x and x<4, Ret(x^2))) ,
  Bind(Gaussian(0,1),`x`,piecewise(`x` < 0,Ret(0),4 < `x`,Ret(16),0 <= `x` and `x` <= 4,Ret(`x`^2)))
 }:

TestHakaru(model4, model4s, label = "piecewise test");
sliced :=
  Weight(exp(-lu),
  Bind(Gaussian(0,1), x,
  Weight(1/exp(-x^2/2),
  piecewise(exp(-lu) < exp(-x^2/2), piecewise(0<exp(-lu), Ret(x)))))):
sliceds :=
  Weight(2*sqrt(lu)*exp(-lu)/sqrt(Pi),
  Uniform(-sqrt(2*lu),sqrt(2*lu))):
TestHakaru(sliced, sliceds, label = "slice sampling", ctx = [lu>0]);
module()
  local d, m, uMax, kb, result;
  for d in [Gaussian(0,1), GammaD(1,1)] do
    m := Bind(d, x,
              piecewise(And(0<u, u<density[op(0,d)](op(d))(x)),
              Weight(1/density[op(0,d)](op(d))(x),
              Ret(x)))):
    uMax := maximize(density[op(0,d)](op(d))(x), x):
    kb := KB:-assert(And(0<u, u<uMax), KB:-empty):
    result := fromLO(improve(toLO(m), _ctx=kb), _ctx=kb):
    CodeTools[Test](type(result, 'Weight(anything, Uniform(anything, anything))'),
                    true,
                    label = sprintf("slice sampling %a", d)):
  end do;
end module:

# test with uniform.  No change without simplifier, eliminates it with
# call to value.
introLO_opt := {
  Bind(Uniform(0,1),x,
  Bind(Uniform(0,1),y,
  piecewise(x<y,Ret(true),x>=y,Ret(false)))),
  Bind(Uniform(0,1),x,
  Bind(Uniform(0,1),y,
  piecewise(x<y,Ret(true),Ret(false)))) }:
introLO := op(1,introLO_opt):
introLOs := Msum(Weight(1/2, Ret(false)), Weight(1/2, Ret(true))):

TestHakaru(introLO, introLO_opt, simp = (x -> x), label = "2 uniform - no change");
TestHakaru(introLO, introLOs, simp = ((x,y) -> value(x)), label = "2 uniform + value = elimination");
TestHakaru(introLO, introLOs, label = "2 uniform + simplifier  elimination");

# a variety of other tests
TestHakaru(LO(h,(x*applyintegrand(h,5)+applyintegrand(h,3))/x), Weight(1/x, Msum(Weight(x, Ret(5)), Ret(3))));
TestHakaru(Bind(Gaussian(0,1),x,Weight(x,Msum())), Msum());
TestHakaru(Bind(Uniform(0,1),x,Weight(x^alpha,Ret(x))), Weight(1/(1+alpha),BetaD(1+alpha,1)));
TestHakaru(Bind(Uniform(0,1),x,Weight(x^alpha*(1-x)^beta,Ret(x))), Weight(Beta(1+alpha,1+beta),BetaD(1+alpha,1+beta)));
TestHakaru(Bind(Uniform(0,1),x,Weight((1-x)^beta,Ret(x))), Weight(1/(1+beta),BetaD(1,1+beta)));

# tests that basic densities are properly recognized
# continuous
TestHakaru(Bind(Uniform(0,1),x,Weight(x*2,Ret(x))), BetaD(2,1),
  label="BetaD(2,1) recog");
TestHakaru(BetaD(alpha,beta), label="BetaD recog.");
TestHakaru(GammaD(a,b), label="GammaD recog.", ctx = [a>0,b>0]);
TestHakaru(GammaD(1/2,2), label="GammaD(1/2,2) recog.");
TestHakaru(LO(h, int(exp(-x/2)*applyintegrand(h,x),x=0..infinity)), Weight(2,GammaD(1,2)), label="GammaD(1,2) recog.");
TestHakaru(LO(h, int(x*exp(-x/2)*applyintegrand(h,x),x=0..infinity)), Weight(4,GammaD(2,2)), label="GammaD(2,2) recog.");
TestHakaru(Bind(Lebesgue(-infinity,infinity), x, Weight(1/x^2, Ret(x))));
TestHakaru(Cauchy(loc,scale), ctx = [scale>0], label="Cauchy recog.");
TestHakaru(StudentT(nu,loc,scale), ctx=[nu>0,scale>0], label = "StudentT recog.");
TestHakaru(StudentT(1,loc,scale),Cauchy(loc,scale), ctx = [scale>0],
  label = "StudentT(1,loc,scale) recog.");
# discrete
TestHakaru(Weight(1/26,Counting(17,43)), label="Discrete uniform recog.");
TestHakaru(NegativeBinomial(r,1-p-q), label="NegativeBinomial recog.");
TestHakaru(Poisson(foo/exp(bar)), label="Poisson recog.");

# how far does myint get us?
TestHakaru(
  Bind(Uniform(0,1),x,Weight(x,Ret(Unit))),
  Weight(1/2,Ret(Unit)),
  label = "eliminate Uniform");

# just the front-end is already enough to get this
TestHakaru(
  Bind(Weight(1/2,Ret(Unit)),x,Ret(Unit)),
  Weight(1/2, Ret(Unit)),
  label = "integrate at work");

# and more various
model_exp := Bind(Uniform(-1,1),x,Ret(exp(x))):
TestHakaru(model_exp, model_exp, label = "uniform -1..1 into exp");
TestHakaru(LO(h, IntegrationTools[Expand](Int((1+y)*applyintegrand(h,y),y=0..1))), Msum(Uniform(0,1), Weight(1/2,BetaD(2,1))));
TestHakaru(Bind(Uniform(0,1),x,Bind(LO(h, IntegrationTools[Expand](Int((1+y)*applyintegrand(h,y),y=0..1))),y,Ret([x,y]))), Weight(3/2,Bind(Uniform(0,1),x,Msum(Weight(2/3,Bind(Uniform(0,1),y,Ret([x,y]))),Weight(1/3,Bind(BetaD(2,1),y,Ret([x,y])))))));

# easy-easy-HMM
eeHMM := Bind(GammaD(1,1),t,
                 Weight(NewSLO:-density[Gaussian](0,1/sqrt(t))(a),
                 Ret(t))):
ees := Weight(1/(a^2+2)^(3/2), GammaD(3/2, 1/((1/2)*a^2+1))):

TestHakaru(eeHMM, ees, ctx = [a::real], label = "easy-easy-HMM");

# from an email conversation on Sept. 11
model6 := Bind(Gaussian(0,1),x, piecewise(x>4,Ret(4),Ret(x))):
TestHakaru(model6, model6, label = "clamped Gaussian");

# and now models (then tests) taken from Tests.RoundTrip
t1 := Bind(Uniform(0, 1), a0, Msum(Weight(a0, Ret(Unit)))):
t2 := BetaD(1,1):
t2s := Uniform(0, 1):
t3 := Gaussian(0,10):
t4 := Bind(BetaD(1, 1), a0,
      Bind(Msum(Weight(a0, Ret(true)),
                Weight((1-a0), Ret(false))), a1,
      Ret(Pair(a0, a1)))):
t4s := Bind(Uniform(0, 1), a0,
       Msum(Weight(a0, Ret(Pair(a0, true))),
            Weight((1+(a0*(-1))), Ret(Pair(a0, false))))):
t5 := Bind(Msum(Weight((1/2), Ret(Unit))), a0, Ret(Unit)):
t5s := Weight((1/2), Ret(Unit)):
t6 := Ret(5):
t7 := Bind(Uniform(0, 1), a0,
  Bind(Msum(Weight((a0+1), Ret(Unit))), a1, Ret((a0*a0)))):
t7s := Bind(Uniform(0,1),a3,Weight(a3+1,Ret(a3^2))):
t7n := Bind(Uniform((-1), 0), a0,
  Bind(Msum(Weight((a0+1), Ret(Unit))), a1, Ret((a0*a0)))):
t7ns := Bind(Uniform(-1,0),a3,Weight(a3+1,Ret(a3^2))):
t8 := Bind(Gaussian(0, 10), a0, Bind(Gaussian(a0, 20), a1, Ret(Pair(a0, a1)))):
t9 := Bind(Lebesgue(-infinity,infinity), a0,
  Bind(Msum(Weight(piecewise(And((3<a0), (a0<7)), (1/2), 0), Ret(Unit))), a1, Ret(a0))):
t9a := Bind(Lebesgue(-infinity,infinity), a0,
  piecewise(3>=a0, Msum(), a0>=7, Msum(), Weight(1/2, Ret(a0)))):
t9s := Weight(2, Uniform(3,7)):

#t23, "bayesNet", to show exact inference.  Original used bern, which
# is here expanded in terms of MSum.
t23 :=
  Bind(Msum(Weight((1/2), Ret(true)), Weight((1-(1/2)), Ret(false))), a0,
  Bind(Msum(Weight(piecewise(a0 = true, (9/10), (1/10)), Ret(true)),
            Weight((1-piecewise(a0 = true, (9/10), (1/10))), Ret(false))), a1,
  Bind(Msum(Weight(piecewise(a0 = true, (9/10), (1/10)), Ret(true)),
            Weight((1-piecewise(a0 = true, (9/10), (1/10))), Ret(false))), a2,
  Ret(Pair(a1, a2))))):
t23s := Msum(Weight(41/100,Ret(Pair(true,true))),
             Weight(9/100,Ret(Pair(true,false))),
             Weight(9/100,Ret(Pair(false,true))),
             Weight(41/100,Ret(Pair(false,false)))):

# to exercise myint_pw
model_pw := Bind(Uniform(0,4), x,
  piecewise(x<1, Ret(x), x<2, Ret(2*x), x<3, Ret(3*x), Ret(5*x))):
model_pw1 := { Bind(Uniform(0,4), x,
  piecewise(x<1, Ret(x), x<2, Ret(2*x), x<3, Ret(3*x), 3<x, Ret(5*x))),
               Bind(Uniform(0,4),`x`,
  piecewise(x<1, Ret(x), 1<=x and x<2, Ret(2*x), 2<=x and x<3, Ret(3*x),3 <= x,Ret(5*x)))
}:
model_pw2 := Bind(Uniform(0,4), x, Weight(piecewise(x<1, 1, x<2, 2, x<3, 3, 5),Ret(x))):
model_pw3 := Bind(Uniform(0,4), x,
  piecewise(x<1, Ret(x), x<2, Weight(2,Ret(x)), x<3, Weight(3,Ret(x)), x>=3, Weight(5,Ret(x)))):
model_pw3_r := { Bind(Uniform(0,4), x,
  piecewise(x<1, Ret(x), x<2, Weight(2,Ret(x)), x<3, Weight(3,Ret(x)), x>=3, Weight(5,Ret(x)))) ,
               Bind(Uniform(0,4),x,
  piecewise(x < 1,Ret(x),1 <= x and x < 2,Weight(2,Ret(x)),2 <= x and x < 3,Weight(3,Ret(x)),3 <= x,Weight(5,Ret(x))))
}:
model_pw5 := Bind(Uniform(0,4), x, Weight(piecewise(x<1, 1, x<2, 2, x<3, 3, x>=3, 5),Ret(x))):
TestHakaru(model_pw , model_pw1, label = "multi-branch choice");
TestHakaru(model_pw3, model_pw3_r, label = "fake multi-branch weight");
TestHakaru(model_pw2, model_pw5, label = "proper multi-branch weight");

fake_pw1 := Bind(Uniform(1,3),x, piecewise(0 < x, Msum(), Ret(x))):
fake_pw2 := Bind(Uniform(1,3),x, piecewise(x < 0, Ret(x), Msum())):
TestHakaru(fake_pw1, Msum(), label = "fake piecewise 1"):
TestHakaru(fake_pw2, Msum(), label = "fake piecewise 2"):

# t43 without the explicit lam
t43 := piecewise(x0=true, Uniform(0, 1), Bind(BetaD(1, 1), a1, Ret(a1))):
t43s := Uniform(0, 1):

t80 := Bind(GammaD(1, 1), a0, Gaussian(0, a0)):

TestHakaru(t1, t5s, label = "t1");
TestHakaru(t2, t2s, label = "t2");
TestHakaru(t3, t3, label = "t3");
TestHakaru(t4, t4s, label = "t4");
TestHakaru(t5, t5s, label = "t5");
TestHakaru(t6, t6, label = "t6");
TestHakaru(t7, t7s, label = "t7");
TestHakaru(t7n, t7ns, label = "t7n");
TestHakaru(t8, t8, label = "t8");
TestHakaru(t9, t9s, label = "t9");
TestHakaru(t9a, t9s, label = "t9a");
TestHakaru(t23, t23s, label = "t23");
TestHakaru(t43, t43s, label = "t43"):

TestHakaru(t80, t80, label = "t80");

## based on https://github.com/hakaru-dev/hakaru/issues/60
module()
  local kb, clamp_t, clamp_r;
  kb := [y::real]:
  clamp_t :=
    Bind(Gaussian(0,1),x,
         Partition(And(0<y,y<density[Gaussian](0,1)(x)),
                   Weight(density[Uniform](0,density[Gaussian](0,1)(x))(y),
                          Ret(x)),
                   Msum())):
  clamp_r :=
    Partition(And(0<y, y<density[Gaussian](0,1)(0)),
              Weight(2 * sqrt(-ln(2)-ln(Pi)-2*ln(y)),
                     Uniform(-sqrt(-ln(2)-ln(Pi)-2*ln(y)),
                             sqrt(-ln(2)-ln(Pi)-2*ln(y)))),
              Msum()):
  TestHakaru( clamp_t, clamp_r, ctx=kb, label="clamp condition to move it out of integral" ):
end module:

###
# From disintegration paper
disint1 :=
Bind(Lebesgue(-infinity,infinity),y, Weight(piecewise(0<y and y<1, 1, 0), Weight(y/2, Ret(y)))):

TestHakaru(disint1, Weight(1/4,BetaD(2,1)), label="minor miracle");

ind1  := Bind(Uniform(0,1),x, Weight(piecewise(x>0,1,0), Weight(piecewise(x>1/2,0,1), Weight(piecewise(0<x,1,0), Ret(x))))):
ind1s := Weight(1/2, Uniform(0,1/2)):
ind2  := Bind(Lebesgue(-infinity,infinity),x, Weight(piecewise(x<0,0,x<1,x,0), Ret(x))):
ind2s := Weight(1/2, BetaD(2,1)):
ind3  := Bind(Uniform(0,1),x, Weight(piecewise(1<x and x<0,1,0), Ret(x))):
ind3s := Msum():
TestHakaru(ind1, ind1s, label="exponentiated indicator");
TestHakaru(ind2, ind2s, label="negated and conjoined indicator");
TestHakaru(ind3, ind3s, label="bounds ordering");
TestHakaru(Msum(ind1,ind2), Msum(ind1s,ind2s), label="simplify under sum");
TestHakaru(piecewise(c>0,ind1,ind2), piecewise(c>0,ind1s,ind2s), label="simplify under piecewise");

# test how banish handles piecewise
m1 := Uniform(0,1):
m2 := BetaD(2,1):
m3 := BetaD(1,2):
bp := proc() Bind(Gaussian(0,1), x,
             Bind(Gaussian(0,1), y,
             piecewise(_passed))) end proc:
TestHakaru(bp(x>y, m1, m2         ), Msum(Weight(1/2, m1), Weight(1/2, m2)                 ), label="banish piecewise 1");
TestHakaru(bp(x>0, m1, m2         ), Msum(Weight(1/2, m1), Weight(1/2, m2)                 ), label="banish piecewise 2");
TestHakaru(bp(y>0, m1, m2         ), Msum(Weight(1/2, m1), Weight(1/2, m2)                 ), label="banish piecewise 3");
TestHakaru(bp(x>y, m1, x>0, m2, m3), Msum(Weight(1/2, m1), Weight(1/8, m2), Weight(3/8, m3)), label="banish piecewise 4");
TestHakaru(bp(x>0, m1, x>y, m2, m3), Msum(Weight(1/2, m1), Weight(1/8, m2), Weight(3/8, m3)), label="banish piecewise 5");
TestHakaru(bp(y>x, m1, y>0, m2, m3), Msum(Weight(1/2, m1), Weight(1/8, m2), Weight(3/8, m3)), label="banish piecewise 6");
TestHakaru(bp(y>0, m1, y>x, m2, m3), Msum(Weight(1/2, m1), Weight(1/8, m2), Weight(3/8, m3)), label="banish piecewise 7");

# Simplify is not yet idempotent
TestHakaru(Bind(Uniform(0,1), x, Weight(x, Uniform(0,x))), Weight(1/2, BetaD(1, 2)));

# Test for change of variables; see Tests/Relationships.hs
# t1
# cv1 := Bind(Gaussian(mu, sigma), x, Ret((x-mu)/sigma)):
# cv1s := Gaussian(0,1):
# TestHakaru(cv1, cv1s, ctx = [sigma>0], label = "renormalize Gaussian");

# t28
# cv2 := Bind(BetaD(a,b), x, Ret(1-x)):
# cv2s := BetaD(b,a):
# TestHakaru(cv2, cv2s, ctx=[a>0,b>0], label = "swap BetaD");

unk_pw := Bind(m, y, Bind(Gaussian(0,1), x, piecewise(x<0, Ret(-x), Ret(x)))):
unk1   := Bind(Gaussian(0,1), x, Bind(m, y, Bind(Gaussian(x,1), z, Ret([y,z])))):
unk1s  := Bind(m, y, Bind(Gaussian(0,sqrt(2)), z, Ret([y,z]))):
unk2   := Bind(Gaussian(0,1), x, Bind(Gaussian(x,1), z, Bind(m, y, Ret([y,z])))):
unk2s  := Bind(Gaussian(0,sqrt(2)), z, Bind(m, y, Ret([y,z]))):
unk3   := Bind(Gaussian(0,1), x, Bind(m(x), y, Bind(Gaussian(x,1), z, Ret([y,z])))):
unk4   := Bind(Gaussian(0,1), x, Bind(Gaussian(x,1), z, Bind(m(x), y, Ret([y,z])))):
TestHakaru(unk_pw, unk_pw, label="Don't simplify Integrand willy-nilly");
TestHakaru(unk1, unk1s, label="Banish into Integrand 1");
TestHakaru(unk2, unk2s, label="Banish into Integrand 2");
TestHakaru(unk3, unk3, label="Banish into Integrand 3");
TestHakaru(unk4, unk4, label="Banish into Integrand 4");

# Disintegration of easierRoadmapProg1 -- variables to be integrated out
rmProg1 := Msum(Weight(1, Msum(Weight(1, Msum(Weight(1, Msum(Weight(1, Bind(Lebesgue(-infinity,infinity), a4, Msum(Weight(1, Msum(Weight(1, Bind(Lebesgue(-infinity,infinity), a5, Msum(Weight(1, Bind(Lebesgue(-infinity,infinity), a6, Msum(Weight(((exp((-(((p3-a6)*(p3-a6))*(1/(2*exp((ln(a5)*2)))))))*(1/a5))*(1/exp((ln((2*Pi))*(1/2))))), Msum(Weight(1, Bind(Lebesgue(-infinity,infinity), a7, Msum(Weight(((exp((-(((a6-a7)*(a6-a7))*(1/(2*exp((ln(a4)*2)))))))*(1/a4))*(1/exp((ln((2*Pi))*(1/2))))), Msum(Weight(((exp((-(((p2-a7)*(p2-a7))*(1/(2*exp((ln(a5)*2)))))))*(1/a5))*(1/exp((ln((2*Pi))*(1/2))))), Msum(Weight(((exp((-((a7*a7)*(1/(2*exp((ln(a4)*2)))))))*(1/a4))*(1/exp((ln((2*Pi))*(1/2))))), Msum(Weight((1/3), Msum(Weight(1, piecewise((a5<4), piecewise((1<a5), Msum(Weight((1/5), Msum(Weight(1, piecewise((a4<8), piecewise((3<a4), Ret(Pair(a4, a5)), Msum()), Msum())), Weight(1, Msum())))), Msum()), Msum())), Weight(1, Msum())))))))))))))))))))))), Weight(1, Msum())))))), Weight(1, Msum())))))):
rmProg1_w := Weight(exp(-(1/2)*(2*p2^2*a4^2+p2^2*a5^2-2*p2*p3*a4^2+p3^2*a4^2+p3^2*a5^2)/(a4^4+3*a4^2*a5^2+a5^4))/sqrt(a4^4+3*a4^2*a5^2+a5^4), Ret(Pair(a4, a5))):
rmProg1_r := {
Weight(1/(2*Pi), Bind(Uniform(1, 4), a5, Bind(Uniform(3, 8), a4, rmProg1_w))),
Weight(1/(2*Pi), Bind(Uniform(3, 8), a4, Bind(Uniform(1, 4), a5, rmProg1_w)))}:
TestHakaru(rmProg1, rmProg1_r, label="Tests.RoundTrip.rmProg1");
# easierRoadmapProg4 -- MH transition kernel with unclamped acceptance ratio
module()
  local unpair, fst, snd, rmProg4;
  unpair := proc(f,p)
    f(fst(p),snd(p))
  end proc;
  fst := proc(p, $)
    if p :: Pair(anything,anything) then op(1,p)
    elif p :: 't_piecewiselike' then map_piecewiselike(procname, _passed)
    else 'procname(_passed)' end if
  end proc:
  snd := proc(p, $)
    if p :: Pair(anything,anything) then op(2,p)
    elif p :: 't_piecewiselike' then map_piecewiselike(procname, _passed)
    else 'procname(_passed)' end if
  end proc:
  rmProg4 := lam(x0,blah, app(lam(x1,blah, lam(x2,blah, Bind(unpair(( (p3 , p4) -> Msum(Weight((1/2), Bind(Uniform(3, 8), a5, Ret(Pair(a5, p4)))), Weight((1/2), Bind(Uniform(1, 4), a6, Ret(Pair(p3, a6)))))), x2), a7, Ret(Pair(a7, (app(x1,a7)/app(x1,x2))))))),lam(x8,blah, unpair(( (p151 , p152) -> app(p152,lam(x153,blah, 1))), app(app(lam(x9,blah, lam(x10,blah, unpair(( (p11 , p12) -> app(lam(x13,blah, app(lam(x14,blah, Pair(Msum(Weight(x13, unpair(( (p15 , p16) -> p15), x14))), lam(x17,blah, (0+(x13*app(unpair(( (p18 , p19) -> p19), x14),x17)))))),app(lam(x20,blah, app(lam(x21,blah, Pair(Msum(Weight(x20, unpair(( (p22 , p23) -> p22), x21))), lam(x24,blah, (0+(x20*app(unpair(( (p25 , p26) -> p26), x21),x24)))))),app(lam(x27,blah, app(lam(x28,blah, app(lam(x29,blah, app(lam(x30,blah, Pair(Msum(Weight(x27, unpair(( (p31 , p32) -> p31), x28)), Weight(x29, unpair(( (p33 , p34) -> p33), x30))), lam(x35,blah, ((0+(x27*app(unpair(( (p36 , p37) -> p37), x28),x35)))+(x29*app(unpair(( (p38 , p39) -> p39), x30),x35)))))),Pair(Msum(), lam(x40,blah, 0)))),1)),app(lam(x41,blah, app(lam(x42,blah, Pair(Msum(Weight(x41, unpair(( (p43 , p44) -> p43), x42))), lam(x45,blah, (0+(x41*app(unpair(( (p46 , p47) -> p47), x42),x45)))))),app(lam(x48,blah, app(lam(x49,blah, app(lam(x50,blah, app(lam(x51,blah, Pair(Msum(Weight(x48, unpair(( (p52 , p53) -> p52), x49)), Weight(x50, unpair(( (p54 , p55) -> p54), x51))), lam(x56,blah, ((0+(x48*app(unpair(( (p57 , p58) -> p58), x49),x56)))+(x50*app(unpair(( (p59 , p60) -> p60), x51),x56)))))),Pair(Msum(), lam(x61,blah, 0)))),1)),app(lam(x62,blah, app(lam(x63,blah, Pair(Msum(Weight(x62, unpair(( (p64 , p65) -> p64), x63))), lam(x66,blah, (0+(x62*app(unpair(( (p67 , p68) -> p68), x63),x66)))))),unpair(( (p69 , p70) -> unpair(( (p71 , p72) -> unpair(( (p73 , p74) -> unpair(( (p75 , p76) -> unpair(( (p77 , p78) -> unpair(( (p79 , p80) -> unpair(( (p81 , p82) -> unpair(( (p83 , p84) -> unpair(( (p85 , p86) -> unpair(( (p87 , p88) -> app(lam(x89,blah, app(lam(x90,blah, Pair(Msum(Weight(x89, unpair(( (p91 , p92) -> p91), x90))), lam(x93,blah, (0+(x89*app(unpair(( (p94 , p95) -> p95), x90),x93)))))),app(lam(x96,blah, app(lam(x97,blah, Pair(Msum(Weight(x96, unpair(( (p98 , p99) -> p98), x97))), lam(x100,blah, (0+(x96*app(unpair(( (p101 , p102) -> p102), x97),x100)))))),app(lam(x103,blah, app(lam(x104,blah, app(lam(x105,blah, app(lam(x106,blah, Pair(Msum(Weight(x103, unpair(( (p107 , p108) -> p107), x104)), Weight(x105, unpair(( (p109 , p110) -> p109), x106))), lam(x111,blah, ((0+(x103*app(unpair(( (p112 , p113) -> p113), x104),x111)))+(x105*app(unpair(( (p114 , p115) -> p115), x106),x111)))))),Pair(Msum(), lam(x116,blah, 0)))),1)),piecewise((p12<4), piecewise((1<p12), app(lam(x117,blah, app(lam(x118,blah, Pair(Msum(Weight(x117, unpair(( (p119 , p120) -> p119), x118))), lam(x121,blah, (0+(x117*app(unpair(( (p122 , p123) -> p123), x118),x121)))))),app(lam(x124,blah, app(lam(x125,blah, app(lam(x126,blah, app(lam(x127,blah, Pair(Msum(Weight(x124, unpair(( (p128 , p129) -> p128), x125)), Weight(x126, unpair(( (p130 , p131) -> p130), x127))), lam(x132,blah, ((0+(x124*app(unpair(( (p133 , p134) -> p134), x125),x132)))+(x126*app(unpair(( (p135 , p136) -> p136), x127),x132)))))),Pair(Msum(), lam(x137,blah, 0)))),1)),piecewise((p11<8), piecewise((3<p11), app(lam(x138,blah, app(lam(x139,blah, Pair(Msum(Weight(x138, unpair(( (p140 , p141) -> p140), x139))), lam(x142,blah, (0+(x138*app(unpair(( (p143 , p144) -> p144), x139),x142)))))),app(lam(x145,blah, Pair(Ret(x145), lam(x146,blah, app(x146,x145)))),Pair(p11, p12)))),5), Pair(Msum(), lam(x147,blah, 0))), Pair(Msum(), lam(x148,blah, 0))))),1))),(1/5)), Pair(Msum(), lam(x149,blah, 0))), Pair(Msum(), lam(x150,blah, 0))))),1))),(1/3)))),((((1/Pi)*exp((((((((((p69*p71)*(p11*p11))*2)+((((p11*p11)*p73)*p76)*(-2)))+((p78*p80)*(p11*p11)))+((p12*p12)*(p81*p83)))+((p12*p12)*(p86*p88)))*(1/((((p11*p11)*(p11*p11))+(((p12*p12)*(p11*p11))*3))+((p12*p12)*(p12*p12)))))*(-(1/2)))))*exp((ln(((exp(((x -> piecewise(x<0, -337, ln(x)))(p11)*4))+((exp(((x -> piecewise(x<0, -337, ln(x)))(p12)*2))*exp(((x -> piecewise(x<0, -337, ln(x)))(p11)*2)))*3))+exp(((x -> piecewise(x<0, -337, ln(x)))(p12)*4))))*(-(1/2)))))*(1/10)))), x9)), x9)), x9)), x9)), x9)), x9)), x9)), x9)), x9)), x9))),1))),1))),1))),1))),1))),1)), x10))),x0),x8))))):
  TestHakaru(app(app(rmProg4,Pair(r1,r2)),Pair(p1,p2)), Msum(Weight(1/2, Bind(Uniform(3, 8), a5, Ret(Pair(Pair(a5, p2), exp((1/2)*(-p1+a5)*(a5+p1)*(p1^2*p2^2*r1^2+p1^2*p2^2*r2^2+2*p1^2*r1^2*a5^2-2*p1^2*r1*r2*a5^2+p1^2*r2^2*a5^2+p2^4*r1^2+2*p2^4*r1*r2+2*p2^4*r2^2+p2^2*r1^2*a5^2+p2^2*r2^2*a5^2)/((p2^4+3*p2^2*a5^2+a5^4)*(p1^4+3*p1^2*p2^2+p2^4)))*sqrt(p1^4+3*p1^2*p2^2+p2^4)/sqrt(p2^4+3*p2^2*a5^2+a5^4))))), Weight(1/2, Bind(Uniform(1, 4), a6, Ret(Pair(Pair(p1, a6), exp((1/2)*(-p2+a6)*(a6+p2)*(5*p1^4*r1^2-6*p1^4*r1*r2+2*p1^4*r2^2+2*p1^2*p2^2*r1^2-2*p1^2*p2^2*r1*r2+p1^2*p2^2*r2^2+2*p1^2*r1^2*a6^2-2*p1^2*r1*r2*a6^2+p1^2*r2^2*a6^2+p2^2*r1^2*a6^2+p2^2*r2^2*a6^2)/((p1^4+3*p1^2*a6^2+a6^4)*(p1^4+3*p1^2*p2^2+p2^4)))*sqrt(p1^4+3*p1^2*p2^2+p2^4)/sqrt(p1^4+3*p1^2*a6^2+a6^4)))))),
  label="rmProg4",
  ctx= [3<p1, p1<8, 1<p2, p2<4]);
end module:

#####################################################################
#
# conjugacy tests
#
#####################################################################
gaussian_gaussian   := Bind(Gaussian(mu0,sigma0),mu, Weight(NewSLO:-density[Gaussian](mu,sigma1)(x), Ret(mu))):
gaussian_gaussian_s := Weight((1/2)*sqrt(2)*exp(-(1/2)*(mu0-x)^2/(sigma0^2+sigma1^2))/(sqrt(Pi)*sqrt(sigma0^2+sigma1^2)), Gaussian((mu0*sigma1^2+sigma0^2*x)/(sigma0^2+sigma1^2), sigma0*sigma1/sqrt(sigma0^2+sigma1^2))):
TestHakaru(gaussian_gaussian, gaussian_gaussian_s,
  label="gaussian_gaussian conjugacy",
  ctx = [mu0::real, sigma0>0, sigma1>0, x::real]);
invgamma_gaussian   := Bind(GammaD(shape,scale),lambda, Weight(NewSLO:-density[Gaussian](mu,lambda^(-1/2))(x), Ret(lambda))):
invgamma_gaussian_s := Weight(GAMMA(1/2+shape)*sqrt(scale)*((1/2)*scale*mu^2-scale*mu*x+(1/2)*scale*x^2+1)^(-shape)/(GAMMA(shape)*sqrt(scale*mu^2-2*scale*mu*x+scale*x^2+2)*sqrt(Pi)), GammaD(1/2+shape, 2*scale/(scale*mu^2-2*scale*mu*x+scale*x^2+2))):
TestHakaru(invgamma_gaussian, invgamma_gaussian_s,
  label="invgamma_gaussian conjugacy",
  ctx = [mu::real, shape>0, scale>0, x::real]);
gaussian_invgamma_gaussian   := Bind(GammaD(shape,scale),tau, Bind(Gaussian(mu0*sqrt(tau),1/sqrt(nu)),mu, Weight(NewSLO:-density[Gaussian](mu*tau^(-1/2),tau^(-1/2))(x), Ret([mu,tau])))):
gaussian_invgamma_gaussian_s := Weight(GAMMA(1/2+shape)*sqrt(nu)*sqrt(scale)*(scale*mu0^2*nu-2*scale*mu0*nu*x+scale*nu*x^2+2*nu+2)^(-1/2-shape)*(2*nu+2)^shape/(GAMMA(shape)*sqrt(Pi)), Bind(GammaD(1/2+shape, 2*scale*(nu+1)/(scale*mu0^2*nu-2*scale*mu0*nu*x+scale*nu*x^2+2*nu+2)),tau, Bind(Gaussian(sqrt(tau)*(mu0*nu+x)/(nu+1), 1/sqrt(nu+1)),mu, Ret([mu,tau])))):
TestHakaru(gaussian_invgamma_gaussian, gaussian_invgamma_gaussian_s,
  label="gaussian_invgamma_gaussian conjugacy",
  ctx = [mu0::real, nu>0, shape>0, scale>0, x::real]);
gamma_gamma   := Bind(GammaD(alpha0,1/beta0),beta, Weight(NewSLO:-density[GammaD](alpha,1/beta)(x), Ret(beta))):
gamma_gamma_s := Weight(beta0^alpha0*x^(alpha-1)*GAMMA(alpha+alpha0)*(beta0+x)^(-alpha-alpha0)/(GAMMA(alpha0)*GAMMA(alpha)), GammaD(alpha+alpha0, 1/(beta0+x))):
TestHakaru(gamma_gamma, gamma_gamma_s, label="gamma_gamma conjugacy",
  ctx = [alpha0>0, beta0>0, alpha>0, x>0]);
gamma_poisson   := Bind(GammaD(shape,scale),lambda, Weight(NewSLO:-density[PoissonD](lambda)(k), Ret(lambda))):
gamma_poisson_s := Weight(scale^k*(scale+1)^(-k-shape)*GAMMA(k+shape)/(GAMMA(shape)*GAMMA(k+1)), GammaD(k+shape, scale/(scale+1))):
TestHakaru(gamma_poisson, gamma_poisson_s, label="gamma_poisson conjugacy",
           ctx=[shape>0, scale>0, k::nonnegint]);

# For the following test, banishing fails because Maple currently evaluates
#   int(piecewise(x < y, 1/(1-x), 0), x = 0 .. 1) assuming y<1
# to "undefined".  (See ppaml-l discussion "NewSLO giving weird output"
# around 2015-12-10.)  The test is that we handle the failure gracefully.
TestHakaru(Bind(Uniform(0,1), x, Uniform(x,1)),
          {Bind(Uniform(0,1), x, Uniform(x,1)),
           Bind(Uniform(0,1), x, Weight(ln(1/(1-x)), Ret(x)))},
  label="roundtrip despite banishing failure");

TestHakaru(Bind(Ret(ary(n,i,i*2)), v, Ret(idx(v,42))), Ret(84), label="basic array indexing");

#####################################################################
#
# Tests for type-directed simplification
#
#####################################################################

module()
  local hpair, hbool,         heither,       hreal, hprob,
        ppair, ptrue, pfalse, pleft, pright,
        dpair, dtrue, dfalse, dleft, dright,
        unpair;
  uses CodeTools;
  hpair   := (t1,t2) -> HData(DatumStruct(pair,[Konst(t1),Konst(t2)])):
  hbool   := HData(DatumStruct(true,[]),DatumStruct(false,[])):
  heither := (t1,t2) -> HData(DatumStruct(left,[Konst(t1)]),DatumStruct(right,[Konst(t2)])):
  hreal   := HReal():
  hprob   := HReal(Bound(`>=`,0)):
  ppair   := (p1,p2) -> PDatum(pair,PInl(PEt(p1,PEt(p2,PDone)))):
  ptrue   := PDatum(true,PInl(PDone)):
  pfalse  := PDatum(false,PInr(PInl(PDone))):
  pleft   := p -> PDatum(left,PInl(PEt(p,PDone))):
  pright  := p -> PDatum(right,PInr(PInl(PEt(p,PDone)))):
  dpair   := (d1,d2) -> Datum(pair,Inl(Et(d1,Et(d2,Done)))):
  dtrue   := Datum(true,Inl(Done)):
  dfalse  := Datum(false,Inr(Inl(Done))):
  dleft   := d -> Datum(left,Inl(Et(d,Done))):
  dright  := d -> Datum(right,Inr(Inl(Et(d,Done)))):
  unpair  := (e1,x,y,e0) -> case(e1,Branches(Branch(ppair(PVar(x),PVar(y)),e0))):
  TestSimplify(lam(pr,hpair(hprob,hreal),unpair(pr,x,y,sqrt(x^2)+y)),
               HFunction(hpair(hprob,hreal),hreal),
               lam(pr,hpair(hprob,hreal),unpair(pr,x,y,x+y)),
               label="Transfer function argument type to KB");
  TestSimplify(lam(b,hbool,case(b,Branches(Branch(pfalse,3),Branch(ptrue,4)))),
               HFunction(hbool,HInt()),
               lam(b,hbool,case(b,Branches(Branch(ptrue,4),Branch(pfalse,3)))),
               label="Eta-expand pure function from boolean");
  TestSimplify(lam(pr,hpair(hbool,hbool),unpair(pr,x,y,dpair(y,x))),
               HFunction(hpair(hbool,hbool), hpair(hbool,hbool)),
               lam(pr,hpair(hbool,hbool),case(pr,Branches(
                 Branch(ppair(ptrue ,ptrue ), dpair(dtrue ,dtrue )),
                 Branch(ppair(ptrue ,pfalse), dpair(dfalse,dtrue )),
                 Branch(ppair(pfalse,ptrue ), dpair(dtrue ,dfalse)),
                 Branch(ppair(pfalse,pfalse), dpair(dfalse,dfalse))))),
               label="Eta-expand pure function from pair of booleans");
  TestSimplify(lam(b,hbool,
                     Bind(Uniform(0,1),p,
                     Weight(case(b,Branches(Branch(pfalse,1-p),
                                            Branch(ptrue,p))),
                     Ret(p)))),
               HFunction(hbool,HMeasure(hreal)),
               lam(b,hbool,
                     case(b,Branches(Branch(ptrue , Weight(1/2,BetaD(2,1))),
                                     Branch(pfalse, Weight(1/2,BetaD(1,2)))))),
               label="Eta-expand function from boolean to measure");
end module:

###########
#
# RoundTrip tests; due to bug in CodeTest, use lower-level function directly

TestTools:-SetRecord(false):
module()
  local kb, rt1, rt1r;
  kb := [x::real]:
  rt1 := Context(kb, Bind(Gaussian(0,1), x, Ret(0))):
  rt1r := sprintf("%a", eval(ToInert(Context(kb, Ret(0))), _Inert_ATTRIBUTE=NULL)):
  TestTools:-Try("simple case of RoundTrip", RoundTrip(rt1), rt1r);
end module:

# Categorical distribution
TestHakaru(Counting(0,11),
           label="Counting roundtrip");
TestHakaru(Bind(Counting(0,2),x,Weight(p^x*(1-p)^(1-x),Ret(x))),
           Categorical([1-p,p]),
           label="Categorical introduction",
           ctx=[0<=p, p<=1]);
TestHakaru(Bind(Categorical([p,1-p]),x,Weight(p^x*(1-p)^(1-x),Ret(x))),
           Weight(p*(1-p),Counting(0,2)),
           label="Categorical elimination",
           ctx=[0<=p, p<=1]);
module()
  local bern, burglary;
  bern := p -> Bind(Categorical(ary(2,i,piecewise(i=0,p,1-p))), i, Ret(i=0)):
  burglary := Bind(bern(1/10000), b,
              Bind(bern(piecewise(b=true,19/20,1/100)), a,
              Ret(a))):
  TestHakaru(burglary, bern(5047/500000), label="burglary");
  TestHakaru(Bind(BetaD(alpha,beta), p, bern(p)),
             Weight(1/(alpha+beta), Bind(Categorical([alpha,beta]), i, Ret(i=0))),
             label="integrate BetaD out of BetaD-Bernoulli"):
end module:


text:
# The Boys embedding, as a map from the Riemann sphere to R^3

complex_boys_embedding := proc(z)
 local a,V,M,B,rM,i;

 if z = 0 then return [0,0,2]; fi;

 a := 1/(z^3-1/z^3+sqrt(5)):
 V := [I*(z^2-1/z^2), z^2+1/z^2, 2/3*I*(z^3+1/z^3)]:
 M := map(Re, a *~ V +~ [0, 0, 1/2]):
 rM := add(M[i]^2,i=1..3);
 B := M /~ rM;
 return B;
end:

sphere_reduce0 := (x) -> (u) -> rem(u,x[1]^2+x[2]^2+x[3]^2-1,x[3]);

sphere_reduce1 := (x) -> proc(u)
  local n,d,n0,d0,n1,d1,f;
  if type(u,list) then return map(sphere_reduce1(x),u); fi;

  n := sphere_reduce0(x)(numer(u));
  d := sphere_reduce0(x)(denom(u));
  n0 := coeff(n,x[3],0);
  n1 := coeff(n,x[3],1);
  d0 := coeff(d,x[3],0);
  d1 := coeff(d,x[3],1);
  n := sphere_reduce0(x)(expand(n * (d0 - d1*x[3])));
  d := sphere_reduce0(x)(expand(d * (d0 - d1*x[3])));
  return factor(n/d);
end:

# The map boys_embedding() is the Boys embedding, as a map from the
# round sphere S^2 to R^3.
make_boys_embedding := proc()
 local x,Z,B,BN,dB,i;
 global boys_embedding0,boys_normal0;

 assume(x[1] > 0, x[2] > 0, x[3] > 0);
 Z := (x[1] + I * x[2])/(1 - x[3]);
 B := complex_boys_embedding(Z);
 B := map(sphere_reduce1(x),B);
 dB := denom(B[1]);
 B := factor(dB *~ B);
 B := [op(B),dB];
 boys_embedding0 := unapply(B,x);

 BN := subs({seq(x[i] = t*x[i],i=1..3)},
            [B[1]/B[4],B[2]/B[4],B[3]/B[4]]);
 BN := subs(t = 1, map(diff,BN,t));
 BN := sphere_reduce1(x)(BN);
 BN := expand(factor(BN *~ denom(BN[1]) /~ 480)):

 boys_normal0 := unapply(BN,x);
end:

boys_embedding := proc(x)
 local B;
 B := boys_embedding0(x);
 return [B[1]/B[4],B[2]/B[4],B[3]/B[4]];
end:

boys_normal := proc(x)
 local n,M,i;
 if x = [0,0,1] then return [0,0,1]; fi;
 M := evalf(boys_jacobian_matrix(x));
 n := convert((1/Transpose(Matrix(M))) . Vector(x),list);
 n := n /~ sqrt(add(n[i]^2,i=1..3));
 return n;
end:

make_boys_jacobian := proc()
 local x,p,xred,B,BJM,BJN,BJNP,DBJNP,p_rule,x0,rot,i,j;
 global boys_jacobian_matrix,boys_jacobian_square_norm,boys_poly,boys_p0,boys_corners;

 B := boys_embedding(x);

 # The embedding map extends naturally over a neighbourhood of S2 in R3,
 # and this is the Jacobian of that extension.
 BJM := [seq([seq(sphere_reduce1(x)(diff(B[i],x[j])),j=1..3)],i=1..3)];
 boys_jacobian_matrix := unapply(BJM,x);

 # This is the square norm of the Jacobian restricted to S2
 BJN := sphere_reduce1(x)(expand(add(add(BJM[i][j]^2,j=1..3),i=1..3) -
                          add(add(BJM[i][j]*x[j],j=1..3)^2,i=1..3))):
 boys_jacobian_square_norm := unapply(BJN,x);

 p_rule := {x[1] = (1-p^2)/(1+p^2)/sqrt(2),
            x[2] = (1-p^2)/(1+p^2)/sqrt(2),
	    x[3] = 2 * p/(1+p^2)};

 BJNP := factor(expand(simplify(subs(p_rule,BJN))));
 DBJNP := factor(expand(diff(BJNP,p)));
 boys_poly := unapply(select(g -> degree(g,p) = 18,{op(numer(DBJNP))})[1],p);
 boys_p0 := fsolve(boys_poly(p),p=0.15..0.16);

 x0 := evalf(subs(p = boys_p0,subs(p_rule,[x[1],x[2],x[3]])));

 rot := (x) -> [(-x[1]-sqrt(3)*x[2])/2,(sqrt(3)*x[1]-x[2])/2,x[3]];
 boys_corners :=
  evalf([x0,rot(x0),rot(rot(x0)),
         -~ x0, -~ rot(x0), -~ rot(rot(x0))]);
end:

# We now construct a map boys_embedding_alt() which is obtained by precomposing
# boys_embedding() with a linear map and then postcomposing with an affine map.
# This has the effect that boys_embedding_alt sends [1,1,1]/sqrt(3) to itself
# and is equivariant for the third-twist about that axis.  Moreover, the
# "corner points" for boys_embedding_alt are at the other points of the form
# [+/- 1, +/- 1, +/- 1]/sqrt(3), where the three signs are not all the same.

make_boys_embedding_alt := proc()
 global boys_M0,boys_M1,boys_a1,boys_embedding_alt;
 local u,a,x,M,M0,be0,F,F0,rels,sol,i;

 if not(type(boys_corners,list(list))) then
  make_boys_jacobian():
 fi;

 u := table():
 u[ 0] := [ 1, 1, 1] /~ sqrt(3):
 u[ 1] := [ 1, 1,-1] /~ sqrt(3):
 u[ 2] := [ 1,-1, 1] /~ sqrt(3):
 u[ 3] := [-1, 1, 1] /~ sqrt(3):

 M := Matrix(3,3,[seq(a[i],i=1..9)]):

 M0 := subs(solve(evalf(map(op,[seq(convert(M . Vector(u[i]),list) -~ boys_corners[i],i=1..3)]))),M):
 boys_M0 := M0;
 be0 := unapply(simplify(expand(evalf(boys_embedding(M0 . <x[1],x[2],x[3]>)))),x):

 F := (y) -> convert(M . Vector(y),list) +~ [a[10],a[11],a[12]]:
 rels := map(op,evalf([
  F(be0(u[0])) -~ u[0], F(be0(u[1])) +~ u[1], F(be0(u[2])) +~ u[2], F(be0(u[3])) +~ u[3]
 ])):
 sol := solve(rels);
 boys_M1 := subs(sol,M);
 boys_a1 := subs(sol,[a[10],a[11],a[12]]);

 boys_embedding_alt := (x) ->
  evalf(convert(boys_M1 . Vector(boys_embedding(boys_M0 . <x[1],x[2],x[3]>)),list) +~
   boys_a1);

end:

ext: mpl


text:
# Collection of small procs to works with CRN
# No documentation provided, they are self explanatory.

with(LinearAlgebra);
with(ListTools);
with(ArrayTools);

## Suit of useful procs.
equfy := proc(u)
local i, v;
    v := copy(u);
    for i in seq(1..Size(v,1))
    do v[i] := v[i] = 0;
    end do;
    return v;
end proc:

## Angelica
SignCoeffs := proc(p, var)
    return ListTools[MakeUnique](map(signum, [coeffs(collect(p, convert(var,set), distributed), var)]));
end proc:

## Angelica
PolyToVector := proc(p, var)
local i, j, thecoeffs, terms;
    thecoeffs := [coeffs(collect(p, var, distributed), var, 'terms')];
    return [seq([thecoeffs[j], seq(degree(terms[j], var[i]), i = 1 .. numelems(var))], j = 1 .. numelems([terms]))];
end proc:

SetFirstToOne := proc(l::list)::list;
local i, t;
    t := Array(1 .. nops(l));
    for i to nops(l) do
        t[i] := l[i];
        t[i][1] := 1;
    end do;
    return convert(t, list);
end proc:

SearchAllPredicate := proc(f, L)
local i, e, index;
    i := 0;
    index := Vector[row]();
    for e in L do
        i := i + 1;
        if f(e) then
            index := <index | i>;
        end if;
    end do;
    return index;
end proc:

PotentiallyNegCoeffs := proc(poly, parameters)
local signcffs, index, negcffs, negexps;
    # polyvec := PolToVector(poly, vars);
    signcffs := map(X -> SignCoeffs(X, parameters), [seq(poly[i][1], i=1..nops(poly))]);
    print(signcffs);
    index := SearchAllPredicate(X -> (-1) in X, signcffs);
    negcffs := [seq(factor(poly[i][1]), i in index)];
    negexps := [seq(poly[i][2..], i in index)];
    print(Vector[column](negcffs));
    return (negcffs, negexps, index);
end proc:

krelations := proc(poly, ks)
local lc, facts, krel, fact;
    facts := factors(poly);
    krel := 1;
    for fact in facts[2] do
        if type(fact[2], odd) and (-1) in SignCoeffs(fact[1], {op(ks)}) then
            krel := krel*fact[1];
        end if;
    end do;
    return facts[1]*krel;
end proc:

SearchingMultistationarity := proc(sys, DJ, depenvars, xs, ks)
local param, freevars, DJparam, polyA, vecpolyA, negcffA, negexpA, negA, krels;
    param := solve(convert(equfy(sys), set), convert(depenvars, set));
    freevars := [op(convert(xs, set) minus convert(depenvars, set))];
    DJparam := subs(param, DJ);
    polyA := numer(DJparam);
    print("Denominator of DJparam:(check that its sign just depends on the numerator)");
    print(denom(DJparam));
    vecpolyA := PolyToVector(polyA, freevars);
    negcffA, negexpA, negA := PotentiallyNegCoeffs(vecpolyA, convert(ks, set));
    krels := map(X -> krelations(X, ks), negcffA);
    return (krels, negcffA, negexpA, negA, collect(polyA, freevars, 'distributed', factor), vecpolyA, param);
end proc:

## Angelica, Elisenda
StabilityMatrix := proc(M::Matrix)
local p, i, j, d, pol, m, k;
global H;
    d := LinearAlgebra[Dimension](M)[1] - LinearAlgebra[Rank](M);
    p := simplify(LinearAlgebra[CharacteristicPolynomial](M, y)/y^d);
    H := HurwitzDet(degree(p, y));
    m := numelems(H);
    for i to m do H[i] := subs(seq(a[k] = coeff(p, y, k), k = 0 .. degree(p, y)), H[i]);
    end do;
    print(m*'expressions*will*be*studied');
    for j to m - 1 do print(['Hurwitz*determinant'*'H'[j], SignCoeffs(numer(H[j]), indets(numer(H[j]))), SignCoeffs(denom(H[j]), indets(denom(H[j])))]);
    end do;
    print(["Lowest degree term", SignCoeffs(numer(H[m]), indets(numer(H[m]))), SignCoeffs(denom(H[m]), indets(denom(H[m])))]);
end proc:

## Angelica, Elisenda
HurwitzDet := proc(n::integer)::list;
local s, t, H, M, i, j, k;
    M := Matrix(n);
    for i to n do
        for j to n do
            if (n - 2*i + j) in [seq(t, t = 0 .. n)] then M[i, j] := a[n - 2*i + j];
            end if;
        end do;
    end do;
    H := [];
    for i to n - 1 do H := [op(H), factor(LinearAlgebra[Determinant](LinearAlgebra[SubMatrix](M, [seq(k, k = 1 .. i)], [seq(k, k = 1 .. i)])))];
    end do;
    H := [op(H), a[0]];
    return H;
end proc:

# splitkrel := proc(krel, ks)
# local poskrel, negkrel, terms, coefflist, i;
#     poskrel := 0;
#     negkrel := 0;
#     coefflist := coeffs(collect(krel, ks, distributed), ks, 'terms');
#     for i from 1 to numelems(coefflist) do
#         if coefflist[i] < 0 then
#             negkrel := negkrel - coefflist[i]*terms[i];
#         else
#             poskrel := poskrel + coefflist[i]*terms[i];
#         end if;
#     end do;
#     return [factor(negkrel), factor(poskrel)]
# end proc:

splitkrel := proc(krel, ks)
local poskrel, negkrel, terms, coefflist, ct;
    poskrel := 0;
    negkrel := 0;
    coefflist := coeffs(collect(krel, ks, distributed), ks, 'terms');
    for ct in zip(`[]`, [coefflist], [terms]) do
        if ct[1] < 0 then
            negkrel := negkrel - ct[1]*ct[2];
        else
            poskrel := poskrel + ct[1]*ct[2];
        end if;
    end do;
    return [factor(negkrel), factor(poskrel)]
end proc:

nicekrelstoLaTeX := proc(krel, ks)
local skrel;
    skrel := splitkrel(krel, ks);
    return latex(skrel[2]<skrel[1]);
end proc:

Velocities := proc(Y, vars)
local i, j, v;
    v := Vector[column](LinearAlgebra[ColumnDimension](Y));
    for i to LinearAlgebra[ColumnDimension](Y) do
        v[i] := 1;
        for j to LinearAlgebra[RowDimension](Y) do
            v[i] := v[i]*vars[j]^Y[j, i];
        end do;
    end do;
    return v;
end proc:

max_stars_repo_name: LauraBMo/CRNTscripts

max_issues_repo_name: LauraBMo/CRNTscripts

max_forks_repo_name: LauraBMo/CRNTscripts


text:
(*
 Copyright (C) 2020 Susi Lehtola

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

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

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

(* The m_min avoids divisions by zero *)
lkt_f0 := s -> 1/cosh(params_a_a * m_min(200, s)) + 5*s^2/3 :
lkt_f := x -> lkt_f0(X2S*x):

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


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

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

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

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

kt_fx := (rs, z, xs) ->
   1 - params_a_gamma/X_FACTOR_C*n_spin(rs, z)^(4/3)*xs^2/(n_spin(rs, z)^(4/3) + params_a_delta):

(* we want energy per particle *)
f := (rs, zeta, xt, xs0, xs1) -> gga_exchange_nsp(kt_fx, rs, zeta, xs0, xs1):


text:
# Display contents of the packages since
# we might use them latter
with(LinearAlgebra):
with(Groebner):

writeto("output.txt"):

#F := [x^3 + 12, 2*y^2 - x]:
#a, b := Basis(F, plex(x, y), output=extended):

#simplify(Multiply(convert(b, Matrix), Transpose(convert(F, Matrix)))):
#simplify(Transpose(convert(a, Matrix))):

#_m_order_1 = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]:
#order_1 = 'matrix'(_m_order_1, [x, y, z]):

#f_example_monomial_order := proc (poly, order_1::{MonomialOrder, ShortMonomialOrder})
#return sort([op(poly)], (b, a) -> TestOrder(a, b, order_1))
#end proc:

#F := 1/2*x + 2*y*z:

#print(f_example_monomial_order(F, 'matrix'([[1,0,0],[0,1,0],[0,0,1]], [x, y, z])        )):
#print(f_example_monomial_order(F, plex(x, y, z)   )):
#print(f_example_monomial_order(F, grlex(x, y, z)  )):

# ---------------------------------------------------------------------------
truncatePolynomial := proc (poly, order_1, order_2)
local leading_coeff_1, leading_mon_1, leading_term, curr_index, polys:

    leading_coeff_1, leading_mon_1 := LeadingTerm(poly, order_1):
    leading_term                   := leading_coeff_1 * leading_mon_1:
    curr_index                     := 1:

    polys := sort([op(poly)], (b, a) -> TestOrder(a, b, order_2)):

    while evalb(polys[curr_index] <> leading_term) do
        curr_index := curr_index + 1:
    end do:

    if curr_index = 1 then
        return polys[1]:
    else
        return add(x, x in polys[1..curr_index]):
    end if:
end proc:
# ---------------------------------------------------------------------------

# ---------------------------------------------------------------------------
# Test truncatePolynomial
print( truncatePolynomial(x^2 - y*z + 1, grlex(x, y, z), plex(x, y, z)) ):
print( truncatePolynomial(y^2 - x, grlex(x, y, z), plex(x, y, z)) ):
# ---------------------------------------------------------------------------

basisConversion := proc (basis, order_1, order_2)
local F, F_t, num_iter, F_t_gb, multipliers, G, repeat, curr_index, num_elements:

    print("Step 1"):
    F   := Basis(basis, order_1):
    F_t := map(v -> truncatePolynomial(v, order_1, order_2), F):
    print("Current F", F):
    print("Current F_t", F_t):

    num_iter := 1:

    while true do
        print("Iteration ", num_iter):
        num_iter := num_iter + 1:

        print("Step 3 & 4"):
        F_t_gb, multipliers := Basis(F_t, order_2, output=extended):
        print("Matrix of multipliers M': ", multipliers):
        print("H (= M' F_t): ", F_t_gb):

        print("Step 5"):
        G := convert(simplify(Multiply(convert(multipliers, Matrix), Transpose(convert(F, Matrix)))), list):
        print("G (= M' * F): ", G):

        print("Step 6 and 7"):
        repeat := false:
        curr_index := 1:
        num_elements := numelems(F_t_gb):
        while curr_index < num_elements do
            if evalb(LeadingMonomial(F_t_gb[curr_index], order_2) <> LeadingMonomial(G[curr_index], order_2)) then
                print("There are the witness polynomials that prove G_s is not empty (g_j, h_j) ", F_t_gb[curr_index], G[curr_index]):
                repeat := true:
                break:
            end if:
            curr_index := curr_index + 1:
        end do:

        if repeat = false then
            print("Done"):
            return InterReduce(G, order_2):
        else
            F := G:
            F_t := map(v -> truncatePolynomial(v, order_1, order_2), F):
        end if:
    end do:

end proc:

# ---------------------------------------------------------------------------
# Test basisConversion
print( basisConversion([y^2-x,x^2-y*z-1,z^2-x], grlex(x, y, z), plex(x, y, z)) ):
# ---------------------------------------------------------------------------

max_stars_repo_name: typesAreSpaces/Basis-Conversion

max_issues_repo_name: typesAreSpaces/Basis-Conversion

max_forks_repo_name: typesAreSpaces/Basis-Conversion


text:
# 求解的成立条件

$ifndef _CONDITION_
$define _CONDITION_

$include "Utils.mpl"

(*
* 求 解的 定义域
* 返回的条件中求解了单变量约束
* 如果无解则返回 x=undefined 的约束用作标记
* 保留了多变量约束
*)
findSolutionDomain:=proc(s)
    local con;
    con:=`union`(findDomain~(rhs~({s[]}))[]);
    return remove(x->type(x,`=`) and (lhs(x)=rhs(x)),classifySolve(con))
        union select(x->type(rhs(x),numeric),{s[]});
end proc:

(*
* 求解表达式的定义域
* 只考虑 + * ^ ln
*)
findDomain:=proc(ee)
    local S,r;
    S:={};
    findDomainCondtions(ee,S);
    return S;
end proc:

(*
* 求解表达式的定义域约束条件
* 只考虑 + * ^ ln
*)
findDomainCondtions:=proc(e,S::evaln(set))
    local _e;
    if type(e,`^`) and (not type(op(1,e),numeric)) then
        if (op(2,e)<0) then
            if type(op(2,e),fraction) and type(denom(op(2,e)),even) then
                S:=eval(S) union {op(1,e)>0};
            else
                S:=eval(S) union {op(1,e)<>0};
            end if;
        else
            if type(op(2,e),fraction) and type(denom(op(2,e)),even) then
                S:=eval(S) union {op(1,e)>=0};
            end if;
        end if;
    elif type(e,`+`) or type(e,`*`) then
        for _e in e do
            findDomainCondtions(_e,S);
        end do;
    elif (op(0,e)='ln') and (not type(op(1,e),numeric)) then
        S:=eval(S) union {op(e)>0};
    end if;
end proc:

# 不等式约束分类求解
# 合并单变量约束
# 保留多变量约束
classifySolve:=proc(con::set)
    local t,sd,ns,c,ind,x,sol;
    t:=table();
    ns,sd:=selectremove(x->(numelems(indets(x,name))=1),con);
    for c in ns do
        tappend(t,indets(c,name)[],c);
    end do;
    ind:=[indices(t,nolist)];
    for x in ind do
        sol:=RealDomain:-solve(t[x],{x});
        if (sol=NULL) then
            t[x]:={x=undefined};
        else
            t[x]:=sol;
        end if;
    end do;
    ns:=`union`(entries(t,nolist));
    return sd union ns;
end proc:

$endif


text:
// Delete all portals specified by the params and kill the params
impulse process ACV_remote_deletePortals {
/scoreboard players set @e[name=ACV_remote_deletePortals_PARAM,tag=ACV_Blue] ACV_COLOR 0
/scoreboard players set @e[name=ACV_remote_deletePortals_PARAM,tag=ACV_Orange] ACV_COLOR 1
/scoreboard players set @e[name=ACV_remote_deletePortals_PARAM,tag=ACV_Purple] ACV_COLOR 2
/scoreboard players set @e[name=ACV_remote_deletePortals_PARAM,tag=ACV_Red] ACV_COLOR 3
/entitydata @e[name=ACV_remote_deletePortals_PARAM] {CustomName:"ACV_deletePortals_PARAM"}
ACV_deletePortals()
}


text:
# A cactus lobe is a map f : R/Z -> R/Z of degree one that has
# constant speed on some finite collection of half-open intervals,
# and speed zero on the complement of those intervals.  We encode
# f as a pair [c,T], where c = f(0) is in R/Z (as defined by
# `is_element/RZ`) and T is the half-open interval where the
# speed is nonzero (as defined by `is_element/RZ_set`).

`is_element/cactus_lobes` := proc(cT)
 global reason;
 local c,T;

 if not(type(cT,list) and nops(cT) = 2) then
  reason := ["is_element/cactus_lobes","cT is not a list of length 2",cT];
  return false;
 fi;

 c,T := op(cT);

 if not `is_element/RZ`(c) then
  reason := ["is_element/cactus_lobes","c is not in R/Z",c,cT];
  return false;
 fi;

 if not `is_element/RZ_set`(T) then
  reason := ["is_element/cactus_lobes","T is not an RZ_set",T,cT];
  return false;
 fi;

 if T = `empty/RZ` then
  reason := ["is_element/cactus_lobes","T is empty",T,cT];
  return false;
 fi;

 return true;
end:

`is_equal/cactus_lobes` := proc(cT0,cT1)
 return `is_equal/RZ`(cT0[1],cT1[1]) and
        `is_equal/RZ_set`(cT0[2],cT1[2]);
end:

`random_element/cactus_lobes` := proc()
 local T;
 T := `empty/RZ_set`;
 while T = `empty/RZ_set` do T := `random_element/RZ_set`(); od:
 return [`random_element/RZ`(),T];
end:

`is_leq/cactus_lobes` := NULL:
`list_elements/cactus_lobes` := NULL:
`count_elements/cactus_lobes` := NULL:

`start/cactus_lobes` := (cT) -> cT[1];
`motion_set/cactus_lobes` := (cT) -> cT[2];
`speed/cactus_lobes` := (cT) -> 1/`length/RZ_set`(cT[2]);

`eval/cactus_lobes` := (cT) -> proc(t)
 local t0,c,T,n,i,u,s;

 c,T := op(cT);
 n := nops(T)/2;
 i := 1;
 u := c;
 s := 1/`length/RZ_set`(T);
 t0 := t - floor(t);
 while i <= n and T[2*i-1] <= t0 do
  u := u + s*(min(t0,T[2*i]) - T[2*i-1]);
  i := i+1;
 od:

 u := u - floor(u);
 return u;
end:

`id/cactus_lobes` := [0,[0,1]];

`o/cactus_lobes` := proc(cT,dU)
 local c,T,d,U,n,m,sT,sU,T1,T2,n1,V,a,b,y,z,i,j,aa,bb;
 c,T := op(cT);
 d,U := op(dU);
 n := nops(T)/2;
 m := nops(U)/2;
 sT := 1/`length/RZ_set`(T);
 sU := 1/`length/RZ_set`(U);
 T1 := `unroll/RZ_set`(T);
 n1 := nops(T1)/2;
 V := NULL;
 y := d;
 for i from 1 to m do
  a := U[2*i-1];
  b := U[2*i];
  z := y + sU * (b - a);
  T2 := (T1 -~ y)/~sU +~ a;
  for j from 1 to n1 do
   if T2[2*j] >= a and T2[2*j-1] <= b then
    aa := max(T2[2*j-1],a);
    bb := min(T2[2*j],b);
    if aa < bb then
     V := V,aa,bb;
    fi;
   fi;
  od;
  y := z;
 od;

 return [`eval/cactus_lobes`(cT)(d),[V]];
end:

`source_shift/cactus_lobes` := (t) -> proc(cT)
 local c,d,T,U;
 c,T := op(cT);
 U := `shift/RZ_set`(t)(T);
 d := `eval/cactus_lobes`(cT)(-t);
 return [d,U];
end:

`target_shift/cactus_lobes` := (t) -> proc(cT)
 local c,d,T;
 c,T := op(cT);
 d := c+t;
 d := d - floor(d);
 return [d,T];
end:

`plot/cactus_lobes` := proc(cT)
 local c,T,L,x,y,z,a,b,n,i,s,A;
 c,T := op(cT);
 L := NULL;
 if T[1] > 0 then
  L := L,[[0,c],[T[1],c]];
 fi;
 y := c;
 n := nops(T)/2;
 s := 1/`length/RZ_set`(T);
 for i from 1 to n do
  a := T[2*i-1];
  b := T[2*i];
  z := y + s * (b-a);
  if z <= 1 then
   L := L,[[a,y],[b,z]];
   y := z;
  else
   x := a + (1 - y)/s;
   L := L,[[a,y],[x,1]],[[x,0],[b,z-1]];
   y := z-1;
  fi;
  x := `if`(i < n,T[2*i+1],1);
  if b < x then
   L := L,[[b,y],[x,y]];
  fi;
 od;

 A := args[2..-1];
 return
  display(map(u -> line(u[1],u[2],A),[L]),scaling=constrained,axes=none);
end:

`full_plot/cactus_lobes` := proc(cT)
 local x;

 display(
  `plot/cactus_lobes`(cT),
  seq(line([x,0],[x,1],colour=red),x in [0,1]),
  seq(line([x,0],[x,1],colour=red,linestyle=dot),x in cT[2]),
  line([0,0],[1,0],colour=blue),
  line([0,1],[1,1],colour=blue),
  line([0,cT[1]],[1,cT[1]],colour=green,linestyle=dash)
 );
end:


text:
# Subspaces of (Z/p)^d are represented as pairs [U,V], where
# U and V are lists of lists, which correspond to RREF matrices.
# The rows of U are a basis for the subspace, and the rows
# of V are a basis for the annihilator.

`dim/subspaces`   := (p,d) -> U -> nops(U[1]);
`codim/subspaces` := (p,d) -> U -> nops(U[2]);

`span/subspaces` := (p,d) -> proc(U)
 local U0,r,c,i;

 if U = [] then return [[],convert(IdentityMatrix(d),listlist)]; fi;

 U0 := map(modp,U,p);
 if type(U0,list) then U0 := Matrix(U0); fi;
 r,c := Dimension(U0);
 if c <> d then error("Dimension mismatch"); fi;
 r := LinearAlgebra[Modular][MatBasis](p,U0,r,false);
 U0 := Matrix([op(1..r,convert(U0,listlist)),
               seq([0$c],i=1..c-r)]);
 LinearAlgebra[Modular][MatBasis](p,U0,r,true);
 U0 := convert(U0,listlist);
 return [[op(1..r,U0)],[op(r+1..c,U0)]];
end:

`ann/subspaces` := (p,d) -> proc(U)
 local U0;
 U0 := U;
 if type(U,list(list(integer))) then
  U0 := `span/subspaces`(p,d)(U0);
 fi;
 return `span/subspaces`(p,d)(U0[2]);
end:

`bot/subspaces`   := (p,d) -> [[],convert(IdentityMatrix(d),listlist)];

`top/subspaces`   := (p,d) -> [convert(IdentityMatrix(d),listlist),[]];

`sum/subspaces`   := (p,d) -> (U,V) -> `span/subspaces`(p,d)([op(U[1]),op(V[1])]);

`inter/subspaces` := (p,d) -> (U,V) -> `ann/subspaces` (p,d)([op(U[2]),op(V[2])]);

`is_member/subspaces` := (p,d) -> proc(u,V)
 if V[2] = [] then return true; fi;
 return is_all_zero(map(modp, Matrix(V[2]) . Vector(u),p));
end:

`is_le/subspaces` := (p,d) -> proc(U,V)
  if U[1] = [] or V[2] = [] then return true; fi;
  return is_all_zero(map(modp, Matrix(V[2]) . Transpose(Matrix(U[1])),p));
end:

`random_element/subspaces` := (p,d) -> proc()
  local r,U,i,j;
  r := rand(0..d)();
  U := [seq([seq(rand(0..p-1)(),j=1..d)],i=1..r)];
  return `span/subspaces`(p,d)(U);
end:


text:
$ifndef _TEQ_SOL_
$define _TEQ_SOL_

$include "Utils.mpl"

TeqSol:=module()
    option object;
    export
            rsol,           # 对应特解
            teq,            # 变换方程
            teqInd,
            tsols,          # 变换方程的解
            tcons,          # 解对应的条件

            printTsol::static,
            getNSols::static,
            getRep::static,
            ModuleApply::static,    # 构造函数
            hasSol::static,         # 判断是否有解
            getCon::static;         # 给出最简约束条件

    ModuleApply:=proc(spec,ESC)
        local s,k;
        s:=Object(TeqSol);
        s:-rsol:=spec;
        s:-teq:=[ESC[1][1],ESC[2][1]];
        s:-teqInd:=[1$numelems(ESC[1][2]),2$numelems(ESC[2][2])];
        s:-tsols:=[ESC[1][2][],ESC[2][2][]];
        s:-tcons:=[ESC[1][3][],ESC[2][3][]];
        return s;
    end proc:

    hasSol:=proc(s::TeqSol)
        return ormap(x->evalb(x<>[]),s:-tsols);
    end proc:

    printTsol:=proc(s::TeqSol)
        local i,n:=numelems(s:-tsols);
        printf("=======================================\n");
        printf("特解\n");
        print(s:-rsol);
        printf("变换方程的解\n");
        for i from 1 to n do
            if s:-teqInd[i]=1 then
                printf("正向\n");
            else
                printf("逆向\n");
            end if;
            print(s:-tsols[i]);
            print(s:-tcons[i]);
        end do;
        printf("=======================================\n");
        return;
    end proc:

    getNSols:=proc(s::TeqSol)
        return numelems(s:-tsols);
    end proc:

    getRep:=proc(s::TeqSol)
        return add(v[i]*s:-rsol[i],i=1..numelems(s:-rsol));
    end proc:

end module:

$endif


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

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

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

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

ms_f := (x, u, t) -> ms_f0(X2S^2*x^2, 0) + \
  ms_fa(ms_alpha(t,x))*(ms_f0(X2S^2*x^2, params_a_c) - ms_f0(X2S^2*x^2, 0)):

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


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

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

(* type: mgga_exc *)

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

rmggac_gamma1 := 0.08:
rmggac_gamma2 := 0.3:
rmggac_g := (alpha, s) ->
  (1 + rmggac_gamma1)*alpha/(rmggac_gamma1 + alpha + rmggac_gamma2*s^2):

rmggac_f2 := (alpha, s) ->
  3*rmggac_g(alpha, s)^3/(1 + rmggac_g(alpha, s)^3 + rmggac_g(alpha, s)^6):
rmggac_f1 := (alpha, s) ->
  1 - rmggac_f2(alpha, s):

rmggac_gamma := 0.031091:

(* from gga_c_regtpss *)
beta_a := 0.066724550603149220:
beta_b := 0.1:
beta_c := 0.1778:
mbeta := (rs, t) -> beta_a*(1 + beta_b*rs)/(1 + beta_c*rs):
(* from mgga_c_r2scan *)
rmggac_w1 := (rs, z) -> exp(-f_pw(rs, z)/(rmggac_gamma*mphi(z)^3)) - 1:
(* from gga_c_pbe *)
A := (rs, z, t) -> mbeta(rs, t)/(rmggac_gamma * rmggac_w1(rs, z)):
(* from gga_c_scan_e0 *)
scan_e0_g := (rs, z, t) -> (1 + 4*A(rs, z, t)*t^2)^(-1/4):
(* from mgga_c_r2scan *)
rmggac_H1 := (rs, z, t) -> rmggac_gamma*mphi(z)^3*log(1 + rmggac_w1(rs, z) * (1 - scan_e0_g(rs, z, t))):

(* from mgga_c_scan *)
scan_alpha := (z, xt, ts0, ts1) ->
  (t_total(z, ts0, ts1) - xt^2/8)/(K_FACTOR_C*t_total(z, 1, 1)):
scan_b1c := 0.0285764:
scan_b2c := 0.0889:
scan_b3c := 0.125541:
scan_eclda0 := rs -> -scan_b1c/(1 + scan_b2c*sqrt(rs) + scan_b3c*rs):

scan_chi_infty := 0.12802585262625815:
scan_g_infty := s -> 1/(1 + 4*scan_chi_infty*s^2)^(1/4):
scan_G_cnst := 2.3631:
scan_Gc := z -> (1 - scan_G_cnst*(2^(1/3) - 1)*f_zeta(z))*(1 - z^12):

scan_H0 := (rs, s) ->
  scan_b1c*log(1 + (exp(-scan_eclda0(rs)/scan_b1c) - 1)*(1 - scan_g_infty(s))):
scan_e0 := (rs, z, s) ->
  (scan_eclda0(rs) + scan_H0(rs, s))*scan_Gc(z):

(* define the functional *)
rmggac_eps1 := (rs, z, t) ->
  (f_pw(rs, z) +  rmggac_H1(rs, z, t)):

rmggac_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  + scan_e0(rs, z, X2S*2^(1/3)*xt)
      * rmggac_f1(scan_alpha(z, xt, ts0, ts1), X2S*2^(1/3)*xt)
  + rmggac_eps1(rs, z, tt(rs, z, xt))
      * rmggac_f2(scan_alpha(z, xt, ts0, ts1), X2S*2^(1/3)*xt):

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


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

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

# 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_L := (p,n,r) -> (k) -> p^(n*(k+r-1));

# 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 conjecturally 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)
 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:

# 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[i,k] and ee[i,k], with tridegrees as follows:

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

# (This gives all the odd generators, but not all the even ones.)

# The generators bb[i,k] and ee[i,k] can be expressed in terms of
# the original generators b[i] and e[i]. as follows:

tanabe_bbb := (p,n,r) -> (i,k) -> b[i+tanabe_N_bar(p,n,r)(k-1)]^(p^k);

tanabe_eee := (p,n,r) -> (i,k) ->
 u^(p^(n*(r+k))-1) * e[i] * mul(b[i+tanabe_N_bar(p,n,r)(j)]^(p^j*(p-1)),j=0..k-1);

tanabe_AHSS_d_rule := (p,n,r) -> proc(k,m)
 local R,Nb,i,j,n0,n1;
 R := NULL;
 Nb := tanabe_N_bar(p,n,r);
 for j from 0 to tanabe_s(p,n,r)(m) do
  n0 := `if`(j = 0,0,Nb(j-1));
  n1 := Nb(j);
  if j <= k then
   for i from n0 to min(n1-1,m) do
    R := R,b[i]^(p^j) = 0;
   od:
  else
   for i from n0 to min(n1-1,m) do
    R := R,b[i]^(p^k) = tanabe_eee(p,n,r)(i-Nb(k),k);
   od;
  fi;
 od;
 return {R};
end:

# 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) -> (k) ->
 (1-t)^(p^(n*r)) * mul((1 - t^(p^j))^(p^(n*(r+j))*(1-1/p^n)),j=1..k);

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_trideg0 := (p,n,r) -> proc(x)
 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);
 else
  return [0,0,0];
 fi;
end:

tanabe_trideg := (p,n,r) ->
 apply_deg(tanabe_trideg0(p,n,r));

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

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

 L      := tanabe_L(p,n,r);
 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(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(p^k*(N_bar(k) - N_bar(k-1)) - (p^n-1)*L(k)) = 0,
  "p^k*(N_bar(k) - N_bar(k-1)) = (p^n-1)*L(k)"
 );

 _ASSERT(
  simplify(N_bar(k-1) - L(k)/p^k - N_star(k)) = 0,
  "N_bar(k-1) = L(k)/p^k + P(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:


text:
program parseoutputformat;
begin
    writeln(char('M') : 10)
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


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

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

(* type: gga_exc *)

$define gga_x_pbe_params
$include "gga_x_pbe.mpl"
$include "gga_x_sfat.mpl"

ityh_enhancement := xs -> pbe_f(xs):


text:
// Param & Return Types:
// AS = armor_stand
// SCV = Scoreboard Value
// SC = Scoreboard
// ASwS = armor_stands with Score
// ASwT = armor_stands with Tag
// IF = ItemFrame

// ACV_DIRECTION:
// 5 = Floor Portal
// 4 = Ceiling Portal
// 3 = East Portal
// 2 = North Portal
// 1 = West Portal
// 0 = South Portal

// ACV_Color:
// 0 = Blue
// 1 = Orange
// 2 = Purple
// 3 = Red

Summon blue Portal:
/summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block

Summon orange Portal:
/summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Orange"],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block

Summon purple Portal:
/summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Purple"],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block

Summon red Portal:
/summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Red"],Duration:2147483647}
/execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block

Summon Repulsion Gel:
/summon item ~ ~ ~ {CustomName:"ACV_RepulsionGel",PickupDelay:32767,Item:{Count:1,Damage:11,id:"minecraft:wool"}}

Summon Propulsion Gel:
/summon item ~ ~ ~ {CustomName:"ACV_PropulsionGel",PickupDelay:32767,Item:{Count:1,Damage:1,id:"minecraft:wool"}}

Summon Cleansing Gel
/summon item ~ ~ ~ {CustomName:"ACV_CleansingGel",PickupDelay:32767,Item:{Count:1,id:"minecraft:glass"}}

Summon Material Emancipation Grill:
/summon area_effect_cloud ~ ~ ~ {CustomName:ACV_MEG,Duration:2147483647,Passengers:[{id:creeper,CustomName:ACV_MEG,powered:1,ActiveEffects:[{Amplifier:1,Duration:1000000,Id:14,ShowParticles:0b}],PersistenceRequired:1,Invulnerable:1,NoAI:1}]}

Kill Material Emancipation Grill:
/tp @e[name=ACV_MEG,dy=1] ~ -100 ~

Control_Panel {

Special_Thanks {
/setblock ~ ~ ~ standing_sign 0 replace {Text1:"{\"text\":\"Special Thanks\",\"bold\":\"true\"}",Text2:"[{\"text\":\"click on the \"},{\"text\":\"red\",\"color\":\"red\"}]",Text3:"{\"text\":\"names for more\"}",Text4:"{\"text\":\"information\"}"}
/setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"Adrodoc55\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"Adrodoc55's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" YouTube-Channel by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"https://www.youtube.com/user/Adrodoc55\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"{\"text\":\"because I wanted\"}",Text4:"{\"text\":\"to thank myself\"}"}
/setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"The Farlanders\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"The Farlander's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" Video by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"http://youtu.be/jSiP51XKK1Q\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"{\"text\":\"for the whole\"}",Text4:"{\"text\":\"Idea\"}"}
/setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"Dragnoz\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"Dragnoz's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" Video by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"http://youtu.be/J0ulXBtxQvo\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"{\"text\":\"for floor and\"}",Text4:"{\"text\":\"ceiling Portals\"}"}
/setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"SimplySarc\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"SimplySarc's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" Video by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"http://youtu.be/O8dv9P49cKk\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"{\"text\":\"for the\"}",Text4:"{\"text\":\"Chunk-Loader\"}"}
/setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"NanoRex\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"NanoRex's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" Video by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"http://youtu.be/lyROYbzLCZ0\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"{\"text\":\"for the Material\"}",Text4:"{\"text\":\"Emancipation Grid\"}"}
/setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"Ersatz_\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"Ersatz_'s\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" YouTube-Channel by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"https://www.youtube.com/channel/UClfq20CeYnZ2_Rm72VgUPXg\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"{\"text\":\"for his\"}",Text4:"{\"text\":\"Contribution\"}"}
Special_Thanks }

Equipment {
/setblock ~ ~ ~ standing_sign 0 replace {Text1:"{\"text\":\"\"}",Text2:"{\"text\":\"Equipment\",\"color\":\"dark_purple\",\"bold\":\"true\"}",Text3:"{\"text\":\"\"}",Text4:"{\"text\":\"\"}"}

/give @p minecraft:bow 1 0 {Unbreakable:1,ench:[{id:48,lvl:6},{id:51,lvl:1}],HideFlags:5,display:{Name:"Blue Portalgun",Lore:["Makes blue Portals!"]}}
/give @p minecraft:bow 1 0 {Unbreakable:1,ench:[{id:48,lvl:7},{id:51,lvl:1}],HideFlags:5,display:{Name:"Orange Portalgun",Lore:["Doesn't make blue Portals!"]}}
/give @p minecraft:bow 1 0 {Unbreakable:1,ench:[{id:48,lvl:8},{id:51,lvl:1}],HideFlags:5,display:{Name:"Purple Portalgun",Lore:["Makes purple Portals!"]}}
/give @p minecraft:bow 1 0 {Unbreakable:1,ench:[{id:48,lvl:9},{id:51,lvl:1}],HideFlags:5,display:{Name:"Red Portalgun",Lore:["Makes red Portals!"]}}
/give @p minecraft:arrow 1 0 {display:{Name:"Picking up Cubes",Lore:["Hold this to pick up the closest Cube"]}}
/replaceitem entity @p slot.armor.feet minecraft:chainmail_boots 1 0 {Unbreakable:1,ench:[{id:2,lvl:1}],HideFlags:5,display:{Name:"Long Fall Boots",Lore:["No more Falldamage!"]}}
Equipment }

Kill all Cubes {
/kill @e[name=ACV_Cube]
Kill all Cubes }

Kill closest Cubes {
/kill @e[name=ACV_Cube,c=1]
Kill closest Cubes }

Summon a Weighted Storage Cube {
/summon armor_stand ~ ~-0.8 ~2 {CustomName:"ACV_Cube",DisabledSlots:2035728,ArmorItems:[{},{},{},{Damage:6,id:"stone",Count:1}],NoGravity:1b,Invisible:1b,Invulnerable:1b}
Summon a Weighted Storage Cube }

Summon a Discourgement Redirection Cube {
/summon armor_stand ~ ~-0.8 ~2 {CustomName:"ACV_Cube",DisabledSlots:2035728,ArmorItems:[{},{},{},{id:"dispenser",Count:1}],Tags:["ACV_Redirecting"],NoGravity:1b,Invisible:1b,Invulnerable:1b}
Summon a Discourgement Redirection Cube }

Delete all Portals {
/execute @e[name=ACV_remote_deleteAllPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
Delete all Portals }

Delete blue Portal {
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
/execute @e[name=ACV_remote_deletePortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
Delete blue Portal }

Delete orange Portal {
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Orange"],Duration:2147483647}
/execute @e[name=ACV_remote_deletePortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
Delete orange Portal }

Delete purple Portal {
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Purple"],Duration:2147483647}
/execute @e[name=ACV_remote_deletePortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
Delete purple Portal }

Delete red Portal {
/summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Red"],Duration:2147483647}
/execute @e[name=ACV_remote_deletePortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
Delete red Portal }

Block_Option {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"{\"text\":\"\"}",Text2:"{\"text\":\"Block-Option\",\"color\":\"black\",\"bold\":\"true\"}",Text3:"{\"text\":\"\"}",Text4:"{\"text\":\"\"}"}

/execute @e[name=ACV_confirmBlockOption] ~ ~ ~ setblock ~ ~ ~ redstone_block
Block_Option }

Mapmaker_Helper {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"{\"text\":\"When building\"}",Text2:"{\"text\":\"a map please\"}",Text3:"{\"text\":\"consider the \"}",Text4:"{\"text\":\"Mapmaker Helper\",\"color\":\"dark_red\"}"}

/tellraw @a [{"text":"Download the ","color":"yellow"},{"text":"Mapmaker Helper","color":"gold"},{"text":" by clicking ","color":"yellow"},{"text":"here","color":"aqua","hoverEvent":{"action":"show_text","value":{"text":"Click to open the Download Page","color":"red"}},"clickEvent":{"action":"open_url","value":"https://github.com/Adrodoc55/MapMakerHelper/releases"}}]
Mapmaker_Helper }

Portalsystem {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"[{\"text\":\"Portal\",\"color\":\"dark_blue\",\"bold\":true},{\"text\":\"system\",\"color\":\"gold\",\"bold\":true}]",Text2:"{\"text\":\"off\",\"color\":\"dark_red\"}",Text3:"{\"text\":\"|\",\"color\":\"black\"}",Text4:"{\"text\":\"on\",\"color\":\"dark_green\"}"}

/execute @e[name=ACV_portals] ~ ~ ~ setblock ~ ~ ~ redstone_block
/execute @e[name=ACV_teleportation] ~ ~ ~ setblock ~ ~ ~ redstone_block
/tellraw @a [{"text":"Welcome to","color":"gold"},{"text":" ApertureCraft Vanilla\n","color":"dark_aqua"},{"text":"   Version: ","color":"gold"},{"text":"9.0.1\n","color":"dark_aqua","italic":"true"},{"text":"   made by ","color":"yellow","italic":"true"},{"text":"Adrodoc55\n","color":"aqua","hoverEvent":{"action":"show_text","value":{"text":"Click to open YouTube Channel","color":"red"}},"clickEvent":{"action":"open_url","value":"https://www.youtube.com/user/Adrodoc55"}},{"text":"\nVisit me on ","color":"yellow"},{"text":"[Github]","color":"aqua","hoverEvent":{"action":"show_text","value":{"text":"Click to open Github","color":"red"}},"clickEvent":{"action":"open_url","value":"https://github.com/Adrodoc55/ApertureCraftVanilla"}},{"text":" to report Bugs and check for Updates","color":"yellow"}]

/execute @e[name=ACV_portals] ~ ~ ~ setblock ~ ~ ~ stone
/execute @e[name=ACV_teleportation] ~ ~ ~ setblock ~ ~ ~ stone
Portalsystem }

Smooth_Teleport {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"{\"text\":\"Smooth Teleport\"}",Text2:"{\"text\":\"off\",\"color\":\"dark_red\"}",Text3:"{\"text\":\"|\",\"color\":\"black\"}",Text4:"{\"text\":\"on\",\"color\":\"dark_green\"}"}

/execute @e[name=ACV_teleportation] ~ ~ ~ setblock ~ ~ ~ stone
/entitydata @e[name=ACV_teleportation] {CustomName:"ACV_simpleTeleportation"}
/entitydata @e[name=ACV_smoothTeleportation] {CustomName:"ACV_teleportation"}

/execute @e[name=ACV_teleportation] ~ ~ ~ setblock ~ ~ ~ stone
/entitydata @e[name=ACV_teleportation] {CustomName:"ACV_smoothTeleportation"}
/entitydata @e[name=ACV_simpleTeleportation] {CustomName:"ACV_teleportation"}

/execute @e[name=ACV_portals] ~ ~ ~ detect ~ ~ ~ redstone_block * execute @e[name=ACV_teleportation] ~ ~ ~ setblock ~ ~ ~ redstone_block
Smooth_Teleport }

Cubes {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"{\"text\":\"Cubes\",\"color\":\"dark_gray\",\"bold\":true}",Text2:"{\"text\":\"off\",\"color\":\"dark_red\"}",Text3:"{\"text\":\"|\",\"color\":\"black\"}",Text4:"{\"text\":\"on\",\"color\":\"dark_green\"}"}

/execute @e[name=ACV_cubes] ~ ~ ~ setblock ~ ~ ~ redstone_block

/execute @e[name=ACV_cubes] ~ ~ ~ setblock ~ ~ ~ stone
Cubes }

Lasers {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"{\"text\":\"Lasers\",\"color\":\"dark_red\",\"bold\":true}",Text2:"{\"text\":\"off\",\"color\":\"dark_red\"}",Text3:"{\"text\":\"|\",\"color\":\"black\"}",Text4:"{\"text\":\"on\",\"color\":\"dark_green\"}"}

/execute @e[name=ACV_lasers] ~ ~ ~ setblock ~ ~ ~ redstone_block

/execute @e[name=ACV_lasers] ~ ~ ~ setblock ~ ~ ~ stone
Lasers }

Material Emancipation Grill {
/setblock ~ ~ ~ minecraft:standing_sign 8 replace {Text1:"{\"text\":\"Material\"}",Text2:"{\"text\":\"Emancipation Grill\"}",Text3:"{\"text\":\"off\",\"color\":\"dark_red\"}",Text4:"{\"text\":\"on\",\"color\":\"dark_green\"}"}

/execute @e[name=ACV_megs] ~ ~ ~ setblock ~ ~ ~ redstone_block

/execute @e[name=ACV_megs] ~ ~ ~ setblock ~ ~ ~ stone
Material Emancipation Grill }

Gels {
/setblock ~ ~ ~ standing_sign 8 replace {Text1:"{\"text\":\"Gels\",\"color\":\"blue\",\"bold\":true}",Text2:"{\"text\":\"off\",\"color\":\"dark_red\"}",Text3:"{\"text\":\"|\",\"color\":\"black\"}",Text4:"{\"text\":\"on\",\"color\":\"dark_green\"}"}

/execute @e[name=ACV_gels] ~ ~ ~ /setblock ~ ~ ~ redstone_block

/execute @e[name=ACV_gels] ~ ~ ~ /setblock ~ ~ ~ stone
Gels }

Control_Panel }


text:
StrongWolfe := module()

  description "StrongWolfe linesearch";

  # Module defined as a package (i.e.) collection of procedures
  option package, load = ModuleLoad, unload = ModuleUnLoad;

  export lines, linesearch, interpolate, zoom, esatta, AB_list, doplot, doplot2;

  local ModuleLoad,
        ModuleUnLoad,
        _debug,
        c_1_vec, c_2_vec, rho, max_iterations,
        Dphi_0, phi_0, A_list, B_list;

  uses LinearAlgebra;

  ModuleUnLoad := proc()
    printf("StrongWolfe unload\n");
    NULL;
  end proc:

  ModuleLoad := proc()
    printf("StrongWolfe load\n");
    _debug         := true;
    c_1_vec        := [0.1,0.01,1e-4]; # Armijo condition
    c_2_vec        := [0.1,0.5,0.9];   # second (strong) Wolfe condition
    rho            := 2;  # bracket growth
    max_iterations := [5,5,10];
    NULL;
  end proc;

  # Explicitly call ModuleLoad here so the type is registered when this
  # code is cut&pasted in.  ModuleLoad gets called when the module is
  # read from the repository, but the code used to define a module
  # (like the command below) is not called at library read time.
  ModuleLoad();

  AB_list := proc()
    A_list, B_list;
  end proc;

  lines := proc( f0, df0, L )
    [ [0,f0], [L,f0+L*rho*df0] ], [ [0,f0], [L,f0+L*sigma*df0] ];
  end proc;

  # a_lo = a_{i - 1}
  # a_hi = a_{i}
  interpolate := proc ( a_lo, a_hi, phi_lo, phi_hi, Dphi_lo, Dphi_hi )
    local d1, d2;
    d1 := Dphi_lo + Dphi_hi - 3 * (phi_lo - phi_hi) / (a_lo - a_hi);
    d2 := sqrt(d1 * d1 - Dphi_lo * Dphi_hi);
    return a_hi - (a_hi - a_lo) * ((Dphi_hi + d2 - d1) / (Dphi_hi - Dphi_lo + 2 * d2))
  end:

  zoom := proc( a_lo_in, a_hi_in, phi, Dphi )
    local kkk, iteration,
          c1_Dphi_0, c2_Dphi_0,
          a_lo, a_hi, a_j,
          phi_lo, phi_hi, phi_j,
          Dphi_lo, Dphi_hi, Dphi_j;

    a_lo := a_lo_in; phi_lo := phi(a_lo); Dphi_lo := Dphi(a_lo);
    a_hi := a_hi_in; phi_hi := phi(a_hi); Dphi_hi := Dphi(a_hi);

    # Shrink bracket
    for kkk from 1 to 3 do
      c1_Dphi_0 := c_1_vec[kkk]*Dphi_0;
      c2_Dphi_0 := c_2_vec[kkk]*Dphi_0;
      for iteration from 1 to max_iterations[kkk] do

        A_list := [ op(A_list), a_lo ];
        B_list := [ op(B_list), a_hi ];

        # Interpolate a_j
        if a_lo < a_hi then
          a_j := interpolate( a_lo, a_hi, phi_lo, phi_hi, Dphi_lo, Dphi_hi );
        else
          # TODO: Check if this is needed
          a_j := interpolate( a_hi, a_lo, phi_hi, phi_lo, Dphi_hi, Dphi_lo );
        end;

        # Evaluate phi(a_j)
        phi_j  := phi(a_j);
        Dphi_j := Dphi(a_j);

        # Check Armijo
        if phi_j > phi_0 + a_j * c1_Dphi_0 or phi_j > phi_lo then
          a_hi    := a_j;
          phi_hi  := phi_j;
          Dphi_hi := Dphi_j;
        else

          if abs(Dphi_j) <= -c2_Dphi_0 then
            return a_j;
          end;

          if Dphi_j * (a_hi - a_lo) >= 0 then
            a_hi    := a_lo;
            phi_hi  := phi_lo;
            Dphi_hi := Dphi_lo;
          end;

          a_lo    := a_j;
          phi_lo  := phi_j;
          Dphi_lo := Dphi_j;
        end;
      end;
    end;

    # Quasi-error response
    return a_j;
  end;

  # `StrongWolfe`: This linesearch algorithm guarantees that
  # the step length satisfies the (strong) Wolfe conditions.
  # See Nocedal and Wright - Algorithms 3.5 and 3.6
  # This algorithm is mostly of theoretical interest,
  # users should most likely
  # `MoreThuente`, `HagerZhang` or `BackTracking`.

  linesearch := proc( phi, Dphi, alpha0 )
    local kkk, iteration, c1_Dphi_0, c2_Dphi_0,
          i, a_0, a_i, a_im1, phi_i, Dphi_i, phi_im1, a_star;

    phi_0  := phi(0);
    Dphi_0 := Dphi(0);

    A_list := [];
    B_list := [];

    # Step-sizes
    a_0     := 0;
    a_im1   := a_0;
    a_i     := alpha0;
    phi_im1 := phi_0;

    for kkk from 1 to 3 do
      c1_Dphi_0 := c_1_vec[kkk]*Dphi_0;
      c2_Dphi_0 := c_2_vec[kkk]*Dphi_0;
      for iteration from 1 to max_iterations[kkk] do

        phi_i := phi(a_i);

        # Test Wolfe conditions
        if (phi_i > phi_0 + a_i * c1_Dphi_0) or
           (phi_i >= phi_im1 and (iteration > 1 or kkk > 1) ) then
          a_star := zoom( a_im1, a_i, phi, Dphi );
          return a_star, phi(a_star);
        end;

        Dphi_i := Dphi(a_i);

        # Check condition 2
        if abs(Dphi_i) <= -c2_Dphi_0 then
          return a_i, phi_i;
        end;

        # Check condition 3
        if Dphi_i >= 0 then # FIXME untested!
          a_star := zoom( a_i, a_im1, phi, Dphi );
          return a_star, phi(a_star);
        end;

        A_list := [ op(A_list), a_im1 ];
        B_list := [ op(B_list), a_i ];

        # Choose a_iplus1 from the interval (a_i, a_max)
        a_im1 := a_i;
        a_i   := a_i * rho;

        # Update phi_im1
        phi_im1 := phi_i;
      end;
    end;

    # Quasi-error response TODO make this error instead
    printf("\n\nlinesearch failed!\n\n");
    return a_im1, phi(a_im1);
  end:

  doplot := proc( f_in, df_in, alpha, amax)
    local f0, df0, x, AA, BB, CC, DD;
    f0  := evalf(f_in(0));
    df0 := evalf(df_in(0));

    AA := plot(f_in(x),x=0..amax);
    BB := plot([[0,f0],[amax,f0+amax*rho*df0]],color="LimeGreen");
    CC := plot([[0,f0],[amax,f0+amax*sigma*df0]],color="blue");
    DD := plot([[alpha,f_in(alpha)]],color="red",style = point);
    display(AA,BB,CC,DD);
  end;

  doplot2 := proc( f_in, df_in, alpha, amax, A_list, B_list )
    local i, f0, df0, x, AA, BB, CC, DD, EE;
    f0  := evalf(f_in(0));
    df0 := evalf(df_in(0));

    AA := plot(f_in(x),x=0..amax);
    BB := plot([[0,f0],[amax,f0+amax*rho*df0]],color="LimeGreen");
    CC := plot([[0,f0],[amax,f0+amax*sigma*df0]],color="blue");
    DD := plot([[alpha,f_in(alpha)]],color="red",style = point);
    EE := plot([seq([[A_list[i],-i/5],[B_list[i],-i/5]],i=1..nops(A_list))],color="black");
    display(AA,BB,CC,DD,EE);
  end;

end module:

max_stars_repo_name: ebertolazzi/NLtoolbox

max_issues_repo_name: ebertolazzi/NLtoolbox

max_forks_repo_name: ebertolazzi/NLtoolbox


text:
program primefactor;	{ program18 }
var factor : array[200] of integer;
	i, n, can : integer;
begin
	i := 0;
	while i < 200 do begin
		factor[i] := 0;
		i := i + 1
	end;
	writeln('Input positive integer');
	readln(n);
	can := 2;
	while (can <= n div 2) and (can < 200) do begin
		if (n - (n div can) * can) <> 0 then can := can + 1
		else begin
			factor[can] := factor[can] + 1;
			n := n div can
		end
	end;
	i := 0;
	if n < 200 then begin
		factor[n] := factor[n]+1;
		n := 1
	end;
	while i < 200 do begin
		if factor[i] <> 0 then begin
			writeln('    ', i, ' ** ', factor[i]);
		end;
		i := i + 1
	end;
	if n > 1 then writeln('    ', n, ' ** ', 1)
end.

max_issues_repo_name: shouth/LanguageProcessing


text:
# Berechnungen für Systeme mit impliziten kinematischen Zwangsbedingungen
# Beschreibung
# Allgemeine Berechnungen zu den impliziten Zwangsbedingungen
# Die Bedingungen müssen vorher in einer Roboterspezifischen Datei berechnet werden (händische Erstellung der kinematischen Schleifen und sonstigen Zwangsbedingungen
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# kinematic_constraints_calculations_implicit -> Berechnungen bezüglich implizit definierter kinematischer Zwangsbedingungen
# Quellen
# [Docquier2013] Docquier, Nicolas and Poncelet, Antoine and Fisette, Paul: ROBOTRAN: a powerful symbolic gnerator of multibody models (2013)
# [DoThanhKotHeiOrt2009b] Do Thanh et al.: On the inverse dynamics problem of general parallel robots (2009)
# [ParkChoPlo1999] Park, FC and Choi, Jihyeon and Ploen, SR: Symbolic formulation of closed chain dynamics in independent coordinates
# Autor
# Moritz Schappler, schappler@imes.uni-hannover.de, 2018-02
# Institut fuer Mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
interface(rtablesize=30):
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 2: # Hoher Optimierungsgrad.
;
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
# Definition und Zwangsbedingungen
read "../robot_codegen_definitions/robot_env_IC":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name_OL):
read "../robot_codegen_definitions/robot_env_IC":
# Term-Vereinfachungen: 0=keine; 1=fast alle; 2=auch Zeitableitung der Projektionsmatrix
if not assigned(simplify_options) or simplify_options(1)=-1 then
  use_simplify := 1: # standardmäßig simplify-Befehle anwenden.
else
  use_simplify := simplify_options(1): # erster Eintrag ist für Zwangsbedingungen
end if:

# Lesen der Zwangsbedingungen
kin_constraints_exist := false:
constrfile := sprintf("../codeexport/%s/tmp/kinematic_constraints_implicit_maple.m", robot_name):
if FileTools[Exists](constrfile) then
  read constrfile:
  implconstr_s := implconstr_s: # Für Sichtbarkeit der Variablen
  implconstr_t := implconstr_t: # Für Sichtbarkeit der Variablen
  kin_constraints_exist := kin_constraints_exist: # Für Sichtbarkeit der Variablen
end if:

if kin_constraints_exist = true then:
  printf("Kinematische Zwangsbedingungen in impliziter Form gelesen.\n"):
else
  printf("Es gibt keine impliziten Zwangsbedingungen. Offene Struktur oder nur explizit definierte Bedingungen. Keine weiteren Berechnungen notwendig."):
  restart: # Funktioniert nicht. Gedacht, damit im worksheet-Modus nichts mehr passieren kann.
  robot_name := "": # Damit werden die Export-Befehle ungültig
  codegen_act := false:
  quit:
end if:
NIZB := RowDimension(implconstr_s):
# Gelenkdefinitionen
# Setze aktuierte Gelenke als Minimalkoordinaten voraus
NAJ := add(mu(k), k=1..NJ):
NPJ := NQJ-NAJ:
# Bestimme Indizes der aktiven und passiven Gelenke
# Entspricht Partitionierung der Gelenkwinkel in q1 und q2 in [Docquier2013]. Hier kann die Reihenfolge von aktiven und passiven Winkeln auch durchmischt sein.
IndAct := Matrix(NAJ, 1):
IndPass:= Matrix(NPJ, 1):
ka := 1: kp := 1:
for i from 1 to NQJ do
  if mu(i) = 1 then
    IndAct(ka) := i:
    ka := ka + 1:
  else
    IndPass(kp) := i:
    kp := kp + 1:
  end if:
end do:
printf("Indizes der %d aktiven Gelenke:\n", NAJ);
Transpose(IndAct);
printf("Indizes der %d passiven Gelenke:\n", NPJ);
Transpose(IndPass);
printf("%d Zwangsbedingungsgleichungen\n", NIZB);
# Positionsvektor der Minimalkoordinaten zur Umrechnung OL zu TE/DE
# Zum Vergleich des Ergebnisses aus den expliziten Zwangsbedingungen und den des offenen Systems
posNQJ := Matrix(NJ,1):
for i from 1 to NJ do
  if SearchText("qJ",convert(theta(i),string)) = 1 or SearchText("qJ",convert(d(i),string)) = 1 then
    posNQJ(i) := 1:
  end if:
end do:
# Jacobi-Matrix der Impliziten Zwangsbedingungen in Abhängigkeit der unabhängigen Koordinaten
#
# Entspricht J1 in [Docquier2013], A in [ParkChoPlo1999]
Phia_s := Matrix(NIZB, NAJ):
# Jacobi-Matrix berechnen
for i from 1 to NIZB do
  for j from 1 to NAJ do
    k := IndAct(j):
    Phia_s(i,j) := diff(implconstr_s(i,1), qJ_s(k,1)):
  end do:
end do:
# Jacobi-Matrix der Impliziten Zwangsbedingungen in Abhängigkeit der abhängigen Koordinaten
#
# Entspricht J2 in [Docquier2013], P in [ParkChoPlo1999]
Phip_s := Matrix(NIZB, NPJ):
# Jacobi-Matrix berechnen
for i from 1 to NIZB do
  for j from 1 to NPJ do
    k := IndPass(j):
    Phip_s(i,j) := diff(implconstr_s(i,1), qJ_s(k,1)):
  end do:
end do:
# Vereinfache die Ausdrücke
#
if use_simplify >= 1 then
  tmp_t11 := time():
  tmp_l11 := length(Phia_s):
  printf("%s: Beginne Term-Vereinfachung der %dx%d Zwangsbedingungs-Jacobi für aktive Gelenke. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(Phia_s), ColumnDimension(Phia_s), tmp_l11):
  Phia_s := simplify2(Phia_s):
  tmp_l21 := length(Phia_s):
  tmp_t21 := time():
  tmp_l12 := length(Phip_s):
  printf("%s: Beginne Term-Vereinfachung der %dx%d Zwangsbedingungs-Jacobi für passive Gelenke. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(Phip_s), ColumnDimension(Phip_s), tmp_l12):
  Phip_s := simplify2(Phip_s):
  tmp_l22 := length(Phip_s):
  tmp_t22 := time():
  printf("%s: Terme für Zwangsbedingungs-Jacobis vereinfacht. Länge (aktiv/passiv): %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_t21):
end if:
# Zeitableitung der Jacobi-Matrix
# Wird benötigt, um die Beschleunigung der Gelenkkoordinaten der offenen Struktur zu berechnen.
Phia_t := convert_s_t(Phia_s):
PhiaD_t := diff~(Phia_t, t):
PhiaD_s := convert_t_s(PhiaD_t):

Phip_t := convert_s_t(Phip_s):
PhipD_t := diff~(Phip_t, t):
PhipD_s := convert_t_s(PhipD_t):

# Vereinfache die Ausdrücke
#

if use_simplify >= 1 then
  tmp_t11 := time():
  tmp_l11 := length(PhiaD_s):
  printf("%s: Beginne Term-Vereinfachung der Zeitableitung der %dx%d Zwangsbedingungs-Jacobi für aktive Gelenke. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(PhiaD_s), ColumnDimension(PhiaD_s), tmp_l11):
  PhiaD_s := simplify2(PhiaD_s):
  tmp_l21 := length(PhiaD_s):
  tmp_t21 := time():
  tmp_l12 := length(PhipD_s):
  printf("%s: Beginne Term-Vereinfachung der Zeitableitung der %dx%d Zwangsbedingungs-Jacobi für passive Gelenke. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(PhipD_s), ColumnDimension(PhipD_s), tmp_l12):
  PhipD_s := simplify2(PhipD_s):
  tmp_l22 := length(PhipD_s):
  tmp_t22 := time():
  printf("%s: Terme für Zwangsbedingungs-Jacobis (Zeitableitung) vereinfacht. Länge (aktiv/passiv): %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_t21):
end if:
# Exportiere Code für folgende Skripte
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save Phia_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_active_jacobian_maple.m", robot_name):
save PhiaD_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_active_jacobian_time_derivative_maple.m", robot_name):
save Phip_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_passive_jacobian_maple.m", robot_name):
save PhipD_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_passive_jacobian_time_derivative_maple.m", robot_name):
save posNQJ, IndAct, IndPass, sprintf("../codeexport/%s/tmp/positionVector_NQJ_maple.m", robot_name):
printf("Ausdrücke für Kinematische ZB gespeichert (Maple)\n"):
if codegen_act then
  MatlabExport(implconstr_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_matlab.m", robot_name), codegen_opt):
  MatlabExport(Phia_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_active_jacobian_matlab.m", robot_name), codegen_opt):
  MatlabExport(PhiaD_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_active_jacobianD_matlab.m", robot_name), codegen_opt):
  MatlabExport(Phip_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_passive_jacobian_matlab.m", robot_name), codegen_opt):
  MatlabExport(PhipD_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_passive_jacobianD_matlab.m", robot_name), codegen_opt):
  MatlabExport(posNQJ, sprintf("../codeexport/%s/tmp/positionVector_NQJ_matlab.m", robot_name), codegen_opt):
  printf("Ausdrücke für Kinematische ZB gespeichert (Matlab)\n"):
end if:

# Invertierung der Gradientenmatrix
# Invertiere die Gradientenmatrix bezüglich der passiven Gelenke
n1 := RowDimension(Phip_s):
n2 := ColumnDimension(Phip_s):
# Erstelle Matrix ppf mit der Form der Gradientmatrix Phip_s
ppf := Matrix(n1,n2,symbol=uu):
for i from 1 to n1 do
  for j from 1 to n2 do
    if Phip_s(i,j) = 0 then
      ppf(i,j) := 0:
    end if:
  end do:
end do:
# Invertiere die allgemeine Form der Matrix und setze dann die Einträge ein.
InvPhip_sf := MatrixInverse(ppf):
InvPhip_s := InvPhip_sf:

for i from 1 to n1 do
  for j from 1 to n2 do
    if not Phip_s(i,j) = 0 then
      InvPhip_s := subs({ppf(i,j)=Phip_s(i,j)}, InvPhip_s):
    end if:
  end do:
end do:

# Vereinfache den Ausdruck
if use_simplify >= 1 then
  tmp_t1 := time(): tmp_l1 := length(InvPhip_s):
  printf("%s: Beginne Term-Vereinfachung der inversen %dx%d Zwangsbedingungs-Jacobi für passive Gelenke. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(InvPhip_s), ColumnDimension(InvPhip_s), tmp_l1):
  InvPhip_s := simplify2(InvPhip_s):
  tmp_t2 := time(): tmp_l2 := length(InvPhip_s):
  printf("%s. Inverse Passiv-Gelenk-Jacobi der ZB vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t11):
end if:
save InvPhip_s, sprintf("../codeexport/%s/tmp/kinematic_constraints_explicit_passive_jacobian_inv_maple", robot_name):
if codegen_act then
  MatlabExport(InvPhip_sf, sprintf("../codeexport/%s/tmp/kinconstr_impl_passive_jacobian_inv_matlab.m", robot_name), codegen_opt):
end if:

# Projektionsmatrix aus impliziter Gradientenmatrix
#
# Form der Projektionsmatrix anschauen
B21f := -InvPhip_sf.Phia_s:
# Projektionsmatrix berechnen. [Docquier2013], Text nach Gl. 12
B21 := -InvPhip_s.Phia_s:
# Vereinfache den Ausdruck
if use_simplify >= 1 then
  tmp_t1 := time(): tmp_l1 := length(B21):
  printf("%s: Beginne Term-Vereinfachung der %dx%d Zwangsbedingungs-Projektionsmatrix. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(B21), ColumnDimension(B21), tmp_l1):
  B21 := simplify2(B21):
  tmp_t2 := time(): tmp_l2 := length(B21):
  printf("%s. ZB-Projektionsmatrix vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t11):
end if:
save B21, sprintf("../codeexport/%s/tmp/kinconstr_impl_projection_jacobian_maple", robot_name):
if codegen_act then
  MatlabExport(B21, sprintf("../codeexport/%s/tmp/kinconstr_impl_projection_jacobian_matlab.m", robot_name), codegen_opt):
end if:

# Zeitableitung der Projektionsmatrix
# [ParkChoPlo1999], Text nach Gl. 73
# [DoThanhKotHeiOrt2009b], Gl. 20
B21D := (InvPhip_s . PhipD_s . InvPhip_s . Phia_s ) + (-InvPhip_s.PhiaD_s):
# Vereinfache den Ausdruck
if use_simplify >= 2 then
  tmp_t1 := time(): tmp_l1 := length(B21D):
  printf("%s: Beginne Term-Vereinfachung der Zeitableitung der %dx%d Zwangsbedingungs-Projektionsmatrix. Länge: %d.\n", \
     FormatTime("%Y-%m-%d %H:%M:%S"), RowDimension(B21D), ColumnDimension(B21D), tmp_l1):
  B21D := simplify2(B21D):
  tmp_t2 := time(): tmp_l2 := length(B21D):
  printf("%s. ZB-Projektionsmatrix-Zeitableitung vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t11):
end if:
save B21D, sprintf("../codeexport/%s/tmp/kinconstr_impl_projection_jacobian_derivative_maple", robot_name):
if codegen_act then
  MatlabExport(B21D, sprintf("../codeexport/%s/tmp/kinconstr_impl_projection_jacobian_derivative_matlab.m", robot_name), codegen_opt):
end if:

printf("Fertig\n"):

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

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

$define gga_x_pbe_params
$include "gga_x_pbe.mpl"

f := (rs, x, t, u) ->
  + params_a_csi_HF*f_pbe(x)*mgga_series_w(params_a_a, 12, t):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


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

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

(* type: work_gga_c *)
(* prefix:
  gga_x_kt_params *params;

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

fx := (rs, z, xs) ->
   params_a_gamma*n_spin(rs, z)^(4/3)*xs^2/(n_spin(rs, z)^(4/3) + params_a_delta):

f0 := (rs, zeta, xt, xs0, xs1) ->
  - (X_FACTOR_C - fx(rs,  zeta, xs0))*n_spin(rs,  zeta)^(4/3)
  - (X_FACTOR_C - fx(rs, -zeta, xs1))*n_spin(rs, -zeta)^(4/3):

(* we want energy per particle *)
f := (rs, zeta, xt, xs0, xs1) -> f0(rs, zeta, xt, xs0, xs1)/n_total(rs):