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

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

# Number of links = 6
# Synonyms of /Series/playbacknet/mrs_bool/mute =
# Number of synonyms = 1
# /Series/playbacknet/Gain/gt/mrs_bool/mute
# Synonyms of /Series/playbacknet/mrs_bool/hasData =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/mrs_bool/hasData
# Synonyms of /Series/playbacknet/mrs_natural/nChannels =
# Number of synonyms = 2
# /Series/playbacknet/SoundFileSource/src/mrs_natural/nChannels
# /Series/playbacknet/AudioSink/dest/mrs_natural/nChannels
# Synonyms of /Series/playbacknet/mrs_natural/pos =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/mrs_natural/pos
# Synonyms of /Series/playbacknet/mrs_real/israte =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/mrs_real/israte
# Synonyms of /Series/playbacknet/mrs_string/filename =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/mrs_string/filename

# nComponents = 3

# MarSystem
# Type = SoundFileSource
# Name = src

# MarControls = 17
# /SoundFileSource/src/mrs_bool/debug = 0
# /SoundFileSource/src/mrs_bool/mute = 0
# /SoundFileSource/src/mrs_bool/hasData = 0
# /SoundFileSource/src/mrs_bool/noteon = 0
# /SoundFileSource/src/mrs_natural/inObservations = 1
# /SoundFileSource/src/mrs_natural/inSamples = 512
# /SoundFileSource/src/mrs_natural/nChannels = 1
# /SoundFileSource/src/mrs_natural/onObservations = 1
# /SoundFileSource/src/mrs_natural/onSamples = 512
# /SoundFileSource/src/mrs_natural/pos = 0
# /SoundFileSource/src/mrs_natural/size = 661500
# /SoundFileSource/src/mrs_real/duration = -1
# /SoundFileSource/src/mrs_real/frequency = 0
# /SoundFileSource/src/mrs_real/israte = 22050
# /SoundFileSource/src/mrs_real/osrate = 22050
# /SoundFileSource/src/mrs_real/repetitions = 1
# /SoundFileSource/src/mrs_string/filename = /home/gtzan/data/sound/music_speech/music/gravity.au

# Number of links = 0

# MarSystem
# Type = Gain
# Name = gt

# MarControls = 9
# /Gain/gt/mrs_bool/debug = 0
# /Gain/gt/mrs_bool/mute = 0
# /Gain/gt/mrs_natural/inObservations = 1
# /Gain/gt/mrs_natural/inSamples = 512
# /Gain/gt/mrs_natural/onObservations = 1
# /Gain/gt/mrs_natural/onSamples = 512
# /Gain/gt/mrs_real/gain = 1
# /Gain/gt/mrs_real/israte = 22050
# /Gain/gt/mrs_real/osrate = 22050

# Number of links = 0

# MarSystem
# Type = AudioSink
# Name = dest

# MarControls = 10
# /AudioSink/dest/mrs_bool/debug = 0
# /AudioSink/dest/mrs_bool/init = 0
# /AudioSink/dest/mrs_bool/mute = 0
# /AudioSink/dest/mrs_natural/inObservations = 1
# /AudioSink/dest/mrs_natural/inSamples = 512
# /AudioSink/dest/mrs_natural/nChannels = 1
# /AudioSink/dest/mrs_natural/onObservations = 1
# /AudioSink/dest/mrs_natural/onSamples = 512
# /AudioSink/dest/mrs_real/israte = 22050
# /AudioSink/dest/mrs_real/osrate = 22050

# Number of links = 0

max_stars_repo_name: jaouahbi/VampPlugins

max_issues_repo_name: jaouahbi/VampPlugins

max_forks_repo_name: jaouahbi/VampPlugins


text:
program gcd2;

function gcd(n1: integer, n2: integer): integer;
begin
    if n2 > n1 then
    begin
        var new: integer;
        new := n2 % n1;
        n2 := new;
    end
    else
    begin
        var new: integer;
        new := n1 % n2;
        n1 := n2;
        n2 := new;
    end;

    if n2 > 0 then
    begin
        return gcd(n1, n2);
    end
    else
    begin
        return n1;
    end;
end;

begin
    var n1, n2: integer;
    read(n1);
    read(n2);
    writeln(gcd(n1, n2));
end.

ext: mpl

max_stars_repo_name: sebschmi/codegen-course

max_issues_repo_name: sebschmi/codegen-course

max_forks_repo_name: sebschmi/codegen-course


text:
use DocumentTools in
if parse(GetProperty("RunSIAN", value)) then

	SetProperty("printSolutions", enabled, true):
	SetProperty("bypass", enabled, true):
	SetProperty("p", enabled, true):
	SetProperty("params", enabled, true):
	SetProperty("replicas", enabled, true):
else
	SetProperty("bypass", enabled, false):
	SetProperty("printSolutions", enabled, false):
	SetProperty("p", enabled, false):
	SetProperty("params", enabled, false):
	SetProperty("replicas", enabled, false):
fi:
end use;

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
kernelopts(assertlevel=2): # be strict on all assertions while testing
with(Hakaru):
with(KB):
with(Summary):
with(NewSLO, simplify_factor_assuming):

module()
  local kb, b, mr, f, bkt;
  kb := empty;
  b, kb := genType(b, HInt(Bound(`>=`,0),Bound(`<=`,size(as)-1)), kb);
  mr, f := summarize(piecewise(b=idx(z,i), idx(t,i), 0),
                     kb, i, summary);
  CodeTools[Test](indets(mr, 'Add(anything)'), {Add(idx(t,i))},
                  label="Simple example from Summary.txt");
  f := simplify_factor_assuming(f, kb);
  bkt := bucket(mr, i=0..size(t)-1);
  CodeTools[Test](eval(f, summary=bkt),
                  sum(piecewise(b=idx(z,i), idx(t,i), 0), i=0..size(t)-1),
                  simplify);
end module:

module()
  local kb, zNew, b, mr, f, bkt;
  kb := empty;
  zNew, kb := genType(zNew, HInt(Bound(`>=`,0),Bound(`<=`,size(as)-1)), kb);
  b, kb := genType(b, HInt(Bound(`>=`,0),Bound(`<=`,size(as)-1)), kb);
  mr, f := summarize(piecewise(b=piecewise(i=docUpdate,zNew,idx(z,i)),
                               idx(t,i)),
                     kb, i, summary);
  CodeTools[Test](indets(mr, 'Add(anything)'), {Add(idx(t,i))},
                  label="First \"summate i\" in offshore under gmm_gibbs.hk");
  f := simplify_factor_assuming(f, kb);
  bkt := bucket(mr, i=0..size(t)-1);
  CodeTools[Test](eval(f, summary=bkt),
                  piecewise(b=zNew,
                            sum(piecewise(i=docUpdate, idx(t,i), 0),
                                i=0..size(t)-1),
                            0) +
                  sum(piecewise(And(b=idx(z,i), Not(i=docUpdate)),
                                idx(t,i),
                                0),
                      i=0..size(t)-1),
                  simplify);
end module:

module()
  local kb, docUpdate, zNew, k, i, mr, f, bkt;
  kb := empty;
  docUpdate, kb := genType(docUpdate, HInt(Bound(`>=`,0),Bound(`<=`,size(z)-1)), kb);
  zNew, kb := genType(zNew, HInt(Bound(`>=`,0),Bound(`<=`,size(topic_prior)-1)), kb);
  k, kb := genType(k, HInt(Bound(`>=`,0),Bound(`<=`,size(topic_prior)-1)), kb);
  i, kb := genType(i, HInt(Bound(`>=`,0),Bound(`<=`,size(word_prior)-1)), kb);
  mr, f := summarize(piecewise(idx(doc,j)=docUpdate,
                               piecewise(k=zNew,
                                         piecewise(i=idx(w,j),1,0),
                                         0),
                               0),
                     kb, j, summary);
  CodeTools[Test](indets(mr, 'Add(anything)'), {Add(1)},
                  label="First \"summate j\" in offshore under naive_bayes_gibbs.hk");
  f := simplify_factor_assuming(f, kb);
  bkt := bucket(mr, j=0..size(w)-1);
  CodeTools[Test](eval(f, summary=bkt),
                  piecewise(k=zNew,
                            sum(piecewise(And(i=idx(w,j),
                                              docUpdate=idx(doc,j)),
                                          1,
                                          0),
                                j=0..size(w)-1),
                            0),
                  simplify);
end module:

module()
  local kb, i18, mr, f, bkt;
  kb := empty;
  i18, kb := genType(i18, HInt(Bound(`>=`,0),Bound(`<=`,word_prior_size-1)), kb);
  mr, f := summarize(piecewise(docUpdate=idx(doc,i13),
                               piecewise(And(i=zNew8, i18=idx(w,i13)), 1, 0),
                               0),
                     kb, i13, summary);
  CodeTools[Test](indets(mr, 'Add(anything)'), {Add(1)},
                  label="Spenser's example from 2017-01-23 email");
  f := simplify_factor_assuming(f, kb);
  bkt := bucket(mr, i13=0..word_prior_size-1);
  CodeTools[Test](eval(f, summary=bkt),
                  piecewise(i=zNew8,
                            sum(piecewise(And(i18=idx(w,i13),
                                              docUpdate=idx(doc,i13)),
                                          1,
                                          0),
                                i13=0..word_prior_size-1), 0),
                  simplify);
end module:

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
# gpgcd-bezout-test 2nd variation
# GPGCD test routines for Bezout resultant
# Akira Terui, Boming Chi, October 2018

with (LinearAlgebra);
with (PolynomialTools);
with (ArrayTools);
with (VectorCalculus);

# The update function
bezout_update := proc (Btemplate, J, degH, F0, G0, V0, F1, G1, V1, alpha)
local i, j, dimF, dimG, dimV, dimVV1, F1poly, G1poly, B1, VV1, Q1, P1, R1, dim, BB, vlist, varlist, F0grad, F1grad, subslistB, subslistV, sublist, J1, Jrow, Jcol, A1, D1, DD1, DD2, VV2, F2, G2, V2, Normfg, st;
    # Jrow, Jcol, A1, R1, D1, f, k, B2, V2;

    # inputs:

    # B:            Bezout matrix (template)
    # J:             Jacobi matrix (template)
    # degH:     degree of the approximate GCD
    # F0, G0:    initial coefficient vector of F and G
    # F1, G1:    current coefficient vector of F and G
    # V0:          initial singular vector of the Bezout matrix
    # V1:          current singular vector of the Bezout matrix
    # alpha:      step width for an iteration

    # outputs:

    # F2, G2:    new (updated) coefficient vector of F and G
    # V2:          new (updated) singular vector of the Bezout matrix
    # Normfg:  the distance between F1, G1 and F2, G2

    # find dimension of initial coefficient vector( F0, G0 and V0)

    dimF := Dimension(F0);
    dimG := Dimension(G0);
    dimV := Dimension(V0);
    userinfo(2, bezout_update, `[dimF, dimG, dimV] =`,
             print([dimF, dimG, dimV]));

    # construct polynomials from coefficient vector( F1 and G1)

    st := time();
    F1poly := FromCoefficientVector(F1, x);
    G1poly := FromCoefficientVector(G1, x);
    st := time() - st;
    userinfo(1, bezout_update, `stpoly=`, print(st));
    userinfo(2, bezout_update, `F1poly=`, print(F1poly));
    userinfo(2, bezout_update, `G1poly=`, print(G1poly));

    # construct substitution table for B1 & J1

    st := time();
    sublist := [seq(f[j] = F1[j+1], j=0..(dimF-1)), seq(g[j] = G1[j+1], j=0..(dimG-1)), seq(v[j] = V1[j], j=1..dimV)];
    st := time() - st;
    userinfo(1, bezout_update, `stsublist=`, print(st));
    userinfo(2, bezout_update, `sublist = `, print(sublist, whattype(sublist)));

    # construct Bezout Matrix of F1 and G1 -> B1

    st := time();
    B1 := BezMatrix(F1poly, G1poly, x);
    userinfo(2, bezout_update, `B1 =`, print(B1));
    st := time() - st;
    userinfo(1, bezout_update, `stB1=`, print(st));

    # construct vector of objective function

    st := time();
    VV1 := <F1, G1, V1>;
    userinfo(2, bezout_update, `VV1 =`, print(VV1, whattype(VV1)));
    dimVV1 := dimF + dimG + dimV;
    st := time() - st;
    userinfo(1, bezout_update, `stVV1=`, print(st));

    # evaluate the constraint value -> Q1

    st := time();
    Q1 := B1 . V1[ 1 .. dimF-1 ];
    st := time() - st;
    userinfo(1, bezout_update, `stQ1=`, print(st));
    userinfo(2, bezout_update, `Q1 =`, print(Q1, whattype(Q1)));

    st := time();
    P1 := < (F1 - F0), (G1 - G0), Vector(dimV, datatype=double) >;
    st := time() - st;
    userinfo(1, bezout_update, `stP1=`, print(st));
    userinfo(2, bezout_update, `P1 =`, print(P1, whattype(P1)));

    # construct right-hand-side vector -> R1

    st := time();
    R1 := - <P1, Q1>;
    st := time() - st;
    userinfo(1, bezout_update, `stR1=`, print(st));
    userinfo(2, bezout_update, `R1 =`, print(R1, whattype(R1)));

    # construct Jacobian matrix -> J1

    st := time();
    J1 := subs(sublist, J);
    st := time() - st;
    userinfo(1, bezout_update, `stJ1=`, print(st));
    userinfo(2, bezout_update, `J1 = `, print(J1, whattype(J1)));

    # construct coefficient matrix for the linear system -> A1

    st := time();
    Jrow, Jcol := Dimension(J);
    A1 := Matrix(Jrow+Jcol);
    for i from 1 to Jcol do A1(i,i):=1 end do:
    for i from 1 to Jrow do for j from 1 to Jcol do A1(i+Jcol,j):=J1(i,j) end do end do:
    for i from 1 to Jcol do for j from 1 to Jrow do A1(i,j+Jcol):=-J1(j,i) end do end do:
    st := time() - st;
    userinfo(1, bezout_update, `stA1=`, print(st));
    userinfo(2, bezout_update, `A1 =`, print(A1));

    # solve the linear system: A1 . D1 == R1

    st := time();
    D1 := LinearSolve(A1, R1);
    st := time() - st;
    userinfo(1, bezout_update, `stD1=`, print(st));
    userinfo(2, bezout_update, `D1 =`, print(D1));

    DD1 := D1[1..dimVV1];
    userinfo(2, bezout_update, `DD1 =`, print(DD1));

    DD2 := D1[1..(dimF+dimG)];
    userinfo(2, bezout_update, `DD2 =`, print(DD2));

    st := time();
    VV2 := (VV1 + alpha * DD1)[..,1];
    st := time() - st;
    userinfo(1, bezout_update, `stVV2=`, print(st));
    userinfo(2, bezout_update, `VV2 =`, print(VV2, whattype(VV2)));

    F2 := VV2[1..dimF];
    G2 := VV2[(dimF + 1)..(dimF + dimG)];
    V2 := VV2[(dimF + dimG + 1)..(dimVV1)];
    Normfg := LinearAlgebra:-Norm(DD2,2);
    userinfo(2, bezout_update, `Normfg =`, print(Normfg));

    return F2, G2, V2, Normfg;

