text:
# Base Parameter Energy Regressor for Robot based on MDH frames
# Einleitung
# Erstellung einer parameterlinearen und -minimalen Regressorform
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# chain -> Berechnung für eine serielle Struktur (nicht: Baumstruktur)
# floatb -> floating base. Zusätzliche Betrachtung der Dynamikparameter der Basis.
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# energy -> Berechnung der Energie
# regressor -> Regressorform (parameterlinear)
# Authors
# Alexander Tödtheide, toedtheide@irt.uni-hannover.de
# Moritz Schappler, schappler@irt.uni-hannover.de, 2017-02
#
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Quellen
# [HRL_IDR] Skript Humanoid Robotics Lab - Serial Chain Robot Identification
# [GautierKhalil1988] A Direct Determination of Minimum Inertial Parameters of Robots
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [Gautier1990] M. Gautier: Numerical calculation of the base inertial parameters of robots
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
interface(rtablesize=100): # Zur Anzeige von größeren Vektoren
;
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 "../robot_codegen_definitions/robot_env":
printf("Generiere Minimalparameterregressor (Floating Base %s) der Energie für %s\n", base_method_name, robot_name):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name, base_method_name):
# Ergebnisse der Energie laden
read sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_worldframe_par2_maple.m", robot_name, base_method_name):
read sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_linkframe_par2_maple.m", robot_name, base_method_name):
T_floatb := T:
U_floatb := U_grav:
# Die kinetischen und potentiellen Energien aus (2) und (3) stehen ab hier durch T_floatb und U_floatb zur Verfügung.
# Der Parametervektor 'PV2_vec' aus (13) wurde in 'robot_tree_floatb_twist_definitions.mw' aufgestellt.
# Parameterlinearisierung
# Parameterlinearisierung auf Basis von [HRL_IDR] (14) und (15)
# Linearisierung
t_ges := Matrix(1, 10*NL):
u_ges := Matrix(1, 10*NL):
for i to 10*NL do
  t_ges[1,i] := diff(T_floatb,PV2_vec[i,1]);
  u_ges[1,i] := diff(U_floatb,PV2_vec[i,1]);
