text:
HagerZhang := module()

  description "HagerZhang";

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

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

  local ModuleLoad,
        ModuleUnLoad,
        _debug,
        FINAL, ITER, PARAMETERS, GRADIENT, SEARCHDIR,
        ALPHA, BETA, BRACKET, LINESEARCH, UPDATE, SECANT2,
        BISECT, BARRIERCOEF, DEFAULTDELTA, DEFAULTSIGMA,
        max_iter, A_list, B_list;

  uses LinearAlgebra;

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

  ModuleLoad := proc()
    printf("HagerZhang load\n");
    _debug   := true;
    # Display flags are represented as a bitfield
    FINAL       := 1;
    ITER        := 2;
    PARAMETERS  := 8;
    GRADIENT    := 16;
    SEARCHDIR   := 32;
    ALPHA       := 64;
    BETA        := 128;
    BRACKET     := 256;
    LINESEARCH  := 512;
    UPDATE      := 1024;
    SECANT2     := 2048;
    BISECT      := 4096;
    BARRIERCOEF := 8192;

    DEFAULTDELTA := 0.1; # Values taken from HZ paper (Nocedal & Wright recommends 0.01?)
    DEFAULTSIGMA := 0.9; # Values taken from HZ paper (Nocedal & Wright recommends 0.1 for GradientDescent)

    max_iter    := 40;
    NULL;
  end proc;

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

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

  min_quadratic := proc ( a, fa, Dfa, b, fb )
    local t, delta; # A+B*(x-a)+C*(x-a)^2
    # p  = fa + Dfa*(x-a)+t*(x-a)^2
    # p' = Dfa+2*t*(x-a)
    ##t     := ((fb-fa)/(b-a)-Dfa)/(b-a);
    ##delta := -Dfa/(2*t);
    t     := (fb-fa)/(b-a)-Dfa;
    delta := max(tau1,min(1-tau1,-Dfa/(2*t)));
    a+(b-a)*delta;
  end proc;

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

  #
  # Conjugate gradient line search implementation from:
  #   W. W. Hager and H. Zhang (2006) Algorithm 851: CG_DESCENT, a
  #     conjugate gradient method with guaranteed descent. ACM
  #     Transactions on Mathematical Software 32: 113–137.
  #
  # Code comments such as "HZ, stage X" or "HZ, eqs Y" are with
  # reference to a particular point in this paper.
  #
  # There are some modifications and/or extensions from what's in the
  # paper (these may or may not be extensions of the cg_descent code
  # that can be downloaded from Hager's site; his code has undergone
  # numerous revisions since publication of the paper):
  # linesearch: the Wolfe conditions are checked only after alpha is
  #   generated either by quadratic interpolation or secant
  #   interpolation, not when alpha is generated by bisection or
  #   expansion. This increases the likelihood that alpha will be a
  #   good approximation of the minimum.
  #
  # linesearch: In step I2, we multiply by psi2 only if the convexity
  #   test failed, not if the function-value test failed. This
  #   prevents one from going uphill further when you already know
  #   you're already higher than the point at alpha=0.
  #
  # both: checks for Inf/NaN function values
  #
  # both: support maximum value of alpha (equivalently, c). This
  #   facilitates using these routines for constrained minimization
  #   when you can calculate the distance along the path to the
  #   disallowed region. (When you can't easily calculate that
  #   distance, it can still be handled by returning Inf/NaN for
  #   exterior points. It's just more efficient if you know the
  #   maximum, because you don't have to test values that won't
  #   work.) The maximum should be specified as the largest value for
  #   which a finite value will be returned.  See, e.g., limits_box
  #   below.  The default value for alphamax is Inf. See alphamaxfunc
  #   for cgdescent and alphamax for HagerZhang.

  #const DEFAULTDELTA = 0.1 # Values taken from HZ paper (Nocedal & Wright recommends 0.01?)
  #const DEFAULTSIGMA = 0.9 # Values taken from HZ paper (Nocedal & Wright recommends 0.1 for GradientDescent)

  # NOTE:
  #   [1] The type `T` in the `HagerZhang{T}` need not be the same `T` as in
  #       `hagerzhang!{T}`; in the latter, `T` comes from the input vector `x`.
  #   [2] the only method parameter that is not included in the
  #       type is `iterfinitemax` since this value needs to be
  #       inferred from the input vector `x` and not from the type information
  #       on the parameters
  #
  # Conjugate gradient line search implementation from:
  # W. W. Hager and H. Zhang (2006) Algorithm 851: CG_DESCENT, a
  # conjugate gradient method with guaranteed descent. ACM
  # Transactions on Mathematical Software 32: 113–137.

  # TODO: Should we deprecate the interface that only uses the phi and phid\phi arguments?
  linesearch := proc( phi, dphi, c, phi_0, dphi_0 )
    local delta, sigma, alphamax, rho, epsilon, gamma,
          linesearchmax, psi3, display, mayterminate;

    if not (isfinite(phi_0) and isfinite(dphi_0)) then
       error("Value and slope at step length = 0 must be finite.");
    end;
    if dphi_0 >= eps * abs(phi_0) then
      error("Search direction is not a descent direction.");
    elif dphi_0 >= 0 then
      return 0, phi_0
    end

    # Prevent values of x_new = x+αs that are likely to make
    # phi(x_new) infinite
    iterfinitemax::Int = ceil(Int, -log2(eps(T)))
    alphas = [0] # for bisection
    values = [phi_0]
    slopes = [dphi_0]
    if display and LINESEARCH > 0 then
      println("New linesearch\n");
    end

    phi_lim = phi_0 + epsilon * abs(phi_0);
    @assert c >= 0
    c <= eps(T) && return zeroT, phi_0
    @assert isfinite(c) && c <= alphamax
    phi_c, dphi_c = phidphi(c)
    iterfinite = 1
    while !(isfinite(phi_c) && isfinite(dphi_c)) && iterfinite < iterfinitemax
        mayterminate[] = false
        iterfinite += 1
        c *= psi3
        phi_c, dphi_c = phidphi(c)
    end
    if !(isfinite(phi_c) && isfinite(dphi_c))
        @warn("Failed to achieve finite new evaluation point, using alpha=0")
        mayterminate[] = false # reset in case another initial guess is used next
        return zeroT, phi_0
    end
    push!(alphas, c)
    push!(values, phi_c)
    push!(slopes, dphi_c)

    # If c was generated by quadratic interpolation, check whether it
    # satisfies the Wolfe conditions
    if mayterminate[] &&
          satisfies_wolfe(c, phi_c, dphi_c, phi_0, dphi_0, phi_lim, delta, sigma)
        if display & LINESEARCH > 0
            println("Wolfe condition satisfied on point alpha = ", c)
        end
        mayterminate[] = false # reset in case another initial guess is used next
        return c, phi_c # phi_c
    end
    # Initial bracketing step (HZ, stages B0-B3)
    isbracketed = false
    ia = 1
    ib = 2
    @assert length(alphas) == 2
    iter = 1
    cold = -one(T)
    while !isbracketed && iter < linesearchmax
        if display & BRACKET > 0
            println("bracketing: ia = ", ia,
                    ", ib = ", ib,
                    ", c = ", c,
                    ", phi_c = ", phi_c,
                    ", dphi_c = ", dphi_c)
        end
        if dphi_c >= zeroT
            # We've reached the upward slope, so we have b; examine
            # previous values to find a
            ib = length(alphas)
            for i = (ib - 1):-1:1
                if values[i] <= phi_lim
                    ia = i
                    break
                end
            end
            isbracketed = true
        elseif values[end] > phi_lim
            # The value is higher, but the slope is downward, so we must
            # have crested over the peak. Use bisection.
            ib = length(alphas)
            ia = 1
            if c ≉  alphas[ib] || slopes[ib] >= zeroT
                error("c = ", c)
            end
            # ia, ib = bisect(phi, lsr, ia, ib, phi_lim) # TODO: Pass options
            ia, ib = bisect!(phidphi, alphas, values, slopes, ia, ib, phi_lim, display)
            isbracketed = true
        else
            # We'll still going downhill, expand the interval and try again.
            # Reaching this branch means that dphi_c < 0 and phi_c <= phi_0 + ϵ_k
            # So cold = c has a lower objective than phi_0 up to epsilon.
            # This makes it a viable step to return if bracketing fails.

            # Bracketing can fail if no cold < c <= alphamax can be found with finite phi_c and dphi_c.
            # Going back to the loop with c = cold will only result in infinite cycling.
            # So returning (cold, phi_cold) and exiting the line search is the best move.
            cold = c
            phi_cold = phi_c
            if nextfloat(cold) >= alphamax
                mayterminate[] = false # reset in case another initial guess is used next
                return cold, phi_cold
            end
            c *= rho
            if c > alphamax
                c = alphamax
                if display & BRACKET > 0
                    println("bracket: exceeding alphamax, using c = alphamax = ", alphamax,
                    ", cold = ", cold)
                end
            end
            phi_c, dphi_c = phidphi(c)
            iterfinite = 1
            while !(isfinite(phi_c) && isfinite(dphi_c)) && c > nextfloat(cold) && iterfinite < iterfinitemax
                alphamax = c # shrinks alphamax, assumes that steps >= c can never have finite phi_c and dphi_c
                iterfinite += 1
                if display & BRACKET > 0
                    println("bracket: non-finite value, bisection")
                end
                c = (cold + c) / 2
                phi_c, dphi_c = phidphi(c)
            end
            if !(isfinite(phi_c) && isfinite(dphi_c))
                if display & BRACKET > 0
                    println("Warning: failed to expand interval to bracket with finite values. If this happens frequently, check your function and gradient.")
                    println("c = ", c,
                            ", alphamax = ", alphamax,
                            ", phi_c = ", phi_c,
                            ", dphi_c = ", dphi_c)
                end
                return cold, phi_cold
            end
            push!(alphas, c)
            push!(values, phi_c)
            push!(slopes, dphi_c)
        end
        iter += 1
    end
    while iter < linesearchmax
        a = alphas[ia]
        b = alphas[ib]
        @assert b > a
        if display & LINESEARCH > 0
            println("linesearch: ia = ", ia,
                    ", ib = ", ib,
                    ", a = ", a,
                    ", b = ", b,
                    ", phi(a) = ", values[ia],
                    ", phi(b) = ", values[ib])
        end
        if b - a <= eps(b)
            mayterminate[] = false # reset in case another initial guess is used next
            return a, values[ia] # lsr.value[ia]
        end
        iswolfe, iA, iB = secant2!(phidphi, alphas, values, slopes, ia, ib, phi_lim, delta, sigma, display)
        if iswolfe
            mayterminate[] = false # reset in case another initial guess is used next
            return alphas[iA], values[iA] # lsr.value[iA]
        end
        A = alphas[iA]
        B = alphas[iB]
        @assert B > A
        if B - A < gamma * (b - a)
            if display & LINESEARCH > 0
                println("Linesearch: secant succeeded")
            end
            if nextfloat(values[ia]) >= values[ib] && nextfloat(values[iA]) >= values[iB]
                # It's so flat, secant didn't do anything useful, time to quit
                if display & LINESEARCH > 0
                    println("Linesearch: secant suggests it's flat")
                end
                mayterminate[] = false # reset in case another initial guess is used next
                return A, values[iA]
            end
            ia = iA
            ib = iB
        else
            # Secant is converging too slowly, use bisection
            if display & LINESEARCH > 0
                println("Linesearch: secant failed, using bisection")
            end
            c = (A + B) / convert(T, 2)

            phi_c, dphi_c = phidphi(c)
            @assert isfinite(phi_c) && isfinite(dphi_c)
            push!(alphas, c)
            push!(values, phi_c)
            push!(slopes, dphi_c)

            ia, ib = update!(phidphi, alphas, values, slopes, iA, iB, length(alphas), phi_lim, display)
        end
        iter += 1
    end

    throw(LineSearchException("Linesearch failed to converge, reached maximum iterations $(linesearchmax).",
                              alphas[ia]))