end proc:

# The initial function
bezout_gpgcd := proc(F, G, x, degH, stopcriterion, alpha)
local degF, degG, hh, Ftemplate, Gtemplate, Btemplate, V, BV, i, j, m, n, J, Fcoef, Gcoef, B, V0, VV0, v0, st, st0, st1, st2, ff, gg, DF, DG, DD, numofiteration, sublist, B0;

    # inputs:

    # F: the initial polynomial F
    # G: the initial polynomial G
    # degH: the degree of the approximate GCD
    # alpha: step width for an iteration

    # outputs:

    # ff, gg: the polynomials finally we get
    # hh: the approximate GCD finally we get
    # DF, DG: the perturbation of ff and gg
    # DD: the whole perturbation
    # numofiteration: number of iterations
    # B0: the Bezout Matrix of ff and gg

    # Construct the Bezout matrix of the initial polynomials

    st0 := time();
    degF := degree(F,x);
    degG := degree(G,x);
    # if degF < degG then degF, degG := degG, degF; F, G := G, F end if;
    userinfo(2, bezout_update, `degG =`, print(degG));
    Ftemplate := add(f[i]*x^i, i=0..degF);
    Gtemplate := add(g[i]*x^i, i=0..degG);
    userinfo(2, bezout_update, `Gtemplate =`, print(Gtemplate));    st := time();
    Btemplate := BezMatrix(Ftemplate, Gtemplate, x);

    # Construct the Jacobian matrix

    st := time() - st;
    userinfo(1, bezout_update, `stBtemplate=`, print(st));
    V := Vector(degF,symbol=v);
    BV := Vector(degF);
    for i from 1 to degF do BV[i] := add(Btemplate[i,j]*V[j], j=1..degF) end do;
    userinfo(2, bezout_update, `BV =`, print(BV));
    st := time();
    J := bezout_jacobian(degF, degG, Btemplate);
    st := time() - st;
    userinfo(1, bezout_update, `stJtemplate=`, print(st));
    userinfo(1, bezout_update, `J =`, print(J));

    # Construct coefficient vectors

    Fcoef := PolynomialTools:-CoefficientVector(F, x);
    Gcoef := PolynomialTools:-CoefficientVector(G, x);

    # Construct singular vector of the d-th from the right side

    B := BezMatrix(F, G, x);
    V0 := LinearAlgebra:-Transpose(LinearAlgebra:-SingularValues(B, output='Vt'));
    v0 := V0[..,-degH];
    userinfo(2, bezout_update, `v0 =`, print(v0));
    st1 := time() - st0;
    userinfo(1, bezout_update, `stbeforeiteration=`, print(st1));
    ff, gg, hh, DF, DG, DD, numofiteration, B0 := bezout_iteration(Btemplate, J, degF, degG, degH, Fcoef, Gcoef, v0, Fcoef, Gcoef, v0, stopcriterion, 1.0, 0);
    st2 := time() - st0;
    userinfo(1, bezout_update, `stwholeiteration=`, print(st2));

    return ff, gg, hh, DF, DG, DD, numofiteration, B0;
    # return B, J, V0, v0
end proc:

# The Jacobian Matrix
bezout_jacobian := proc (m, n, Btemplate)
local J, i, j;

    # inputs:
    # degF, degG： the degree of polynomials F and G

    J := Matrix(m, 2*m + n + 2);

    for i from 1 to m do for j from i+1 to m+1 do J(i,j):= -g[i-1]*v[j-1] end do end do;
    for i from 2 to m do for j from i+1 to m+1 do J(i,j):= J(i,j)+J(i-1,j-1) end do end do;

    for j from 1 to n do for i from j to n do J(i,j):= g[i]*v[j] end do end do;
    for j from n-1 by -1 to 1 do for i from n-1 by -1 to j do J(i,j):= J(i,j)+J(i+1,j+1) end do end do;

    for i from 1 to n+1 do for j from i+m+2 to m+n+2 do J(i,j):= f[i-1]*v[j-m-2] end do end do;
    for i from 2 to n+1 do for j from i+m+2 to m+n+2 do J(i,j):= J(i,j)+J(i-1,j-1) end do end do;

    for j from m+2 to 2*m+1 do for i from j-m-1 to m do J(i,j):= -f[i]*v[j-m-1] end do end do;
    for j from 2*m by -1 to m+2 do for i from m-1 by -1 to j-m-1 do J(i,j):= J(i,j)+J(i+1,j+1) end do end do;

    for j from m+n+3 to 2*m+n+2 do for i from 1 to m do J(i,j):= Btemplate(i,j-m-n-2) end do end do;
    return J
end proc:

# The iteration function
bezout_iteration := proc (Btemplate, J, degF, degG, degH, F0, G0, V0, F1, G1, V1, stopcriterion, alpha, numofiteration)
local FF, GG, VV, Normfg, B, ff, gg, f0, g0, DF, DG, DD, hh, st, numofiteration0, sublist;

    # inputs:

    # Btemplate:       the template of the Bezout Matrix
    # J:                      the template of the Jacobian Matrix
    # degF, degG:     the degree of the initial polynomials
    # degH:               the degree of the approximate GCD
    # F0, G0:              initial coefficient vector of F and G
    # F1, G1:              current coefficient vector of F and G
    # V0:                    initial singular vector of the Bezout matrix
    # V1:                    current singular vector of the Bezout matrix
    # alpha:                step width for an iteration
    # stopcriterion:     the stop criterion
    # numofiteration: number of current iteration( initial number of iteration: 0)

    # outputs:

    # ff, gg: the polynomials finally we get
    # hh: the approximate GCD finally we get
    # DF, DG: the perturbation of ff and gg
    # DD: the whole perturbation
    # numofiteration: number of current iteration
    # B0: the Bezout Matrix of ff and gg

    FF, GG, VV, Normfg := bezout_update(Btemplate, J, degH, F0, G0, V0, F1, G1, V1, alpha);
    if Normfg < stopcriterion then
    userinfo(2, bezout_update, `FF=`, print(FF));
    st := time();
    f0 := FromCoefficientVector(F0,x);
    g0 := FromCoefficientVector(G0,x);
    ff := FromCoefficientVector(FF,x);
    gg := FromCoefficientVector(GG,x);
    DF := Norm(PolynomialTools:-CoefficientVector(ff-f0,x),2);
    DG := Norm(PolynomialTools:-CoefficientVector(gg-g0,x),2);
    DD := sqrt(DF^2+DG^2);
    sublist := [seq(f[j] = FF[j+1], j=0..degF), seq(g[j] = GG[j+1], j=0..degG)];
    B := BezMatrix(ff, gg, x);
    userinfo(2, bezout_update, `B=`, print(B));
    hh := bezout_gcd(B, max(degF, degG), degH);
    st := time() - st;
    userinfo(1, bezout_update, `stafteriteration=`, print(st));
        return ff, gg, hh, DF, DG, DD, numofiteration, B;
    else
    numofiteration0 := numofiteration + 1;
    userinfo(2, bezout_update, `numofiteration0 =`, print(numofiteration0));
    return bezout_iteration(Btemplate, J, degF, degG, degH, F0, G0, V0, FF, GG, VV, stopcriterion, alpha, numofiteration0); end if
end proc:

    # Find Bezout Matrix from coefficient vector

BezMatrix := proc(F, G, x)
local degF, degG, Bez, Bez1, i, j, fv, gv, gv0;

    # inputs:
    # F:    the polynomial F
    # G:    the polynomial G

    # outputs:

    # Bez:    the Bezout Matrix of F and G

    degF := degree(F, x);
    degG := degree(G, x);
    fv := CoefficientVector(F, x);
    gv0 := CoefficientVector(G, x);
    gv := Vector(degF+1);
    for i from 1 to degG+1 do gv(i):=gv0(i) end do;

    # let the degree of F be the larger one

    #if degF < degG then F, G := G, F ; degF, degG := degG, degF end if;

    # construct the Bezout Matrix

    Bez := Matrix(degF);

    for j from 1 to degF do
       for i from 1 to j do
          Bez(i,j) := fv[i]*gv[j+1]-gv[i]*fv[j+1]
       end do
    end do;

    for i from 2 to degF-1 do
       for j from i to degF-1 do
          Bez(i,j) := Bez(i,j)+Bez(i-1,j+1)
        end do
    end do;

    for j from 1 to degF-1 do
       for i from j+1 to degF do
          Bez(i,j) := Bez(j,i)
       end do
    end do;

    return Bez
end proc:

# Find norm between two polynomials

normpoly := proc(f, g, x)
local n;

    n := Norm(CoefficientVector(f-g, x),2);
    return n
end proc:

# Find GCD of polynomials from their Bezout Matrix

bezout_gcd := proc (B, n, h)
local p, l, u, b, i, j, X, k, y, CoeffVector, gcdB;

    # inputs:
    # B: the Bezout Matrix
    # n: the degree of the Bezout Matrix
    # k: the degree of GCD

    # outputs:

    # gcdB: the GCD calculated by the Bezout Matrix

    p,l,u:=LUDecomposition(B[..,h+1..n]);
    CoeffVector := Vector(h+1);
    CoeffVector(h+1) := 1;
    for k from 1 to h do
        b:=Transpose(p).B[..,k];
        y:=Vector(1..n);
        for i from 1 to n do y(i):=b(i)-add(l[i,j]*y[j],j=1..i-1) end do:
        X:=Vector(1..n-h);
        for i from 1 to n-h do X(n-h+1-i):=(y(n-h+1-i)-add(u[n-h+1-i,j]*X[j],j=n-h+2-i..n-h))/u(n-h+1-i,n-h+1-i) end do:
        CoeffVector(k) := X(1);
    end do:
    #p,l,u:=LUDecomposition(B[1..n-h,h+1..n]);
    #CoeffVector := Vector(h+1);
    #CoeffVector(h+1) := 1;
    #for k from 1 to h do
    #    LinearSolve([p,l,u],B[..,k],method=LU);
    #end do:
    gcdB := FromCoefficientVector(CoeffVector,x);
    return gcdB
end proc:


text:
//impulse: /give @p spawn_egg 1 0 {display:{Name:"Propulsion Gel Spawner"},EntityTag:{id:"minecraft:chicken",CustomName:"MMH_PropulsionGelSpawner",Silent:1,NoGravity:1}}
//MMH_loadChunks()
///clone ~1 ~1 ~1 ~3 ~5 ~3 20 1 24

#MMH
repeat process MMH_propulsionGelSpawner {
if: /testfor @e[type=Chicken,name=MMH_PropulsionGelSpawner]
then {
  MMH_loadChunks()
  /execute @e[type=Chicken,name=MMH_PropulsionGelSpawner] ~ ~ ~ clone 20 1 24 22 5 26 ~-1 ~ ~-1 masked
  /execute @e[type=Chicken,name=MMH_PropulsionGelSpawner] ~ ~ ~ fill ~ ~ ~ ~ ~ ~ air 0 replace pumpkin
  /tp @e[type=Chicken,name=MMH_PropulsionGelSpawner] ~ -100 ~
}
}

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


text:
# Suppose that u is in Lie(A), and that a is in A.  Then there is a unique
# way to express u as a linear combination of fully right associated words
# with a as the last term.  The function below calculates this.

`comb/lie_operad` := (a) -> proc(u)
 local x,y,z,n,ub,uc,t;
 if u = 0 then
  return 0;
 elif type(u,`+`) or type(u,list) or type(u,set) then
  return map(`comb/lie_operad`(a),u);
 elif type(u,`*`) then
  ub,uc := selectremove(type,[op(u)],specfunc(bracket));
  if nops(ub) <> 1 then
   return FAIL;
  fi;
  ub := ub[1];
  uc := `*`(op(uc));
  return(uc * `comb/lie_operad`(a)(ub));
 elif type(u,specfunc(bracket)) then
  n := nops(u);
  if n = 0 then
   return FAIL;
  fi;
  if n = 1 then
   if op(u) = a then
    return u;
   else
    return FAIL;
   fi;
  else
   x := op(1,u);
   y := op(2,u);
   z := op(3..n,u);
   if x = a then
    if n = 2 then
     return -bracket(y,x);
    else
     return `comb/lie_operad`(a)(bracket(y,bracket(a,z))) -
      eval(subs(t = bracket(y,a),`comb/lie_operad`(t)(bracket(t,z))));
    fi;
   elif has(x,a) then
    return eval(subs(t = `comb/lie_operad`(a)(x),
                     `comb/lie_operad`(t)(bracket(t,y,z))));
   else
    return bracket(x,`comb/lie_operad`(a)(bracket(y,z)));
   fi;
  fi;
 fi;