end do:
# Export
save t_ges, sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_regressor_maple.m", robot_name, base_method_name):
save u_ges, sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_regressor_maple.m", robot_name, base_method_name):
if codegen_act then
  MatlabExport(convert_t_s(t_ges), sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_regressor_matlab.m", robot_name, base_method_name), codegen_opt):
end if:
if codegen_act then
  MatlabExport(convert_t_s(u_ges), sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_regressor_matlab.m", robot_name, base_method_name), codegen_opt):
end if:
# Parameterminimierung
# Minimalparametervekor
# Definiere Parametermatrix
XX := Matrix(NL, 1, PV2_mat[1 .. NL, 1]):
XY := Matrix(NL, 1, PV2_mat[1 .. NL, 2]):
XZ := Matrix(NL, 1, PV2_mat[1 .. NL, 3]):
YY := Matrix(NL, 1, PV2_mat[1 .. NL, 4]):
YZ := Matrix(NL, 1, PV2_mat[1 .. NL, 5]):
ZZ := Matrix(NL, 1, PV2_mat[1 .. NL, 6]):
mX := Matrix(NL, 1, PV2_mat[1 .. NL, 7]):
mY := Matrix(NL, 1, PV2_mat[1 .. NL, 8]):
mZ := Matrix(NL, 1, PV2_mat[1 .. NL, 9]):
m :=  Matrix(NL, 1, PV2_mat[1 .. NL, 10]):
# Rekursive Berechnung der Minimalparameter
# Rekursive Berechnung der Minimalparameter mit [HRL_IDR] (23), siehe Abbildung 1.
# [GautierKhalil1990] (eq. 15, 16)
# Die Sonderfälle (eq. 19, 20) für parallele oder senkrechte Schubachsen scheinen für Floating-Base-Modelle nicht zu gelten, da die dadurch gruppierten Dynamikparameter auf die Basis einen voneinander unabhängigen Einfluss haben.
# Festgestellt durch Vergleich mit numerischer Berechnung der Minimalparameter nach [Gautier1990]
for i from NL by -1 to 2 do
  j := i-1: # Laufvariable für die MDH-Kinematikparameter (Inertialparameter von Segment 0 fangen auch bei 1 an).
  # printf("i=%d, j=%d, sigma=%d\n", i, j, sigma[j,1]):
  if sigma[j,1] = 0 then # Drehgelenk (eq. 15)
    XX[i,1] := XX[i,1]-YY[i,1]:
    XX[i-1,1] := XX[i-1,1] + YY[i,1]+d[j,1]^2*m[i,1]+2*d[j,1]*mZ[i,1]:
    XY[i-1,1] := XY[i-1,1] + a[j,1]*sin(alpha[j,1])*mZ[i,1]+a[j,1]*d[j,1]*sin(alpha[j,1])*m[i,1]:
    XZ[i-1,1] := XZ[i-1,1] - a[j,1]*cos(alpha[j,1])*mZ[i,1]-a[j,1]*d[j,1]*cos(alpha[j,1])*m[i,1]:
    YY[i-1,1] := YY[i-1,1] + cos(alpha[j,1])^2*YY[i,1]+2*d[j,1]*cos(alpha[j,1])^2*mZ[i,1]+(a[j,1]^2+d[j,1]^2*cos(alpha[j,1])^2)*m[i,1]:
    YZ[i-1,1] := YZ[i-1,1] + cos(alpha[j,1])*sin(alpha[j,1])*YY[i,1]+2*d[j,1]*cos(alpha[j,1])*sin(alpha[j,1])*mZ[i,1]+d[j,1]^2*cos(alpha[j,1])*sin(alpha[j,1])*m[i,1]:
    ZZ[i-1,1] := ZZ[i-1,1]+sin(alpha[j,1])^2*YY[i,1]+2*d[j,1]*sin(alpha[j,1])^2*mZ[i,1]+(a[j,1]^2+d[j,1]^2*sin(alpha[j,1])^2)*m[i,1]:
    mX[i-1,1] := a[j,1]*m[i,1]+mX[i-1,1]:
    mY[i-1,1] := mY[i-1,1]-sin(alpha[j,1])*mZ[i,1]-d[j,1]*sin(alpha[j,1])*m[i,1]:
    mZ[i-1,1] := mZ[i-1,1]+cos(alpha[j,1])*mZ[i,1]+d[j,1]*cos(alpha[j,1])*m[i,1]:
    m[i-1,1] := m[i-1,1]+m[i,1]:
  else: # Schubgelenk (eq. 16)
    XX[i-1,1] := XX[i-1,1] + cos(theta[j,1])^2*XX[i,1]
                           - 2*cos(theta[j,1])*sin(theta[j,1])*XY[i,1]
                           + sin(theta[j,1])^2*YY[i,1]:
    XY[i-1,1] := XY[i-1,1] + cos(theta[j,1])*sin(theta[j,1])*cos(alpha[j,1])*XX[i,1]
                           + (cos(theta[j,1])^2-sin(theta[j,1])^2)*cos(alpha[j,1])*XY[i,1]
                           - cos(theta[j,1])*sin(alpha[j,1])*XZ[i,1]
                           - cos(theta[j,1])*sin(theta[j,1])*cos(alpha[j,1])*YY[i,1]
                           + sin(theta[j,1])*sin(alpha[j,1])*YZ[i,1]:
    XZ[i-1,1] := XZ[i-1,1] + cos(theta[j,1])*sin(theta[j,1])*sin(alpha[j,1])*XX[i,1]
                           + (cos(theta[j,1])^2-sin(theta[j,1])^2)*sin(alpha[j,1])*XY[i,1]
                           + cos(theta[j,1])*cos(alpha[j,1])*XZ[i,1]
                           - cos(theta[j,1])*sin(theta[j,1])*sin(alpha[j,1])*YY[i,1]
                           - sin(theta[j,1])*cos(alpha[j,1])*YZ[i,1]:
    YY[i-1,1] := YY[i-1,1] + sin(theta[j,1])^2*cos(alpha[j,1])^2*XX[i,1]
                           + 2*cos(theta[j,1])*sin(theta[j,1])*cos(alpha[j,1])^2*XY[i,1]
                           - 2*sin(theta[j,1])*cos(alpha[j,1])*sin(alpha[j,1])*XZ[i,1]
                           + cos(theta[j,1])^2*cos(alpha[j,1])^2*YY[i,1]
                           - 2*cos(theta[j,1])*cos(alpha[j,1])*sin(alpha[j,1])*YZ[i,1]
                           + sin(alpha[j,1])^2*ZZ[i,1]:
    YZ[i-1,1] := YZ[i-1,1] + sin(theta[j,1])^2*cos(alpha[j,1])*sin(alpha[j,1])*XX[i,1]
                           + 2*cos(theta[j,1])*sin(theta[j,1])*cos(alpha[j,1])*sin(alpha[j,1])*XY[i,1]
                           + sin(theta[j,1])*(cos(alpha[j,1])^2-sin(alpha[j,1])^2)*XZ[i,1]
                           + cos(theta[j,1])^2*cos(alpha[j,1])*sin(alpha[j,1])*YY[i,1]
                           + cos(theta[j,1])*(cos(alpha[j,1])^2-sin(alpha[j,1])^2)*YZ[i,1]
                           - cos(alpha[j,1])*sin(alpha[j,1])*ZZ[i,1]:
    ZZ[i-1,1] := ZZ[i-1,1] + sin(theta[j,1])^2*sin(alpha[j,1])^2*XX[i,1]
                           + 2*cos(theta[j,1])*sin(theta[j,1])*sin(alpha[j,1])^2*XY[i,1]
                           + 2*sin(theta[j,1])*cos(alpha[j,1])*sin(alpha[j,1])*XZ[i,1]
                           + cos(theta[j,1])^2*sin(alpha[j,1])^2*YY[i,1]
                           + 2*cos(theta[j,1])*cos(alpha[j,1])*sin(alpha[j,1])*YZ[i,1]
                           + cos(alpha[j,1])^2*ZZ[i,1]:
    # Keine Sonderregeln (eq. 19,20)
  end if:
end do:
# Parameter ohne Einfluss "Markieren"
# Es kann keine Parameter mehr ohne Einfluss geben, da die Basis-Bewegung die Trägheit aller Segmente anregt.
# Parametervektor in Richtiger Reihenfolge aufstellen und Entfernen von mZ, YY, m
# Anwendung von [HRL_IDR] Regel 3 (S. 5) auf den Parametervektor
# [GautierKhalil1990] (eq. 18)
nt:=add(sigma[i,1],i=1..NJ): # Anzahl Translatorische Gelenke
nr:=NJ-nt: # Anzahl rotatorische Gelenke
;
MPV_n_max := 10+7*nr+4*nt: # ignoriere Abzüge in der Formel
Paramvec := Matrix(MPV_n_max, 1):
# Basis-Parameter direkt übernehmen
Paramvec(1..10,1) := <XX[1,1],XY[1,1],XZ[1,1],YY[1,1],YZ[1,1],ZZ[1,1],mX[1,1],mY[1,1],mZ[1,1],m[1,1]>:
# Formel für Parameter, die durch Gelenke bewegt werden, anwenden
ii := 11: # Indizes 1..10 für Basis s.o.
for i from 1 to NJ do
  if sigma[i,1] = 0 then # Drehgelenk (eq. 15)
    Paramvec[ii, 1] := XX[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := XY[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := XZ[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := YZ[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := ZZ[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := mX[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := mY[i+1, 1]: ii:= ii+1:
  else: # Schubgelenk (eq. 16)
    # Keine Sonderregeln (eq. 19,20)
    Paramvec[ii, 1] := mX[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := mY[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] := mZ[i+1, 1]: ii:= ii+1:
    Paramvec[ii, 1] :=  m[i+1, 1]: ii:= ii+1:
  end if:
end do:
# Entfernung von Nullelementen
# Durch Entfernung der Nullelemente erfolgt 'beta_b' aus [HRL_IDR] (22)
p := 0:
for i to MPV_n_max do
  if Paramvec[i, 1] = 0 then
    p := p+1
  end if:
end do:
# Determine size of matrix dependent on number of Zeroes
Paramvec_size := MPV_n_max-p:
Paramvec2 := Matrix(Paramvec_size, 1):
p := 0:
for i to MPV_n_max do
  # count Zeroes
  if Paramvec[i, 1] = 0 then
    p := p+1:
  else
    # Shift by number of currently found Zeroes
    Paramvec2[i-p, 1] := Paramvec[i, 1]:
  end if:
end do:
printf("Dimension des Minimalparametervektors: %dx1\n", Paramvec_size):
Paramvec2;
# Export - Minimalparametervektor
save Paramvec2, sprintf("../codeexport/%s/tmp/minimal_parameter_vector_floatb_%s_maple", robot_name, base_method_name):
if codegen_act then
   MatlabExport(Paramvec2, sprintf("../codeexport/%s/tmp/minimal_parameter_vector_floatb_%s_matlab.m", robot_name, base_method_name), codegen_opt):
end if;
# Minimal geometrievektor t_i und u_i
# Markieren von t_mZ, t_YY, t_m, u_mZ, u_YY, u_m_j
# Anwendung von [HRL_IDR] Regel 3 (S. 5) auf die Geometrievektoren t und u der Energien T und U
# Die Regel bezieht sich auf  Körper, die durch Drehgelenke bewegt werden. Einträge für den Basiskörper werden daher nicht geändert.
# Die Parameter YY,MZ und M sind immer linear abhängig. Die entsprechenden Spalten der Regressormatrizen werden daher zum Entfernen markiert.
# [GautierKhalil1990], p.369
for i from 1 to NL-1 do
  if sigma[i,1] = 0 then: # Drehgelenk
    t_ges[1,i*10+ 4]:=REMOVE; #YY
    t_ges[1,i*10+ 9]:=REMOVE; #MZ
    t_ges[1,i*10+10]:=REMOVE; #M
    u_ges[1,i*10+ 4]:=REMOVE; #YY
    u_ges[1,i*10+ 9]:=REMOVE; #MZ
    u_ges[1,i*10+10]:=REMOVE; #M
  else:
    t_ges[1,i*10+ 1]:=REMOVE; #XX
    t_ges[1,i*10+ 2]:=REMOVE; #XY
    t_ges[1,i*10+ 3]:=REMOVE; #XZ
    t_ges[1,i*10+ 4]:=REMOVE; #YY
    t_ges[1,i*10+ 5]:=REMOVE; #YZ
    t_ges[1,i*10+ 6]:=REMOVE; #ZZ

    u_ges[1,i*10+ 1]:=REMOVE; #XX
    u_ges[1,i*10+ 2]:=REMOVE; #XY
    u_ges[1,i*10+ 3]:=REMOVE; #XZ
    u_ges[1,i*10+ 4]:=REMOVE; #YY
    u_ges[1,i*10+ 5]:=REMOVE; #YZ
    u_ges[1,i*10+ 6]:=REMOVE; #ZZ
  end if:
end do:
# Entfernungen von markierten Elementen
# Durch Entfernung der markierten Elemente, äquivalent zum Minimalparametervektor beta_b, haben t und u die gleiche Spaltenanzahl wie 'C_b'.
p:=0:
for i from 1 to 10*NL do       #Nullen Zählen um die Matrixgröße von t_ges zu bestimmen
  if t_ges[1,i]=REMOVE and u_ges[1,i]=REMOVE then
    p:=p+1;
  end if
end do:
size_Matrix:=10*NL-p:
t_ges_minpar:=Matrix(1,size_Matrix):
u_ges_minpar:=Matrix(1,size_Matrix):
printf("Dimension der Minimalparameter-Regressormatrix: 1x%d\n", size_Matrix):
p:=0:
for i from 1 to 10*NL do       #Nullen Zählen
   if t_ges[1,i]=REMOVE and u_ges[1,i]=REMOVE then
      p:=p+1;
   else
      t_ges_minpar[1,i-p]:=t_ges[1,i]; #Um die Anzahl der im Iterationsschritt gezählten Nullen verschieben
      u_ges_minpar[1,i-p]:=u_ges[1,i]; #Um die Anzahl der im Iterationsschritt gezählten Nullen verschieben
   end if
end do:
save t_ges_minpar, sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_regressor_minpar_maple.m", robot_name, base_method_name):
save u_ges_minpar, sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_regressor_minpar_maple.m", robot_name, base_method_name):
# Export
if codegen_act then
  MatlabExport(convert_t_s(t_ges_minpar), sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_regressor_minpar_matlab.m", robot_name, base_method_name), codegen_opt):
end if:
if codegen_act then
  MatlabExport(convert_t_s(u_ges_minpar), sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_regressor_minpar_matlab.m", robot_name, base_method_name), codegen_opt):
end if:


text:
Fun2 := module()

  description "fun2";

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

  export fun, grad, hess, F, JF;

  local ModuleLoad,
        ModuleUnLoad,
        f1, f2, f12,
        CST;

  uses LinearAlgebra;

  ModuleLoad := proc()
    CST := 53/27;
    NULL;
  end proc;

  ModuleUnLoad := proc()
    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();

  f1  := unapply( sin(3*x+5*y)+x+exp(3*(x-y))-1, (x,y) );
  f2  := unapply( (CST*x^2+2*y^2+cos(x-y)^2), (x,y) );
  f12 := unapply( f1(x,y)^2+f2(x,y)^2, (x,y) );

  fun := proc( X::Vector )
    local x, y;
    x := X[1];
    y := X[2];
    f12(x,y);
  end proc;

  grad := proc( X::Vector )
    local x, y;
    x := X[1];
    y := X[2];
    eval(<D[1](f12)(x,y),D[2](f12)(x,y)>);
  end proc;

  hess := proc( X::Vector )
    local x, y, d11, d12, d22;
    x   := X[1];
    y   := X[2];
    d11 := eval(D[1,1](f12)(x,y));
    d12 := eval(D[1,2](f12)(x,y));
    d22 := eval(D[2,2](f12)(x,y));
    <<d11,d12>|<d12,d22>>;
  end proc;

  F := proc( X::Vector )
    local x, y;
    x := X[1];
    y := X[2];
    <f1(x,y),f2(x,y)>;
  end proc;

  JF := proc( X::Vector )
    local x, y;
    x := X[1];
    y := X[2];
    <<D[1](f1)(x,y),D[1](f2)(x,y)>|
     <D[2](f1)(x,y),D[2](f2)(x,y)>>;
  end proc;

end module:

ext: mpl

max_stars_repo_name: ebertolazzi/NLtoolbox

max_issues_repo_name: ebertolazzi/NLtoolbox

max_forks_repo_name: ebertolazzi/NLtoolbox


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

a  := [0.4581652932831429, 2.217058676663745,  0.7405551735357053, 0.01968227878617998 ]:
ap := [0.119086804055547,  0.6157402568883345, 0.1574201515892867, 0.003532336663397157]:
b  := [1.0000000000000000, 4.504130959426697,  1.110667363742916,  0.02359291751427506 ]:
bp := [0.000000000000000,  0.2673612973836267, 0.2052004607777787, 0.004200005045691381]:

f  := (rs, zeta) ->
  - add((a[i] + f_zeta(zeta)*ap[i])*rs^(i-1), i=1..4) /
    add((b[i] + f_zeta(zeta)*bp[i])*rs^i,     i=1..4):


text:
#sample
tinput := "5483143223
2745854711
5264556173
6141336146
6357385478
4167524645
2176841721
6882881134
4846848554
5283751526";

#toy
tinput := "11111
19991
19191
19991
11111";

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

ogrid := map((parse~),StringTools:-Explode~(StringTools:-Split(input)));
gridw := nops(ogrid[1]); gridl := nops(ogrid);

neighbors := proc(pt) # include diagonals
local x, y, out;
global gridw, gridl;
    out := NULL;
    (x,y) := pt[];
    if x > 1 then
        out := out, [x-1, y];
        if y > 1 then
            out := out, [x-1, y-1];
        end if;
    end if;
    if y > 1 then
        out := out, [x, y-1];
        if x< gridw then
           out := out, [x+1, y-1];
        end if;
     end if;
    if x < gridw then
        out := out, [x+1, y];
        if y < gridl then
           out := out, [x+1, y+1];
        end if;
    end if;
    if y < gridl then
        out := out, [x, y+1];
        if x > 1 then
           out := out, [x-1, y+1];
        end if;
     end if;
    return [out];
end proc:

flash := proc(x,y)
local n, nb := neighbors([x,y]);
global grid, flashed;
    if flashed[x,y] <> 0 then
    # flashed is a table to track that we flashed only once
        return;
    end if;
    flashed[x,y] := 1;
    for n in nb do
    # increment neighbors and maybe trigger their flashes
        grid[n[]] += 1;
        if grid[n[]] > 9 then
            flash(n[]);
        end if;
     end do;
end proc: # flash
totalflashes := 0;
grid := Array(ogrid);

for dd from 1 to 1000 while add(grid) <> 0 do
    grid := grid +~ 1; # increment grid
    flashed := table(sparse=0); # clear flash tracking table
    # do the flashes - flash may trigger flashes earlier too
    for i from 1 to gridw do
        for j from 1 to gridl do
            if grid[i,j] > 9 then
                flash(i,j);
            end if;
        end do;
    end do;
    # count and reset flashed
    for i from 1 to gridw do
        for j from 1 to gridl do
            if grid[i,j] > 9 then
                grid[i,j] := 0;
                if i <= 100 then
                    totalflashes+=1;
                end if;
            end if;
        end do;
    end do;

end do: # days

answer1 := totalflashes;
answer2 := (dd-1);


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

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

(* type: work_mgga_c *)
(* prefix:
  mgga_c_m05_params *params;

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

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

$include "b97.mpl"

m05_comp := (rs, z, spin, xs, ts) ->
  + lda_stoll_par(f_pw, rs,  z,  1)
  * b97_g(params_a_gamma_ss, params_a_css, xs)
  * Fermi_D(xs, ts):

(* The parallel and perpendicular components of the energy *)
m05_fpar  := (rs, z, xs0, xs1, ts0, ts1) ->
  + m05_comp(rs,  z,  1, xs0, ts0)
  + m05_comp(rs, -z, -1, xs1, ts1):

m05_fperp := (rs, z, xs0, xs1, ts0, ts1) ->
  + lda_stoll_perp(f_pw, rs,  z)
  * b97_g(params_a_gamma_ab, params_a_cab, sqrt(xs0^2 + xs1^2)):

f_m05 := (rs, z, xs0, xs1, ts0, ts1) ->
  + m05_fpar (rs, z, xs0, xs1, ts0, ts1)
  + m05_fperp(rs, z, xs0, xs1, ts0, ts1):

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

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

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

pbeint_mu := s -> params_a_muGE + (params_a_muPBE - params_a_muGE)* \
   params_a_alpha*s^2/(1 + params_a_alpha * s^2):

(* this is the gga_x_pbe expression *)
pbeint_f0 := s -> 1 + params_a_kappa * (1 - params_a_kappa/(params_a_kappa + pbeint_mu(s)*s^2)):
pbeint_f  := x -> pbeint_f0(X2S * x):

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


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

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

(* type: mgga_exc *)

(* prefix:
  mgga_x_mvsb_params *params;

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

$include "mgga_x_mvs.mpl"

mvsb_beta := (t, x) -> mvs_alpha(t, x)*K_FACTOR_C/(t - K_FACTOR_C):

mvsb_f := (x, u, t) -> (1 + params_a_k0*mvs_fa(mvsb_beta(t,x)))
       / (1 + params_a_b*(X2S*x)^4)^(1/8):

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


text:
# # Computation for the equivariant cohomology of G/B
# # by Shizuo Kaji

# it requires coxeter and weyl package by J.Stembridge
# http://www.math.lsa.umich.edu/~jrs/maple.html

# symbols
# e(i), a[i]: elements in H^*(BT) in abstract coordinate and root coordinate
# x[i], z[i]: elements in H^*(G/T) in abstract coordinate and root coordinate

# look at the demo worksheet "BGGT_example.mw" for the details.

with(coxeter); with(weyl);with(WeylOps);

# decompose a word into products of smaller words
# used in "double"
decomp := proc(weyl::list) local L, i, temp;
 if nops(weyl)=1 then return {[weyl]}; end if;
 temp:={};
 for i in descent(weyl) do
   for L in procname(weyl.[i]) do
     if nops(L)=1 then
	   temp:={op(temp), [op(L),[i]], [L[1].[i]]};
     else
       temp:={op(temp), [op(L),[i]], [op(1..-2,L),L[-1].[i]]};
     end if;
   end do;
 end do;
 return temp;
end proc:

# Enumerate subwords of "longweyl" which is equal to "weyl"
subwords := proc(weyl::list, longweyl::list) local i,L,W; global R;
  W := {};
  for L in combinat[choose]([seq(i,i=1..nops(longweyl))],nops(weyl)) do
    if reduce([seq(longweyl[i],i=L)], R)=weyl then
	   W:=W union {L};
    end if;
  end do;
  return W;
end proc;

## Utility Functions
# convert between e-poly and z-poly
e2x := proc(f::polynom);
	return eval(f,[seq(e(j)=x[j], j = 1 .. dimr(R))]);
end proc:
x2e := proc(f::polynom);
	return eval(f,[seq(x[j]=e(j) , j = 1 .. dimr(R))]);
end proc:
e2u := proc(f::polynom);
	return eval(f,[seq(e(j)=u[j], j = 1 .. dimr(R))]);
end proc:

# convert between a-poly and z-poly
a2z := proc(f::polynom);
	return eval(f,[seq(a[j]=z[j], j = 1 .. rank(R))]);
end proc:

# BGG operator on GKM poly
BGG_GKM := proc(rt::integer, f::list(polynom), LR::symbol := 'R') local j,k,ff,v,u; global reduw, R, B, optype; option remember;
  ff := [seq(0,j=1..nops(reduw))];
  if LR='R' then
	for j from 1 to nops(reduw) do
		u := reflect(seq(B[reduw[j][k]],k=1..nops(reduw[j])),B[rt]);
		member(reduw[j].[rt],reduw,'v');
        if optype = 'e' then
            ff[j] := -simplify((f[j] - f[v])/u);
        else
            ff[j] := -simplify((f[j] - f[v])/e2a(u));
        end if;
	end do;
  elif LR='L' then
	for j from 1 to nops(reduw) do
		member([rt].reduw[j], reduw,'v');
        if optype = 'e' then
            ff[j]:=simplify((f[j]-subs({seq(e(i)=reflect(B[rt],e(i))
                                            ,i=1..dimr(R))},f[v]))/B[rt]);
        else
            ff[j]:=simplify((f[j]-subs({seq(a[i]=e2a(reflect(B[rt],B[i]))
                                            ,i=1..nops(B))},f[v]))/a[rt]);
        end if;
	end do;
  end if;
  return ff;
end proc:

# BGG on double polynomials
BGGT_e := proc(rt::integer, f::polynom, LR::symbol := 'R')
  if LR='R' then
    return -simplify((f - act([rt], f))/e2x(B[rt]));
  else
    return simplify((f - act([rt], f, 'L'))/B[rt]);
  end if;
  return -infinity;
end proc;
#
BGGT_a := proc(rt::integer, f::polynom, LR::symbol := 'R')
  if LR='R' then
      return -simplify((f - act([rt], f))/z[rt]);
  else
      return simplify((f - act([rt], f, 'L'))/a[rt]);
  end if;
  return -infinity;
end proc;

# BGG operator for poly in "X[weyl]"
BGGT_X := proc(rt::integer, f::polynom, LR::symbol := 'R') local ff,w,term;
    if degree(f) < 1 then
        return 0;
    elif LR = 'L' then
        if optype = 'e' then
            return (f-act([rt],f,'L'))/B[rt];
        else
            return (f-act([rt],f,'L'))/a[rt];
        end if;
    elif LR ='R' then
        if type(f, 'indexed') and op(0,f)=`X` then
            w:=op(1,f).[rt];
            if nops(w)=0 then return 1;
            elif nops(w) < nops(op(1,f)) then return(X[w]);
            else return 0;
            end if;
        elif op(0, f) = `+` then
            return expand(map2(procname, rt, f));
        elif op(0, f) = `*` then
            if type(op(1,f), 'rational') then
                return op(1,f)*procname(rt, f/op(1,f));
            end if;
            term[1] := op(1, f);
            term[2] := mul(op(i, f), i = 2 .. nops(f));
            return expand(procname(rt, term[1])*term[2]
                          +ref_rt(rt, term[1])*procname(rt, term[2]));
        elif op(0, f) = `^` then
            term[1] := op(1, f);
            return expand(procname(rt, term[1])*add(term[1]^(i-1)*act([rt], term[1])^(op(2,f)-i),i=1..op(2,f)));
        end if;
        print("error in BGGT_X:",f);
    end if;
end proc;

# BGG operator for weyl words
BGGweylT := proc(weyl::list, f::{list(polynom),polynom}, LR::symbol := 'R')
                 global optype; local i,temp,func;
	temp := f;
	if type(f, list(polynom)) then
	   func := 'BGG_GKM';
    elif has(f,`X`) then
        func := 'BGGT_X';
    elif type(f, polynom) then
	   func := cat('BGGT_',optype);
    else
        print("error in BGGWeylT:",f);
        return -infinity;
    end if;
	for i from nops(weyl) to 1 by -1 do
		temp := func(weyl[i], temp, LR);
	end do;
	return temp;
end proc:

# Weyl group action on polynomials
act := proc(weyl::list, f::ratpoly, LR::symbol := 'R') local ff,t,rt; global R,optype;
    if nops(weyl)=0 then
        return f;
    elif nops(weyl)>1 then
        return procname(weyl[1..-2],procname([weyl[-1]],f,LR),LR);
    elif nops(weyl)=1 then
        rt := weyl[-1];
        if LR='R' then
            ff:=eval(f, [seq(x[i] = subs([seq(e(j) = x[j], j = 1 .. dimr(R))],
                               reflect(B[rt], e(i))), i = 1 .. dimr(R))]);
            ff:=eval(ff, [seq(z[i] = add(root_coords(reflect(B[rt],B[i]),R)[j]*z[j],
                                               j = 1 .. rank(R)), i = 1 .. rank(R))]);
            return eval(ff,[seq(t=ref_weyl(rt,op(1,t),'R'), t=indetsX(ff))]);
       elif LR='L' then
            ff:= eval(f, [seq(e(i)=reflect(B[rt], e(i)),
                              i = 1 .. dimr(R))]);
            ff := eval(ff, [seq(a[i] = add(root_coords(reflect(B[rt], B[i]),R)[j]*a[j],
                                            j = 1 .. rank(R)), i = 1 .. rank(R))]);
            return eval(ff,[seq(t=ref_weyl(rt,op(1,t),'L'), t=indetsX(ff))]);
        end if;
    end if;
    print("error in act:",weyl,f);
    return -infinity;
end proc;

# enumerate X in a polynomial
indetsX := proc(f) local t,V;
    V:=[];
    for t in indets(f) do
        if type(t,'indexed') and op(0,t)=`X` then
            V:=[op(V),t];
        end if;
    end do;
    return V;
end proc:

# Weyl group action on X
ref_weyl := proc(rt::integer,weyl::list, LR::symbol := 'R') local b,j,v,u,w,temp; global B;
    if LR = 'R' then
        w := weyl.[rt];
        if nops(w) > nops(weyl) then
            return X[weyl];
        end if;
        temp := X[weyl]+act(w,B[rt],'L')*X[w];
        for b in pos_roots(B) do
            v := reflect(b,interior_pt(B));
            v := reflect(seq(B[j],j=w),v);
            vec2fc(v,R,'u');    # u = w.s_b
            if nops(u)=nops(weyl) then
                temp := temp - iprod(2/iprod(b,b)*b,B[rt])*X[u];
            end if;
        end do;
        return temp;
    elif LR = 'L' then
        w := [rt].weyl;
        if nops(w) > nops(weyl) then
            return X[weyl]
        else
            return -B[rt]*X[w]+X[weyl];
        end if;
    end if;
    print("error in ref_weyl:",rt,weyl);
    return -infinity;
end proc:

# evaluate poly at "w" by putting x[i]=w^{-1}e[i]
ev := proc(weyl::list, f::polynom) local ff; global R,optype;
      ff := eval(act(weyl,f), [seq(x[i] = e(i), i = 1 .. dimr(R))]);
      return eval(ff, [seq(z[i] = a[i], i = 1 .. rank(R))]);
end proc:

# orbit sum
orbit_sum := proc(f::polynom) global reduw,R;
  return expand(add(ev(weyl,f), weyl=reduw)/size(R));
end proc:

# push-forward in cohomology
push := proc(f::polynom) local d,weyl; global reduw,R,P;
    d := mul(e2x(i), i = P);
    return expand(ev([],(-1)^(nops(P))*simplify(add(act(weyl,f/d), weyl=reduw)))):
end proc:

# convert double polynomials to Schubert basis
P2X := proc(f::polynom) local weyl, temp, i; global R, redu;
    temp := ev([],f);
    for i to min(degree(f),nops(redu)) do
      temp := temp+add(
    	  ev([],BGGweylT(weyl, f))*X[weyl], weyl=redu[i]);
    end do;
	return col(temp);
end proc:

# Convert GKM polynomials to Schubert basis using upper triangularity
GKM2X := proc (f::list(polynom)) local i, j, it, lc, ff; global reduw;
   for it to nops(f) do
     if f[it] <> 0 then break; end if;
   end do;
   if nops(f) < it then return 0; end if;
   lc := normal(f[it]/GKM(reduw[it])[it]);
   ff := normal(zip(`-`, f, zip(`*`, lc, GKM(reduw[it]))));
   return col(lc*X[reduw[it]]+procname(ff));
end proc;

# Convert GKM polynomials to Schubert basis using right BGG operators
GKM2X2 := proc (f::list(polynom)) local i, deg, temp, weyl; global reduw;
    temp := f[1];
    deg := max(map(degree,f));
    for i to min(deg,nops(redu)) do
        temp := temp+add(BGGweylT(weyl, f)[1]*X[weyl], weyl=redu[i]);
    end do;
    return col(temp);
end proc;

# convert poly in "X" to Schubert basis
X2X := proc(f::polynom);
	return GKM2X(X2GKM(f));
end proc:

# Convert poly in "X" into GKM-poly
X2GKM := proc (f::polynom) local i, weyl; global reduw;
    return [seq(eval(f, [seq(X[weyl] = GKM(weyl)[i], weyl = reduw)]),i=1..nops(reduw))] ;
end proc;

# convert double poly into GKM-poly
P2GKM := proc(f::polynom);
    return ([seq(simplify(ev(w, f)), w = reduw)])
end proc:

# GKM polys (Billey's formula)
GKM :=proc(weyl::list, upto:=-1) option remember; global reduw;
 return [seq(GKM_internal(weyl,v),v=reduw[1..upto])];
end proc:

GKM_internal := proc(weyl::list, longweyl::list) local i,j,k,u,L,K,M; global R,B,optype;
  K := 0;
  for L in subwords(weyl, longweyl) do
    M := 1;
    for i from 1 to nops(weyl) do
      u := reflect(seq(B[longweyl[j]],j=1..L[i]-1), B[longweyl[L[i]]]);
	  if optype = 'e' then
	  	  M := M * u;
	  else
		  M := M * e2a(u);
	  end if;
	end do;
    K := K + M;
  end do;
  return simplify(K);
end proc:

# initialize for Lie type
setupT := proc (lie_type, roots::set :={}, upto::integer := 0) local i;
     global R, B, P, S, ST, redu, udim, reduw, optype, optypeK, w0;
  forget(GKM);
  forget(BGG_GKM);
  forget(refK); forget(strong_bruhat);
  unassign('S','ST');

  R := lie_type;
  B := base(R);
  P := pos_roots(R);
    optype := 'e';
    if R = cat('A',rank(R)) then
        optypeK := 'x'
    else
        optypeK := 'u'
    end if;
  udim := `if`(upto>0,upto,nops(P));
  redu := PReducedWd(R, roots, udim);
  reduw := [[],seq(op(i),i=redu)];
  w0 := longest_elt(R);
  return R;
end proc;

# TopSchubert poly
top := proc () local k, j, i, temp; global R,P;
   temp := 1;
   if optype='e' then
	if R = G2 then
	    return (-3*e2x(w[1])^6+2*e2x(w[1])^5*e2x(w[2]))/6;
    elif R = cat('A',rank(R)) then
	     return mul(x[i]^(rank(R)-i+1), i = 1 .. rank(R));
	elif R = cat('B',rank(R)) then
    	return mul(x[i]^(rank(R)-i), i = 1 .. rank(R))*mul(c(i,`x`)/2, i = 1 .. rank(R));
    elif lie_type = cat('C',rank(R)) then
		 return mul(x[i]^(2*rank(R)-2*i+1), i = 1 .. rank(R));
    elif R = E7 then
		temp := mul(x[i]-c(1,`x`), i = 1 .. 7);
		for i to 6 do temp := temp*mul(x[i]-x[j], j = i+1 .. 7) end do;
	    for i to 5 do
		 for j from i+1 to 6 do
			temp := temp*mul(x[i]+x[j]+x[k]-c(1,`x`), k = j+1 .. 7);
	     end do;
		end do;
		return temp;
    elif R = E8 then
		temp := mul(t[i], i = 1 .. 8);
		for i to 7 do temp := temp*mul(t[i]-t[j], j = i+1 .. 8) end do;
		for i to 7 do temp := temp*mul(t[i]+t[j]-t[0], j = i+1 .. 8) end do;
		for i to 6 do
			for j from i+1 to 7 do
				temp := temp*mul(t[i]+t[j]+t[k]-t[0], k = j+1 .. 8)
			end do;
		end do;
		return temp;
    else
		return mul(e2x(i), i = P)/mul(i, i=degrees(R));
    end if;
   else
	if R = G2 then
   		return z[1]^5*z[2]/2;
	else
		return mul(a2z(-e2a(i)), i = P)/mul(i, i=degrees(R));
    end if;
   end if;
end proc;

# Top double poly
topT := proc() global R,optype,S,w0;  local u, i, j, k, D, temp;
  if optype='e' then
     if R = cat('A',rank(R)) then
	     return mul(mul(x[i]-e(j), j = 1 .. rank(R)-i+1), i = 1 .. rank(R)+1);
	 elif R = cat('B',rank(R)) then
	    #  D:=Matrix(rank(R));
		#  for i from 1 to rank(R) do
		#     for j from 1 to rank(R) do
		#          k:=rank(R)+1+j-2*i;
		# 		 if k<=rank(R) and k>=0 then
		# 		   D[i,j]:=(c(k,`x`)+eval(c(k,`t`),[seq(t[h]=e(h),h=1..rank(R))]))/2;
		# 		 end if;
		# 	end do;
		# end do;
		# return act(longest_elt(R),mul(mul(x[i]-e(j), j = 1 .. i-1), i = 2 .. rank(R))
		#   *LinearAlgebra[Determinant](D));
   return  mul(mul(x[i]-e(j), j = 1 .. i), i = 1 .. rank(R))
                  *mul(mul(x[i]+e(j), j = 1 .. i-1), i = 1 .. rank(R))/(-2)^rank(R);

	 elif R = cat('C',rank(R)) then
	    #  D:=Matrix(rank(R));
		#  for i from 1 to rank(R) do
		#     for j from 1 to rank(R) do
		#          k:=rank(R)+1+j-2*i;
		# 		 if k<=rank(R) and k>=0 then
		# 		   D[i,j]:=(c(k,`x`)+eval(c(k,`t`),[seq(t[h]=e(h),h=1..rank(R))]));
		# 		 end if;
		# 	end do;
		# end do;
		# return act(longest_elt(R),mul(mul(x[i]-e(j), j = 1 .. i-1), i = 2 .. rank(R))
		#   *LinearAlgebra[Determinant](D));
        return   mul(mul(x[i]-e(j), j = 1 .. i), i = 1 .. rank(R))
                  *mul(mul(x[i]+e(j), j = 1 .. i-1), i = 1 .. rank(R));
	 elif R = cat('D',rank(R)) then
	     D:=Matrix(rank(R));
		 for i from 1 to rank(R) do
		    for j from 1 to rank(R) do
		         k:=rank(R)+j-2*i;
				 if k<=rank(R) and k>=0 then
				     D[i,j]:=(c(k,`x`)+eval(c(k,`t`),[seq(t[h]=e(h),h=1..rank(R))]))/2;
				 end if;
			end do;
		end do;
		return act(w0,mul(mul(x[i]-e(j), j = 1 .. i-1), i = 2 .. rank(R))
                   *LinearAlgebra[Determinant](D));
     elif R = G2 then
         return ((x[1]-x[2]+e2-e1)/2)*(x[2]-x[1]-e2+e3)*(x[2]-x[1]-e1+e3)
         *(x[1]-x[2]-e2+e3)*(x[1]-x[2]-e1+e3)*(-x[1]+2*x[2]-x[3]+e2-2*e1+e3);
#	     return -(x[3]-(2/3)*e3+(1/3)*e1-(2/3)*e2)*((3/2)*x[1]-(3/2)*e1)*(3*x[1]-3*e2)
#         *(3*x[1]-3*e3)*(x[2]-(2/3)*e1-(2/3)*e2+(1/3)*e3)*(x[2]-(2/3)*e3+(1/3)*e1-(2/3)*e2);
	 elif type(name_of(R),'indexed') then       # type I2
	      D:=[]; temp:=x[2]-e2;
	      for i from 1 to trunc(op(1,name_of(R))/2-1) do
		     D:=[1,2,op(D)];
		     temp:=temp*(x[2]-act(D,e2,'L'));
		  end do;
		  D:=[];
		  for i from 1 to trunc((op(1,name_of(R))-1)/2) do
		     D:=[2,1,op(D)];
		     temp:=temp*(x[2]-act(D,e2,'L'));
		  end do;
	      return temp*(x[1]+act(w0,e1,'L'));
	 else
	    return double(w0);
	 end if;
  elif optype='a' then
     if R = G2 then
         return ((-z[1]+a[1])/2)*(z[1]+a[1]+a[2])*(z[1]+2*a[1]+a[2])
         *(-z[1]+a[1]+a[2])*(-z[1]+2*a[1]+a[2])*(-z[2]+3*a[1]+a[2]);
     elif R = cat('A',rank(R)) then
        for i from 1 to rank(R)+1 do
	       u[i] := e2a(-add(e(j), j = 1 .. rank(R)+1)/(rank(R)+1) + e(i));
	    end do;
        return mul(mul(eval(u[i],[seq(a[k]=z[k],k=1..rank(R))])-u[j], j = 1 .. rank(R)-i+1),
	        i = 1 .. rank(R)+1);
     else
	    return double(w0);
	 end if;
  end if:
end proc:

# enumerate Schubert polynomials
schubertpolys := proc () global S,reduw,w0; local w,f;
  f:=top();
  for w in reduw do
	 S[w] := simplify(BGGweylT(inv(w).w0,f));
  end do;
  print("The top schubert poly is",f);
end proc:

# enumerate Double Schubert polynomials
schubertpolysT := proc () global optype,ST,reduw,w0; local w,f;
  f:=topT();
  for w in reduw do
	 ST[w] := simplify(BGGweylT(inv(w).w0,f));
  end do;
  print("The top double schubert poly is",f);
end proc:

# single to double formula
# compute double Schubert poly "ST[weyl]" from single Schubert polys "S[]"
double:=proc(weyl::list) global S; local temp,L;
  if optype = 'e' then
      return add((-1)^(nops(L))*mul(x2e(S[L[i]]),i=1..nops(L)-1)
                 *(x2e(S[L[-1]])-S[L[-1]]),L=decomp(weyl));
  else
      return add((-1)^(nops(L))*mul(z2a(S[L[i]]),i=1..nops(L)-1)
                 *(z2a(S[L[-1]])-S[L[-1]]),L=decomp(weyl));
  end if;
end proc:

# double to single
single:=proc(f::polynom) local ff; global R;
    ff := eval(f,[seq(e(j) = 0, j = 1 .. dimr(R))]);
    return eval(ff, [seq(a[j] = 0, j = 1 .. rank(R))]);
end proc:

# Cauchy formula for type A
Cauchy:=proc() global reduw,w0 ;local weyl;
 return add((-1)^(nops(weyl))*x2e(S[weyl])*S[weyl.w0],weyl=reduw);
end proc:

# Chevalley formula computes "X[weyl]X[rt]"
Chev := proc(weyl::list,rt::integer) local j,v,u,temp; global w,P;
	 temp := 0;
     for v in P do
       	 vec2fc(reflect(seq(B[j], j = weyl),v,interior_pt(R)),R,'u');
         if nops(u)=nops(weyl)+1 then
		   temp := temp + iprod(2/iprod(v,v)*v,w[rt])*X[u];
         end if;
     end do;
     return (w[rt]-reflect(seq(B[j], j = weyl),w[rt]))*X[weyl]+temp;
end proc:

# collect terms in Schubert presentation
col := proc(f::algebraic) local g; global reduw;
    g := collect(simplify(f),[seq(X[weyl],weyl=reduw)]);
    return collect(g,[seq(XK[weyl],weyl=reduw)]);
end proc:

max_stars_repo_name: shizuo-kaji/Maple-flag-cohomology

max_issues_repo_name: shizuo-kaji/Maple-flag-cohomology

max_forks_repo_name: shizuo-kaji/Maple-flag-cohomology


text:
program assignboolean;	{sample22}
var x y : boolean;
begin
 x := true;
 y := false;
end.

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: shouth/LanguageProcessing

max_forks_repo_name: Hiroya-W/lang-processing


text:
# Domain types which are registered/unregistered with TypeTools.
# Note that the types have to be quoted (additionally to the quote one would
# normally place on certain types) to work properly
local DomainTypes := table(
       # Domain bounds
       [(DomBoundVar = 'DomBoundVar_type')
       ,(DomBoundRange = 'DomBoundRange_type')
       ,(DomBoundBinder = ''DInto(DomBoundVar, DomBoundRange, DomBoundKind)'' )
       ,(DomBoundKind   = 'And(name, satisfies(nm->assigned(Domain:-ExtBound[nm])))' )
       ,(DomBound       = ''And(specfunc(DBound)
                               ,{anyfunc(list(DomBoundBinder))
                                ,anyfunc(list(DomBoundBinder),DomCtx)
                                ,anyfunc(list(DomBoundBinder),DomCtx,table) })'' )
       # Domain shape
       ,(DomConstrain = 'specfunc(relation, `DConstrain`)' )
       ,(DomSum       = 'specfunc(DomShape, `DSum`)' )
       ,(DomSplit     = ''DSplit(Partition(DomShape))'' )
       ,(DomInto      = ''Or(DInto(DomBoundVar, DomBoundRange, DomShape)
                            ,DInto(DomBoundVar, DomShape) )'' )
       ,(DomShape     = 'Or( DomConstrain, DomSum, DomSplit, DomInto )' )
       # Domain
       ,(DomCtx = ''set({relation, `::`})'')
       ,(HDomain = ''DOMAIN(DomBound, DomShape)'' )
       # Maybe domain
       ,(DomNoSol  = ''Not(freeof(`DNoSol`))'' )
       ,(HDomain_mb = ''Or(HDomain, DOMAIN(DomBound, DomNoSol))'' )
       ] );

local DomBoundVar_type := proc(nm)
    local ixs := [indices(Domain:-ExtBound, nolist)], vs;
    ormap(i->type(nm,Domain:-ExtBound[i]:-VarType), ixs);
end proc;

local DomBoundRange_type := proc(nm)
    local ixs := [indices(Domain:-ExtBound, nolist)], vs;
    ormap(i->type(nm,Domain:-ExtBound[i]:-RangeType), ixs);
end proc;

export DOMAIN := proc()
    local errs := indets([args[2..-1]], DomNoSol);
    if errs <> {} then
        'procname'(args[1], 'DNoSol'(map(op,errs)[]));
    end if;
    'procname'(args);
end proc;

export DBound := proc(a)
    local b,cs;
    b := `if`(nargs>=2,[args[2]],[{}])[];
    cs := `if`(nargs>=3,[args[3..-1]],[])[];
    'procname'(a,b,cs) ;
end proc;

export DConstrain := proc()
    local as := {args};
    if false in as then return DSum() end if;
    as := remove(x->x::identical(false),as);
    'procname'(op(as));
end proc;

export DSum := proc()
    local as := [args];
    as := subsindets(as, specfunc(`DSum`), xs->
                     map(x->if op(0,x)=`DSum` then op(x) else x end if,
                         xs));
    if nops(as) = 1 then return op(1,as) end if;
    'procname'(op(as));
end proc;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


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

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

(* type: gga_exc *)

(* constants from text in section 2 *)
pbetrans_kappa_pbe := 0.814: (* This is probably a misprint in the manuscript as KAPPA_PBE is 0.8040 *)
pbetrans_kappa_revpbe := 1.227:
pbetrans_mu := 0.219:

(* parameters from section 4 *)
pbetrans_alpha := 2*(3*Pi^2)^(1/3):
pbetrans_beta := 3:

(* eq 3 *)
pbetrans_fermi := s -> 1/(1+exp(-pbetrans_alpha*(s-pbetrans_beta))):
(* eq 5 *)
pbetrans_kappa := s -> (1-pbetrans_fermi(s))*pbetrans_kappa_revpbe + pbetrans_fermi(s)*pbetrans_kappa_pbe:
(* eq 4 *)
pbetrans_f0 := s -> 1 + pbetrans_kappa(s)*(1 - pbetrans_kappa(s)/(pbetrans_kappa(s) + pbetrans_mu*s^2)):

pbetrans_f  := x -> pbetrans_f0(X2S*x):

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


text:
func $foo1 (var %p1 i8, var %p2 i16, var %p3 i32, var %p4 i64, ...)i32 {
   return (
       constval i32 33)}
func $foo2 (var %p1 u8, var %p2 u16, var %p3 u32, var %p4 u64,...) i32 {
   return (
       constval i32 33)}
func $foo4 (var %p1 f32, var %p2 f64, ... ) i32 {
   return (
       constval i32 33)}

func $foo5 (var %p1 a32, var %p2 a64,...) i32
func $foo6 (var %p1 c64, var %p2 c128, ...) i32

func $foo7 (var %p1 <* i8>, var %p2 <* i16>, var %p3 <* i32>, var %p4 <* i64>) i32 {
   return (
       constval i32 33)}
func $foo8 (var %p1 <* u8>, var %p2 <* u16>, var %p3 <* u32>, var %p4 <* u64>) i32 {
   return (
       constval i32 33)}
func $foo9 (var %p1 <* u1>) i32 {
   return (
       constval i32 33)}
func $foo10 (var %p1 <* f32>, var %p2 <* f64>) i32 {
   return (
       constval i32 33)}
func $foo11 (var %p1 <* void>)void

func $foo12 (var %p1 <[10] i8>) i32 {
   return (
       constval i32 33)}
func $foo13 (var %p1 <[10][32] i32>) i32 {
   return (
       constval i32 33)}
func $foo14 (var %p1 <[10][32][43] i32>) i32 {
   return (
       constval i32 33)}
func $foo15 (var %p1 <[10][32][43][54][56][78][90][122][223][334] i32>) i32 {
   return (
       constval i32 33)}

type $SS <struct {@f1 i32, @f2:3 i8, @f3:5 i8}>
func $foo16 (var %x <$SS>) i32 {
   return (
       constval i32 33)}

 # 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:
toy_input := "start-A
    start-b
    A-c
    A-b
    b-d
    A-end
    b-end":

toy2_input:="dc-end
    HN-start
    start-kj
    dc-start
    dc-HN
    LN-dc
    HN-end
    kj-sa
    kj-HN
    kj-dc":

sample_input := "fs-end
    he-DX
    fs-he
    start-DX
    pj-DX
    end-zg
    zg-sl
    zg-pj
    pj-he
    RW-he
    fs-DX
    pj-RW
    zg-RW
    start-pj
    he-WI
    zg-he
    pj-fs
    start-RW":

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

nodes := map(s->{StringTools:-Split(s,"-")[]}, StringTools:-Split(input)):
G := GraphTheory:-Graph( {nodes[]} );
small := select(StringTools:-IsLower,{map(op, nodes)[]}) minus {"end"};

numpaths := proc(a, s2:=FAIL, visited_s2:=false,
                 excluding:={"start"}, path:=[] )
global G, small;
local newex, newvs2;

    if a = "end" then
       if s2<>FAIL and numboccur(path, s2) <> 2 then
          return 0;
       end if;
#       print([path[],"end"]);
       return 1;
    end if;

    local nb := {GraphTheory:-Neighborhood(G, a)[]} minus excluding;
    if a = s2 and visited_s2 = false then
         newex := excluding;
         newvs2 := true;
    elif a in small then
         newex := excluding union {a};
         newvs2 := visited_s2;
    else
         newex := excluding;
         newvs2 := visited_s2;
    end if;

    return add(numpaths(i, s2, newvs2, newex, [path[],a]), i in nb);

 end proc:

answer1 := numpaths("start");
answer2 := answer1 + add(numpaths("start", pp), pp in small minus {"start"});


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

`is_element/WFbar` := (N::posint) -> (A::set) -> proc(TTxl)
 local TT,x,l,TT1,T,n,u,children;
 global reason;

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

 TT,x,l := op(TTxl);

 if not(`is_element/trees`(A)(TT)) then
  reason := [convert(procname,string),"TT is not a tree on A",TT,A];
  return false;
 fi;

 TT1 := select(U -> nops(U) > 1,TT);

 if not(type(x,table)) then
  reason := [convert(procname,string),"x is not a table",x];
  return false;
 fi;

 if map(op,{indices(x)}) <> TT1 then
  reason := [convert(procname,string),"x is not indexed by TT1",x,TT1];
  return false;
 fi;

 children := `children_map`(A)(TT);
 for T in TT1 do
  if not(`is_element/Fbar`(N)(children[T])(x[T])) then
   reason := [convert(procname,string),"x[T] has invalid type",T,x[T],reason];
   return false;
  fi;
 od;

 if not(type(l,table)) then
  reason := [convert(procname,string),"l is not a table",l];
  return false;
 fi;

 if map(op,{indices(l)}) <> TT then
  reason := [convert(procname,string),"l is not indexed by TT",l,TT];
  return false;
 fi;

 n := nops(A);

 for T in TT do
  if not(`is_element/RR`(l[T]) and l[T] >= 0 and l[T] <= 1) then
   reason := [convert(procname,string),"l[T] is not in the unit interval",T,l[T]];
   return false;
  fi;

  if nops(T) = 1 and l[T] <> 1 then
   reason := [convert(procname,string),"l[{a}] <> 1",op(T),l[T]];
   return false;
  fi;

  if nops(T) = n and l[T] <> 1 then
   reason := [convert(procname,string),"l[A] <> 1",T,l[T]];
   return false;
  fi;
 od;

 return true;
end;

`is_interior/WFbar` := (N::posint) -> (A::set) -> proc(TTxl)
 local TT,x,l,n,T;
 global reason;

 TT,x,l := op(TTxl);
 n := nops(A);

 for T in TT do
  if nops(T) > 1 and nops(T) < n and l[T] = 1 then
   reason := [convert(procname,string),"not interior, l[T] = 1",T];
   return false;
  fi;
 od:

 return true;
end;

`is_reduced/WFbar` := (N::posint) -> (A::set) -> proc(TTxl)
 local TT,x,l,T;
 global reason;

 TT,x,l := op(TTxl);

 for T in TT do
  if l[T] = 0 then
   reason := [convert(procname,string),"not reduced, l[T] = 0",T];
   return false;
  fi;
 od:

 return true;
end;

`reduce_once/WFbar` := (N::posint) -> (A::set) -> proc(TTxl,T)
 local TT,x,l,parent,children,n,TT1,x1,l1,U,B,C,X,p,y,z;

 TT,x,l := op(TTxl);

 parent := `parent_map`(A)(TT);
 children := `children_map`(A)(TT);

 TT1 := TT minus {T};

 U := parent[T];
 C := children[U];
 B := C minus {T} union children[T];

 p := table();
 for X in C minus {T} do p[X] := X; od;
 for X in children[T] do p[X] := T; od;

 y := table();
 for X in C minus {T} do y[X] := `eta/Fbar`(N)({X}); od;
 y[T] := x[T];

 z := `gamma/Fbar`(N)(B,C)(p)(x[U],y);

 x1 := table();
 l1 := table();

 for X in TT1 do
  l1[X] := l[X];
  if nops(X) > 1 then
   x1[X] := `if`(X = U,eval(z),eval(x[X]));
  fi;
 od;

 return([TT1,eval(x1),eval(l1)]);
end;

`reduce/WFbar` := (N::posint) -> (A::set) -> proc(TTxl)
 local is_reduced,X,TT,x,l,T;

 is_reduced := false;
 X := TTxl;

 while not(is_reduced) do
  TT,x,l := op(X);
  is_reduced := true;

  for T in TT do
   if l[T] = 0 then
    is_reduced := false;
    X := `reduce_once/WFbar`(N)(A)(X,T);
    break;
   fi;
  od:
 od;

 return X;
end;

`is_equal/WFbar` := (N::posint) -> (A::set) -> proc(X0,X1)
 local Y0,Y1,TT0,TT1,x0,x1,l0,l1,children,T;
 global reason;

 Y0 := `reduce/WFbar`(N)(A)(X0);
 Y1 := `reduce/WFbar`(N)(A)(X1);
 TT0,x0,l0 := op(Y0);
 TT1,x1,l1 := op(Y1);

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

 children := `children_map`(A)(TT0);

 for T in TT0 do
  if l0[T] <> l1[T] then
   reason := [convert(procname,string),"l0[T] <> l1[T]",T,l0[T],l1[T]];
   return false;
  fi;

  if nops(T) > 1 then
   if not(`is_equal/Fbar`(N)(children[T])(x0[T],x1[T])) then
    reason := [convert(procname,string),"x0[T] <> x1[T]",T,x0[T],x1[T],reason];
    return false;
   fi;
  fi;
 od;

 return true;
end;

`is_leq/WFbar` := NULL;

`random_element/WFbar` := (N::posint) -> (A::set) -> proc()
 local TT,x,l,TT1,T,n,children,d;

 d := 3;
 TT := `random_element/full_trees`(A)();

 TT1 := select(U -> nops(U) > 1,TT);

 x := table();
 l := table();

 children := `children_map`(A)(TT);
 for T in TT1 do
  x[T] := `random_element/Fbar`(N)(children[T])();
 od;

 n := nops(A);

 for T in TT do
  if nops(T) = 1 or nops(T) = n then
   l[T] := 1;
  else
   l[T] := rand(0..d)()/d;
  fi;
 od;

 return [TT,eval(x),eval(l)];
end;

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


text:
(*
 Copyright (C) 2020 Susi Lehtola

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

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

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

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

$define mgga_x_br89_params
$include "mgga_x_br89.mpl"

(* This is a fake parameter in libxc *)
params_a_at := 0:

(* Equation 9, same-spin correlation *)
b94_css := (rs, z, xs, us, ts) ->
  - 0.01 * (1 + z)^(8/3) * 2^(-8/3) * n_total(rs)^(5/3) * (2*ts - xs^2/4)
  * b94_zss(params_a_css, br89_f, rs, z, xs, us, ts)^4 * (
    1 - 2*log(1 + b94_zss(params_a_css, br89_f, rs, z, xs, us, ts)/2)
      / b94_zss(params_a_css, br89_f, rs, z, xs, us, ts)
  ):

(* Same-spin correlation overall *)
b94_par := (rs, z, xs0, xs1, us0, us1, ts0, ts1) ->
  + my_piecewise3(screen_dens(rs,  z), 0, b94_css(rs, z_thr( z), xs0, us0, ts0))
  + my_piecewise3(screen_dens(rs, -z), 0, b94_css(rs, z_thr(-z), xs1, us1, ts1)):

(* Equation 8, opposite-spin correlation *)
b94_cab := (rs, z, xs0, xs1, us0, us1, ts0, ts1) ->
  - 0.8 * (1 - z^2)/4 * n_total(rs)
  * b94_zab(params_a_cab, br89_f, rs, z, xs0, xs1, us0, us1, ts0, ts1) * (
     b94_zab(params_a_cab, br89_f, rs, z, xs0, xs1, us0, us1, ts0, ts1)
     - log(1 + b94_zab(params_a_cab, br89_f, rs, z, xs0, xs1, us0, us1, ts0, ts1))
  ):

(* Whole functional *)
b94_c_f := (rs, z, xs0, xs1, us0, us1, ts0, ts1) ->
  + b94_cab(rs,  z, xs0, xs1, us0, us1, ts0, ts1)
  + b94_par(rs,  z, xs0, xs1, us0, us1, ts0, ts1):

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


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

htbs_s1 := 0.6:
htbs_s2 := 2.6:

(* The equations to solve in order to obtain the coeficients cc are
  G(s1) = 0
  G(s2) = 1
 G'(s1) = 0
 G'(s2) = 0
G''(s1) = 0
G''(s2) = 0
*)

htbs_cc0 :=  htbs_s1^3*(htbs_s1^2 - 5*htbs_s1*htbs_s2 + 10*htbs_s2^2)/(htbs_s1 - htbs_s2)^5:
htbs_cc1 := -30*htbs_s1^2*htbs_s2^2/(htbs_s1 - htbs_s2)^5:
htbs_cc2 :=  30*htbs_s1*htbs_s2*(htbs_s1 + htbs_s2)/(htbs_s1 - htbs_s2)^5:
htbs_cc3 := -10*(htbs_s1^2 + 4*htbs_s1*htbs_s2 + htbs_s2^2)/(htbs_s1 - htbs_s2)^5:
htbs_cc4 :=  15*(htbs_s1 + htbs_s2)/(htbs_s1 - htbs_s2)^5:
htbs_cc5 := -6/(htbs_s1 - htbs_s2)^5:

$define gga_x_rpbe_params
$include "gga_x_rpbe.mpl"
$include "gga_x_wc.mpl"

htbs_g := s -> htbs_cc0 + htbs_cc1*s + htbs_cc2*s^2 + htbs_cc3*s^3 + htbs_cc4*s^4 + htbs_cc5*s^5:

htbs_f0 := s -> my_piecewise3(
   s < htbs_s1, wc_f0(s),
   my_piecewise3(s > htbs_s2, rpbe_f0(s), htbs_g(s)*rpbe_f0(s) + (1 - htbs_g(s))*wc_f0(s))
):

htbs_f  := x -> htbs_f0(X2S*x):

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


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

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

 global reason;

 if not `is_element/ACP`(N)(A)(Q) then
  reason := [convert(procname,string),"Q in ACP(N)(A)",reason];
  return false;
 fi;

 if not `is_separated/preord`(A)(Q[N]) then
  reason := [convert(procname,string),"Q[N] is not separated",Q[N]];
  return false;
 fi;

 return true;
end:

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

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

 return true;
end:

`is_leq/ICP` := (N::posint) -> (A::set) -> proc(Q1,Q2)
 local i;

 for i from 1 to N do
  if Q2[i] minus Q1[i] <> {} then
   return false;
  fi;
 od;

 return true;
end:

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

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

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

 R := `random_element/total_preord`(A)();
 Q := [R];
 for i from 2 to N do
  pi := `block_partition/preord`(A)(R);
  R := {};
  for B in pi do
   if i < N then
    S := `random_element/total_preord`(B)();
   else
    C := `random_element/ord`(B)();
    n := nops(C);
    S := {seq(seq([C[i],C[j]],j=i..n),i=1..n)};
   fi;

   R := R union S;
  od;
  Q := [op(Q),R];
 od;

 return Q;
end:

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

`build/ICP` := (N::posint) -> (A::set) -> proc(Ru)
 local R,u,n,i,j,p,a,b,aa,bb,Q;
 R,u := op(Ru);
 n := nops(A);
 aa := table():
 bb := table():
 for i from 1 to n do
  for p from 1 to N do
   a := i;
   while a > 1 and u[a-1] > p do a := a - 1; od;
   aa[i,p] := a;
   b := i;
   while b < n and u[b] >= p do b := b + 1; od;
   bb[i,p] := b;
  od:
 od:
 Q := [seq(
  {seq(seq([R[i],R[j]],j=aa[i,p]..bb[i,p]),i=1..n)}
 ,p=1..N)];

 return Q;
end:

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

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

 U := [[]];
 n := nops(A);
 for i from 1 to n-1 do
  U := [seq(seq([op(u),j],j=1..N),u in U)];
 od:

 RR := `list_elements/ord`(A);

 [seq(seq(`build/ICP`(N)(A)([R,u]),u in U),R in RR)];
end:

`count_elements/ICP` := (N::posint) -> (A::set) ->
 nops(A)! * N^(nops(A) - 1);

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

`list_ordered_elements/ICP` := (N::posint) -> proc(A::{set,list})
 local U,u,n,i,j,R,A0;

 U := [[]];
 n := nops(A);
 for i from 1 to n-1 do
  U := [seq(seq([op(u),j],j=1..N),u in U)];
 od:

 R := [op(A)];
 A0 := {op(A)};

 [seq(`build/ICP`(N)(A0)([R,u]),u in U)];
end:

`count_ordered_elements/ICP` := (N::posint) -> (A::set) ->
 N^(nops(A) - 1);

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

# Note that we omit the rank of Q[N], because it is always equal to N here.

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

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

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

`phi/SEM/ICP` := (N::posint) -> (A::set) -> proc(eta)
 local A2,a,b,i;

 A2 := {seq(seq([a,b],b in A),a in A)};

 return [seq(select(ab -> `is_preceq/E`(N)(eta[op(ab)],epsilon^i),A2),i=0..N-1)];
end;

`psi/ICP/SEM` := (N::posint) -> (A::set) -> proc(Q)
 local eta,a,b,i;

 eta := table();
 for a in A do
  for b in A do
   if a = b then
    eta[a,b] := 0;
   else
    i := 1;
    while (i < N and member([a,b],Q[i]) and member([b,a],Q[i])) do
     i := i+1;
    od;
    if member([a,b],Q[i]) then
     eta[a,b] :=  epsilon^(i-1);
    else
     eta[a,b] := -epsilon^(i-1);
    fi;
   fi;
  od;
 od;

 return eval(eta);
end:

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

`totalise/ICP/ord` := (N::posint) -> (A::set) -> proc(Q)
 local R,T,r,i,a;

 R := `id/autorel`(A);
 for i from 1 to N do
  R := R union (Q[i] minus `op/autorel`(A)(Q[i]));
 od:

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

 return [seq(T[i],i=1..nops(A))];
end:

`res/ICP` := (N::posint) -> (A::set,B::set) -> proc(Q)
 local BB;

 BB := `top/autorel`(B);

 return(map(`intersect`,Q,BB));
end;

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

`is_fibre_sphere/ICP` := (N::posint) -> (A::set) -> (a,b,Q) -> proc(P)
 local B,i,U,V,L;

 if not(member(a,A)) then error("a not in A"); fi;
 B := A minus {a};

 if not(`is_equal/ICP`(N)(B)(`res/ICP`(N)(A,B)(P),Q)) then
  return false;
 fi;

 i := 1;
 while i < N and member([a,b],P[i]) and member([b,a],P[i]) do
  i := i+1;
 od;

 U := select(x -> member([x,a],P[i]),B);
 V := select(x -> member([a,x],P[i]),B);

 if U intersect V <> {} then return false; fi;
 if member(b,U) then
  L := select(x -> not(member([x,b],P[i])),U);
  if L <> {} then return false; fi;
 fi;
 if member(b,V) then
  L := select(x -> not(member([b,x],P[i])),V);
  if L <> {} then return false; fi;
 fi;

 return true;
end;

######################################################################
# Each of the subsets S_b(Q) is isomorphic to ICP_N({a,b})

`f/fibre_sphere/ICP` := (N::posint) -> (A::set) -> (a,b,Q) -> proc(R)
 local i,e;

 i := 1;
 while member([a,b],R[i]) and member([b,a],R[i]) do
  i := i+1;
 od;
 if member([a,b],R[i]) then
  e := -1;
 else
  e := +1;
 fi;

 return `f_alt/fibre_sphere/ICP`(N)(A)(a,b,Q)([i,e]);
end:

`f_alt/fibre_sphere/ICP` := (N::posint) -> (A::set) -> (a,b,Q) -> proc(ie)
 local i,e,PT,j,x,L,U;

 i,e := op(ie);

 PT := table();
 for j from 1 to i-1 do
  PT[j] := {[a,a]} union Q[j];
  for x in Q[j] do
   if x[1] = b then PT[j] := {op(PT[j]),[a,x[2]]}; fi;
   if x[2] = b then PT[j] := {op(PT[j]),[x[1],a]}; fi;
  od;
 od;
 PT[i] := Q[i] union {[a,a]};
 L := select(x -> member([x,b],Q[i]),A);
 U := select(x -> member([b,x],Q[i]),A);
 if e = -1 then
  L := L minus U;
 else
  U := U minus L;
 fi;
 PT[i] := PT[i] union {seq([x,a],x in L)}
                union {seq([a,x],x in U)};
 for j from i+1 to N do
  PT[j] := Q[j] union {[a,a]};
 od;

 return [seq(PT[j],j=1..N)];
end:

`g/fibre_sphere/ICP` := (N::posint) -> (A::set) -> (a,b,Q) -> proc(P)
 if not(`is_fibre_sphere/ICP`(N)(A)(a,b,Q)(P)) then
  return FAIL;
 fi;
 return `res/ICP`(N)(A,{a,b})(P);
end:

`g_alt/fibre_sphere/ICP` := (N::posint) -> (A::set) -> (a,b,Q) -> proc(P)
 local R,i,e;

 R := `g/fibre_sphere/ICP`(N)(A)(a,b,Q)(P);
 if R = FAIL then return FAIL; fi;
 i := 1;
 while member([a,b],R[i]) and member([b,a],R[i]) do
  i := i+1;
 od;
 if member([a,b],R[i]) then
  e := -1;
 else
  e := +1;
 fi;
 return [i,e];
end:

######################################################################
# This function provides a point of intesection between
# S_b(Q) and S_b1(Q), where b1 is the successor of b in the total
# order induced by Q.

`m/fibre_sphere/ICP` := (N::posint) -> (A::set) -> proc(a,b,Q)
 local B,L,n,r,b1,i,R;

 B := A minus {a};
 L := `totalise/ICP/ord`(N)(B)(Q);
 n := nops(B);
 r := table();
 for i from 1 to n do
  r[L[i]] := i;
 od;
 if not(r[b] < n) then
  return FAIL;
 fi;
 b1 := L[r[b]+1];
 i := 0;
 while i < N and member([b1,b],Q[i+1]) do
  i := i+1;
 od;
 R := [{[a,a],[a,b],[b,a],[b,b]}$i,
       {[a,a],[b,a],[b,b]},
       {[a,a],[b,b]}$(N-i-1)];
 return `f/fibre_sphere/ICP`(N)(A)(a,b,Q)(R);
end:

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

`bump/ICP` := (N::posint) -> (A::set) -> proc(Q)
 local r,n,i,j,R,E,Q1;

 r := `rank_vector/ICP`(N)(A)(Q);
 n := nops(A);
 i := N-1;
 while i > 0 and r[i] = n-1 do
  i := i - 1;
 od;
 if i = 0 then
  return FAIL;
 fi;
 R := `bump/striped_preord`(A)(Q[i]);
 E := R intersect `op/autorel`(A)(R);
 Q1 := [seq(Q[j],j=1..i-1),
        R,
        seq(Q[j] intersect E,j=i+1..N)];
 return Q1;
end:

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

`describe/ICP` := eval(`describe/ACP`):


text:
include "MMH_megs.mpl"
include "MMH_portalSpawner.mpl"
include "MMH_lightBridgeSpawner.mpl"
include "MMH_laserSpawner.mpl"
include "MMH_weightedStorageCubeSpawner.mpl"
include "MMH_discouragementRedirectionCubeSpawner.mpl"
include "MMH_repulsionGelSpawner.mpl"
include "MMH_propulsionGelSpawner.mpl"
include "MMH_cleansingGelSpawner.mpl"

project MMH {
  orientation "zy-x"
}

///setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"Adrodoc55\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"Adrodoc55's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" YouTube-Channel by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"https://www.youtube.com/user/Adrodoc55\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"because I wanted",Text4:"to thank myself"}
///setblock ~ ~ ~ standing_sign 0 replace {Text1:"{\"text\":\"Special Thanks\",\"bold\":\"true\"}",Text2:"{\"text\":\"click on the \",\"extra\":[{\"text\":\"red\",\"color\":\"red\"}]}",Text3:"names for more",Text4:"information"}
///setblock ~ ~ ~ wall_sign 3 replace {Text1:"{\"text\":\"SQRTdude\",\"color\":\"red\",\"clickEvent\":{\"action\":\"run_command\",\"value\":\"/tellraw @a {\\\"text\\\":\\\"\\\",\\\"extra\\\":[{\\\"text\\\":\\\"Check out \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"SQRTdude's\\\",\\\"color\\\":\\\"gold\\\"},{\\\"text\\\":\\\" Video by clicking \\\",\\\"color\\\":\\\"yellow\\\"},{\\\"text\\\":\\\"here\\\",\\\"color\\\":\\\"aqua\\\",\\\"clickEvent\\\":{\\\"action\\\":\\\"open_url\\\",\\\"value\\\":\\\"http://youtu.be/g7xfCBRjGRQ\\\"}}]}\"}}",Text2:"{\"text\":\"\"}",Text3:"for the",Text4:"Copy Helper"}

install {
  /scoreboard objectives add MMH_Rotation dummy
  start @e[tag=MMH]
  /tellraw @a [{"text":"Welcome to","color":"dark_green"},{"text":" Mapmaker Helper\n","color":"green"},{"text":"   Version: ","color":"aqua"},{"text":"2.1.0\n","color":"dark_aqua","italic":"true"},{"text":"   made by ","color":"aqua","italic":"true"},{"text":"Adrodoc55\n","color":"dark_aqua","hoverEvent":{"action":"show_text","value":{"text":"Click to open YouTube Channel","color":"red"}},"clickEvent":{"action":"open_url","value":"https://www.youtube.com/user/Adrodoc55"}}]
}

uninstall {
  /scoreboard objectives remove MMH_Rotation
  stop @e[tag=MMH]
}

#MMH
repeat process MMH_rotation {
  /scoreboard players set @a[rym=45,ry=135] MMH_Rotation 0
  /scoreboard players set @a[rym=135,ry=-135] MMH_Rotation 1
  /scoreboard players set @a[rym=-135,ry=-45] MMH_Rotation 2
  /scoreboard players set @a[rym=-45,ry=45] MMH_Rotation 3
}


text:
`bideg0/Lambda` := proc(u)
 if type(u,specfunc(nonnegint,lambda)) then
  return [nops(u),`+`(op(u)) + nops(u)];
 else
  return [0,0];
 fi;
end:

`bideg/Lambda` := apply_deg(`bideg0/Lambda`);

`mu0/Lambda` := (u,v) -> `reduce/Lambda`(op(u),op(v));

`mu/Lambda` := apply_linear_assoc_mod(`mu0/Lambda`,lambda(),2);

`is_admissible/Lambda` := proc()
 if nargs = 1 and type(args[1],specfunc(nonnegint,lambda)) then
  return `is_admissible/Lambda`(op(args[1]));
 elif type([args],list(nonnegint)) then
  return evalb(min(op((2 *~ [infinity,args]) -~ [args,0])) >= 0);
 else
  return FAIL;
 fi;
end:

`reduce/Lambda` := proc()
 local K,t,left,right,i,j,n,u,c;

 option remember;

 if `is_admissible/Lambda`(args) then
  return lambda(args);
 fi;

 K := [args];
 t := 1;
 while t < nargs and K[t+1] <= 2*K[t] do t := t+1; od;
 left := lambda(seq(K[i],i=1..t-1));
 right := `reduce/Lambda`(seq(K[i],i=t+2..nops(K)));
 i := K[t];
 n := K[t+1] - 2*i - 1;
 u := 0;
 for j from 0 to floor((n-1)/2) do
  c := modp(binomial(n-1-j,j),2);
  if c <> 0 then
   u := modp(u + c * `mu/Lambda`(left,`reduce/Lambda`(i+n-j,2*i+1+j),right),2);
  fi;
 od;
 return u;
end:

`d0/Lambda` := proc(n::nonnegint)
 local u,j,c;
 option remember;

 u := 0;
 for j from 1 to floor(n/2) do
  c := modp(binomial(n-j,j),2);
  if c <> 0 then
   u := u + `reduce/Lambda`(n-j,j-1);
  fi;
 od;
 u := modp(u,2);
 return u;
end:

`d1/Lambda` := proc(u)
 local K,n,v,i,j,x,y,z;
 option remember;

 if not(type(u,specfunc(nonnegint,lambda))) then
  return FAIL;
 fi;

 K := [op(u)];
 n := nops(K);
 v := 0;
 for i from 1 to n do
  x := lambda(seq(K[j],j=1..i-1));
  y := `d0/Lambda`(K[i]);
  z := lambda(seq(K[j],j=i+1..n));
  v := modp(v + `mu/Lambda`(x,y,z),2);
 od;
 return v;
end:

`d/Lambda` := apply_linear_mod(`d1/Lambda`,2);

`basis/Lambda` := proc(s::integer,t::integer)
 local B,C,r;
 option remember;

 if s < 0 then
  return [];
 elif s = 0 then
  if t = 0 then
   return [lambda()];
  else
   return [];
  fi;
 elif s = 1 then
  if t >= 1 then
   return [lambda(t-1)];
  else
   return [];
  fi;
 else
  B := NULL;
  for r from 0 to t-1 do
   C := `basis/Lambda`(s-1,r);
   C := select(v -> (op(-1,v) >= t-r-1),C);
   C := map(v -> lambda(op(v),t-r-1),C);
   B := B,op(C);
  od;
  return [B];
 fi;
end:

`vec/Lambda` := (s,t) -> proc(u)
 local m;

 [seq(modp(coeff(u,m,1),2),m in `basis/Lambda`(s,t))];
end:

`d_matrix/Lambda` := proc(s,t)
 local M;
 option remember;

 M := map(m -> `vec/Lambda`(s+1,t)(`d/Lambda`(m)),`basis/Lambda`(s,t));
 M := LinearAlgebra[Modular][Mod](2,Transpose(Matrix(M)),integer[]);
 return M;
end:


text:
`is_element/vector_functions` := (N::posint) -> (A::set) -> proc(x)
 local a;
 global reason;

 if not(is_table_on(A)(x)) then
  reason := [convert(procname,string),"x is not a table on A",x,A];
  return false;
 fi;

 for a in A do
  if not `is_element/R`(N)(x[a]) then
   reason := [convert(procname,string),"x[a] is not in R^N",a,x[a],N];
   return false;
  fi;
 od;

 return true;
end;

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

`is_equal/vector_functions` := (N::posint) -> (A::set) -> proc(x,y)
 local a;
 global reason;

 for a in A do
  if not(`is_equal/R`(N)(x[a],y[a])) then
   reason := [convert(procname,string),"x[a] <> y[a]",a,x,y];
   return false;
  fi;
 od;

 return true;
end;

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

`is_nonnegative/vector_functions` := (N::posint) -> (A::set) -> proc(x)
 local a,i;

 for a in A do
  for i from 1 to N do
   if x[a][i] < 0 then
    return false;
   fi;
  od;
 od;

 return true;
end:

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

`is_zero/vector_functions` := (N::posint) -> (A::set) -> proc(x)
 local a,i;

 for a in A do
  for i from 1 to N do
   if x[a][i] <> 0 then
    return false;
   fi;
  od;
 od;

 return true;
end:

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

`is_leq/vector_functions` := (N::posint) -> (A::set) -> proc(x,y)
 local a,i;

 for a in A do
  for i from 1 to N do
   if x[a][i] > y[a][i] then
    return false;
   fi;
  od;
 od;

 return true;
end:

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

`plus/vector_functions` := (N) -> (A::set) -> proc(x,y)
 local z,a;

 z := table();
 for a in A do
  z[a] := x[a] +~ y[a];
 od:

 return eval(z):
end:

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

`times/vector_functions` := (N) -> (A::set) -> proc(t,x)
 local z,a;

 z := table();
 for a in A do
  z[a] := t *~ x[a];
 od:

 return eval(z):
end:

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

`norm/vector_functions` := (N::posint) -> (A::set) -> proc(x)
 local a,n;

 n := 0;
 for a in A do n := n + `norm_2/R`(N)(x[a])^2; od;

 return sqrt(n);
end:

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

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

 u := [0$N];
 for a in A do u := u +~ x[a]; od;

 return u;
end:

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

`average/vector_functions` := (N::posint) -> (A::set) -> proc(x)
 if A = {} then
  return FAIL;
 fi;

 return `sum/vector_functions`(N)(A)(x) /~ nops(A);
end:

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

`dist/vector_functions` := (N::posint) -> (A::set) -> proc(x,y)
 local a,n;

 n := 0;
 for a in A do n := n + `d_2/R`(N)(x[a],y[a])^2; od;

 return sqrt(n);
end:

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

`dot/vector_functions` := (N::posint) -> (A::set) -> proc(x,y)
 local a,d;

 d := 0;
 for a in A do d := d + `dot/R`(N)(x[a],y[a]); od;

 return d;
end:

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

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

 x := table();
 for a in A do
  x[a] := `random_element/R`(N)();
 od:

 return eval(x);
end;

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

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


text:
# Kinematik-Berechnung Modell PalH2m1 (Hybrider Palettierer, Variante 2, Modellierung 1)
# Beschreibung
# Ein Palettierroboter mit fünf Achsen als Näherung für einen Palettierroboter mit Viergelenkketten
# Die vierte Achse des Ersatzmodells wird als abhängiges Gelenk in Abhängigkeit der Achsen 2 und 3 beschrieben.
# Dadurch wird der Effekt der Parallelogramme kinematisch nachgebildet.
# Die Gelenkmomente sollten genauso wie durch die Parallelogramme übertragen werden (als ob deren Dynamikparameter Null sind).
# Quellen
# Autor
# Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-12
# (C) Institut für Mechatronische Systeme, Universität Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 2: # Hoher Optimierungsgrad.
;
read "../helper/proc_MatlabExport":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Variable mit Winkeln der Zwangsbedingungen nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
kintmp_qt := Matrix(RowDimension(kintmp_s),1):
# Ersetzungsausdrücke definieren.
# Variable zum speichern des Sinus und Cosinus der Winkel. Für dieses System ist das eigentlich nicht notwendig. Erstelle Variable, da sie von den anderen Skripten erwartet wird
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
# Zwangsbedingung für abhängigen Gelenkwinkel
# Der vierte Winkel ist der passive Winkel nahe des Endeffektors. Er wird durch die Parallelstruktur (aus zwei Parallelogrammen) bei hybriden Palettierrobotern bestimmt.
# Durch die Parallelogramme ergibt sich ein konstantes Verhältnis zu den Gelenkwinkeln der Achsen 2 und 3, die auch in derselben Ebene liegen.
rho4_qt := -theta(2) -theta(3):
kintmp_qt(1,1) := rho4_qt:
kintmp_qs := convert_t_s(kintmp_qt):
# Exportiere Code für folgende Skripte
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
printf("Ausdrücke für kintmp_subsexp gespeichert (Maple)\n"):
save kintmp_qs, kintmp_qt, kin_constraints_exist, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kintmp_qs, kintmp_qt, kin_constraints_exist, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_qs )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m", robot_name), 2):
printf("Fertig\n"):

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


text:
TITLE
   KandW3Scen;

STOCHASTIC

STAGES
   t := 1..3;

SCENARIO
   sc := 1..9;

TREE
   TSTAGE(1,3,3,2,3,3,2,3,3);
!   NWAY(3);

INDEX
   rmat := 1..2; !material
   prod := 1..2; !prod

PROBABILITIES
   Prob[sc] := (0.06,0.15,0.09,0.12,0.16,0.12,0.12,0.12,0.06);

RANDOM DATA
   Demand[t,prod,sc] := SPARSEFILE("KandW3.dat");

DATA
   RawMatCost[rmat]       := (2.0, 3.0);

   RawMatUsage[prod,rmat] := (2.0, 6.0,
                              3.0, 3.4);

   OutSourceCost[t>=2,prod] := ( 7.0, 10.0,
                                12.0, 15.0);

   InvtCap := 50;

VARIABLE
   RawMatPurch[rmat, t<=2] -> x;
   OutSourced[prod, t>=2] -> y;

MODEL
   MIN TCosts = SUM(rmat, t<=2:  RawMatCost * RawMatPurch)
              + SUM(prod, t>=2:  OutSourceCost * OutSourced);

SUBJECT TO
   Dummy[t=1]:
       SUM(rmat: RawMatPurch) >=0;

   MaxPurchase -> INV:

       SUM(rmat, t<=2: RawMatPurch) <= InvtCap;

   MeetDemand[prod, t>=2] -> DMD:

       SUM(rmat: RawMatUsage * RawMatPurch[t-1]) + OutSourced >= Demand;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


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

`eta/simplex_interior` :=
 (A::set) -> `if`(nops(A)=1,table({op(A)=1}),FAIL);

`gamma/simplex_interior` := (A::set,B::set) -> (p) -> proc(y,x)
 local a;
 return table({seq(a = x[p[a]][a] * y[p[a]],a in A)});
end:


text:
euler_angle_matrix := (t) ->
<<cos(t[2])|
  -cos(t[3])*sin(t[2])|
  sin(t[2])*sin(t[3])>,
 <cos(t[1])*sin(t[2])|
  cos(t[1])*cos(t[2])*cos(t[3])-sin(t[1])*sin(t[3])|
  -cos(t[3])*sin(t[1])-cos(t[1])*cos(t[2])*sin(t[3])>,
 <sin(t[1])*sin(t[2])|
  cos(t[1])*sin(t[3])+cos(t[2])*cos(t[3])*sin(t[1])|
  cos(t[1])*cos(t[3])-cos(t[2])*sin(t[1])*sin(t[3])>
>;

mobius_embedding := (R,r) -> (t,u) ->
 (R + u*r*cos(2*Pi*t)) *~ [cos(4*Pi*t),sin(4*Pi*t),0] +~ [0,0,u *r* sin(2*Pi*t)];

torus_embedding := (R,r) -> (t,u) ->
 (R + r * cos(2*Pi*u)) *~ [cos(2*Pi*t),sin(2*Pi*t),0] +~ [0,0,(r * sin(2*Pi*u))];

simplex_embedding := proc(t)
 if nops(t) = 1 then
  return 0;
 elif nops(t) = 2 then
  return 2*t[1]-1;
 elif nops(t) = 3 then
  return [t[1]-t[2]/2-t[3]/2,sqrt(3)/2*(t[2]-t[3])];
 elif nops(t) = 4 then
  return [sqrt(2)*(2*t[2]-t[3]-t[4])/3,
          sqrt(2/3)*(t[3]-t[4]),
	  t[1]-(t[2]+t[3]+t[4])/3];
 else
  error("Simplex embedding only defined in dim <= 3");
 fi;
end:

mobius_arc := (R,r) -> proc(t1,u1,t2,u2)
 local s1,s2;
 if u1 = 0 then
  s2 := t2;
  s1 := t1 - round(2*(t1 - t2))/2;
 elif u2 = 0 then
  s1 := t1;
  s2 := t2 - round(2*(t2 - t1))/2;
 else
  s1 := t1;
  s2 := t2 - round(t2-t1);
 fi;
 return spacecurve(mobius_embedding(R,r)(s1 + p*(s2-s1),u1 + p*(u2-u1)),p=0..1,args[5..-1]);
end:

sphere_arc := proc(a,b)
 local r,c,d,theta,i;
 r := evalf(sqrt(add(a[i]^2,i=1..3)));
 c := b -~ a;
 d := b +~ a;
 c := evalf(c *~ (r/sqrt(add(c[i]^2,i=1..3))));
 d := evalf(d *~ (r/sqrt(add(d[i]^2,i=1..3))));
 theta := evalf(arccos(add(a[i] * b[i],i=1..3)/r^2)/2);
 return spacecurve(d *~ cos(t) +~ c *~ sin(t),t = -theta..theta,args[3..-1]);
end:

spherical_angle := proc(a,b,c)
 local dp,ab,ac;
 dp := `dot/R`(3);
 ab := b -~ dp(b,a) *~ a;
 ac := c -~ dp(c,a) *~ a;
 return arccos(dp(ab,ac) /~ sqrt(dp(ab,ab) * dp(ac,ac)));
end:

spherical_triangle_area := proc(a,b,c)
 return spherical_angle(a,b,c) +
        spherical_angle(b,c,a) +
        spherical_angle(c,a,b) - Pi;
end:

spherical_triangle_area_alt := proc(a,b,c)
 local dp,f,x,y,z,F;
 dp := `dot/R`(3);
 f := (c) -> sqrt(1-c)/(sqrt(2) + sqrt(1+c));
 x,y,z := f(dp(b,c)),f(dp(c,a)),f(dp(a,b));
 F := (p,q,r) -> (p + q + r - p*q*r)/(1 - p*q - q*r - r*p);
 return 4 * arctan(sqrt(F( x, y, z) * F(-x, y, z) * F( x,-y, z) * F( x, y,-z)));
end:

simplex_outline := proc(d)
 local e,i,j,opts;
 if d <> 2 and d <> 3 then
  error("Simplex outline is only defined in dimensions 2 and 3");
 fi;

 e := table():
 for i from 0 to d do
  e[i] := simplex_embedding([0$i,1,0$(d-i)]);
 od:

 opts := [];
 if nargs > 1 then opts := [args[2..-1]]; fi;
 if d = 3 and opts = [] then opts := [colour=black]; fi;

 return display(seq(seq(line(e[i],e[j],op(opts)),j=i+1..d),i=0..d),
                scaling=constrained,axes=none);
end:

trefoil_embedding := (R,r) -> proc(t,u)
 local x,y,z,i;
 x := [sin(t) + 2*sin(2*t),cos(t) - 2 * cos(2*t),-sin(3*t)];
 y := [72*sin(2*t)+3*sin(8*t)-13*sin(4*t)+3*sin(7*t)-14*sin(5*t)+3*sin(t),
       3*cos(t)-3*cos(8*t)+3*cos(7*t)-72*cos(2*t)+14*cos(5*t)-13*cos(4*t),
       10*sin(6*t)-34*sin(3*t)];
 z := [-391*cos(t)+2*cos(8*t)-29*cos(7*t)+85*cos(2*t)-99*cos(5*t)+24*cos(4*t)+9*cos(10*t)-9*cos(11*t),
       -9*sin(11*t)+29*sin(7*t)+85*sin(2*t)-9*sin(10*t)+2*sin(8*t)-24*sin(4*t)-99*sin(5*t)+391*sin(t),
       -570-34*cos(3*t)-94*cos(6*t)+18*cos(9*t)];
 y := y /~ sqrt(add(y[i]^2,i=1..3));
 z := z /~ sqrt(add(z[i]^2,i=1..3));
 return R *~ x +~ (r * cos(u)) *~ y +~ (r * sin(u)) *~ z;
end:

trefoil_b := (t) -> [(2+cos(3*t))*cos(2*t),(2+cos(3*t))*sin(2*t),sin(3*t)];

klein_embedding := (t,u) ->
[(0.1*sin(3*Pi*t)+0.1*sin(Pi*t)+0.4*cos(Pi*t))*sin(2*Pi*u)-0.5*sin(4*Pi*t)+sin(2*Pi*t),
  0.2*sin(2*Pi*u)*sin(3*Pi*t)+2.*cos(2*Pi*t)+0.5,
  0.25*cos(2*Pi*u)*sin(2*Pi*t)+0.4*cos(2*Pi*u)];

klein_seam_a := (t) -> [0.8621-0.0015*cos(4*Pi*t)+0.0171*sin(2*Pi*t)-0.0002*sin(6*Pi*t),t];
klein_seam_b := (t) -> [0.1179+0.0031*cos(4*Pi*t)-0.0385*sin(2*Pi*t)+0.0005*sin(6*Pi*t),
                        -0.25+0.0597*cos(2*Pi*t)-0.0004*cos(6*Pi*t)+0.0049*sin(4*Pi*t)];
klein_seam := (t) -> [-0.2811+0.0051*cos(4*Pi*t)-0.1401*sin(2*Pi*t)+0.0005*sin(6*Pi*t),
                      1.7871-0.0061*cos(4*Pi*t)+0.3532*sin(2*Pi*t)-0.0005*sin(6*Pi*t),
                      0.2090*cos(2*Pi*t)-0.0010*cos(6*Pi*t)+0.0086*sin(4*Pi*t)];

orient_face := proc(ii,v)
 local dp,u0,u1,u2,vv,i;
 dp := (x,y) -> add(x[i] * y[i],i=1..3);
 u0 := [0,0,0];
 for i in ii do u0 := u0 +~ v[i]; od;
 u0 := evalf(u0 /~ sqrt(dp(u0,u0)));
 u1 := evalf(v[ii[1]]);
 u1 := u1 -~ dp(u1,u0) *~ u0;
 u1 := evalf(u1 /~ sqrt(dp(u1,u1)));
 u2 := [u0[2] * u1[3] - u0[3] * u1[2],
        u0[3] * u1[1] - u0[1] * u1[3],
        u0[1] * u1[2] - u0[2] * u1[1]];
 vv := map(i -> [evalf(arctan(dp(u2,v[i]),dp(u1,v[i]))),i],ii);
 return map(x -> x[2],sort(vv));
end:


text:
program assignchar;	{sample23}
var x,y : char;
begin
x := 'X';y := 'Y'  end.


text:
# We use the definition below for generalised Narayana numbers because it is
# good for m > 0 and for n = 0.  In the most common case where m = 0 and
# n > 0 it is equivalent to the function narayana_number_basic(n,j).

narayana_number := proc(n,j,m := 0)
 (m+1)/(n+1)*binomial(n+1,j) * binomial(n-m-1,j-1);
end:

narayana_number(0,0) := 1;

narayana_number_basic := (n,j) ->
 binomial(n,j) * binomial(n,j-1) / n;

# This is the sum of narayana_number(n,j,0) x^n y^j for n,j >= 0
narayana_series := (x,y) ->
 (1 + x*(1-y) - sqrt(1-2*x*(1+y)+x^2*(1-y)^2))/(2*x);

# This is the series used by Stanley which has slightly different
# conventions for n = 0 and so ends up as 1 + narayana_series(x,y)

narayana_series_stanley := (x,y) ->
 (1-x-x*y - sqrt((1-x-x*y)^2-4*x^2*y))/(2*x);

`is_element/narayana_paths` := (len::posint,peaks::posint) -> proc(u)
 local v,w,y,i;

 if not(type(u,list(integer))) then
  return false;
 fi;

 if {op(u),1,-1} <> {1,-1} then
  return false;
 fi;

 if nops(u) <> 2*len then
  return false;
 fi;

 v := 0;
 y := 0;
 for i from 1 to 2*len do
  y := y + u[i];
  v := v,y;
 od:
 v := [v];

 if v[2*len+1] <> 0 then
  return false;
 fi;

 if min(op(v)) < 0 then
  return false;
 fi;

 w := select(i -> (u[i] = 1 and u[i+1] = -1),[seq(i,i=1..2*len-1)]);
 if nops(w) <> peaks then
  return false;
 fi;

 return true;
end:

`list_elements/narayana_paths` := proc(len,peaks)
 local A,B,L,a,b,u,ok,y,i;
 A := map(a -> [0,op(a),len+1],combinat[choose](len  ,peaks-1));
 B := map(b -> [0,op(b),len  ],combinat[choose](len-1,peaks-1));
 L := NULL;
 for a in A do
  for b in B do
   u := map(op,[seq([(+1)$(a[i+1]-a[i]),(-1)$(b[i+1]-b[i])],i=1..peaks)]);
   ok := true;
   y := 0;
   for i from 1 to 2*len+1 do
    y := y + u[i];
    if y = 0 then
     ok := false;
     break;
    fi;
   od:
   if ok then
    L := L,u[2..-1];
   fi;
  od:
 od:
 return [L];
end:

# This is very unintelligent
`list_elements_basic/narayana_paths` := proc(len,peaks)
 local L,M,i,k,u,h;
 L := [[1]]:
 for i from 1 to 2*len-1 do
  M := NULL;
  for u in L do
   h := `+`(op(u));
   if h > 0 then
    M := M,[op(u),-1],[op(u),1];
   else
    M := M,[op(u),1];
   fi;
  od:
  L := [M];
 od:
 M := NULL;
 for u in L do
  k := nops(select(i -> (u[i]=1 and u[i+1]=-1),[seq(i,i=1..2*len-1)]));
  if k = peaks and `+`(op(u)) = 0 then
   M := M,u;
  fi;
 od:
 return [M];
end:

`count_elements/narayana_paths` := (len,peaks) ->
 narayana_number(len,peaks);

`plot/narayana_path` := proc(u)
 local P,y,i;
 P := NULL;
 y := 0;
 for i from 1 to nops(u) do
  P := P,line([i-1,y],[i,y+u[i]]);
  y := y + u[i];
 od:
 display(P,scaling=constrained,axes=none);
end:

`digraph/narayana_path` := proc(u)
 local n,y,v,V,E,i,j,h;
 n := nops(u);
 y := 0;
 v := NULL;
 for i from 1 to n do
  y := y + u[i];
  v := v,y;
 od:
 v := [v];
 V := [0,op(select(i -> u[i] = 1,[seq(i,i=1..n)]))];
 E := NULL;
 for i in V do
  if i = 0 then h := 0; else h := v[i]; fi;
  j := i+1;
  while j <= n and v[j] >= h do
   if v[j] = h + 1 and u[j] = 1 then
    E := E,[i,j];
   fi;
   j := j+1;
  od;
 od:
 return [V,[E]];
end:

`extended_stasheff_tree/narayana_path` := proc(u)
 local n,V,E,L,ix,c,v,e,m,i,T,T0,t,UU;
 n := nops(u);
 V,E := op(`digraph/narayana_path`(u));
 L := select(i -> nops(select(e -> e[1]=i,E))=0,V);

 ix := table():
 for i from 1 to nops(L) do ix[L[i]] := i; od:

 c := table();
 for v in V do
  c[v] := {};
 od:
 c[0] := {};
 for e in E do
  c[e[1]] := {op(c[e[1]]),e[2]};
 od:
 m := nops(V);
 for i from m to 1 by -1 do
  v := V[i];
  if member(v,L) then
   T[v] := {v};
  else
   T[v] := map(j -> op(T[j]),c[v]);
  fi;
 od:

 t := table():
 for v in V do
  T0 := T[v];
  if type(t[T0],integer) then
   t[T0] := t[T0] + 1;
  else
   t[T0] := 0;
  fi;
 od:

 UU := NULL;
 for T0 in map(op,[indices(t)]) do
  UU := UU,[map(v -> ix[v],T0),t[T0]];
 od:
 UU := {UU};

 return UU;
end:

check_narayana := proc()
 local x,y,d,j,n,m,F0,F1,F2,F3;

 _ASSERT(
  `and`(seq(seq(evalb(narayana_number(n,j) = narayana_number_basic(n,j)),j=0..20),n=1..20)),
  "narayana_number_basic(n,j) = narayana_number(n,j,0) for n > 0"
 );

 F0 := narayana_series(x,y);
 _ASSERT(
  simplify(x * F0^2 + (x*y-x-1) * F0 + 1) = 0,
  "Quadratic relation for Narayana series"
 );

 F1 := narayana_series_stanley(x,y);

 _ASSERT(
  simplify(F0 - F1 - 1) = 0,
  "narayana_series(x,y) = narayana_series_stanley(x,y) + 1"
 );

 d := 9;
 m := 4;
 F2 := expand(convert(expand(series(
  (e*x)^m * narayana_series(e*x,e*y)^(m+1),
  e=0,d+2)),polynom,e));
 F3 := add(add(narayana_number(n,j,m) * x^n * y^j * e^(n+j),j=0..min(n-m,d-n)),n=m..d);
 _ASSERT(simplify(F3 - F2) = 0,"generating function");
end:


text:
(*
 Copyright (C) 2020 Susi Lehtola

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

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

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

lk_p := x -> X2S^2*x^2:
lk_q := u -> X2S^2*u:

(* Equation (10) *)
lk_delta := (p, q) -> 8/81*q^2 - 1/9*p*q + 8/243*p^2:
(* Equation (15) *)
lk_f0 := (x1, x2) -> 1 + params_a_kappa*(2 - 1/(1+x1/params_a_kappa) - 1/(1+x2/params_a_kappa)):
(* Equation (16) *)
lk_x1 := (p, q) -> 5/27*p + lk_delta(p,q) + (5/27*p)^2/params_a_kappa:
(* Equation (17) *)
lk_x2 := (p, q) -> 2*(5/27*p)*lk_delta(p,q)/params_a_kappa + (5/27*p)^3/params_a_kappa^2:

(* Full functional *)
lk_f := (x, u) ->
  lk_f0(lk_x1(lk_p(x),lk_q(u)), lk_x2(lk_p(x),lk_q(u))):

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


text:
# Example from Section 5.1 of the paper

read "../ComputeIdentifiableFunctions.mpl";

model := [
    diff(x1(t), t) = (a * x1(t) - b * x1(t) * x2(t)),
    diff(x2(t), t) = (-c * x2(t) + d * x1(t) * x2(t) + e * u(t)),
    y(t) = x1(t)
]:

# Computing function identifiable from a single experiment
printf("Single-experiment identifiable functions: %a\n", SingleExperimentIdentifiableFunctions(model)):

me := MultiExperimentIdentifiableFunctions(model, simplified_generators=true):
printf("The bound for the number of experiments is %a, this means that the field of multi-experiment identifiable functions is the same, %a\n", me[1], me[3]):

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


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

(* replace: "BesselI\(0, " -> "xc_bessel_I0(" *)
(* replace: "BesselI\(1, " -> "xc_bessel_I1(" *)

$define xc_dimensions_2d

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

(* This is the solution of solve((y-1)*exp(y) = x/Pi) *)
prhg07_y := x -> LambertW(m_max(x/Pi, -0.9999999999) * exp(-1)) + 1:

prhg07_v := y -> Pi/X_FACTOR_2D_C * BesselI(0, y/2):

prhg07_f := (x, u, t) ->
  prhg07_v(prhg07_y(prhg07_C(x, u, t)))/2:

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


text:
# Berechne kinematische Zwangsbedingungen für den Palettierroboter MPL800-YASKAWA aus 3 Viergelenkketten
# Einleitung

# Die kinematischen Zwangsbedingungen werden als Ersetzungsausdruck für die abhängigen Winkel aufgestellt.
#
# palh1m1TE -> MPL800 II-Yaskawa, modellierung der Zwangsbedingungen mit Ausdrücken für trigonometrische Elimination
# kinematic_constraint -> Kinematische Zwangsbedingungen
# Datenblatt des Roboters unter :   https://www.motoman.com/industrial-robots/mpl800-ii
# Voraussetzung
# Viergelenkkette muss mit der Toolbox berechnet worden sein (Arbeitsblatt "fourbar1TE_kinematic_constraints.mw")
# Quelle
# SA Bejaoui: Bejaoui2018_S749; "Modellierung kinematischer Zwangsbedingungen für hybride serielle Roboter mit planaren Parallelmechanismen"
# Autor
# Abderahman Bejaoui
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2018-08
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
# Import
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
# Konstante Winkel bereits hineinschreiben
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Variablen definieren für die Hilfswinkel
#
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
# Ähnliches Vorgehen wie in [1].
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:

# Sicherung der aktuellen Daten
backup_kintmp_qs := kintmp_qs:
backup_kintmp_qt := kintmp_qt:
backup_kintmp_subsexp := kintmp_subsexp:

read sprintf("../codeexport/fourbar1TE/tmp/kinematic_constraints_maple_inert.m"): # mit .m datei wird es schneller
;
kintmp_subsexp_fourbar := kintmp_subsexp:
kintmp_qs_fourbar := kintmp_qs:
kintmp_qt_fourbar := kintmp_qt:
kintmp_subsexp_fourbar:
winkel(6,2):
kintmp_subsexp_fourbar(6,2):
kintmp_subsexp_fourbar(3,2):
winkel_neu := Matrix(6,2):
winkel_neu(1,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(3,2)):
winkel_neu(2,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(4,2)):
winkel_neu(3,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(1,2)):
winkel_neu(4,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(2,2)):
winkel_neu(5,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(5,2)):
winkel_neu(6,2) := subs({qJ1s=phi_s}, kintmp_subsexp_fourbar(6,2)):
# Vorherige Werte für dieses System wieder herstellen
kintmp_qs := backup_kintmp_qs:
kintmp_qt := backup_kintmp_qt:
kintmp_subsexp := backup_kintmp_subsexp:
winkel_neu(1..6,1):
# Viergelenkkette 1: A-M-L-B
winkel1:=<<sin_rho28|sin(rho28_s)>;<cos_rho28|cos(rho28_s)>;<sin_rho312|sin(rho312_s)>;<cos_rho312|cos(rho312_s)>;<sin_rho89|sin(rho89_s)>;<cos_rho89|cos(rho89_s)>>:
for i from 1 to 6 do
	winkel1(i,2):=kintmp_subsexp_fourbar(i,2):
     cos_phi_s :=-(sin(qJ3s)*cos(phi312)+cos(qJ3s)*sin(phi312)):
	sin_phi_s := -(cos(qJ3s)*cos(phi312)-sin(qJ3s)*sin(phi312)):
	winkel1(i,2):=subs({sin(qJ1s)=sin_phi_s},winkel1(i,2)):
	winkel1(i,2):=subs({cos(qJ1s)=cos_phi_s},winkel1(i,2)):
     winkel1(i,2) := subs({l1=BL}, winkel1(i,2)):
     winkel1(i,2) := subs({l2=AB}, winkel1(i,2)):
     winkel1(i,2) := subs({l3=AM}, winkel1(i,2)):
     winkel1(i,2) := subs({l4=ML}, winkel1(i,2)):
end do:
winkel1:
# Umrechnung aus Zwangsbedingugen für Winkel aus MDH-Tabelle
sin_rho28_s:=winkel1(1,2): # schon in MDH
;
cos_rho28_s:=winkel1(2,2):  # schon in MDH
;
sin_rho312_s:=winkel1(3,2):  # schon in MDH
;
cos_rho312_s:=winkel1(4,2):  # schon in MDH
;
sin_rho89_s:=winkel1(5,2):    # für MDH , am 20.11 der neuen fourbar1 angepasst
;
cos_rho89_s:=winkel1(6,2):  # für MDH , am 20.11 der neuen fourbar1 angepasst
;
#sin_rho312_s:=sin_rho3_s:  # schon in MDH
;
#cos_rho312_s:=-cos_rho3_s:  # schon in MDH
;
# Viergelenkkette 2: A-B-C-D
winkel2:=<<sin_eta1|sin(eta1_s)>;<cos_eta1|cos(eta1_s)>;<sin_eta3|sin(eta3_s)>;<cos_eta3|cos(eta3_s)>;<sin_eta4|sin(eta4_s)>;<cos_eta4|cos(eta4_s)>>:
for i from 1 to 6 do
	winkel2(i,2):=kintmp_subsexp_fourbar(i,2):  #
	cos_phi_s :=sin(qJ2s)*cos(phi2)-cos(qJ2s)*sin(phi2):
	sin_phi_s :=cos(qJ2s)*cos(phi2)+sin(qJ2s)*sin(phi2):
	winkel2(i,2):=subs({sin(qJ1s)=sin_phi_s},winkel2(i,2)):
	winkel2(i,2):=subs({cos(qJ1s)=cos_phi_s},winkel2(i,2)):
     winkel2(i,2) := subs({l1=AB}, winkel2(i,2)):
     winkel2(i,2) := subs({l2=DA}, winkel2(i,2)):
     winkel2(i,2) := subs({l3=DC}, winkel2(i,2)):
     winkel2(i,2) := subs({l4=BC}, winkel2(i,2)):
 end do:
winkel2(1..6,1):
sin_eta1_s:=winkel2(1,2):
cos_eta1_s:=winkel2(2,2):
sin_eta3_s:=winkel2(3,2):
cos_eta3_s:=winkel2(4,2):
sin_eta4_s:=winkel2(5,2):
cos_eta4_s:=winkel2(6,2):
cos_eta16_s:=cos(phi1)*cos_eta1_s-sin(phi1)*sin_eta1_s:
sin_eta16_s:=cos(phi1)*sin_eta1_s+cos_eta1_s*sin(phi1):
cos_eta711_s:=-cos_eta4_s:
sin_eta711_s:=sin_eta4_s:
cos_eta27_s:=cos(phi711)*cos_eta3_s+sin(phi711)*sin_eta3_s: # am 20.11 der neuen fourbar1 angepasst
;
sin_eta27_s:=cos(phi711)*sin_eta3_s-cos_eta3_s*sin(phi711): # am 20.11 der neuen fourbar1 angepasst
;
# Viergelenkkette 3: B-E-P-G
winkel3:=<<sin_xi1016|sin(xi1016_s)>;<cos_xi1016|cos(xi1016_s)>;<sin_xi3|sin(xi3_s)>;<cos_xi3|cos(xi3_s)>;<sin_xi2|sin(xi2_s)>;<cos_xi2|cos(xi2_s)>>:
for i from 1 to 6 do
	winkel3(i,2):=kintmp_subsexp_fourbar(i,2):
	cos_phi_s :=(cos(phi711+phi710)*(sin(qJ3s)*cos_eta3_s-cos(qJ3s)*sin_eta3_s)+sin(phi711+phi710)*(cos(qJ3s)*cos_eta3_s+sin(qJ3s)*sin_eta3_s)):
	sin_phi_s := (cos(phi711+phi710)*(cos(qJ3s)*cos_eta3_s+sin(qJ3s)*sin_eta3_s)-sin(phi711+phi710)*(sin(qJ3s)*cos_eta3_s-cos(qJ3s)*sin_eta3_s)):
	winkel3(i,2):=subs({sin(qJ1s)=sin_phi_s},winkel3(i,2)):
	winkel3(i,2):=subs({cos(qJ1s)=cos_phi_s},winkel3(i,2)):
     winkel3(i,2) := subs({l1=BG}, winkel3(i,2)):
     winkel3(i,2) := subs({l2=BE}, winkel3(i,2)):
     winkel3(i,2) := subs({l3=EP}, winkel3(i,2)):
     winkel3(i,2) := subs({l4=GP}, winkel3(i,2)):
 end do:
winkel3(1..6,1):
sin_xi1016_s:=winkel3(1,2):
cos_xi1016_s:=winkel3(2,2):
sin_xi3_s:=winkel3(3,2):  # schon in MDH
;
cos_xi3_s:=winkel3(4,2): # schon in MDH
;
sin_xi2_s:=winkel3(5,2):
cos_xi2_s:=winkel3(6,2):
cos_xi34_s:=cos_xi3_s*cos(phi413)+sin_xi3_s*sin(phi413): # am 20.11 der neuen fourbar1 angepasst
;
sin_xi34_s:=-cos_xi3_s*sin(phi413)+sin_xi3_s*cos(phi413): # am 20.11 der neuen fourbar1 angepasst
;
cos_xi413_s:=-cos_xi2_s:
sin_xi413_s:=sin_xi2_s:
# Kintmp_subsexp
kintmp_subsexp(1,2):=sin_xi34_s:
kintmp_subsexp(2,2):=cos_xi34_s:
kintmp_subsexp(3,2):=sin_eta16_s:
kintmp_subsexp(4,2):=cos_eta16_s:
kintmp_subsexp(5,2):=sin_eta27_s:
kintmp_subsexp(6,2):=cos_eta27_s:
kintmp_subsexp(7,2):=sin_rho28_s:
kintmp_subsexp(8,2):=cos_rho28_s:
kintmp_subsexp(9,2):=sin_rho89_s:
kintmp_subsexp(10,2):=cos_rho89_s:
kintmp_subsexp(11,2):=sin_xi1016_s:
kintmp_subsexp(12,2):=cos_xi1016_s:
kintmp_subsexp(13,2):=sin_eta711_s:
kintmp_subsexp(14,2):=cos_eta711_s:
kintmp_subsexp(15,2):=sin_rho312_s:
kintmp_subsexp(16,2):=cos_rho312_s:
kintmp_subsexp(17,2):=sin_xi413_s:
kintmp_subsexp(18,2):=cos_xi413_s:
kintmp_subsexp(11..18,1):
# Kintmp_qs
kintmp_qs(1,1):=%arctan(sin_xi34_s,cos_xi34_s):
kintmp_qs(2,1):=%arctan(sin_eta16_s,cos_eta16_s):
kintmp_qs(3,1):=%arctan(sin_eta27_s,cos_eta27_s):
kintmp_qs(4,1):=%arctan(sin_rho28_s,cos_rho28_s):
kintmp_qs(5,1):=%arctan(sin_rho89_s,cos_rho89_s):
kintmp_qs(6,1):=%arctan(sin_xi1016_s,cos_xi1016_s):
kintmp_qs(7,1):=%arctan(sin_eta711_s,cos_eta711_s):
kintmp_qs(8,1):=%arctan(sin_rho312_s,cos_rho312_s):
kintmp_qs(9,1):=%arctan(sin_xi413_s,cos_xi413_s):
# Export
kintmp_qt := convert_s_t(kintmp_qs):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
#printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d",robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert" ,robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
#printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
#kc_symbols :=Transpose(kc_symbols);
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);
#printf("Fertig. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):