end

# Check Wolfe & approximate Wolfe
function satisfies_wolfe(c::T,
                         phi_c::Real,
                         dphi_c::Real,
                         phi_0::Real,
                         dphi_0::Real,
                         phi_lim::Real,
                         delta::Real,
                         sigma::Real) where T<:Number
    wolfe1 = delta * dphi_0 >= (phi_c - phi_0) / c &&
               dphi_c >= sigma * dphi_0
    wolfe2 = (2 * delta - 1) * dphi_0 >= dphi_c >= sigma * dphi_0 &&
               phi_c <= phi_lim
    return wolfe1 || wolfe2
end

# HZ, stages S1-S4
function secant(a::Real, b::Real, dphi_a::Real, dphi_b::Real)
    return (a * dphi_b - b * dphi_a) / (dphi_b - dphi_a)
end
function secant(alphas, values, slopes, ia::Integer, ib::Integer)
    return secant(alphas[ia], alphas[ib], slopes[ia], slopes[ib])
end
# phi
function secant2!(phidphi,
                  alphas,
                  values,
                  slopes,
                  ia::Integer,
                  ib::Integer,
                  phi_lim::Real,
                  delta::Real = DEFAULTDELTA,
                  sigma::Real = DEFAULTSIGMA,
                  display::Integer = 0)
    phi_0 = values[1]
    dphi_0 = slopes[1]
    a = alphas[ia]
    b = alphas[ib]
    dphi_a = slopes[ia]
    dphi_b = slopes[ib]
    T = eltype(slopes)
    zeroT = convert(T, 0)
    if !(dphi_a < zeroT && dphi_b >= zeroT)
        error(string("Search direction is not a direction of descent; ",
                     "this error may indicate that user-provided derivatives are inaccurate. ",
                      @sprintf "(dphi_a = %f; dphi_b = %f)" dphi_a dphi_b))
    end
    c = secant(a, b, dphi_a, dphi_b)
    if display & SECANT2 > 0
        println("secant2: a = ", a, ", b = ", b, ", c = ", c)
    end
    @assert isfinite(c)
    # phi_c = phi(tmpc, c) # Replace
    phi_c, dphi_c = phidphi(c)
    @assert isfinite(phi_c) && isfinite(dphi_c)

    push!(alphas, c)
    push!(values, phi_c)
    push!(slopes, dphi_c)

    ic = length(alphas)
    if satisfies_wolfe(c, phi_c, dphi_c, phi_0, dphi_0, phi_lim, delta, sigma)
        if display & SECANT2 > 0
            println("secant2: first c satisfied Wolfe conditions")
        end
        return true, ic, ic
    end

    iA, iB = update!(phidphi, alphas, values, slopes, ia, ib, ic, phi_lim, display)
    if display & SECANT2 > 0
        println("secant2: iA = ", iA, ", iB = ", iB, ", ic = ", ic)
    end
    a = alphas[iA]
    b = alphas[iB]
    doupdate = false
    if iB == ic
        # we updated b, make sure we also update a
        c = secant(alphas, values, slopes, ib, iB)
    elseif iA == ic
        # we updated a, do it for b too
        c = secant(alphas, values, slopes, ia, iA)
    end
    if (iA == ic || iB == ic) && a <= c <= b
        if display & SECANT2 > 0
            println("secant2: second c = ", c)
        end
        # phi_c = phi(tmpc, c) # TODO: Replace
        phi_c, dphi_c = phidphi(c)
        @assert isfinite(phi_c) && isfinite(dphi_c)

        push!(alphas, c)
        push!(values, phi_c)
        push!(slopes, dphi_c)

        ic = length(alphas)
        # Check arguments here
        if satisfies_wolfe(c, phi_c, dphi_c, phi_0, dphi_0, phi_lim, delta, sigma)
            if display & SECANT2 > 0
                println("secant2: second c satisfied Wolfe conditions")
            end
            return true, ic, ic
        end
        iA, iB = update!(phidphi, alphas, values, slopes, iA, iB, ic, phi_lim, display)
    end
    if display & SECANT2 > 0
        println("secant2 output: a = ", alphas[iA], ", b = ", alphas[iB])
    end
    return false, iA, iB
end

# HZ, stages U0-U3
# Given a third point, pick the best two that retain the bracket
# around the minimum (as defined by HZ, eq. 29)
# b will be the upper bound, and a the lower bound
function update!(phidphi,
                 alphas,
                 values,
                 slopes,
                 ia::Integer,
                 ib::Integer,
                 ic::Integer,
                 phi_lim::Real,
                 display::Integer = 0)
    a = alphas[ia]
    b = alphas[ib]
    T = eltype(slopes)
    zeroT = convert(T, 0)
    # Debugging (HZ, eq. 4.4):
    @assert slopes[ia] < zeroT
    @assert values[ia] <= phi_lim
    @assert slopes[ib] >= zeroT
    @assert b > a
    c = alphas[ic]
    phi_c = values[ic]
    dphi_c = slopes[ic]
    if display & UPDATE > 0
        println("update: ia = ", ia,
                ", a = ", a,
                ", ib = ", ib,
                ", b = ", b,
                ", c = ", c,
                ", phi_c = ", phi_c,
                ", dphi_c = ", dphi_c)
    end
    if c < a || c > b
        return ia, ib #, 0, 0  # it's out of the bracketing interval
    end
    if dphi_c >= zeroT
        return ia, ic #, 0, 0  # replace b with a closer point
    end
    # We know dphi_c < 0. However, phi may not be monotonic between a
    # and c, so check that the value is also smaller than phi_0.  (It's
    # more dangerous to replace a than b, since we're leaving the
    # secure environment of alpha=0; that's why we didn't check this
    # above.)
    if phi_c <= phi_lim
        return ic, ib#, 0, 0  # replace a
    end
    # phi_c is bigger than phi_0, which implies that the minimum
    # lies between a and c. Find it via bisection.
    return bisect!(phidphi, alphas, values, slopes, ia, ic, phi_lim, display)
end

# HZ, stage U3 (with theta=0.5)
function bisect!(phidphi,
                 alphas::AbstractArray{T},
                 values,
                 slopes,
                 ia::Integer,
                 ib::Integer,
                 phi_lim::Real,
                 display::Integer = 0) where T
    gphi = convert(T, NaN)
    a = alphas[ia]
    b = alphas[ib]
    # Debugging (HZ, conditions shown following U3)
    zeroT = convert(T, 0)
    @assert slopes[ia] < zeroT
    @assert values[ia] <= phi_lim
    @assert slopes[ib] < zeroT       # otherwise we wouldn't be here
    @assert values[ib] > phi_lim
    @assert b > a
    while b - a > eps(b)
        if display & BISECT > 0
            println("bisect: a = ", a, ", b = ", b, ", b - a = ", b - a)
        end
        d = (a + b) / convert(T, 2)
        phi_d, gphi = phidphi(d)
        @assert isfinite(phi_d) && isfinite(gphi)

        push!(alphas, d)
        push!(values, phi_d)
        push!(slopes, gphi)

        id = length(alphas)
        if gphi >= zeroT
            return ia, id # replace b, return
        end
        if phi_d <= phi_lim
            a = d # replace a, but keep bisecting until dphi_b > 0
            ia = id
        else
            b = d
            ib = id
        end
    end
    return ia, ib