end:

# This returns the basis consisting of fully right associated words
# ending with a.

`basis/lie_operad` := proc(A::set,a_)
 local a,B;

 a := `if`(nargs > 1,a_,A[nops(A)]);
 if not member(a,A) then
  error "a not in A";
 fi;

 B := A minus {a};
 return map(u -> bracket(op(u),a),`list_elements/ord`(B));
end:

`vec0/lie_operad` := (A::set,a) -> proc(u)
 local n,m,ub,uc,T,B,i;
 global `vec_table/lie_operad`;

 n := nops(A);
 m := (n-1)!;
 if u = 0 then
  return [0$m];
 elif type(u,`+`) then
  return map(`vec0/lie_operad`(A,a),u);
 else
  if type(u,specfunc(bracket)) then
   ub := u;
   uc := 1;
  elif type(u,`*`) then
   ub,uc := selectremove(type,[op(u)],specfunc(bracket));
   if nops(ub) <> 1 then
    return FAIL;
   fi;
   ub := ub[1];
   uc := `*`(op(uc));
  else
   return FAIL;
  fi;
  if not(type(`vec_table/lie_operad`[A,a],table)) then
   T := table():
   B := `basis/lie_operad`(A,a);
   for i from 1 to m do
    T[B[i]] := [0$(i-1),1,0$(m-i)];
   od:
   `vec_table/lie_operad`[A,a] := eval(T);
  fi;
  return uc *~ `vec_table/lie_operad`[A,a][ub];
 fi;
end:

`vec/lie_operad` := (A::set,a) -> proc(u)
 `vec0/lie_operad`(A,a)(`comb/lie_operad`(a)(u));
end:

# This returns an alternative basis for Lie(X[n]), where
# X[n] = {x[1],...,x[n]}.  The basis elements for Lie(X[n]) are
# obtained from those for Lie(X[n-1]) by substitutions of the
# form x[i] |-> [x[i],x[n]].  Note that we use the inert Bracket
# operation to prevent automatic expansion by the Jacobi identity.
# One can convert using eval(subs(Bracket = bracket,...))

`alt_basis/lie_operad` := proc(n::nonnegint)
 option remember;
 local i,u;

 if n = 0 then
  return [];
 elif n = 1 then
  return [x[1]];
 else
  return [seq(seq(subs(x[i]=Bracket(x[i],x[n]),u),i=1..n-1),u in `alt_basis/lie_operad`(n-1))];
 fi;
end:


text:
with(Groebner):
with(PolynomialIdeals):

J := PolynomialIdeal({10*x1*x2^2 + 10*x1*x3^2 + 10*x1*x4^2 + 10*x1*x5^2 + 10*x1*x6^2 + 10*x1*x7^2 - 11*x1 + 10, 10*x1^2*x2 + 10*x2*x3^2 + 10*x2*x4^2 + 10*x2*x5^2 + 10*x2*x6^2 + 10*x2*x7^2 - 11*x2 + 10, 10*x1^2*x3 + 10*x2^2*x3 + 10*x3*x4^2 + 10*x3*x5^2 + 10*x3*x6^2 + 10*x3*x7^2 - 11*x3 + 10, 10*x1^2*x4 + 10*x2^2*x4 + 10*x3^2*x4 + 10*x4*x5^2 + 10*x4*x6^2 + 10*x4*x7^2 - 11*x4 + 10, 10*x1^2*x5 + 10*x2^2*x5 + 10*x3^2*x5 + 10*x4^2*x5 + 10*x5*x6^2 + 10*x5*x7^2 - 11*x5 + 10, 10*x1^2*x6 + 10*x2^2*x6 + 10*x3^2*x6 + 10*x4^2*x6 + 10*x5^2*x6 + 10*x6*x7^2 - 11*x6 + 10, 10*x1^2*x7 + 10*x2^2*x7 + 10*x3^2*x7 + 10*x4^2*x7 + 10*x5^2*x7 + 10*x6^2*x7 - 11*x7 + 10}):
st := time[real]():
Basis(J, tdeg(x1,x2,x3,x4,x5,x6,x7), charactesistic=2^31-1):
time[real]() - st;

J := PolynomialIdeal({10*x1*x2^2 + 10*x1*x3^2 + 10*x1*x4^2 + 10*x1*x5^2 + 10*x1*x6^2 + 10*x1*x7^2 + 10*x1*x8^2 - 11*x1 + 10, 10*x1^2*x2 + 10*x2*x3^2 + 10*x2*x4^2 + 10*x2*x5^2 + 10*x2*x6^2 + 10*x2*x7^2 + 10*x2*x8^2 - 11*x2 + 10, 10*x1^2*x3 + 10*x2^2*x3 + 10*x3*x4^2 + 10*x3*x5^2 + 10*x3*x6^2 + 10*x3*x7^2 + 10*x3*x8^2 - 11*x3 + 10, 10*x1^2*x4 + 10*x2^2*x4 + 10*x3^2*x4 + 10*x4*x5^2 + 10*x4*x6^2 + 10*x4*x7^2 + 10*x4*x8^2 - 11*x4 + 10, 10*x1^2*x5 + 10*x2^2*x5 + 10*x3^2*x5 + 10*x4^2*x5 + 10*x5*x6^2 + 10*x5*x7^2 + 10*x5*x8^2 - 11*x5 + 10, 10*x1^2*x6 + 10*x2^2*x6 + 10*x3^2*x6 + 10*x4^2*x6 + 10*x5^2*x6 + 10*x6*x7^2 + 10*x6*x8^2 - 11*x6 + 10, 10*x1^2*x7 + 10*x2^2*x7 + 10*x3^2*x7 + 10*x4^2*x7 + 10*x5^2*x7 + 10*x6^2*x7 + 10*x7*x8^2 - 11*x7 + 10, 10*x1^2*x8 + 10*x2^2*x8 + 10*x3^2*x8 + 10*x4^2*x8 + 10*x5^2*x8 + 10*x6^2*x8 + 10*x7^2*x8 - 11*x8 + 10}):
st := time[real]():
Basis(J, tdeg(x1,x2,x3,x4,x5,x6,x7,x8), characteristic=2^31-1):
time[real]() - st;

J := PolynomialIdeal({10*x1*x2^2 + 10*x1*x3^2 + 10*x1*x4^2 + 10*x1*x5^2 + 10*x1*x6^2 + 10*x1*x7^2 + 10*x1*x8^2 + 10*x1*x9^2 - 11*x1 + 10,
10*x1^2*x2 + 10*x2*x3^2 + 10*x2*x4^2 + 10*x2*x5^2 + 10*x2*x6^2 + 10*x2*x7^2 + 10*x2*x8^2 + 10*x2*x9^2 - 11*x2 + 10,
10*x1^2*x3 + 10*x2^2*x3 + 10*x3*x4^2 + 10*x3*x5^2 + 10*x3*x6^2 + 10*x3*x7^2 + 10*x3*x8^2 + 10*x3*x9^2 - 11*x3 + 10,
10*x1^2*x4 + 10*x2^2*x4 + 10*x3^2*x4 + 10*x4*x5^2 + 10*x4*x6^2 + 10*x4*x7^2 + 10*x4*x8^2 + 10*x4*x9^2 - 11*x4 + 10,
10*x1^2*x5 + 10*x2^2*x5 + 10*x3^2*x5 + 10*x4^2*x5 + 10*x5*x6^2 + 10*x5*x7^2 + 10*x5*x8^2 + 10*x5*x9^2 - 11*x5 + 10,
10*x1^2*x6 + 10*x2^2*x6 + 10*x3^2*x6 + 10*x4^2*x6 + 10*x5^2*x6 + 10*x6*x7^2 + 10*x6*x8^2 + 10*x6*x9^2 - 11*x6 + 10,
10*x1^2*x7 + 10*x2^2*x7 + 10*x3^2*x7 + 10*x4^2*x7 + 10*x5^2*x7 + 10*x6^2*x7 + 10*x7*x8^2 + 10*x7*x9^2 - 11*x7 + 10,
10*x1^2*x8 + 10*x2^2*x8 + 10*x3^2*x8 + 10*x4^2*x8 + 10*x5^2*x8 + 10*x6^2*x8 + 10*x7^2*x8 + 10*x8*x9^2 - 11*x8 + 10,
10*x1^2*x9 + 10*x2^2*x9 + 10*x3^2*x9 + 10*x4^2*x9 + 10*x5^2*x9 + 10*x6^2*x9 + 10*x7^2*x9 + 10*x8^2*x9 - 11*x9 + 10}):
st := time[real]():
Basis(J, tdeg(x1,x2,x3,x4,x5,x6,x7,x8,x9), characteristic=2^31-1):
time[real]() - st;


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

XX := (z, xs) -> xs*opz_pow_n(z,4/3)*2^(-4/3):
YY := (z, xt, xs0, xs1) -> 2*(XX(z, xs0)^2 + XX(-z, xs1)^2) - xt^2:

f_th := (rs, z, xt, xs0, xs1) -> add(params_a_omega[i]
  * (n_spin(rs, z)^params_a_a[i] + n_spin(rs, -z)^params_a_a[i])
  * z^(2*params_a_b[i])
  * 1/2*(XX(z, xs0)^params_a_c[i] + XX(-z, xs1)^params_a_c[i])
  * YY(z, xt, xs0, xs1)^params_a_d[i], i=1..params_a_n)/n_total(rs):


text:
main()inline{
	Trace3(i, j, k);
	do(7){
		Trace3(i, j, k);
		do(3){
			Trace3(i, j, k);
			do(5) Trace3(i, j, k);
		}
	}
	i = j < k ? j : k;
}


text:
# Kinematik-Berechnung 3. Arm KAS5 Modell 7
# Beschreibung
#
# Berechnung der kinematischen Zwangsbedingungen des Systems KAS5m7 in impliziter Form.
# Ansatz: Schließen von Gelenkketten nach [1] und Anwendung der Zahnrad-Zwangsbedingungen.
#
# Quellen
# [1] Khalil and Bennis: Symbolic calculation of the base inertial parameters of closed-loop robots (1995)
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2018-02
# Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Initialisierung
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env_IC":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name_OL):
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name_OL, base_method_name);
read "../robot_codegen_definitions/robot_env_IC":
Trf := Trf:
Trf_c := Trf_c:
# Schleife über B-C-D
# 3-4-12-13-17
T_3_17 := combine( Matrix(Trf(1..4,1..4, 4)) . Matrix(Trf(1..4,1..4,12)) . Matrix(Trf(1..4,1..4,13)) . Matrix(Trf(1..4,1..4,17)) ):
# 3-10-11-20
T_3_20 := combine( Matrix(Trf(1..4,1..4, 10)) . Matrix(Trf(1..4,1..4,11)) . Matrix(Trf(1..4,1..4,20)) ):
h1 := T_3_17(1..3,4) - T_3_20(1..3,4):
# Schleife über 2-F-E
#
# 2-3-10-16
T_2_16 := combine( Matrix(Trf(1..4,1..4, 3)) . Matrix(Trf(1..4,1..4,10)) . Matrix(Trf(1..4,1..4,16)) ):
# 2-8-9-19
T_2_19 := combine( Matrix(Trf(1..4,1..4, 8)) . Matrix(Trf(1..4,1..4,9)) . Matrix(Trf(1..4,1..4,19)) ):
# Vektorkette
h2 := T_2_16(1..3,4) - T_2_19(1..3,4):
# Schleife über 5-B-A