end
© 2021 GitHub, Inc.
Terms
Privacy
Security
Status
Docs
Contact GitHub
Pricing
API
Training
Blog
About
Loading complete

  esatta := proc( fun )
    local res, alpha;
    #res := Optimization[Minimize](
    res := Optimization[NLPSolve](
      fun(alpha),
      {alpha >= 0, alpha <= 3},
      #initialpoint={alpha=1},
      iterationlimit=50
    ):
    subs(res[2],alpha);
  end proc:

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

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

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

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

end module:

max_stars_repo_name: ebertolazzi/NLtoolbox

max_issues_repo_name: ebertolazzi/NLtoolbox

max_forks_repo_name: ebertolazzi/NLtoolbox


text:
fadf:=proc()
    fadsfafadfasd:=module()
    end:
end:

gjhdg:=module()

end:

ext: mpl


text:
#@ Not autoload

# 3-stem
SpherePi(5,2) := [[2,o(eta[2],eta[3],eta[4])]]:

# NB nuprime is usually only defined at p=2.
# Presumably we get a generator from the unstable J-map?
SpherePi(6,3) := [[12,nuprime]]:

SpherePi(7,4) := [[infinity,nu[4]],[12,E(nuprime)]]:

for i from 5 to 20 do
 SpherePi(i+3,i) := [[24,nu[i]]]:
od:

# 4-stem
SpherePi(6,2) := [[12,o(eta[2],nuprime)]]:
SpherePi(7,3) := [[2,o(nuprime,eta[6])]]:
SpherePi(8,4) := [[2,o(nu[4],eta[7])],[2,o(E(nuprime),eta[7])]]:
SpherePi(9,5) := [[2,o(nu[5],eta[8])]]:

for i from 6 to 20 do
 SpherePi(i+4,i) := []:
od:

# 5-stem

SpherePi(7, 2) := [[2,o(eta[2],nuprime,eta[6])]]:
SpherePi(8, 3) := [[2,o(nuprime,eta[6],eta[7])]]:
SpherePi(9, 4) := [[2,o(nu[4],eta[7],eta[8])],
                   [2,o(E(nuprime),eta[7],eta[8])]]:
SpherePi(10, 5) := [[2,o(nu[5],eta[8],eta[9])]]:
SpherePi(11, 6) := [[infinity,w[6]]]:

for i from 7 to 20 do
 SpherePi(i+5,i) := []:
od:

# 6-stem

SpherePi(8, 2) := [[2,o(eta[2],nuprime,eta[6],eta[7])]]:
SpherePi(9, 3) := []:
SpherePi(10, 4) := [[24,o(nu[4],nu[7])],3]:

for i from 5 to 20 do
 SpherePi(i+6,i) := [[2,o(nu[i],nu[i+3])]]:
od:

# 7-stem

SpherePi( 9, 2) := [3]:
SpherePi(10, 3) := [15]:
SpherePi(11, 4) := [15]:
SpherePi(12, 5) := [30]:
SpherePi(13, 6) := [60]:
SpherePi(14, 7) := [120]:
SpherePi(15, 8) := [[infinity,sigma[8]],120]:

for i from 9 to 20 do
 SpherePi(i+7,i) := [[240,sigma[i]]]:
od:

# 8-stem

SpherePi(10, 2) := [15]:
SpherePi(11, 3) := [[2,epsilon[3]]]:
SpherePi(12, 4) := [[2,epsilon[4]]]:
SpherePi(13, 5) := [[2,epsilon[5]]]:
SpherePi(14, 6) := [[2,epsilon[6]],[24,nubar[6]]]:
SpherePi(15, 7) := [[2,epsilon[7]],[2,nubar[7]],[2,o(sigmaprime,eta[14])]]:
SpherePi(16, 8) := [[2,epsilon[8]],
                    [2,nubar[8]],
                    [2,o(E(sigmaprime),eta[15])],
                    [2,o(sigma[8],eta[15])]]:
SpherePi(17, 9) := [[2,epsilon[9]],
                    [2,nubar[9]],
                    [2,o(sigma[9],eta[16])]]:


text:
# Base Parameter Regresor Inverse Dynamics for Robot based on MDH frames
# Einleitung
# Berechnung der inversen Dynamik in Regressorform
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# chain -> Berechnung für eine serielle Struktur (nicht: Baumstruktur)
# floatb -> Floating Base (und auch fixed Base) Modell der Basis.
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# dynamics -> Berechnung der Dynamik
# regressor -> Regressorform (parameterlinear)
#
# Siehe auch: robot_chain_floatb_rotmat_dynamics_regressor.mw.
# Im Gegensatz zu der allgemeinen Dynamik werden hier die fixed-Base-Ausdrücke ohne eine "twist"-Basis-Berechnung übernommen.
# Ursache ist, dass der Parametervektor und damit der Regressor sich zwischen Fixed-Base und Floating-Base unterscheidet und daher ein Ableiten des Fixed-Base falles aus der Floating-Base-Formulierung mit einfacher Basis-Darstellung ("twist") nicht möglich ist.
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
# Einstellungen für Code-Export: Optimierungsgrad (2=höchster) und Aktivierung jedes Terms.
codegen_act := true:
codegen_opt := 2:
codeexport_grav := true:
codeexport_corvec := true:
codeexport_cormat := true:
codeexport_inertia := true:
codeexport_inertiaD := true:
codeexport_invdyn := true:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_LagrangeN":
read "../helper/proc_index_symmat2vector":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
# Term-Vereinfachungen einstellen
if not assigned(simplify_options) or simplify_options(10)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(10): # zehnter Eintrag ist für Dynamik-Regressor
end if:
# Mit diesem Arbeitsblatt werden die Regressor-Ausdrücke für Fixed-Base und Floating-Base Modelle generiert. Erkenne welche Basis-Modellierung aktiv ist
if base_method_name="twist" then # Basis-Methode "twist" wird (hier) nur für fixed Base benutzt
  expstring:="fixb":
elif base_method_name="eulxyz" then
  expstring:="floatb_eulxyz":
else
  printf("Nicht behandelte Basis-Methode: %s\n", base_method_name):
fi:
# Es ist in diesem Arbeitsblatt möglich, zwei verschiedene Regressoren zu generieren und zu exportieren: Basierend auf Minimalparametern und auf vollem Parametersatz (PV2).
# Der Term "regressor" oder "regressor_minpar" ist jeweils in den Dateinamen enthalten.
# Der folgende Befehl muss immer auf "regressor_minpar" gesetzt sein, da diese Zeile durch das Skript robot_codegen_maple_preparation.sh ausgewertet und modifiziert wird.
regressor_modus := "regressor_minpar":
if regressor_modus = "regressor_minpar" then
  ekinfile := sprintf("../codeexport/%s/tmp/energy_kinetic_%s_regressor_minpar_maple.m", robot_name, expstring):
  if FileTools[Exists](ekinfile) then
    read ekinfile:
  else
    printf("%s. Energie-Regressor in Minimalparameterform wurde nicht berechnet. Abbruch der Dynamik-Minimalparameterregressor-Berechnung.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
    quit: # Funktioniert in GUI nicht richtig...
    robot_name := "": # ...Daher auch Löschung des Roboternamens.
  end if:
  read ekinfile:
  read sprintf("../codeexport/%s/tmp/energy_potential_%s_regressor_minpar_maple.m", robot_name, expstring):
  t_ges := t_ges_minpar:
  u_ges := u_ges_minpar:
  printf("%s. Generiere Minimalparameterregressor der Dynamik für %s\n", FormatTime("%Y-%m-%d %H:%M:%S"), robot_name):
end if:
if regressor_modus = "regressor" then
  read sprintf("../codeexport/%s/tmp/energy_kinetic_%s_regressor_maple.m", robot_name, expstring):
  read sprintf("../codeexport/%s/tmp/energy_potential_%s_regressor_maple.m", robot_name, expstring):
  t_ges := t_ges:
  u_ges := u_ges:
  printf("%s. Generiere Regressor der Dynamik für %s (nicht Minimalparameter)\n", FormatTime("%Y-%m-%d %H:%M:%S"), robot_name):
end if:
# Schalter zur Auswahl der unterschiedlichen Terme, die exportiert werden sollen. Für parallele Berechnung interessant.
DynString := "Term:":
if codeexport_grav then
  DynString := sprintf("%s g",DynString):
end if:
if codeexport_corvec then
  DynString := sprintf("%s c",DynString):
end if:
if codeexport_cormat then
  DynString := sprintf("%s C",DynString):
end if:
if codeexport_inertia then
  DynString := sprintf("%s M",DynString):
end if:
if codeexport_inertiaD then
  DynString := sprintf("%s MD",DynString):
end if:
if codeexport_invdyn then
  DynString := sprintf("%s tau",DynString):
end if:
printf("%s. Generiere Dynamik-Regressor (%s) für %s und %s\n", FormatTime("%Y-%m-%d %H:%M:%S"), DynString, robot_name, base_method_name):
# Platzhalter-Vektor der Dynamik-Parameter aufstellen
nDP := ColumnDimension(u_ges):
PV := Matrix(nDP,1):
for i from 1 to nDP do
  if regressor_modus = "regressor_minpar" then
    PV(i,1) := parse(sprintf("MDP%d%", i)):
  elif regressor_modus = "regressor" then
    PV(i,1) := parse(sprintf("DP%d%", i)):
  end if:
end do:
PV:
# Name für Export der Terme als Funktion des Parametervektors
if regressor_modus = "regressor_minpar" then
  regshortname := "mdp":
elif regressor_modus = "regressor" then
  regshortname := "dp":
end if:

# Lagrange Formalismus (mit Funktion)
OutputLagrange := LagrangeN(t_ges, u_ges):
dTdqDdt_s := OutputLagrange[1]:
dTdq_s := OutputLagrange[2]:
dUdq_s := OutputLagrange[3]:
save dUdq_s, sprintf("../codeexport/%s/tmp/floatb_lagrange_dUdq_s_%s_maple.m", robot_name, regressor_modus):
save dTdq_s, sprintf("../codeexport/%s/tmp/floatb_lagrange_dTdq_s_%s_maple.m", robot_name, regressor_modus):
save dTdqDdt_s, sprintf("../codeexport/%s/tmp/floatb_lagrange_dTdqDdt_s_%s_maple.m", robot_name, regressor_modus):
# Extraktion einzelner Terme
# Gravitational Load
# Generate
taug_regressor_s := dUdq_s:
# Terme vereinfachen
if codeexport_grav then
  if use_simplify=1 then
    tmp_t1:=time():
    tmp_l1 := length(taug_regressor_s):
    taug_regressor_s := simplify2(taug_regressor_s):
    tmp_t2:=time():
    tmp_l2 := length(taug_regressor_s):
    printf("%s: Gravitationsmoment-Regressor vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end if:
end if:
save taug_regressor_s, sprintf("../codeexport/%s/tmp/gravload_%s_maple.m", robot_name, regressor_modus):
# Matlab Export
# Belastung der Basis (nur Floating Base)
if codegen_act and codeexport_grav and not(base_method_name="twist") then
  MatlabExport(taug_regressor_s(1..6,..), sprintf("../codeexport/%s/tmp/gravload_base_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Belastung der Gelenke:
# Fixed-Base: Gravitationsvektor im Basis-KS
# Floating Base: Gravitationsvektor im Welt-KS, Basis-Orientierung berücksichtigt
taug := taug_regressor_s . PV:
if codegen_act and codeexport_grav then
  MatlabExport(taug_regressor_s(7..NQ,..), sprintf("../codeexport/%s/tmp/gravload_joint_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(taug(7..NQ,..), sprintf("../codeexport/%s/tmp/gravload_joint_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Kompletter Vektor
if codegen_act and codeexport_grav and not(base_method_name="twist") then
  MatlabExport(taug_regressor_s(1..NQ,..), sprintf("../codeexport/%s/tmp/gravload_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
  MatlabExport(taug(1..NQ,..), sprintf("../codeexport/%s/tmp/gravload_floatb_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:

# Mass Matrix
# Berechnung vollständige Massenmatrix
# Generiere die Massenmatrix erst für Floating-Base-Komplettsystem, dann entnehme die Teilmatrizen für Gelenk-Gelenk, Gelenk-Basis zum separaten Export
Paramvec_size := ColumnDimension(t_ges):
# Berechnung vollständige Massenmatrix
# Initialisiere. Speichere nur den unteren linken Teil der Massenmatrix
# Siehe: https://de.wikipedia.org/wiki/Symmetrische_Matrix
tauMM_regressor_s := dTdqDdt_s(1..NQ,..):
MM_regressor_s := Matrix(NQ*(NQ+1)/2, Paramvec_size):
if codeexport_inertia or codeexport_inertiaD or codeexport_cormat then
i_rr := 0:
for i to NQ do # Zeilenindex der Massenmatrix
  for j to NQ do  # Spaltenindex der Massenmatrix
    if j > i then
      next: # rechte obere Seite der symmetrischen Matrix. Keine neue Information. Nicht berechnen oder speichern.
    end if:
    i_rr := i_rr + 1: # Gehe zeilenweise durch den unteren linken Teil der Massenmatrix (inkl. Diagonale)
    for k to Paramvec_size do # Spaltenindex der Regressormatrix
      MM_regressor_s[i_rr, k] := diff(tauMM_regressor_s[i, k], qDD_s[j, 1]):
    end do:
  end do:
end do:
# Terme vereinfachen
if use_simplify>=1 then
  tmp_t1:=time():
  tmp_l1 := length(MM_regressor_s):
  printf("%s. Beginne Vereinfachung: Massenmatrix-Regressor (%s). Länge: %d.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), regressor_modus, tmp_l1):
  MM_regressor_s := simplify2(MM_regressor_s):
  tmp_t2:=time():
  tmp_l2 := length(MM_regressor_s):
  printf("%s. Massenmatrix-Regressor (%s) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), regressor_modus, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:
end if:
# Export Gesamt-Massenmatrix
if codegen_act and codeexport_inertia and not(base_method_name="twist") then
  MatlabExport(MM_regressor_s, sprintf("../codeexport/%s/tmp/inertia_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Gelenk-Massenmatrix extrahieren
MMjj_regressor_s := Matrix(NQJ*(NQJ+1)/2, Paramvec_size):
itmp:=0:
for i to NQ do
  for j to NQ do  # Spaltenindex der Massenmatrix
    if j > i then
      next: # rechte obere Seite der symmetrischen Matrix. Keine neue Information. Nicht berechnen oder speichern.
    end if:
    if i > 6 and j > 6 then # unterer rechter Teil (Gelenkmoment-Gelenkbeschleunigung-Terme)
      i_MM := index_symmat2vec(NQ,i,j): # Passender Index für zeilenweise ausgewählten symmetrischen Teil (siehe Gesamt-Massenmatrix)
      itmp := itmp + 1:
      MMjj_regressor_s[itmp,..] := MM_regressor_s[i_MM,..]:
    end if:
  end do:
end do:
# Multiplikation der Massenmatrix mit dem Regressorvektor und Export als obere Dreiecksmatrix
MMjj_tmp := MMjj_regressor_s . PV: # Vektor der oberen rechten Dreiecksmatrix
;
# Export Gelenkmassenmatrix
if codegen_act and codeexport_inertia and base_method_name="twist" then
  MatlabExport(MMjj_regressor_s, sprintf("../codeexport/%s/tmp/inertia_joint_joint_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(MMjj_tmp, sprintf("../codeexport/%s/tmp/inertia_joint_joint_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Gelenk-Basis-Massenmatrix
#  (untere linke Teilmatrix der Gesamt-Massenmatrix)
MMjb_regressor_s := Matrix(6*NQJ, Paramvec_size):
# Gehe Schleife über alle Massenmatrix-Elemente durch und entnehme die passenden Elemente für die Teilmatrix)
itmp:=0:
for i to NQ do
  for j to NQ do  # Spaltenindex der Massenmatrix
    if i > 6 and j < 7 then # unterer linker Teil
      i_MM := index_symmat2vec(NQ,i,j): # Passender Index für zeilenweise ausgewählten symmetrischen Teil (siehe Gesamt-Massenmatrix)
      itmp := itmp + 1:
      MMjb_regressor_s[itmp,..] := MM_regressor_s[i_MM,..]:
    end if:
  end do:
end do:
if codegen_act and codeexport_inertia and not(base_method_name="twist") then
  MatlabExport(MMjb_regressor_s, sprintf("../codeexport/%s/tmp/inertia_joint_base_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Basis-Massenmatrix
MMbb_regressor_s := Matrix(6*(6+1)/2, Paramvec_size):
# Gehe Schleife über alle Massenmatrix-Elemente durch und entnehme die passenden Elemente für die Teilmatrix)
itmp:=0:
for i to NQ do
  for j to NQ do  # Spaltenindex der Massenmatrix
    if j > i then
      next: # rechte obere Seite der symmetrischen Matrix. Keine neue Information. Nicht berechnen oder speichern.
    end if:
    if i < 7 and j < 7 then # unterer linker Teil
      i_MM := index_symmat2vec(NQ,i,j): # Passender Index für zeilenweise ausgewählten symmetrischen Teil (siehe Gesamt-Massenmatrix)
      itmp := itmp + 1:
      MMbb_regressor_s[itmp,..] := MM_regressor_s[i_MM,..]:
    end if:
  end do:
end do:
if codegen_act and codeexport_inertia and not(base_method_name="twist") then
  MatlabExport(MMbb_regressor_s, sprintf("../codeexport/%s/tmp/inertia_base_base_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Mass Matrix Time Derivative
# Gesamt Massenmatrix
MM_regressor_t := convert_s_t(MM_regressor_s):
MMD_regressor_t := diff~(MM_regressor_t, t):
MMD_regressor_s := convert_t_s(MMD_regressor_t):
if codegen_act and codeexport_inertiaD and not(base_method_name="twist") then
  MatlabExport(MMD_regressor_s, sprintf("../codeexport/%s/tmp/inertiaD_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Konvertiere Gelenk-Massenmatrix in zeitabhängige Variablen, um Zeitableitung zu berechnen
MMjj_regressor_t := convert_s_t(MMjj_regressor_s):
MMDjj_regressor_t := diff~(MMjj_regressor_t, t):
MMDjj_regressor_s := convert_t_s(MMDjj_regressor_t):
if codegen_act and codeexport_inertiaD and base_method_name="twist" then
  MatlabExport(MMDjj_regressor_s, sprintf("../codeexport/%s/tmp/inertiaD_joint_joint_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
end if:
# Gelenk-Basis-Massenmatrix
MMjb_regressor_t := convert_s_t(MMjb_regressor_s):
MMDjb_regressor_t := diff~(MMjb_regressor_t, t):
MMDjb_regressor_s := convert_t_s(MMDjb_regressor_t):
if codegen_act and codeexport_inertiaD and not(base_method_name="twist") then
  MatlabExport(MMDjb_regressor_s, sprintf("../codeexport/%s/tmp/inertiaD_joint_base_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Coriolis Vector
# Generieren (gleiches Vorgehen für fixed und floating base)
tauC_regressor_s := dTdqDdt_s-dTdq_s:
for i to NQ do
  tauC_regressor_s := subs({qDD_s(i, 1) = 0}, tauC_regressor_s):
end do:
save tauC_regressor_s, sprintf("../codeexport/%s/tmp/coriolisvec_joint_%s_%s_maple.m", robot_name, expstring, regressor_modus):
tauc := tauC_regressor_s . PV:
# Matlab Export
# Belastung der Gelenke
if codegen_act and codeexport_corvec then
  MatlabExport(tauC_regressor_s(7..NQ,..), sprintf("../codeexport/%s/tmp/coriolisvec_joint_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tauc(7..NQ,..), sprintf("../codeexport/%s/tmp/coriolisvec_joint_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Gesamter Vektor für Floating Base
if codegen_act and codeexport_corvec and not(base_method_name="twist") then
  MatlabExport(tauC_regressor_s(1..NQ,..), sprintf("../codeexport/%s/tmp/coriolisvec_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tauc(1..NQ,..), sprintf("../codeexport/%s/tmp/coriolisvec_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Nur Basis-Terme
if codegen_act and codeexport_corvec and not(base_method_name="twist") then
  MatlabExport(tauC_regressor_s(1..6,..), sprintf("../codeexport/%s/tmp/coriolisvec_base_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tauc(1..6,..), sprintf("../codeexport/%s/tmp/coriolisvec_base_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Coriolis Matrix
# Floating-Base Coriolismatrix (Gesamt)
# Calculation with Christoffel Symbol approach
# [KhalilDombre2002], equ. (9.7) (p. 195)
# Funktion zur Anwendung des Algorithmus aus [KhalilDombre2002]:
cijk := proc (i::integer, j::integer, k::integer, A, qs)
  local c:
  c := (1/2)*(diff(A[i, j], qs(k, 1)))+(1/2)*(diff(A[i, k], qs(j, 1)))-(1/2)*(diff(A[j, k], qs(i, 1))):
  return c:
end proc:
# Initialisierung. Speichere die vollständige Coriolismatrix (nicht symmetrisch/schiefsymmetrisch). Unterschied zum Vorgehen bei der Massenmatrix
C_regressor_s := Matrix(NQ*NQ, Paramvec_size):
# Berechnung
if codeexport_cormat then:
i_rr := 0: # Vektor-Index für den Regressor der Coriolismatrix (Ausgabe)
for i to NQ do # Zeilenindex der Coriolismatrix
  for j to NQ do  # Spaltenindex der Coriolismatrix
    i_rr := i_rr + 1: # Gehe zeilenweise durch den unteren linken Teil der Coriolismatrix (inkl. Diagonale)
    for k from 1 to Paramvec_size do # Spaltenindex der Regressormatrix
      # Massenmatrix für Parameter k generieren (für Funktion mit Christoffel-Symbol-Ansatz benötigt)
      MM_k := Matrix(NQ,NQ):
      for ii from 1 to NQ do
        for jj from 1 to NQ do
          MM_k(ii,jj) := MM_regressor_s(index_symmat2vec(NQ, ii, jj), k):
        end do:
      end do:
      for l from 1 to NQ do # Siehe [KhalilDombre2002]
        C_regressor_s[i_rr, k] := C_regressor_s[i_rr,k] + cijk(i,j,l,MM_k,q_s)*qD_s[l,1]:
      end do:
    end do:
  end do:
end do:
# Terme vereinfachen
if use_simplify>=1 then
  tmp_t1:=time():
  tmp_l1 := length(C_regressor_s):
  printf("%s. Beginne Vereinfachung: Coriolis-Matrix-Regressor (%s). Länge: %d.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), regressor_modus, tmp_l1):
  C_regressor_s := simplify2(C_regressor_s):
  tmp_t2:=time():
  tmp_l2 := length(C_regressor_s):
  printf("%s: Coriolis-Matrix-Regressor (%s) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), regressor_modus, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:
end if:
save C_regressor_s, sprintf("../codeexport/%s/tmp/coriolismat_%s_%s_maple.m", robot_name, expstring, regressor_modus):
# Matlab Export: Floating Base Gesamt
if codegen_act and codeexport_cormat and not(base_method_name="twist") then
  MatlabExport(C_regressor_s, sprintf("../codeexport/%s/tmp/coriolismat_floatb_%s_%s_matlab.m", robot_name, base_method_name, regressor_modus), codegen_opt):
end if:
# Gelenk-Coriolismatrix
# Extrahiere Teilmatrix, die für die Gelenkterme (ohne Basiseinfluss) zuständig sind
Cjj_regressor_s := Matrix(NQJ*NQJ, Paramvec_size):
itmp := 0:
for i to NQ do # Zeilenindex der Coriolismatrix
  for j to NQ do  # Spaltenindex der Coriolismatrix
    if i > 6 and j > 6 then # unterer rechter Teil
      i_C := (i-1)*NQ+j: # Passender Index für zeilenweise ausgewähltes Element
      itmp := itmp + 1:
      Cjj_regressor_s[itmp,..] := C_regressor_s[i_C,..]:
    end if:
  end do:
end do:
if codegen_act and codeexport_cormat and base_method_name="twist" then
  MatlabExport(Cjj_regressor_s, sprintf("../codeexport/%s/tmp/coriolismat_joint_fixb_%s_matlab.m", robot_name, regressor_modus), codegen_opt):
end if:
# Inverse Dynamics
# Generate
tau_regressor_s := dTdqDdt_s-dTdq_s+dUdq_s:
save tau_regressor_s, MMjj_regressor_s, tauMM_regressor_s, tauC_regressor_s, taug_regressor_s, sprintf("../codeexport/%s/tmp/invdyn_%s_%s_maple.m", robot_name, expstring, regressor_modus):
tau := tau_regressor_s . PV:
# Gesamter Vektor (floating base)
if codegen_act and codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(tau_regressor_s(1..NQ,..), sprintf("../codeexport/%s/tmp/invdyn_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tau, sprintf("../codeexport/%s/tmp/invdyn_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Belastung der Basis (floating base)
if codegen_act and codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(tau_regressor_s(1..6,..), sprintf("../codeexport/%s/tmp/invdyn_base_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tau(1..6,..), sprintf("../codeexport/%s/tmp/invdyn_base_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Belastung der Gelenke
if codegen_act and codeexport_invdyn then
  MatlabExport(tau_regressor_s(7..NQ,..), sprintf("../codeexport/%s/tmp/invdyn_joint_%s_%s_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tau(7..NQ,..), sprintf("../codeexport/%s/tmp/invdyn_joint_%s_%s_matlab.m", robot_name, expstring, regshortname), codegen_opt):
end if:
# Regressor-Matrix als Vektor hinschreiben. Dabei Null-Einträge herausfinden und nicht weiter berücksichtigen.
# Dadurch wird die obere rechte Dreiecksform der Regressor-Matrix ausgenutzt.
# Platzhalter-Vektor für alle Einträge der Regressor-Matrix
Reg_Vector := Matrix(ColumnDimension(tau_regressor_s)*RowDimension(tau_regressor_s), 1):
# Gelenkmoment in Abhängigkeit des Regressor-Vektors (RV)
tauJ_RV := Matrix(NQJ, 1):
kk := 0: # Zähler für laufende Nummer des Regressor-Vektors
for i from 1 to ColumnDimension(tau_regressor_s) do
  for j from 1 to RowDimension(tau_regressor_s)-6 do
    # Gehe Spalten- und dann Zeilenweise durch die Regressormatrix (nur Gelenk-Teil)
    if tau_regressor_s(j+6,i) <> 0 then # Eintrag ungleich Null
    	 kk := kk + 1:
    	 # Regressor-Vektor-Eintrag als Platzhalter direkt so, wie er in Matlab-Funktion gebraucht wird ("RV" ist Eingabe der Matlab-Funktion).
    	 rv_kk := parse(sprintf("RV(%d)", kk)):
    	 # Eintrag für Regressor-Vektor so schreiben, dass er in einer Matlab-Funktion aus der Regressor-Matrix gebildet werden kann
    	 Reg_Vector(kk,1) := parse(sprintf("RM(%d,%d)", j, i)):
    	 # Gelenkmoment aus Multiplikation mit Parametervektor. Dadurch Berechnung unter Ausnutzung der Nullen.
    	 tauJ_RV(j) := tauJ_RV(j) + rv_kk*PV(i):
    end if:
  end do:
end do:
# Kürzen des Regressor-Vektors auf die tatsächlich benutzten Einträge
Reg_Vector := Reg_Vector(1..kk,..):
save Reg_Vector, sprintf("../codeexport/%s/tmp/invdyn_joint_%s_%s_occupancy_vector_maple", robot_name, expstring, regressor_modus):
if codegen_act and codeexport_invdyn then
  interface(warnlevel=0): # Unterdrücke die folgende Warnung (Ist prinzipbedingt, da RM/RV als Funktionsnamen interpretiert werden. Der Code funktioniert trotzdem.
  MatlabExport(Reg_Vector, sprintf("../codeexport/%s/tmp/invdyn_joint_%s_%s_occupancy_vector_matlab.m", robot_name, expstring, regressor_modus), codegen_opt):
  MatlabExport(tauJ_RV, sprintf("../codeexport/%s/tmp/invdyn_joint_%s_%s_mult_matlab.m", robot_name, expstring, regshortname), codegen_opt):
  interface(warnlevel=3):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
lines.linewidth: 1.0          # line width in points

text.color: 232F3E

axes.facecolor:     FAFAFA    # axes background color
axes.edgecolor:     CACFD0    # axes edge color
axes.linewidth:     1.5       # edge line width
axes.grid:          True      # display grid or not
axes.titlesize:     14        # font size of the axes title
axes.labelsize:     10        # font size of the x and y labels
axes.labelcolor:    FF9900    # color of the x and y labels
axes.titleweight:   bold      # font weight of title
axes.titlecolor:    FF9900    # color of the axes title
#axes.labelcolor:    0073BB
axes.prop_cycle: cycler('color', ['232F3E', '527FFF', '0073BB', '4d27aa', 'b0084d', 'bf0816', 'd45b07', 'FF9900', 'e38b00', '067f68', '16bf9f'])

xtick.color:         0073BB   # color of the ticks
xtick.labelsize:     9        # font size of the tick labels
ytick.color:         0073BB   # color of the ticks
ytick.labelsize:     9        # font size of the tick labels

grid.color:     232F3E        # grid color
grid.alpha:     0.1           # transparency, between 0.0 and 1.0

legend.framealpha:    0.5     # legend patch transparency
legend.facecolor:     fee385  # inherit from axes.facecolor; or color spec
legend.edgecolor:     e38b00  # background patch boundary color
legend.fancybox:      True    # if True, use a rounded box for the

figure.titlesize:   large     # size of the figure title (``Figure.suptitle()``)
figure.titleweight: bold      # weight of the figure title
figure.facecolor:   F2F3F3    # figure face color

max_stars_repo_name: bpatoul/amazon-lookout-for-equipment

max_issues_repo_name: bpatoul/amazon-lookout-for-equipment

max_forks_repo_name: bpatoul/amazon-lookout-for-equipment


text:
input := FileTools:-Text:-ReadFile("AoC-2021-6-input.txt" ):
tally := table(sparse=0,Statistics:-Tally(
    StringTools:-Split(input,",")));
lanternfish := Array(0..8, [tally["0"],tally["1"],tally["2"],tally["3"],
    tally["4"],tally["5"],0,0,0]);

to 256 do
    day0 := lanternfish[0];
    lanternfish[0..7] := lanternfish[1..8];
    lanternfish[8] := day0;
    lanternfish[6] := lanternfish[6] + day0;
end do:

answer2 := add(lanternfish);


text:
importpath "."
importpath ".."
importpath "../ff"
var $x i32
importpath "/a/bb/ccc/"
func $foo (var %y i32) void {
  dassign $x (dread i32 %y)
}
 # 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:
#####################################################################
# this file contains some useful procedures for the computation of atan in
# DDE

restart:
Digits := 200:

with (numapprox):with(orthopoly):
read "double-extended.mpl";
read "common-procedures.mpl";

e := 2^(-6.3):
marge := 2^(-30);
P19 := convert(series(arctan(x),x=0,20),polynom):
P18 := (P19-x)/x;
Q18 := polyExact2Ext(P18, 9);
Q9 := -1/3 +1/5*x - 1/7*x^2 + 1/9*x^3 - 1/11*x^4 + 1/13*x^5 - 1/15*x^6 + 1/17*x^7- 1/19*x^8;
Qprime := poly_exact( -x /3 + x^2/5 - x^3/7 + x^4 /9);
#log2 (infnorm( (arctan(x)-x*(1+Q18))/x, x=0..e));

#function taken from the "old" coef_atan to compute interval parameters.
maxx := bi -> simplify( solve( (x-bi) / (1+x*bi) = e ,x) ):
minx := bi -> simplify(solve( (x-bi) /(1+x*bi)=-e ,x)):
nextbi := proc (x) evalf( max(solve( minx(bi) = x ,bi) )*(1-marge)); end:
allbi := proc (n)
local xi,nbi,x,i,j;
global b,a, nbofai, nbofbi, valuefordicho;
    x := e;
    nbi := 0;
    i := 0;
    while(i<n and nbi < 1/e and nbi >= 0) do
        nbi := nearest ( nextbi(x) );
        b[i] := evalf( nbi );
        a[i] := x;
        x   := evalf(maxx(nbi));
        i:=i+1;
    od;
    j:=0;
    while ( 2^j < i ) do j:=j+1 od:
    nbofai := i;
    nbofbi := i;
    b[i-1] := nearest(1/e+4):
    valuefordicho := 2^j;
    return i,b[0],b[i-1];
end:
allbi(100);
#-------------

#-----------------------
#quick error calc :
#-----Reduction:-------
XredEpsilon := 2^(-64) + 2^(-64) + 2^(-64) + 2^(-64):
Xred2Epsilon := 2*XredEpsilon + 2^(-64);
log[2.](Xred2Epsilon);
errlist:=errlist_quickphase_horner(degree(Qprime),0,0,Xred2Epsilon, 2^(-64));
errorr := compute_horner_rounding_error(Qprime,x,e,errlist,true);
qEpsilon := errorr[1]:
log[2.](qEpsilon);
deltaApprox := e^11/11:
polyDelta := e^3*qEpsilon + 2^(-64)*e+ 3*2^(-64)*e^3:
EpsilonFinal := polyDelta/arctan(e)+2*2^(-64)+deltaApprox/e:
log[2.](EpsilonFinal);
#------no reduction-------
Errlist:=errlist_quickphase_horner(degree(Qprime),0,0,0, 2**(-64)):
qEpsilon:= compute_horner_rounding_error(Qprime,x,e,errlist,true)[1]:
deltaApprox := e^11/5:
EpsilonFinalNoRed := infnorm( (qEpsilon.x^3+2^(-64).x+2^(-64).x^3)/arctan(x),x=2^(-27)..2^(-6) ) + deltaApprox + 2^(-64):
log[2.](EpsilonFinalNoRed);
#-----------------------

# Output :

doth:=proc()
 local filename, fd, i, hi, lo;
   filename := "TEMPATAN/atan_ext.h":
    fd := fopen (filename,WRITE,TEXT):
    fprintf(fd,"/* file generated by atan_ext.mpl*/\n\n"):
    fprintf(fd,"#include \"double_ext.h\"\n"):

    fprintf(fd, "static const double  HALFPI = %1.50e;\n", nearest(Pi/2)):

    fprintf(fd, "#define MIN_REDUCTION_NEEDED %s\n",  printDoubleAsHexInt(e)):
    fprintf(fd,"#define A 0\n"):
    fprintf(fd,"#define B 1\n"):
    fprintf(fd,"#define ATAN_BHI 0\n"):
    fprintf(fd,"#define ATAN_BLO 1\n"):

    fprintf(fd,"#define epsilon %1.50e\n", EpsilonFinal):
    fprintf(fd,"#define epsilon_no_red %1.50e\n",EpsilonFinalNoRed):
    fprintf(fd,"#define TWO_M_64 %1.50e\n",2^(-64)):
    fprintf(fd,"#define TWO_10 %1.50e\n",2^10):

    if 1+1=3 then
        fprintf(fd, "__declspec(align(16)) static const unsigned long long int a_table[%d] = {\n", nbofai );
        for i from 0 to nbofai - 1 do
            fprintf(fd, " /*a[%d] ~= %1.20e   */ ",i,a[i]):
            fprintf(fd, "  %s,\n", printDoubleAsULL(a[i]));
        od:
        fprintf(fd,"\n};\n"):

        fprintf(fd, "static const double b_table[%d] = {\n", nbofai );
        for i from 0 to nbofai - 1 do
            fprintf(fd, "   /*b[%d] = */   %1.50e,\n", i, b[i]):
        od:
        fprintf(fd,"\n};\n"):

    else
        fprintf(fd, "static const struct{long long int a; double b;} ab_table[%d] = {\n", nbofai);
        for i from 0 to nbofai - 1 do
            fprintf(fd, " { /*a[%d] ~= %1.20e   */ ",i,a[i]):
            fprintf(fd, "  %s,\n", printDoubleAsULL(a[i]));
            fprintf(fd, "   /*b[%d] = */   %1.50e},\n", i, b[i]):
        od:
        fprintf(fd,"\n};\n"):

    fi:

    fprintf(fd, "#define atanb_table ((const XC_FLOAT_TYPE *)_atanb_table)\n");
    fprintf(fd, "__declspec(align(16)) static const unsigned short _atanb_table[] = {\n", nbofai );
    for i from 0 to nbofai - 1 do
        fprintf(fd, " /*atan_b[%d] ~= %1.10e*/\n ",i, arctan(b[i])):
        hi:=nearestExt(arctan(b[i]));
        lo:=arctan(b[i])-nearestExt(arctan(b[i]));
        fprintf(fd, " %s, \n", printDoubleExtAsShort(hi));
        fprintf(fd, "   %s, \n", printDoubleExtAsShort(lo));
    od:
    fprintf(fd,"};\n"):

    printPolyExt(fd,Q9,4,"coef_poly"):
    fclose(fd):
end proc;

doth();

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


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

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

(* type: mgga_exc *)

$define gga_x_b86_mgc_params
$include "gga_x_b86.mpl"

cab := 0.63:
css := 0.96:

(* Equation 50, same-spin correlation *)
b88_css := (rs, z, xs, ts) ->
  - 0.01 * (1 + z)^(8/3) * 2^(-8/3) * n_total(rs)^(5/3) * (2*ts - xs^2/4)
  * b88_zss(css, b86_f, rs, z, xs)^4 * (
    1 - 2*log(1 + b88_zss(css, b86_f, rs, z, xs)/2)
      / b88_zss(css, b86_f, rs, z, xs)
  ):

(* Same-spin correlation overall *)
b88_par := (rs, z, xs0, xs1, ts0, ts1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, b88_css(rs, z_thr( z), xs0, ts0))
  + my_piecewise3(screen_dens(rs, -z), 0, b88_css(rs, z_thr(-z), xs1, ts1)):

(* Equation 49, opposite-spin correlation *)
b88_cab := (rs, z, xs0, xs1) ->
  - 0.8 * (1 - z^2)/4 * n_total(rs)
  * b88_zab(cab, b86_f, rs, z, xs0, xs1) * (
      b88_zab(cab, b86_f, rs, z, xs0, xs1) - log(1 + b88_zab(cab, b86_f, rs, z, xs0, xs1))
  ):

(* Whole functional *)
b88_c_f := (rs, z, xs0, xs1, ts0, ts1) ->
  + b88_cab(rs,  z, xs0, xs1)
  + b88_par(rs,  z, xs0, xs1, ts0, ts1):

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


text:
srclang 3
type $SS <struct {@f1 i32, @f2 i8, @f3 <[] u16>}>
func $foo ( var %x <*<$SS>>, var %y <*<$SS>>) void {
  var %i i32
  var %j i32
  iassign <* u16> ( array 1 ptr <*[]u16> (
                iaddrof ptr <* <$SS>> 3 (dread ptr %x), dread i32 %i),
                iread u16 <* u16> ( array 1 ptr <*[]u16> (
                    iaddrof ptr <* <$SS>> 3 (dread ptr %y), dread i32 %j)))
  return()  }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# Center of Mass Jacobian for the Robot based on MDH frames
# Introduction
# Berechnung der Schwerpunkts-Jacobi-Matrix (Veränderung des Schwerpunktes in Abhängigkeit der Gelenkwinkel).
# Diese Matrix wird bei Gleichgewichtsregelungen (z.B. zweibeiniges Laufen) gebraucht.
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb -> floating base wird durch base twist (Geschwindigkeit der Basis) oder vollständige Orientierung (Euler-Winkel) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# kinematics_comjacobi -> Jacobi-Matrix der Körperschwerpunkte
# worldframe -> Berechnung der Geschwindigkeit im Welt-KS (KSW)
# par1 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ)
# Authors
# Moritz Schappler, schappler@irt.uni-hannover.de, 2015-12
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
#
# Sources
# [SugiharaNakIno2002] Sugihara, T. and Nakamura, Y. and Inoue, H.: Real-time humanoid motion generation through ZMP manipulation based on inverted pendulum control (2002)
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_opt := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
read "../robot_codegen_definitions/robot_env":
printf("Generiere Schwerpunktskinematik für %s\n", robot_name):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf := Trf:
Trf_c := Trf_c:
# Schwerpunkt des Roboters laden
read sprintf("../codeexport/%s/tmp/kinematics_com_total_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
# Calculate CoG Jacobian
# CoG Jacobian nach [SugiharaNakIno2002].
# TODO: baseframe in worldframe ändern
J_COG_s := Matrix(3, NQJ):
c_s := convert_t_s(c):
for i from 1 to 3 do
  for j from 1 to NQJ do
    J_COG_s(i,j) := diff(c_s(i,1), qJ_s(j,1)):
  end:
end:
# TODO: Wenn der Ausdruck J_COG_s nicht gespeichert und wieder geladen wird, hängt die Berechnung bei der Optimierung des Ausdrucks mit "tryhard"
save J_COG_s, sprintf("../codeexport/%s/tmp/kinematics_Jcom_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
read sprintf("../codeexport/%s/tmp/kinematics_Jcom_worldframe_floatb_%s_par1_maple.m", robot_name, base_method_name):
if codegen_act then
  MatlabExport(J_COG_s, sprintf("../codeexport/%s/tmp/com_jacobi_baseframe_par1_matlab.m", robot_name), 2):
end if:
# CoG-Jacobian Time derivative
J_COG := convert_s_t(J_COG_s):
JD_COG := diff~(J_COG, t):
JD_COG_s := convert_t_s(JD_COG):
if codegen_act then
  MatlabExport(convert_t_s(JD_COG_s), sprintf("../codeexport/%s/tmp/com_jacobiD_baseframe_par1_matlab.m", robot_name), codegen_opt):
end if:


text:
with(StringTools):
DiagLatex := proc(spart, size:=0)
	local phis, thetas, phithetas, zeds, Lambda, string, cercle, triangle, triangle_cercle, boite, part, endline, k;
	phis:= spart[2];
	phis:= map( x-> [x,"c"], phis);
	thetas:= spart[3];
	thetas:= map(x-> [x,"b"], thetas);
	phithetas:= spart[1];
	phithetas:= map(x-> [x,"d"], phithetas);
	zeds:= spart[4];
	zeds:= map(x-> [x,"a"], zeds);
	Lambda:= sort([op(phithetas),op(phis), op(thetas), op(zeds)]);

	string:="";

	cercle:= "\\yC";
	triangle:= "\\yT";
	triangle_cercle:= "\\yTC";
	boite:= "\\,";
	for k from 1 to nops(Lambda) do
		part:= Lambda[k];
		print(Lambda);
		print(part);
		if k = 1 then endline:=""; else endline:="\\\\ \n"; end if:
		if part[2] = "c" then
			string:= cat(cat(Repeat(cat(boite, "& "), part[1]), triangle, endline), string);
		elif part[2] = "b" then
			string:= cat(cat(Repeat(cat(boite, "& "), part[1]), cercle, endline), string);
		elif part[2] = "a" then
			string:= cat(cat(Repeat(cat(boite, "& "), part[1]-1), boite, endline), string);
		elif part[2] = "d" then
			string:= cat(cat(Repeat(cat(boite, "& "), part[1]), triangle_cercle, endline), string);
		end if;
	end do:
	if size = 0 then
		string:= cat("{\\, \\superYsmall{", string, "}}");
	elif size = 1 then
		string:= cat("{\\, \\superY{", string, "}}");
	elif size = 2 then
		string:= cat("\\superRusse{{\\superY{", string, "}}}");
	end if;
	return string;
end proc:

superLatex := proc(expr,initial_term:=0, wipe:=0)
	local sparts, astring, spart, spart_string, fd, basetype, terms, thecoeff, latexsparts, latex_terms, k, ou;
	if wipe <> 0 then
		fd:= fopen("out.txt", WRITE);
		fprintf(fd, "");
		fclose(fd);
	else
		fd:=fopen("out.txt", APPEND);
		fprintf(fd,"new \n \\begin{gather*} \n");
		fclose(fd);
	end if;
	basetype:= "m";
	terms:= [op(indets(expr, 'superindexed'))];
	terms:= sort(terms, list_sort);
	#print(terms);
	sparts:= map(x-> [op(x)], terms);
	#print(sparts);
	thecoeff:= map(x-> coeff(expr, x), terms);
	#print(thecoeff);
	#latexcoeff:= map(x-> latex(x), thecoeff);

	#Activate following line for diagrams
	#latexsparts:= map(x-> DiagLatex(x,0), sparts);
	#Activate following line for sparts
	latexsparts:= map(x-> SpartLatex(x),sparts);
	latex_terms:= [];
	if initial_term <> 0 then
		fd:=fopen("out.txt", APPEND);
		fprintf(fd,cat(initial_term[1],"\\,",basetype,"_",DiagLatex(initial_term[2], 0),"=","\n+"));
		fclose(fd);
	end if;
	for k from 1 to nops(sparts) do
		fd:= fopen("out.txt", APPEND);
		fprintf(fd, "\t");
		fclose(fd);
		if whattype(thecoeff[k]) = `+` then
			fd:= fopen("out.txt", APPEND);
			fprintf(fd, "(");
			fclose(fd);
		end if;

		latex(thecoeff[k],"out.txt",'append');

		if whattype(thecoeff[k]) = `+` then
			fd:= fopen("out.txt", APPEND);
			fprintf(fd, ")");
			fclose(fd);
		end if;
		fd:= fopen("out.txt", APPEND);
		fprintf(fd, cat(basetype,"_",latexsparts[k],"\n+"));
		fclose(fd);
	end do:
		fd:=fopen("out.txt", APPEND);
		fprintf(fd,"\\end{gather*} \n");
		fclose(fd);
end proc:

SpartLatex:= proc(spart)
	local pts, phis, thetas, zeds, pt_string, phi_string, theta_string, zed_string, SpartString;
	pts:= spart[1];
	phis:= spart[2];
	thetas:= spart[3];
	zeds:= spart[4];

	pt_string:= convert(pts, string);
	phi_string:= convert(phis, string);
	theta_string:= convert(thetas, string);
	zed_string:= convert(zeds, string);

	pt_string:= Substitute(pt_string, "[", "(");
	pt_string:= Substitute(pt_string, "]", ";\\,");
	phi_string:= Substitute(phi_string, "[", "");
	phi_string:= Substitute(phi_string, "]", ";\\,");
	theta_string:= Substitute(theta_string, "[", "");
	theta_string:= Substitute(theta_string, "]", ";\\,");
	zed_string:= Substitute(zed_string, "[", "");
	zed_string:= Substitute(zed_string, "]", ")");

	SpartString:= cat(pt_string, phi_string, theta_string, zed_string);
	SpartString:= cat("{",SpartString, "}");
	return SpartString;
end proc:


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

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

ft97_beta := (rs, z, xs) -> params_a_beta0
  + params_a_beta1*sigma_spin(rs, z, xs)/(params_a_beta2 + sigma_spin(rs, z, xs)):

ft97_fx := (rs, z, xs) -> 1 + ft97_beta(rs, z, xs)*xs^2 /
  (X_FACTOR_C*sqrt(1 + 9*xs^2*ft97_beta(rs, z, xs)^2*arcsinh(xs^2)^2)):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange_nsp(ft97_fx, rs, zeta, xs0, xs1):


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

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

(* type: work_mgga_x *)

mlambda := 0.6866:
mbeta   := 79.873:

(* below Equation (6) *)
y  := p -> (2*mlambda - 1)^2 * p:

(* Equation (7) *)
f0 := x -> (1 + 10*70*y(X2S^2*x^2)/27 + mbeta*y(X2S^2*x^2)^2)^(1/10):

R  := (x, t) -> 1 + 595*(2*mlambda - 1)^2 * X2S^2*x^2/54 \
   - (t - 3*(mlambda^2 - mlambda + 1/2)*(t - K_FACTOR_C - x^2/72))/K_FACTOR_C:

fx_DME := (x, t) -> 1/f0(x)^2 + 7*R(x, t)/(9*f0(x)^4):

malpha := (x, t) -> (t - x^2/8)/K_FACTOR_C:
qtilde := (x, t) -> 9/20*(malpha(x, t) - 1) + 2*X2S^2*x^2/3:

fx_SC_aux := (x, t) -> (1 + 10*( \
       + (MU_GE + 50*X2S^2*x^2/729)*X2S^2*x^2 \
       + 146*qtilde(x, t)^2/2025 \
       - 73*qtilde(x,t)/405*(3*K_FACTOR_C/(5*t))*(1 - K_FACTOR_C/t))
       ):

(* It turns out that the quantity fx_SC_aux can become negative, leading to NaN.
   this appears, for example, for the following densities

   ./xc-get_data 540 1 729.292719352103 0 365940004.708665 0 0 -970050.160877395 0 64600.7149099035 0

   As a turn around, I demand that the function is always larger than 1e-16
*)
fx_SC  := (x, t) -> m_max(fx_SC_aux(x,t), 1e-16)^(0.1):

w := t-> (K_FACTOR_C^2/t^2 + 3*K_FACTOR_C^3/t^3)/(1 + K_FACTOR_C^3/t^3)^2:

f := (rs, x, t, u) -> w(t)*fx_DME(x, t) + (1 - w(t))*fx_SC(x, t):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


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

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

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


text:
#Play with n, deg and s to have larger examples
deg:=2:
for n from 4 to 10 do

    vars:=[seq(cat(x,i),i=1..n+1)]:
    f1:=expand(randpoly(vars,degree=deg, dense));
    f2:=expand(randpoly(vars,degree=deg, dense)):
    f:=resultant(f1, f2, vars[-1]);

    #To decompose
    F:=[f, seq(diff(f,v), v in vars[2..n])]:
    print(nops(F), nops(vars[1..n]), degree(f));
#    FGb:-fgb_matrixn_radical([l*diff(f, vars[1])-1, op(F)], 0, [l, op(vars[1..n])], 0, {"verb"=3});

   str:=cat("./singhypcritpoints_deg_", deg, "_nvars_", n, ".jl");
   fd:=fopen(str, WRITE):
   fprintf(fd, "using SignatureGB\n");
   fprintf(fd, "using Singular\n");

   fprintf(fd, "P, (%a, ", vars[1]);
   for j from 2 to n-1 do
   fprintf(fd, "%a, ", vars[j]);
   od:
   fprintf(fd, "%a)  = Singular.PolynomialRing(Fp(65521), [", vars[n]);
   for j from 1 to n-1 do
   fprintf(fd, "\"%a\", ", vars[j]);
   od:
   fprintf(fd, "\"%a\"]);\n", vars[n]);
   fprintf(fd, "gb = SignatureGB.f5(%a);", vars[1..n]);

   fprintf(fd, "F = %a;\n", F);
   fprintf(fd, "@time gb = SignatureGB.f5(%a, verbose=true);", F);
   fprintf(fd, "[leading_monomial(g) for g in gb];");
   fprintf(fd, "exit()\n");
   fclose(fd):

od:


text:
// Deletes all portals whose color is specified by an entity named ACV_deletePortals_PARAM and then kill the params
// Params: entities with name ACV_deletePortals_PARAM and a score of ACV_COLOR
process ACV_deletePortals {
  triggerAntiBridges()
  closeConnections()
  killPortals()
  lightUpHorizontalPortals()

  // Kill the Params
  /kill @e[name=ACV_deletePortals_PARAM]
}

// Triggers the AntiBridges at every specified Portal
// Params: entities with tag ACV_deletePortals_PARAM and a score of ACV_COLOR
process triggerAntiBridges {
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=0,score_ACV_COLOR=1] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=0,score_ACV_COLOR=1] ~ ~ ~ entitydata @e[type=area_effect_cloud,name=ACV_AntiBridge_Idle,dy=0] {CustomName:"ACV_AntiBridge"}
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=2,score_ACV_COLOR=3] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_PortalOut,score_ACV_COLOR_min=2,score_ACV_COLOR=3] ~ ~ ~ entitydata @e[type=area_effect_cloud,name=ACV_AntiBridge_Idle,dy=0] {CustomName:"ACV_AntiBridge"}
  start ACV_antiBridges
}

// Close the connection between the two portals
// Params: entities with tag ACV_deletePortals_PARAM and a score of ACV_COLOR
process closeConnections {
  ACV_loadChunks()
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=0,score_ACV_COLOR=1] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_PortalIn,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ clone 0 1 0 0 2 0 ~ ~ ~
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=0,score_ACV_COLOR=1] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_PortalIn,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ clone 1 1 0 1 2 0 ~ ~ ~
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=2,score_ACV_COLOR=3] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_PortalIn,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ clone 2 1 0 2 2 0 ~ ~ ~
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=2,score_ACV_COLOR=3] ~ ~ ~ execute @e[type=area_effect_cloud,name=ACV_PortalIn,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ clone 3 1 0 3 2 0 ~ ~ ~
}

// Delete the specified portals
// Params: entities with tag ACV_deletePortals_PARAM and a score of ACV_COLOR
process killPortals {
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ kill @e[tag=ACV_Portal,score_ACV_COLOR_min=0,score_ACV_COLOR=0]
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ kill @e[tag=ACV_Portal,score_ACV_COLOR_min=1,score_ACV_COLOR=1]
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ kill @e[tag=ACV_Portal,score_ACV_COLOR_min=2,score_ACV_COLOR=2]
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ kill @e[tag=ACV_Portal,score_ACV_COLOR_min=3,score_ACV_COLOR=3]
}

// Place barrier blocks at floor and ceiling portals, because otherwise the display turns dark
// Params: entities with tag ACV_deletePortals_PARAM and a score of ACV_COLOR
process lightUpHorizontalPortals {
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ execute @e[type=armor_stand,name=ACV_PortalDisplay,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ setblock ~ ~2 ~ barrier
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ execute @e[type=armor_stand,name=ACV_PortalDisplay,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ setblock ~ ~2 ~ barrier
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ execute @e[type=armor_stand,name=ACV_PortalDisplay,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ setblock ~ ~2 ~ barrier
  /execute @e[name=ACV_deletePortals_PARAM,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ execute @e[type=armor_stand,name=ACV_PortalDisplay,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ setblock ~ ~2 ~ barrier
}


text:
 func $ceilf64toi32 (
  var %i f64
  ) i32 {
   return (
     ceil i32 f64(dread f64 %i))}

 func $ceilf64toi64 (
  var %i f64
  ) i64 {
   return (
     ceil i64 f64(dread f64 %i))}

 func $ceilf32toi32 (
  var %i f32
  ) i32 {
   return (
     ceil i32 f32(dread f32 %i))}

 func $ceilf32toi64 (
  var %i f32
  ) i64 {
   return (
     ceil i64 f32(dread f32 %i))}

# todo float ceil
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
with(Threads:-Task):
randomize():
roll:=rand(10):
nTasks:=10:
TaskExec:=Array(1..nTasks):
fun:=proc()
    local x;
    global TaskExec;
    TaskExec[_passed]:=1;
    x:=roll();
    if x>5 then
        Return(x);
    end if;
    return x;
end proc:
Start(null,Tasks=[fun,seq([i],i=1..10)]);
print(convert(TaskExec,list));


text:
cdir := "/home/scott/Documents/research/software/hqca/hqca/maple/":
with(LinearAlgebra):
with(QuantumChemistry):
loaddata := readdata(cat(cdir, "h4_qc_spatial.rdm"), 6):
Flatten := proc(x) local n, a, i, j, k, l; `local`(a, n, i, j, k, l); `description`("convert chemists to numpy and flatten an array to form a matrix"); n := Size(x); a := Matrix(1 .. n[1]*n[2], 1 .. n[3]*n[4], datatype = float[8]); for i to n[1] do for j to n[2] do for k to n[3] do for l to n[4] do a[(i - 1)*n[1] + j, (k - 1)*n[3] + l] := x[i, k, j, l]; end do; end do; end do; end do; return a; end proc:
New := Array(1 .. 4, 1 .. 4, 1 .. 4, 1 .. 4, datatype = float[8]):
for i in loaddata[3 .. ()] do
    New[round(i[1]), round(i[3]), round(i[2]), round(i[4])] := i[5];
end do;
pure := Purify2RDM(New, spin_free = true, electron_number = 4, conv_tol = 0.10000000):
ExportMatrix(cat(cdir, "h4_qc_spatial_pure.csv"), Flatten(pure[rdm2]));

max_stars_repo_name: damazz/HQCA

max_issues_repo_name: damazz/HQCA

max_forks_repo_name: damazz/HQCA