# 4-5-6-18
T_4_18 := combine( Matrix(Trf(1..4,1..4, 5)) . Matrix(Trf(1..4,1..4,6)) . Matrix(Trf(1..4,1..4,18)) ):
# 4-12-13-14-15-21
T_4_21 := combine( Matrix(Trf(1..4,1..4, 12)) . Matrix(Trf(1..4,1..4,13)) . Matrix(Trf(1..4,1..4,14)) . Matrix(Trf(1..4,1..4,15)) . Matrix(Trf(1..4,1..4,21)) ):
# Vektorkette
h3 := T_4_18(1..3,4) - T_4_21(1..3,4):
# Zahnrad-Kontakt
#
hz1 := -qJ_t(4) + qJ_t(11):
hz2 := -qJ_t(6) + qJ_t(5) - qJ_t(4) + delta18:
h_user := <hz1; hz2>:
# Zusammenstellen aller Zwangsbedingungen
#
implconstr_t := <h1([1, 2],1);h2([1, 3],1);h3([1, 2],1); h_user>:
implconstr_s := convert_t_s(implconstr_t):
# Exportiere Code für folgende Skripte
#
#
kin_constraints_exist:=true:
save implconstr_t, implconstr_s, kin_constraints_exist, sprintf("../codeexport/%s/tmp/kinematic_constraints_implicit_maple.m", robot_name):
# Exportieren des vollständigen Ausdruckes
if codegen_act then
  MatlabExport(implconstr_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_matlab.m", robot_name), 2):
end if:
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( implconstr_s )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_matlab.m", robot_name), 2):
interface(rtablesize=100):
kc_symbols
;


text:
######################################################################
# Here B is assumed to be a subset of A, and SWW(A,B,N) is the space
# denoted by S(NWA,NWB) in the LaTeX document.

`is_element/SWW` := (N::posint) -> (A::set,B::set) -> proc(xy)
 local x,y,x0,y0,a,b,r;
 global reason;

 if not(type(xy,list) and nops(xy) = 2) then
  reason := [convert(procname,string),"xy cannot be split as [x,y]",xy];
  return false;
 fi;

 x,y := op(xy);
 if not(`is_element/SW`(N)(A)(x)) then
  reason := [convert(procname,string),"x is not in SW(N)(A)",x,N,A];
  return false;
 fi;

 if not(`is_element/SW`(N)(B)(y)) then
  reason := [convert(procname,string),"y is not in SW(N)(B)",y,N,B];
  return false;
 fi;

 x0 := table();
 for b in B do x0[b] := x[b]; od;
 x0 := `normalise/W`(N)(B)(x0);
 y0 := `normalise/SW`(N)(B)(y);

 r := 0;
 for b in B do r := max(r,op(map(abs,x0[b]))); od;
 if r = 0 then return true; fi;
 for b in B do
  if simplify(x0[b] -~ r *~ y0[b]) <> [0$N] then
   reason := [convert(procname,string),"x0 and y0 are incompatible at b",x0,y0,b,r];
   return false;
  fi;
 od;

 return true;
end;

`normalise/SWW` := (N::posint) -> (A::set,B::set) -> proc(xy)
 [`normalise/SW`(N)(A)(xy[1]),
  `normalise/SW`(N)(B)(xy[2])];
end;

`is_equal/SWW` := (N::posint) -> (A::set,B::set) -> proc(x,y)
 local p,q,r,s;
 global reason;

 p,q := op(`normalise/SWW`(N)(A,B)(x));
 r,s := op(`normalise/SWW`(N)(A,B)(y));

 if not `is_equal/SW`(N)(A)(p,r) then
  reason := [convert(procname,string),"p <> r",p,r,reason];
  return false;
 fi;

 if not `is_equal/SW`(N)(B)(q,s) then
  reason := [convert(procname,string),"q <> s",q,s,reason];
  return false;
 fi;

 return true;
end;

`is_leq/SWW` := NULL;

`random_element/SWW` := (N::posint) -> (A::set,B::set) -> proc()
 local x,y,b,u,c;

 x := `random_element/SW`(N)(A)();

 if rand(1..3)() = 1 then
  u := `random_element/R`(N)();

  for b in B do
   x[b] := u;
  od;
  y := `random_element/SW`(N)(B)();
 else
  y := `rho/W`(A,B)(x);
  if `is_zero/W`(N)(B)(y) then
   y := `random_element/SW`(N)(B)();
  fi;
 fi;

 return [eval(x),eval(y)];
end;

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


text:
######################################################################
# big_subsets(A) is the set of all subsets B of A with |B| > 1

`is_element/big_subsets` := (A::set) -> proc(B)
 type(B,set) and B minus A = {} and nops(B) > 1;
end;

`is_equal/big_subsets` := (A::set) -> (B,C) -> evalb(B = C):

`is_leq/big_subsets` := (A::set) -> (B,C) -> evalb(B minus C = {}):

`random_element/big_subsets` := (A::set) -> proc()
 local r,B;

 if nops(A) < 2 then return FAIL; fi;

 r := rand(2);
 B := {};

 while nops(B) < 2 do
  B := select(a -> (r() = 1),A);
 od;

 return B;
end;

`list_elements/big_subsets` := (A::set) ->
 select(B -> nops(B) > 1,sort(map(sort,[op(combinat[powerset](A))]))):

`count_elements/big_subsets` := (A::set) ->
 2^nops(A) - nops(A) - 1;


text:
func $switchfunc ( var %n i32) i32 {
 var %alocal i32
 switch (dread i32 %n) @labdft {
   -3: goto @lab0
   1: goto @lab1
   9: goto @lab9 }
@lab0
 return (constval i32 -3)
@labdft
 return (constval i32 100)
@lab9
 return (constval i32 9)
@lab1
 return (constval i32 1) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


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

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

(* Equation 28 squared with Equation 25 built in *)
tpss_xi2 := (z, xt, xs0, xs1) ->
  (1 - z^2)*(t_total(z, xs0^2, xs1^2) - xt^2)/(2*(3*Pi^2)^(1/3))^2:

(* Equation 33 *)
tpss_C00 := (cc, z) ->
  + add(cc[i]*z^(2*(i-1)), i=1..4):

(* Equation 34 *)
(* The series expansion of C0 for z -> +- 1 goes as C + D*(1+-z)^(2/3) + O(1+-z)
   whose first derivative seems to diverge for ferromagnetic densities, leading
   to severe numerical instabilities. I can not correct for the bad design of
   the functional, and this is the possible solution that I found *)
tpss_C0_den := (z, xt, xs0, xs1) ->
  1 + tpss_xi2(z, xt, xs0, xs1)*((1 + z)^(-4/3) + (1 - z)^(-4/3))/2:
tpss_C0 := (cc, z, xt, xs0, xs1) -> my_piecewise3(1 - m_abs(z) <= 1e-12,
  add(cc[i], i=1..4),
  tpss_C00(cc, z) / tpss_C0_den(z_thr(z), xt, xs0, xs1)^4):

(* Equation 11, with tau_W from Equation 12 *)
tpss_aux := (z, xt, ts0, ts1) ->
  m_min(xt^2/(8*t_total(z, ts0, ts1)), 1):

(* n_sigma/n \epsilon^sigma in Equation 25 *)
tpss_par_s0 := (f_gga, rs, z, xt, xs0, xs1) ->
  m_max(f_gga(rs*(2/(1 + z))^(1/3),  1, xs0, xs0, 0), f_gga(rs, z, xt, xs0, xs1))*(1 + z)/2:
tpss_par_s1 := (f_gga, rs, z, xt, xs0, xs1) ->
  m_max(f_gga(rs*(2/(1 - z))^(1/3), -1, xs1, 0, xs1), f_gga(rs, z, xt, xs0, xs1))*(1 - z)/2:

(* Second line of Equation 25 *)
(* The screening of the density is important in order to stabilize this functional
   for ferromagnetic densities *)
tpss_par  := (f_gga, rs, z, xt, xs0, xs1, ts0, ts1) ->
  - (1 + tpss_C0(params_a_C0_c, z, xt, xs0, xs1))*tpss_aux(z, xt, ts0, ts1)^2*(
    + my_piecewise3(screen_dens_zeta(rs,  z),
        f_gga(rs, z_thr(z), xt, xs0, xs1)*(1 + z)/2,
        tpss_par_s0(f_gga, rs, z_thr(z), xt, xs0, xs1)
      )
    + my_piecewise3(screen_dens_zeta(rs, -z),
        f_gga(rs, z_thr(z), xt, xs0, xs1)*(1 - z)/2,
        tpss_par_s1(f_gga, rs, z_thr(z), xt, xs0, xs1)
      )
  ):

(* First line of Equation 25 *)
tpss_perp := (f_gga, rs, z, xt, xs0, xs1, ts0, ts1) ->
  (1 + tpss_C0(params_a_C0_c, z, xt, xs0, xs1)*tpss_aux(z, xt, ts0, ts1)^2)
  * f_gga(rs, z, xt, xs0, xs1):

(* Equation in full 25 *)
tpss_f0 := (f_gga, rs, z, xt, xs0, xs1, ts0, ts1) ->
  + tpss_par (f_gga, rs, z, xt, xs0, xs1, ts0, ts1)
  + tpss_perp(f_gga, rs, z, xt, xs0, xs1, ts0, ts1):

(* Equation 24 *)
tpss_f := (f_gga, rs, z, xt, xs0, xs1, ts0, ts1) ->
  + tpss_f0(f_gga, rs, z, xt, xs0, xs1, ts0, ts1)
  * (1 + params_a_d*tpss_f0(f_gga, rs, z, xt, xs0, xs1, ts0, ts1)*tpss_aux(z, xt, ts0, ts1)^3):


text:
    {    MultPlan.MPL    }

TITLE
    Multistage_Production_Planning ;

INDEX
    time     =  (June,July,August)  -> (Jun,Jul,Aug);
    tire     =  (Nylon,Glass)       -> (Nyl,Gls);
    machine  =  (Wheel,Regular)     -> (W,R);

DATA
    ProdHours[machine,time]  =  (( 700, 300,1000),
                                 (1500, 400, 300));

    Delivery[tire,time]      =  ((4000, 8000, 3000),
                                 (1000, 5000, 5000));

    ProdRates[machine,tire]  =  ((0.15, 0.12),
                                 (0.16, 0.14));

    Prices[tire]             =  (7.00, 9.00)
    PurchasePrice[tire]      =  (3.10, 3.90)

DECISION
    Prod[machine,tire,time]   ->  P
    Invt[tire,time];

MACRO
    TotalRevenues   :=  SUM(tire,time: Prices * Delivery) ;
    ProductionCost  :=  0.50 SUM(machine,tire,time: ProdRates * Prod) ;
    InventoryCost   :=  0.10 SUM(tire,time: Invt) ;

    RawMaterialCost :=  SUM(tire,time: PurchasePrice * Delivery);
    Post&PackCost   :=  0.23 SUM(tire,time: Delivery) ;
    FixedCost       :=  RawMaterialCost + Post&PackCost ;

{----------------------------------------------------------------------------}

  MAX Profit   =  TotalRevenues
                  -  ProductionCost
                  -  InventoryCost
                  -  FixedCost ;

SUBJECT TO

  ProdCap[machine,time] :

     SUM(tire: ProdRates * Prod)  <  ProdHours ;

  Inventory[tire,time] :

     Invt[time-1] + SUM(machine: Prod)  =  Invt + Delivery ;

BOUNDS
     CloseInvt[tire,time=August]  :  Invt = 0 ;
END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
# If X ~ Exponential(λ) then k*exp(X) ~ Pareto(λ, k)

# Proof:

PDF_exponential := (alpha, x) -> exp(-x/alpha)/alpha;

PDF_pareto := (lambda, k, x) -> lambda * k ^ lambda / x ^ (lambda + 1);

simplify(PDF_exponential(1/lambda, log(x/k))*diff(log(x/k), x));

simplify(PDF_pareto(lambda, k, x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


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

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

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

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

$ifdef gga_x_rpbe_params
params_a_rpbe_kappa := KAPPA_PBE:
params_a_rpbe_mu    := MU_PBE:
$endif

rpbe_f0 := s -> 1 + params_a_rpbe_kappa * (
  1 - exp(-params_a_rpbe_mu*s^2/params_a_rpbe_kappa)
):
rpbe_f  := x -> rpbe_f0(X2S*x):

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


text:
 func $foo (
  var %i i32
#var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     depositbits i32 1 23(
       constval i32 0xfff, constval i32 0x11))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
check_koszul := proc()
 local A,n,dA,u,J,i;

 printf("%a()\n",procname);

 A := <<1|5|4|9>,<2|2|7|5>,<3|9|9|6>,<4|3|2|1>>;
 n := Dimension(A)[1];
 dA := Determinant(A);
 u := [4,5,6,11];

 J := Matrix(n,n);
 for i from 1 to n do J[i,n+1-i] := (-1)^(i-1); od;

 _ASSERT(
  Equal(koszul_adjoint(0,A), <<Determinant(A)>>),
  "koszul adjoint(0,A) = det(A)"
 );

 _ASSERT(
  Equal(koszul_adjoint(1,A), Adjoint(A)),
  "koszul_adjoint(1,A) = adj(A)"
 );

 _ASSERT(
  Equal(koszul_adjoint(n-1,A), Transpose((-1)^(n-1) * J . A . J)),
  "koszul_adjoint(n-1,A) is conjugate to A^T"
 );

 _ASSERT(
  Equal(koszul_adjoint(n,A), <<1>>),
  "koszul_adjoint(n,A) = 1"
 );

 _ASSERT(
  `and`(seq(Equal(exterior_power(k,A) . koszul_adjoint(k,A),
                  dA * IdentityMatrix(binomial(n,k))),k=0..n)),
  "exterior_power(k,A) o koszul_adjoint(k,A) = det(A) in all degrees"
 );

 _ASSERT(
  `and`(seq(Equal(koszul_d(k,u) . koszul_d(k+1,u),
                  Matrix(binomial(n,k-1),binomial(n,k+1))),k=1..n-1)),
  "koszul_d ^ 2 = 0"
 );
end:

check_koszul_comparison := proc()
 global v,n,u_vec,v_vec,w_vec,x_vec,x_vars,Iv,Mv,Nv,
        soc_u,soc_v,soc_w,A,B,p;

 v := [x[1]^2+x[2]^3,x[1]*x[2]+x[3],x[1]^4-x[3]^2*x[2],x[4]-x[2]^3];
 koszul_setup(v);

 _ASSERT(
  Equal(map(koszul_pmod,map(expand,u_vec)),
        map(koszul_pmod,map(expand,v_vec . A))),
  "u = v . A"
 );

 _ASSERT(
  Equal(map(koszul_pmod,map(expand,v_vec)),
        map(koszul_pmod,map(expand,w_vec . B))),
  "v = w . B"
 );

 _ASSERT(
   Equal(map(koszul_pmod,map(Nv,soc_v * x_vec)),Transpose(Vector(n))),
   "soc_v is in the socle"
 );

 _ASSERT(
   not(koszul_pmod(Nv(soc_v)) = 0),
   "soc_v is nonzero"
 );

 _ASSERT(
  `and`(seq(
   Equal(map(expand,koszul_d(k,v_vec) . exterior_power(k,A)),
	 map(expand,exterior_power(k-1,A) . koszul_d(k,u_vec))),
   k=1..n
  )),
  "exterior_power(*,A) is a chain map"
 );

 _ASSERT(
  `and`(seq(
   Equal(map(expand,koszul_d(k,w_vec) . exterior_power(k,B)),
	 map(expand,exterior_power(k-1,B) . koszul_d(k,v_vec))),
   k=1..n
  )),
  "exterior_power(*,B) is a chain map"
 );

 _ASSERT(
  `and`(seq(
   Equal(map(expand,koszul_d(k,u_vec) . koszul_adjoint(k,A)),
	 map(expand,koszul_adjoint(k-1,A) . koszul_d(k,v_vec))),
   k=1..n
  )),
  "koszul_adjoint(*,A) is a chain map"
 );

 _ASSERT(
  `and`(seq(
   Equal(map(expand,koszul_d(k,v_vec) . koszul_adjoint(k,B)),
	 map(expand,koszul_adjoint(k-1,B) . koszul_d(k,w_vec))),
   k=1..n
  )),
  "koszul_adjoint(*,B) is a chain map"
 );

end:


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

#####################################################################
#
# plate/array tests
#
#####################################################################

triv := Plate(k, i, Ret(i)):
TestHakaru(triv, Ret(ary(k, i, i)), label="Dirac Plate");

ary0 := Bind(Plate(k, i, Gaussian(0,1)), xs, Ret([xs])):
TestHakaru(ary0, ary0, label="plate roundtripping", ctx = [k::nonnegint]);

ary1  := Bind(Gaussian(0,1), x,
         Bind(Plate(n, i, Weight(density[Gaussian](x,1)(idx(t,i)), Ret(Unit))), ys,
         Ret(x))):
ary1w := 2^(-(1/2)*n+1/2)*exp(sum(idx(t,i),i=0..n-1)^2/(n+1)/2)*exp(-sum(idx(t,i)^2,i=0..n-1)/2)*Pi^(-(1/2)*n)/sqrt(2+2*n):
TestHakaru(ary1,
  Weight(ary1w, Gaussian((sum(idx(t, i), i = 0 .. n-1))/(n+1), 1/sqrt(n+1))),
  label="Wednesday goal", ctx = [n::nonnegint]);
TestHakaru(Bind(ary1, x, Ret(Unit)), Weight(ary1w, Ret(Unit)),
  label="Wednesday goal total", ctx = [n::nonnegint]);
ary2  := Bind(Gaussian(0,1), x,
         Bind(Plate(n, i, Bind(Gaussian(idx(t,i),1),z, Weight(density[Gaussian](x,1)(idx(t,i)), Ret(z+1)))), ys,
         Ret(ys))):
TestHakaru(ary2,
  Weight(ary1w, Bind(Plate(n, i, Gaussian(idx(t,i),1)), zs, Ret(ary(n, i, idx(zs,i)+1)))),
  label="Reason for fission", ctx = [n::nonnegint]);
ary3  := Bind(Gaussian(0,1), x,
         Bind(Plate(n, i, Bind(Gaussian(idx(t,i),1),z, Weight(density[Gaussian](x,1)(idx(t,i)), Ret(z)))), zs,
         Ret(zs))):
TestHakaru(ary3, Weight(ary1w, Plate(n, i, Gaussian(idx(t,i),1))),
  label="Array eta", ctx = [n::nonnegint]);
TestHakaru(Ret(ary(n,i,idx(f(i),i))), label="Don't overdo array eta");

bry1  := Bind(BetaD(alpha,beta), x,
         Bind(Plate(n, i, Weight(x    ^piecewise(idx(y,i)=true ,1) *
                                 (1-x)^piecewise(idx(y,i)=false,1),
                          Ret(Unit))), ys,
         Ret(x))):
bry1s := Weight(Beta(alpha+sum(piecewise(idx(y,i)=true ,1), i=0..n-1),
                     beta +sum(piecewise(idx(y,i)=false,1), i=0..n-1))/Beta(alpha,beta),
         BetaD(alpha+sum(piecewise(idx(y,i)=true ,1), i=0..n-1),
               beta +sum(piecewise(idx(y,i)=false,1), i=0..n-1))):
TestHakaru(bry1, bry1s,
  label="first way to express flipping a biased coin many times",
  ctx = [n::nonnegint]);

bry2  := Bind(BetaD(alpha,beta), x,
         Bind(Plate(n, i, Weight(x    ^(  idx(y,i)) *
                                 (1-x)^(1-idx(y,i)),
                          Ret(Unit))), ys,
         Ret(x))):
bry2s := Weight(Beta(alpha+  sum(idx(y,i),i=0..n-1),
                     beta +n-sum(idx(y,i),i=0..n-1))/Beta(alpha,beta),
         BetaD(alpha+  sum(idx(y,i),i=0..n-1),
               beta +n-sum(idx(y,i),i=0..n-1))):
TestHakaru(bry2, bry2s,
  label="second way to express flipping a biased coin many times",
  ctx = [n::nonnegint]);

fission     := Bind(Plate(k, i, Gaussian(f(i),1)), xs, Plate(k, i, Gaussian(idx(xs,i),1))):
fusion      := Plate(k, i, Bind(Gaussian(f(i),1), x, Gaussian(x,1))):
conjugacies := Plate(k, i, Gaussian(f(i), sqrt(2))):
conjugacies5:= Bind(Gaussian(f(0), sqrt(2)), x0,
               Bind(Gaussian(f(1), sqrt(2)), x1,
               Bind(Gaussian(f(2), sqrt(2)), x2,
               Bind(Gaussian(f(3), sqrt(2)), x3,
               Bind(Gaussian(f(4), sqrt(2)), x4,
               Ret([x0,x1,x2,x3,x4])))))):
TestHakaru(     fission                  ,      conjugacies                   , verify=normal, label="Conjugacy across plates (function)", ctx=[k::nonnegint]);
TestHakaru(     fusion                   ,      conjugacies                   , verify=normal, label="Conjugacy in plate (function)", ctx=[k::nonnegint]);
TestHakaru(eval(fission,{k=5           }),      conjugacies5                  , verify=normal, label="Conjugacy across plates unrolled (function)");
TestHakaru(eval(fusion ,{k=5           }),      conjugacies5                  , verify=normal, label="Conjugacy in plate unrolled (function)");
TestHakaru(eval(fission,{    f=1       }), eval(conjugacies ,{    f=1       }), verify=normal, label="Conjugacy across iid plates", ctx=[k::nonnegint]);
TestHakaru(eval(fusion ,{    f=1       }), eval(conjugacies ,{    f=1       }), verify=normal, label="Conjugacy in iid plate", ctx=[k::nonnegint]);
TestHakaru(eval(fission,{k=5,f=1       }), eval(conjugacies5,{    f=1       }), verify=normal, label="Conjugacy across iid plates unrolled");
TestHakaru(eval(fusion ,{k=5,f=1       }), eval(conjugacies5,{    f=1       }), verify=normal, label="Conjugacy in iid plate unrolled");
TestHakaru(eval(fission,{    f=(i->z^i)}), eval(conjugacies ,{    f=(i->z^i)}), verify=normal, label="Conjugacy across plates", ctx=[k::nonnegint,z>0]);
TestHakaru(eval(fusion ,{    f=(i->z^i)}), eval(conjugacies ,{    f=(i->z^i)}), verify=normal, label="Conjugacy in plate", ctx=[k::nonnegint,z>0]);
TestHakaru(eval(fission,{k=5,f=(i->z^i)}), eval(conjugacies5,{k=5,f=(i->z^i)}), verify=normal, label="Conjugacy across plates unrolled", ctx=[z>0]);
TestHakaru(eval(fusion ,{k=5,f=(i->z^i)}), eval(conjugacies5,{k=5,f=(i->z^i)}), verify=normal, label="Conjugacy in plate unrolled", ctx=[z>0]);

# Don't be confused by extra iterations in plate
extra_1 := Bind(Plate(n+1,i,Gaussian(0,1)),xs, Plate(n,i,Gaussian(idx(xs,i+1),1))):
extra_2 := Bind(Plate(n+1,i,Gaussian(0,1)),xs, Plate(n,i,Gaussian(idx(xs,i  ),1))):
extra_s := Plate(n,i,Gaussian(0,sqrt(2))):
TestHakaru(extra_1, extra_s, label="Don't be confused by extra iterations at beginning of plate", ctx=[n::nonnegint]);
TestHakaru(extra_2, extra_s, label="Don't be confused by extra iterations at end of plate"      , ctx=[n::nonnegint]);

# Simplify by size of array
TestHakaru(Bind(Plate(k,c,Uniform(37,42)),xs,Weight(f(size(xs)),Ret(Unit))),
           Weight(f(k),Ret(Unit)),
           label="plate size", ctx = [k::nonnegint]):

# Simplifying gmm below is a baby step towards index manipulations we need
gmm := Bind(Plate(k, c, Gaussian(0,1)), xs,
       Bind(Plate(n, i, Weight(density[Gaussian](idx(xs,idx(cs,i)),1)(idx(t,i)), Ret(Unit))), ys,
       Ret(xs))):
gmm_s := Weight(2^(-(1/2)*n)*Pi^(-(1/2)*n)*exp(-(1/2)*(sum(idx(t,i)^2, i=0..n-1)))*
  exp((1/2)*(sum((sum(piecewise(c=idx(cs,i), idx(t,i), 0), i=0..n-1))^2/(1+sum(piecewise(c=idx(cs,i), 1, 0), i=0..n-1)), c=0..k-1))) *
  sqrt(1/product(1+sum(piecewise(c=idx(cs,i), 1, 0), i=0..n-1), c=0..k-1)),
    Plate(k, c,
      Gaussian((sum(piecewise(c=idx(cs,i), idx(t,i), 0), i=0..n-1))/
               (sum(piecewise(c=idx(cs,i), 1, 0), i=0..n-1)+1),
               (1/(sum(piecewise(c=idx(cs,i), 1, 0), i=0..n-1)+1))^(1/2)))):
TestHakaru(gmm, gmm_s,
           label="gmm (currently failing -- due to alpha-equivalence)",
           ctx = [k::nonnegint, n::nonnegint]);
TestHakaru(Bind(gmm, xs, Ret(cs)),
           subsop(2=Ret(cs), gmm_s),
           label="gmm total (currently failing -- due to alpha-equivalence)",
           ctx = [k::nonnegint, n::nonnegint]);

# Detecting Dirichlet-multinomial conjugacy when unrolled
dirichlet := proc(as)
  Bind(Plate(size(as)-1, i, BetaD(sum(idx(as,j), j=0..i), idx(as,i+1))), xs,
  Ret(ary(size(as), i,
          product(idx(xs,j), j=i..size(as)-2)
          * piecewise(0=i, 1, 1-idx(xs,i-1)))))
end proc:
categorical := proc(ps)
  Bind(Counting(0, size(ps)), i, Weight(idx(ps,i), Ret(i)))
end proc:
# The next line eta-expands CodeTools[Test] so that its arguments get evaluated
(proc() CodeTools[Test](_passed) end proc)
  (op(2, [unweight(
     fromLO(improve(toLO(
       Plate(n,k,
         Bind(dirichlet(ary(5,i,1)),ps,
           Weight(product(idx(ps,i)^idx(idx(t,k),i),i=0..size(ps)-1),
             Ret(ps))))))))]),
   fromLO(toLO(
     Plate(n,k,
       dirichlet(ary(5,i,1+idx(idx(t,k),i)))))),
   measure(simplify),
   label="Conjugacy between unrolled symmetric Dirichlet and multinomial");
# We'd like the test above to pass even if the count 5 becomes symbolic.
# Below is some progress towards this goal:
TestHakaru(dirichlet(as), label="Dirichlet symbolic prior roundtrip");
TestHakaru(Context(size(as)>0,
             Bind(dirichlet(as),ps,
               Weight(product(idx(ps,i)^idx(t,i),i=0..size(ps)-1),
                 Ret(ps)))),
           Context(size(as)>0,
             Weight(product(Beta(idx(as,i+1)+idx(t,i+1),idx(t,0)+sum(idx(as,k),k=0..i)+sum(idx(t,k),k=1..i)),i=0..size(as)-2)
                    /product(Beta(idx(as,i+1),sum(idx(as,k),k=0..i)),i=0..size(as)-2),
               Bind(Plate(size(as)-1,i,BetaD(idx(t,0)+sum(idx(as,k),k=0..i)+sum(idx(t,k),k=1..i),idx(as,i+1)+idx(t,i+1))),xs,
                 Ret(ary(size(as),i,piecewise(i=0,product(idx(xs,j),j=0..size(as)-2),product(idx(xs,j),j=i..size(as)-2)*(1-idx(xs,i-1)))))))),
           label="Conjugacy between rolled Dirichlet and multinomial");


text:
 func $landi32 (
  var %i i32, var %j i32
  ) i32 {
   return (
     land i32(dread i32 %i, dread i32 %j))}

 func $landi64 (
  var %i i64, var %j i64
  ) i32 {
   return (
     land i64(dread i64 %i, dread i64 %j))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# Berechne kinematische Zwangsbedingungen für Viergelenkkette
# Einleitung

# Die impliziten kinematischen Zwangsbedingungen werden aufgestellt: Annahme : Alle MDH-Winkel seien unabhängig und die Kette sei offen
#
# fivebar1_OL -> Fünfgelenkkette, modellierung impliziten der Zwangsbedingungen mit Ausdrücken einer offenen Kette  und basierend auf der MDH-Tabelle.
# Bei Matlab: Abhängige Winkel durch die Ausdrücke von TE aus der Datei "fourbar1TE_kinconstr_expl_mdh_sym_varpar" ersetzen.
# Quelle
# SA Bejaoui: Bejaoui2018_S749; "Modellierung kinematischer Zwangsbedingungen für hybride serielle Roboter mit planaren Parallelmechanismen"
# Autor
# Abderahman Bejaoui
# Studienarbeit bei: Moritz Schappler, schappler@irt.uni-hannover.de, 2018-08
# (C) Institut fuer Mechatronische Systeme, Leibniz Universitaet Hannover
#
# Initialisierung
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env_IC":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name_OL):
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name_OL, base_method_name);
read "../robot_codegen_definitions/robot_env_IC":
Trf := Trf:
Trf_c := Trf_c:
Trf
;
# Fourbar1 0-3-4-5-6-2-1-0
# Schleife 0-1-2-4
T_0_5 := combine(Matrix(Trf(1..4,1..4, 3)) . Matrix(Trf(1..4,1..4,4)). Matrix(Trf(1..4,1..4, 5)));
# Schleife 0-3-5
T_0_6:= combine(Matrix(Trf(1..4,1..4, 1)) . Matrix(Trf(1..4,1..4,2)) . Matrix(Trf(1..4,1..4,6)));
h1t := T_0_6(1..3,4) - T_0_5(1..3,4);
tmp := Transpose( Matrix(T_0_6(1..3,1..3)) ) . Matrix(T_0_5(1..3,1..3));  # nur anzeigen lassen für h1r
;
 combine(tmp); # nur anzeigen lassen für h1r
;
h1r:=qJ1(t) + qJ2(t) - qJ3(t) - qJ4(t) - qJ5(t);
# Zusammenstellen aller Zwangsbedingungen
implconstr_t := <h1t([1, 2]);h1r>;
implconstr_s := convert_t_s(implconstr_t):
# Exportiere Code für folgende Skripte
kin_constraints_exist:=true:
save implconstr_t, implconstr_s, kin_constraints_exist, sprintf("../codeexport/%s/tmp/kinematic_constraints_implicit_maple.m", robot_name):
# Exportieren des vollständigen Ausdruckes
if codegen_act then
  MatlabExport(implconstr_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_matlab.m", robot_name), 2):
end if:
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( implconstr_s )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_matlab.m", robot_name), 2):

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


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

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

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

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

$define gga_x_pbe_params
$include "gga_x_pbe.mpl"

m05_f := (x, u, t) ->
  + params_a_csi_HF*pbe_f(x)*mgga_series_w(params_a_a, 12, t):

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


text:
$ifndef _CLASSIFY_HOLDER_
$define _CLASSIFY_HOLDER_

$include "Utils.mpl"
ClassifyHolder:=module()
    local   ieqCode,sols,unsolvedSols;
    export  reset,      # 重置状态
            addSol,     # 新增解
            getSols,    # 获取解
            addUnsolvedSol,
            getUnsolvedSols,
            getIeqCode; # 获取不变量方程的编号

    reset:=proc()
        ieqCode:=0;
        sols:={};
        unsolvedSols:={};
        return;
    end proc:

    addSol:=proc(s::InvSol)
        flogf[1]("添加代表元");
        flog[1](s:-rep);
        sols:=sols union {s};
        return;
    end proc:

    getSols:=proc()
        local res:=sols;
        res:=uniqueObj(res,InvSol:-uniqueKey);
        res:=sort(res,'key'=(x->x:-ieqCode));
        return res;
    end proc:

    getIeqCode:=proc()
        ieqCode:=ieqCode+1;
        return ieqCode;
    end proc:

    addUnsolvedSol:=proc(s)
        unsolvedSols:=unsolvedSols union {s};
        return;
    end proc:

    getUnsolvedSols:=proc()
        return unsolvedSols;
    end proc:

end module:

$endif


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

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

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

mbeta  := 15.75592*0.004235: (* the usual value of 0.066726 *)
malpha := mbeta/(16*2^(2/3)):

sogga11_yy := (rs, z, xt) -> -malpha*mphi(z)*xt^2/(rs*f_pw(rs, z)):

sogga11_f0 := (rs, z, xt) -> 1 - 1/(1 + sogga11_yy(rs, z, xt)):
sogga11_f1 := (rs, z, xt) -> 1 - exp(-sogga11_yy(rs, z, xt)):

sogga11_t0 := (rs, z, xt) -> add(params_a_sogga11_a[i]*sogga11_f0(rs, z, xt)^(i-1), i=1..6):
sogga11_t1 := (rs, z, xt) -> add(params_a_sogga11_b[i]*sogga11_f1(rs, z, xt)^(i-1), i=1..6):

sogga11_f := (rs, z, xt, xs0, xs1) ->
  f_pw(rs, z)*(sogga11_t0(rs, z, xt) + sogga11_t1(rs, z, xt)):

f  := (rs, z, xt, xs0, xs1) -> sogga11_f(rs, z, xt, 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: work_mgga_x *)

(* bottom left column on page 687 of article *)
e1 := -1.6665:
c1 :=  0.7438:

(* equation 10 *)
fa := a -> (1 - a) / ((1 + e1*a^2)^2 + c1*a^4)^(1/4):

(* between eqs 8 and 9 *)
k0 := 0.174:
(* after eq 11 *)
b  := 0.0233:

(* eq 7 *)
f := (rs, x, t, u) -> (1 + k0*fa((t-x^2/8)/K_FACTOR_C)) / (1 + b*(X2S*x)^4)^(1/8):

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

$include "gga_c_gapc.mpl"

(* The two parameters were fixed by fitting to the exact correlation
energy per particle of the He atom *)
gaploc_b  := 14.709046:
gaploc_a1 :=  6.54613 + 2:

gaploc_alpha := t  -> (gaploc_a1 + 3*t^3)/(1 + t^3):
gaploc_s     := xt -> X2S*xt*2^(1/3):

(* override definition of gap_G *)
(* The pre-factor is completely different from the one in
   Equation (7). I used the "gfac" from the original code. *)
gap_G := (rs, z, xt, par) -> (9*Pi/4)^(2/3)/2.0
  * gaploc_s(xt)^(gaploc_alpha(gap_t(rs, z, xt)))/rs^2
  * (gaploc_b + gaploc_s(xt)^2)/(1 + gaploc_s(xt)^(gaploc_alpha(gap_t(rs, z, xt)))):


text:
with(DifferentialAlgebra):
ring_diff := DifferentialRing(blocks = [[EGF_EGFR, EGFR, pEGFR, pAkt_S6, pEGFR_Akt, pAkt, S6, pS6, Akt], [y1, y2, y3] , [pro_EGFR]], derivations = [t]):
sys := [
diff(EGF_EGFR(t), t) - (reaction_1_k1*EGF_EGFR(t) - reaction_9_k1*EGF_EGFR(t) - reaction_1_k2*EGF_EGFR(t)),
diff(EGFR(t), t) - (-reaction_1_k1*EGF_EGFR(t) + reaction_1_k2*EGF_EGFR(t) - EGFR(t)*EGFR_turnover + EGFR_turnover*pro_EGFR(t)),
diff(pEGFR(t), t) - (-reaction_4_k1*pEGFR(t) + reaction_9_k1*EGF_EGFR(t) - pEGFR(t)*Akt(t)*reaction_2_k1 + reaction_3_k1*pEGFR_Akt(t) + pEGFR_Akt(t)*reaction_2_k2),
diff(pAkt_S6(t), t) - (pAkt(t)*reaction_5_k1*S6(t) - reaction_6_k1*pAkt_S6(t) - pAkt_S6(t)*reaction_5_k2),
diff(pEGFR_Akt(t), t) - (pEGFR(t)*Akt(t)*reaction_2_k1 - reaction_3_k1*pEGFR_Akt(t) - pEGFR_Akt(t)*reaction_2_k2),
diff(pAkt(t), t) - (-pAkt(t)*reaction_7_k1 - pAkt(t)*reaction_5_k1*S6(t) + reaction_6_k1*pAkt_S6(t) + reaction_3_k1*pEGFR_Akt(t) + pAkt_S6(t)*reaction_5_k2),
diff(S6(t), t) - (pS6(t)*reaction_8_k1 - pAkt(t)*reaction_5_k1*S6(t) + pAkt_S6(t)*reaction_5_k2),
diff(pS6(t), t) - (-pS6(t)*reaction_8_k1 + reaction_6_k1*pAkt_S6(t)),
diff(Akt(t), t) - (pAkt(t)*reaction_7_k1 - pEGFR(t)*Akt(t)*reaction_2_k1 + pEGFR_Akt(t)*reaction_2_k2),
y1(t) - (pEGFR(t)*a1 + a1*pEGFR_Akt(t)),
y2(t) - (a2*pAkt(t) + a2*pAkt_S6(t)),
y3(t) - (pS6(t)*a3)
];
res := CodeTools[CPUTime](RosenfeldGroebner(sys, ring_diff, singsol=none));


text:
//impulse: /give @p spawn_egg 1 0 {display:{Name:"Portal Spawner"},EntityTag:{id:"minecraft:chicken",CustomName:"MMH_PortalSpawner",Silent:1,NoGravity:1}}
//MMH_loadChunks()
///clone ~1 ~1 ~1 ~3 ~3 ~3 16 1 16

#MMH
repeat process MMH_portalSpawner {
if: /testfor @e[type=Chicken,name=MMH_PortalSpawner]
then {
  MMH_loadChunks()
  /execute @e[type=Chicken,name=MMH_PortalSpawner] ~ ~ ~ scoreboard players operation Rotation MMH_Rotation = @p MMH_Rotation
  /scoreboard players test Rotation MMH_Rotation 0 0
  conditional: /execute @e[type=Chicken,name=MMH_PortalSpawner] ~ ~ ~ clone 16 1 17 16 3 17 ~-2 ~ ~ masked
  /scoreboard players test Rotation MMH_Rotation 1 1
  conditional: /execute @e[type=Chicken,name=MMH_PortalSpawner] ~ ~ ~ clone 17 1 16 17 3 16 ~ ~ ~-2 masked
  /scoreboard players test Rotation MMH_Rotation 2 2
  conditional: /execute @e[type=Chicken,name=MMH_PortalSpawner] ~ ~ ~ clone 18 1 17 18 3 17 ~2 ~ ~ masked
  /scoreboard players test Rotation MMH_Rotation 3 3
  conditional: /execute @e[type=Chicken,name=MMH_PortalSpawner] ~ ~ ~ clone 17 1 18 17 3 18 ~ ~ ~2 masked
  /tp @e[type=Chicken,name=MMH_PortalSpawner] ~ -100 ~
}
}

//impulse: /summon area_effect_cloud ~2 ~.1 ~ {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
///execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
//
//impulse: /summon area_effect_cloud ~ ~.1 ~2 {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
///execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
//
//impulse: /summon area_effect_cloud ~-2 ~.1 ~ {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
///execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block
//
//impulse: /summon area_effect_cloud ~ ~.1 ~-2 {CustomName:"ACV_remote_createPortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
///execute @e[name=ACV_remote_createPortals] ~ ~ ~ setblock ~ ~ ~ redstone_block


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

Commands := [ `Simplify`, `Disintegrate`, `Summarize`, `Reparam` ];

RoundTrip := proc(e, t::t_type, {_ret_type := {'print'}
                                ,_postproc := [['convert','Partition','piecewise'],
                                               'stable_Partition_order']
                                ,_gensym_charset := FAIL
                                ,_debug :: truefalse := false
                                ,_command := Simplify})
  local result, ret_type, cs, ifc_opts, syms0, command, pp;
  ret_type := _ret_type; command := _command;
  if not(ret_type::set) then ret_type := {ret_type}; end if;
  if command::string then command := convert(command, name); end if;
  if not (command in Commands) then
    error "unrecognized command: %1", command;
  end if;
  if _gensym_charset <> FAIL then
    syms0 := copy(gensym:-SymbolsToGen);
    gensym:-SymbolsToGen := _gensym_charset;
  else
    syms0 := gensym:-SymbolsToGen;
  end if;

  ifc_opts[0] :=
     screenwidth=9999,
     `if`('print' in ret_type, prettyprint=0, NULL),
     warnlevel=0,
     showassumed=0,
     quiet=true;

  ifc_opts[1] := interface(ifc_opts[0]);
  try
    if not(_debug) then kernelopts(assertlevel=0); end if;
    result := _command(e,t,_rest);
    for pp in _postproc do
      if assigned(RoundTrip_postproc[op(1,pp)]) then
        result := RoundTrip_postproc[op(1,pp)](result, op(2..-1,pp))
      else
        error "invalid postproc option: %1", pp;
      end if;
    end do;
    if not ('expr' in ret_type) then
      result := eval(ToInert(result), _Inert_ATTRIBUTE=NULL);
      if ({'print','string'} intersect ret_type) <> {} then
        result := sprintf("%a",result);
      end if;
      if 'print' in ret_type then
        printf("%s\n",result);
        result := NULL;
      end if;
    end if;
  finally
    interface(zip(`=`, map(lhs,[ifc_opts[0]]), [ifc_opts[1]])[]);
    gensym:-SymbolsToGen := syms0;
  end try;
  return result;
end proc;

RoundTrip_postproc[convert] := proc(r, f::type, t, $)
  subsindets(r, f, x->convert(x,t));
end proc;
RoundTrip_postproc[stable_Partition_order] := proc(r, $)
  subsindets(r, Partition, Simpl:-stable_order)
end proc;

Reparam := proc(e, t::t_type, {ctx :: list := []})
  SimplifyKB_(
    ((e1,_,kb1)->
     fromLO(ReparamDetermined(toLO(e1),kb1),_ctx=kb1)),
    e,t,build_kb(ctx,"Reparam",KB:-empty))
end proc;

Summarize := proc(e, _) Summary:-Summarize(e,_rest) end proc;

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

Rename := proc(e, {_do_rename :: truefalse := false})
  local ns, i;
  if _do_rename then
    ns := [op(bound_names_in(res))];
    subs(zip(`=`,ns,
             [seq(cat(`x`,convert(i,string)),i=1..nops(ns))]),
         e);
  else e
  end if;
end proc;

Disintegrate := proc(e::t_Hakaru,t::t_type, {ctx::list := []})
  local res;
  res := SimplifyKB_((e,t,kb)->disint(e,disint:-htype_to_disint_wrt(t),kb,'do_lam',_rest),
                     e, t, build_kb(ctx,"Disintegrate",KB:-empty));
  Rename(res,_rest);
end proc;

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

SimplifyKB := proc(e,t,kb,$)
  local res;
  res[0] := SimplifyKB_(
  proc(e0,t,kb)
    local e := eval_for_Simplify(e0,kb);
    if t :: HMeasure(anything) then
      fromLO(improve(toLO(e), _ctx=kb), _ctx=kb);
    else
      simplify_assuming(e,kb);
    end if;
  end proc, e, t, kb);

  res[1] := eval(res[0]);
  if res[1] <> res[0] then
    KB:-simplify_assuming(res[1], kb);
  else
    res[0]
  end if;
end proc;

eval_for_Simplify_tbl := table(
  [ `Int`=`int`
  # Some parts (which?) of the simplifier still always need `sum` instead of `Sum`
  , `Sum`=[`sum`,Not(anything)]
  # Some parts (e.g. gmm_gibbs, see #103) of the simplifier still need `product`
  , `Product`=[`product`,Not(satisfies(x->has(x,{Hakaru:-idx,Hakaru:-size})))]
  ]);

eval_for_Simplify := proc(e,kb,$)
  eval_in_ctx(
    proc(x,kb1) local q;
      subsindets(x, {seq(specfunc(q),q=[indices(eval_for_Simplify_tbl, nolist)])},
                 y -> KB:-kb_eval_mb(eval_for_Simplify_tbl[op(0,y)],y,kb1));
    end proc, e, kb);
end proc;

eval_in_ctx_tbl := table(
  [ `PARTITION` = ((e,kb,ev)->KB:-kb_Partition(e,kb,a->a,ev))
  , `piecewise` = ((e,kb,ev)->KB:-kb_piecewise(e,kb,a->a,ev,'no_split_disj'))
  ]);

eval_in_ctx := proc(ev, e, kb, $)
  if assigned(eval_in_ctx_tbl[op(0,e)]) then
    eval_in_ctx_tbl[op(0,e)](e,kb,curry(eval_in_ctx, ev));
  elif has(e, {indices(eval_in_ctx_tbl,nolist)}) then
    ev(map(x->eval_in_ctx(ev,x,kb),e),kb);
  else
    ev(e,kb);
  end if;
end proc;

PrintVersion := proc($)
  local i,v,s;
  s := max(op(map(length@lhs,Hakaru:-Version)));
  for i in map(lhs,Hakaru:-Version) do
    v := table(Hakaru:-Version)[i];
    printf(cat("%-",s,"s : %s\n"), i, `if`(v::string,v,""));
  end do;
end proc;

# Testing

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

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

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

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

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

TestEfficient := proc(e, t::t_type, t_kb := KB:-empty, {label::string := "Test" })
  local out, result;
  out := SimplifyKB(args[1..3]);

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

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

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


text:
  # Like simplify_assuming, but does a lot of hack-y things
  simplify_factor_assuming := module ()

    export ModuleApply;
    local graft_pw, GAMMAratio, wrap, and_info,
          hack_Beta_pw, hack_Beta, hackier_Beta,
          eval_piecewise, bounds_are_simple, eval_loop, eval_factor;

$include "NewSLO/Beta.mpl"
$include "NewSLO/Piecewise.mpl"

    bounds_are_simple := proc(bounds :: range, $)
      local bound, term, i;
      if `-`(op(bounds)) :: integer then return true end if;
      for bound in bounds do
        for term in convert(bound, 'list', `+`) do
          if term :: `*` then term := remove(type, term, integer) end if;
          if term :: 'specfunc(piecewise)' then
            for i from 2 by 2 to nops(term)-1 do
              if not(op(i,term) :: integer) then return false end if
            end do;
            if not(op(-1,term) :: integer) then return false end if
          elif not(term :: integer) then return false end if
        end do
      end do;
      true
    end proc;

    eval_loop := proc(e :: And(specfunc({product,Product,sum,Sum}),
                               anyfunc(anything, name=range)),
                      kb :: t_kb,
                      mode :: identical(`*`,`+`),
                      loops, $)
      local o, body, x, bounds, res, go, y, kb1;
      o := op(0,e);
      body, bounds := op(e);
      x, bounds := op(bounds);
      bounds := map(eval_factor, bounds, kb, `+`, []);
      if bounds_are_simple(bounds) then
        # Expand {product,sum} to {mul,add}
        if o=Product then o:=product elif o=Sum then o:=sum end if;
        bounds := lift_piecewise(bounds,
                    'And(range, Not(range(Not(specfunc(piecewise)))))');
        res := `if`(bounds :: 'specfunc(piecewise)', map_piecewiselike, apply)
                   ((b -> o(go(x), x=b)), bounds);
        # Work around this bug:
        # > product(Sum(f(j),j=0..n-1),j=0..0)
        # Sum(f(0),0=0..n-1)
        res := eval(res, go = (i -> eval(body, x=i)));
        return eval_factor(res, kb, mode, loops);
      end if;
      y, kb1 := genType(x, HInt(closed_bounds(bounds)), kb);
      return eval_factor(subs(x=y,body), kb1, mode, [[o,y=bounds],op(loops)]);
    end proc;

    # eval_factor is a simplifier.  It maintains the following invariants:
    #   eval_factor(e, kb, mode, []) = e
    #   eval_factor(e, kb, `*` , [...,[product,i=lo..hi]])
    #     = product(eval_factor(e, kb, `*`, [...]), i=lo..hi)
    #   eval_factor(e, kb, `+` , [...,[sum    ,i=lo..hi]])
    #     = sum    (eval_factor(e, kb, `+`, [...]), i=lo..hi)
    # It recursively traverses e while "remembering" the loops traversed,
    # as ``context'' information, to help with transformations.
    eval_factor := proc(e, kb :: t_kb, mode :: identical(`*`,`+`), loops, $)
      local res, i, j, k, s, r;
      if not (loops :: 'list'([`if`(mode=`*`, 'identical(product,Product)',
                                              'identical(sum    ,Sum    )'),
                               'name=range'])) then
        error "invalid input: mode=%1, loops=%2", mode, loops;
      end if;
      if e :: mode then
        # Transform product(a*b,...) to product(a,...)*product(b,...)
        # (where e=a*b and loops=[[product,...]])
        return map(eval_factor, e, kb, mode, loops);
      end if;
      if e = mode () then
        return e;
      end if;
      if e :: And(specfunc(`if`(mode=`*`, '{product,Product}', '{sum,Sum}')),
                  'anyfunc(anything, name=range)') then
        return eval_loop(e, kb, mode, loops);
      end if;
      if e :: 'specfunc(piecewise)' then
        return eval_piecewise(e, kb, mode, loops);
      end if;
      if e :: `*` then
        # If we're here, then mode=`+` (else "e :: mode" above would be true)
        s, r := selectremove(depends, e, map2(op,[2,1],loops));
        # Transform sum(a*b(i),i=...) to a*sum(b(i),i=...)
        if r <> 1 then
          return eval_factor(s, kb, `+`, loops)
               * maptype(`*`, eval_factor, r, kb, `+`, []);
        end if;
      end if;
      if mode = `*` then
        i := map2(op,[2,1],loops);
        if e :: '`^`' then
          # Transform product(a(i)^b,i=...) to product(a(i),i=...)^b
          if not depends(op(2,e), i) then
            return eval_factor(op(1,e), kb, `*`, loops)
                 ^ eval_factor(op(2,e), kb, `+`, []);
          end if;
        end if;
        if e :: 'exp(anything)' or e :: '`^`' and not depends(op(1,e), i) then
          # Transform product(a^((b(i)+c(i))^2),i=...)
          #        to a^   sum(b(i)^2   ,i=...)
          #         * a^(2*sum(b(i)*c(i),i=...))
          #         * a^   sum(c(i)^2   ,i=...)
          return mul(subsop(-1=j,e),
                     j in convert(eval_factor(expand(op(-1,e)), kb, `+`,
                                              map2(subsop,1=sum,loops)),
                                  'list', `+`));
        end if;
        # Rewrite ... * idx([p,1-p],i)
        #      to ... * p^idx([1,0],i) * (1-p)^idx([0,1],i)
        # because the latter is easier to integrate and recognize with respect to p
        if e :: 'idx(list, anything)' and not depends(op(1,e), i) then
          return mul(op([1,j],e)
                     ^ eval_factor(idx([seq(`if`(k=j,1,0), k=1..nops(op(1,e)))],
                                       op(2,e)),
                                   kb, `+`, map2(subsop,1=sum,loops)),
                     j=1..nops(op(1,e)));
        end if;
      end if;

      # Try not to use hack_beta for now...

      if mode = `*` and e :: 'specfunc(Beta)' then
        res := hack_Beta(e, kb, loops);
        if res <> FAIL then
          userinfo(3, 'hack_Beta',
                   printf("Beta hack was applied to %a, %a, %a\n"
                          "\tto produce %a\n", e, kb, loops, res));
          return res
        end if;
      end if;

      if nops(loops) > 0 then
        # Emit outermost loop as last resort
        return op([-1,1],loops)(eval_factor(e, kb, mode, loops[1..-2]),
                                op([-1,2],loops));
      end if;
      # In the remainder of this function, assume loops=[] and recur in
      if e :: '{specfunc({GAMMA, Beta, exp, And, Or, Not}), relation, logical}' then
        return map(eval_factor, e, kb, `+`, []);
      end if;
      if e :: `^` then
        return eval_factor(op(1,e), kb, mode, [])
             ^ eval_factor(op(2,e), kb, `+` , []);
      end if;
      if e :: `+` then
        return map(eval_factor, e, kb, mode, []);
      end if;
      return e;
    end proc;

    ModuleApply := proc(e, kb::t_kb, $)
      simplify_assuming(eval_factor(convert(Partition:-PartitionToPW_mb(e), 'Beta'), kb, `*`, []), kb);
    end proc;
  end module; # simplify_factor_assuming


text:
$ifndef _REP_REFINE_
$define _REP_REFINE_

$endif


text:
kernelopts(printbytes=false);
with(RigidMotionsRecoverNMM);
LaunchComputeNMM([a,b,c], "N2", [-2, -1, 0, 1, 2], "__REPLACE__");
done

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
`is_element/E/subdiv` := (p::nonnegint,n::nonnegint,m::nonnegint) -> proc(theta)
 global reason;
 local tag,src,tgt,i,t;
 tag := "is_element/E/subdiv";

 if p = 0 then
  return evalb(n = m and theta = []);
 fi;

 if not(type(theta,list(list(nonnegint)))) then
  reason := [tag,"theta is not a list of lists of natural numbers",theta];
  return false;
 fi;

 if nops(theta) <> p then
  reason := [tag,"theta is not a list of length p",theta,p];
  return false;
 fi;

 tgt := NULL;
 for i from 1 to p do
  if theta[i] = [] then
   reason := [tag,"theta[i] is empty",theta,i];
   return false;
  fi;
  t := max(op(theta[i]));
  if {op(theta[i])} <> {seq(j,j=0..t)} then
   reason := [tag,"theta[i] does not give an epimorphism",theta,i,t];
   return false;
  fi;
  tgt := tgt,t;
 od:
 tgt := [tgt];

 src := map(nops,theta) -~ 1;

 if src[1] <> n then
  reason := [tag,"chain does not start with [n]",theta,src,n];
  return false;
 fi;

 if tgt[p] <> m then
  reason := [tag,"chain does not end with [m]",theta,tgt,m];
  return false;
 fi;

 if [op(src),m] <> [n,op(tgt)] then
  reason := [tag,"sources and targets do not match",theta,src,tgt,n,m];
  return false;
 fi;

 return true;
end:

`is_leq/E/subdiv` := NULL;

`list_elements/E/subdiv/aux` := proc(n::nonnegint,m::nonnegint)
 option remember;
 local N,M,E;

 N := {seq(i,i=0..n)};
 M := {seq(i,i=0..m)};
 E := `list_elements/epi`(N,M);
 E := map(f -> [seq(f[i],i=0..n)],E);
 return E;
end:

`list_elements/E/subdiv` := proc(p::nonnegint,n::nonnegint,m::nonnegint)
 option remember;
 local L,k,E0,E1;

 if p = 0 then
  if n = m then
   return [[]];
  else
   return [];
  fi;
 else
  L := NULL;
  for k from m to n do
   E0 := `list_elements/E/subdiv`(p-1,n,k);
   E1 := `list_elements/E/subdiv/aux`(k,m);
   L := L,seq(seq([op(e0),e1],e1 in E1),e0 in E0);
  od:
 fi;
 return [L];
end:

`random_element/E/subdiv` := (p::nonnegint,n::nonnegint,m::nonnegint) -> proc()
 local k,q,f,g;

 if m > n then
  return FAIL;
 fi;

 if p = 0 then
  if m = n then
   return [];
  else
   return FAIL;
  fi;
 elif p = 1 then
  g := `random_element/epi`({seq(j,j=0..n)},{seq(j,j=0..m)})();
  return [[seq(g[j],j=0..n)]];
 else
  k := rand(m..n)();
  q := rand(1..p-1)();
  f := `random_element/E/subdiv`(q,n,k)();
  g := `random_element/E/subdiv`(p-q,k,m)();
  return [op(f),op(g)];
 fi;
end:

`count_elements/E/subdiv` := proc(p::nonnegint,n::nonnegint,m::nonnegint)
 option remember;
 if p = 0 then
  return `if`(n = m,1,0);
 else
  return add(`count_elements/E/subdiv`(p-1,n,k) *
             Stirling2(k+1,m+1) * (m+1)!,k=m..n);
 fi;
end:

`degrees/E/subdiv` := (p,n,m) -> (theta) ->
 [op(map(nops,theta) -~ 1),m];

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

`is_element/EE/subdiv` := (p::nonnegint,n::nonnegint) -> proc(theta)
 global reason;
 local tag,src,tgt,i,m,t,b;
 tag := "is_element/EE/subdiv";

 if p = 0 then
  if theta = [] then
   return true;
  else
   reason := [tag,"for p=0 the only element is []",theta];
   return phi;
  fi;
 fi;

 if not(type(theta,list(list(nonnegint)))) then
  reason := [tag,"theta is not a list of lists of natural numbers",theta];
  return false;
 fi;

 if nops(theta) <> p then
  reason := [tag,"theta is not a list of length p",theta,p];
  return false;
 fi;

 if theta[p] = [] then
  reason := [tag,"theta[p] = []",theta,p];
  return false;
 fi;

 m := max(op(theta[p]));

 b := `is_element/E/subdiv`(p,n,m)(theta);
 if not(b) then
  reason := [tag,"is_element/E/subdiv failed",m,reason];
  return false;
 fi;

 return true;
end:

`is_leq/EE/subdiv` := NULL;

`list_elements/EE/subdiv` := proc(p::nonnegint,n::nonnegint)
 option remember;
 [seq(op(`list_elements/E/subdiv`(p,n,m)),m=0..n)];
end:

`random_element/EE/subdiv` := (p::nonnegint,n::nonnegint) -> proc()
 local m;
 m := rand(0..n)();
 return `random_element/E/subdiv`(p,n,m)();
end:

`count_elements/EE/subdiv` := proc(p::nonnegint,n::nonnegint)
 add(`count_elements/E/subdiv`(p,n,m),m=0..n);
end:

`find_count_coeffs/EE/subdiv` := proc(n)
 local E,a,p,k;
 E := [seq(`count_elements/EE/subdiv`(p,n) - add((-1)^(n+k)*a[n,k]*((k+1)!)^p,k=0..n),p=1..n+5)];
 return solve(E);
end:

`degrees/EE/subdiv` := (p,n) -> (theta) ->
 [op(map(nops,theta) -~ 1),max(op(theta[p]))];

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

`is_element/threads/subdiv` := (p::nonnegint,n::nonnegint) -> (theta) -> proc(i)
 global reason;
 local tag,k,m;

 tag := "is_element/threads/subdiv";

 if not(type(i,list(nonnegint)) and nops(i) = p+1) then
  reason := [tag,"i is not a list of p+1 natural numbers",i,p];
  return false;
 fi;

 k := `degrees/EE/subdiv`(p,n)(theta);

 for m from 0 to p do
  if i[m+1] > k[m+1] then
   return false;
  fi;

  if m < p and theta[m+1][i[m+1]+1] > i[m+2] then
   return false;
  fi;
 od:

 return true;
end:

`is_leq/threads/subdiv` := NULL;

`list_elements/threads/subdiv` := (p::nonnegint,n::nonnegint) -> proc(theta)
 local i,k,m,phi,L,M;

 if p = 0 then
  return [seq([i],i=0..n)];
 else
  m := max(op(theta[p]));
  phi := [op(1..p-1,theta)];
  L := `list_elements/threads/subdiv`(p-1,n)(phi);
  M := NULL;
  for i in L do
   k := theta[p][i[p]+1];
   M := M,seq([op(i),j],j=k..m);
  od:
  return [M];
 fi;
end:

`random_element/threads/subdiv` := (p::nonnegint,n::nonnegint) -> (theta) -> proc()
 local i,j,k,m0,m1;

 k := `degrees/EE/subdiv`(p,n)(theta);
 i := [rand(0..n)()];
 for j from 1 to p do
  m0 := theta[j][i[j]+1];
  m1 := k[j+1];
  i := [op(i),rand(m0..m1)()];
 od;

 return i;
end:

`count_elements/threads/subdiv` := NULL;

`mu/threads/subdiv` := (p::nonnegint,n::nonnegint) -> (theta) -> proc(i)
 local k,u,m;

 k := `degrees/EE/subdiv`(p,n)(theta);

 u := k[p+1]+1;
 for m from 0 to p-1 do
  u := u * nops(select(j -> (j <= i[m+2]),theta[m+1]));
 od:

 return u;
end:

`xi/EE/subdiv` := (p::nonnegint,n::nonnegint) -> proc(theta)
 local x,II,i;

 if p = 0 then
  return [1$(n+1)]/~(n+1);
 fi;

 x := table([seq(i=0,i=0..n)]);
 II := `list_elements/threads/subdiv`(p,n)(theta);
 for i in II do
  x[i[1]] := x[i[1]] + 1/`mu/threads/subdiv`(p,n)(theta)(i);
 od:
 return [seq(x[i],i=0..n)];
end:

######################################################################
# F(p,n) is the set of vertices in the (p+1)-fold barycentric
# subdivision of the n-simplex

`is_element/FF/subdiv` := (p::nonnegint,n::nonnegint) -> proc(alpha_theta)
 local alpha,theta,m;

 if not(type(alpha_theta,list) and nops(alpha_theta) = 2) then
  return false;
 fi;

 alpha,theta := op(alpha_theta);

 if not(type(alpha,list(nonnegint)) and nops(alpha) > 0) then
  return false;
 fi;

 m := nops(alpha) - 1;

 if not(`is_element/simplicial_mono_alt`(m,n)(alpha)) then
  return false;
 fi;

 if not(`is_element/EE/subdiv`(p,m)(theta)) then
  return false;
 fi;

 return true;
end:

`list_elements/FF/subdiv` := proc(p,n)
 local L,N,m,alpha,theta;

 L := NULL:
 N := [seq(i,i=0..n)];
 for alpha in combinat[powerset](N) do
  m := nops(alpha) - 1;
  if m >= 0 then
   L := L,seq([alpha,theta],theta in `list_elements/EE/subdiv`(p,m));
  fi;
 od:

 return [L];
end:

`random_element/FF/subdiv` := (p,n) -> proc()
 local d,m,alpha,theta;

 d := rand(1..n)();
 alpha := `random_subset_of`({seq(i,i=0..n)},d);
 alpha := sort([op(alpha)]);
 m := nops(alpha) - 1;
 theta := `random_element/EE/subdiv`(p,m)();
 return [alpha,theta];
end:

# This function is not quite right.  It seems to find a correct morphism
# when one exists, but it claims that morphisms exist when that is not the case.

`find_morphism/FF/subdiv` := (p::nonnegint,n::nonnegint) -> proc(alpha_theta0,alpha_theta1)
 local alpha0,alpha1,theta0,theta1,m0,m1,k0,k1,f,h,i,j,u,v,xi,V,W,G,H;

 alpha0,theta0 := op(alpha_theta0);
 alpha1,theta1 := op(alpha_theta1);
 m0 := nops(alpha0) - 1;
 m1 := nops(alpha1) - 1;
 k0 := `degrees/EE/subdiv`(p,m0)(theta0);
 k1 := `degrees/EE/subdiv`(p,m1)(theta1);

 f := table([seq(i=FAIL,i=0..n)]);
 for i from 0 to m1 do f[alpha1[i+1]] := i; od;
 xi := [[seq(f[alpha0[i+1]],i=0..m0)]];
 if member(FAIL,xi[1]) then return FAIL; fi;
 for j from 0 to p-1 do
  V := {seq(w,w=0..k0[j+1])};
  W := {seq(w,w=0..k0[j+1+1])};
  G := [seq(select(v -> theta0[j+1][v+1] <= w,V),w in W)];
  H := map(A -> map(v -> theta1[j+1][xi[j+1][v+1]+1],A),G);
  h := map(max,H);
  if h <> sort([op({op(h)})]) then
   return false;
  fi;
  xi := [op(xi),h];
 od;

 return xi;
end:

`is_leq/FF/subdiv` := (p::nonnegint,n::nonnegint) -> proc(alpha_theta0,alpha_theta1)
 evalb(`find_morphism/FF/subdiv`(p,n)(alpha_theta0,alpha_theta1) <> FAIL);
end:

`xi/FF/subdiv` := (p::nonnegint,n::nonnegint) -> proc(alpha_theta)
 local alpha,theta,i,m,x,y,yt;

 alpha,theta := op(alpha_theta);
 m := nops(alpha) - 1;
 x := `xi/EE/subdiv`(p,m)(theta);
 yt := table([seq(i=0,i=0..n)]);
 for i from 0 to m do
  yt[alpha[i+1]] := x[i+1];
 od:
 y := [seq(yt[i],i=0..n)];
 return y;
end:

# Every vertex of the (p+1)-fold subdivision is the barycentre of
# a d-simplex of the p-fold subdivision for some d.  This function
# returns d.

`rank/FF/subdiv` := (p::nonnegint,n::nonnegint) -> proc(alpha_theta)
 local alpha,theta,m;

 alpha,theta := op(alpha_theta);
 m := nops(alpha) - 1;

 if p = 0 then
  return m;
 else
  return max(op(theta[p]));
 fi;
end:

######################################################################
# Here zeta is a strictly increasing map [j] -> [k] and theta is an
# unordered surjective map [m] -> [k].  The function returns a triple
# [n,beta,phi] where beta : [n] -> [m] is strictly increasing and
# phi : [n] -> [j] is surjective.  In an appropriate context we will
# have alpha^*([x,theta]) = [beta^*(x),phi].

`rewrite/epi/subdiv` := proc(zeta,theta)
 local i,j,m,J,M,k1,N,n,xi,phi;

 j := nops(zeta) - 1;
 m := nops(theta) - 1;
 J := {seq(j0,j0=0..j)};
 M := {seq(m0,m0=0..m)};

 k1 := op(-1,zeta);
 N := select(m0 -> theta[m0+1] <= k1,M);
 n := nops(N) - 1;
 xi := sort([op(N)]);
 phi := NULL;

 for i from 0 to n do
  phi := phi,min(select(j0 -> zeta[j0+1] >= theta[xi[i+1]+1],J));
 od:
 phi := [phi];

 return [n,xi,phi];
end:

`rewrite/FF/subdiv` := (p,n) -> (zeta) -> proc(alpha_theta)
 local alpha,theta,k0,m,xi,phi,alpha1,theta1;
 alpha,theta := op(alpha_theta);

 if p = 0 then
  k0 := nops(zeta) - 1;
  return [[seq(alpha[zeta[i+1]+1],i=0..k0)],[]];
 else
  k0 := nops(zeta) - 1;
  m,xi,phi := op(`rewrite/epi/subdiv`(zeta,theta[p]));
  alpha1,theta1 := op(`rewrite/FF/subdiv`(p-1,n)(xi)([alpha,[op(1..-2,theta)]]));
  return [alpha1,[op(theta1),phi]];
 fi;
end:

`rewrite_aux/FF/subdiv` := (p,n) -> (zeta) -> proc(alpha_theta)
 local alpha,theta,k0,m,xi,phi,xi1,alpha1,theta1;
 alpha,theta := op(alpha_theta);

 if p = 0 then
  k0 := nops(zeta) - 1;
  return [[zeta],[seq(alpha[zeta[i+1]+1],i=0..k0)],[]];
 else
  k0 := nops(zeta) - 1;
  m,xi,phi := op(`rewrite/epi/subdiv`(zeta,theta[p]));
  xi1,alpha1,theta1 := op(`rewrite_aux/FF/subdiv`(p-1,n)(xi)([alpha,[op(1..-2,theta)]]));
  return [[op(xi1),zeta],alpha1,[op(theta1),phi]];
 fi;
end:

`vertices/FF/subdiv` := (p,n) -> proc(alpha_theta)
 local r;
 r := `rank/FF/subdiv`(p,n)(alpha_theta);
 return [seq(`rewrite/FF/subdiv`(p,n)([i])(alpha_theta),i=0..r)];
end:


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

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

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

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

gds08_fs := (rs, z) -> (1 + z)/2 *(
  + params_a_A
  + params_a_B*log(2*n_spin(rs, z))
  + params_a_C*log(2*n_spin(rs, z))^2
):

# Eq. (12)
gds08_f := (rs, z) ->
  + my_piecewise3(screen_dens(rs,  z), 0, gds08_fs(rs, z_thr( z)))
  + my_piecewise3(screen_dens(rs, -z), 0, gds08_fs(rs, z_thr(-z))):

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


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

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

(* type: work_gga_x *)

lambda := y -> 1/2*(1 + (1 - y^2)*log((1 + y)/abs(1 - y))/(2*y)):

f := x -> 1 + lambda(X2S*x/6)*x^2/(8*K_FACTOR_C):


text:
# 字符串编码
str_encode:=s->sprintf("'%s'",HTTP:-URLEncode(s));

# 文件计数器
FileCounter:=module()
    local k:=0;
    export next_file,curr_file;

    curr_file:=proc()
        return sprintf("imgs/fig%03d",k);
    end proc:

    next_file:=proc()
        k:=k+1;
        return sprintf("imgs/fig%03d",k);
    end proc:
end module:

# 打开 MATLAB 并设置路径
Matlab:-openlink();
Matlab:-evalM(
    sprintf("cd(char(java.net.URLDecoder.decode(%s,'UTF8')))",str_encode(currentdir()))
);

# 处理坐标名称
label_str:=x->str_encode(`if`(type(x,string),x,latex(x,output=string)));

# 3D曲面绘图接口
plot3dm:=proc(_eq,{x::range:=0..1,y::range:=0..1,title::string:="",grid::list:=[100,100],labels::list:=[x,y,z],extcmd::string:="",
        view::{range,list(range)}:=[x,y,-infinity..infinity]})
    uses ImageTools;
    local eq,_view,plot3d_args;
    if type(view,range) then
        _view:=[x,y,view];
    else
        _view:=view;
    end if;
    eq:=simplify(_eq);
    plot3d_args:=remove(x->type(x,equation) and lhs(x)=:-extcmd,[_passed]);
    print(old_plot3d(plot3d_args[]));
    Matlab:-evalM(sprintf("plot3dm(%s,{%f,%f,%d},{%f,%f,%d},[%f,%f,%f,%f,%f,%f],%s,{%s,%s,%s},%s,%s)",
        str_encode( sprintf("%a",eq) ),
        op(1,x),op(2,x),grid[1],
        op(1,y),op(2,y),grid[2],
        map(t->(op(1,t),op(2,t)),_view)[],
        str_encode(title),
        label_str(labels[1]),label_str(labels[2]),label_str(labels[3]),
        str_encode(extcmd),
        str_encode(FileCounter:-next_file())
    ));
    Embed(Read(sprintf("%s.png",FileCounter:-curr_file())));
end proc:

# 重载plot3d函数
old_plot3d:=copy(plot3d,deep);
unprotect(plot3d);
plot3d:=plot3dm;