text:
# Definitions for Robot Dynamics Code Generation
# Init
# Erstelle Definitionen für die Maple-Skripte zur Berechnung von Kinematik und Dynamik des Roboters
#
# Quellen:
# [1] Sousa, C. D. and Cortesao, R.: Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach (2014)
# [2] Ayusawa, K. and Venture, G. and Nakamura, Y.: Identifiability and identification of inertial parameters using the underactuated base-link dynamics for legged multibody systems (2013)
# [3] Fujimoto, Y. and Obata, S. and Kawamura, A.: Robust biped walking with active interaction control between foot and ground (1998)
# [4] Khalil, W. and Kleinfinger, J.-F.: Minimum operations and minimum parameters of the dynamic models of tree structure robots (1987)
#
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
read "../helper/proc_MatlabExport":
read "../helper/proc_convert_t_s":
read "../helper/proc_vec2skew":
codegen_act := true:
# Lese Umgebungsvariable für Codegenerierung.
read "../robot_codegen_definitions/robot_env":
printf("Generiere Parameter für %s\n",robot_name):
# Robotics Definitions
# Joint Angles, Velocities and Accelerations in time depending and substitution form.
# Time Depending form is for time differentiation.
# Substitution Form is for partial differentiation and code export.
qJ_t := Matrix(NQJ, 1):
qJ_s := Matrix(NQJ, 1):
qJD_s := Matrix(NQJ, 1):
qJDD_s := Matrix(NQJ, 1):
for i from 1 to NQJ do
  qJ_t(i,1):=parse(sprintf("qJ%d(t)", i)):
  qJ_s(i,1):=parse(sprintf("qJ%ds", i)):
  qJD_s(i,1):=parse(sprintf("qJD%ds", i)):
  qJDD_s(i,1):=parse(sprintf("qJDD%ds", i)):
end do:
qJD_t := diff~(qJ_t, t):
qJDD_t := diff~(qJD_t, t):
# Gravity vector in world frame
if not assigned(g_world) then
  g_world := Matrix(3, 1):
  g_world(1 .. 3, 1) := <g1, g2, g3>:
end if:
# Position und Orientierung der Basis. Die Orientierung ist mit XYZ-Euler-Winkeln definiert, die aber nicht im weiteren Algorithmus verwendet werden (nur Platzhalter).
# Eine Invertierung der Orientierungsdarstellung sollte nicht notwendig werden, von daher kein Problem mit Orientierungsrepräsentationssingularität.
# gem. [2], S. 4 X_base_t SE(3): Position und Orientierung.
NQB := 6:
X_base_t:=Matrix(6, 1, [rx_base(t),ry_base(t),rz_base(t),alphax_base(t), betay_base(t), gammaz_base(t)]):
X_base_s:=Matrix(6, 1, [rxs_base,rys_base,rzs_base,alphaxs_base, betays_base, gammazs_base]):
# Geschwindigkeit der Basis (pelvis)
# Twist: Verallgemeinerte Geschwindigkeit. Diese Darstellung ist nicht konsistent (Ableitung der Position stimmt nicht mit Ausdruck der Geschwindigkeit überein)
V_base_t:=Matrix(6, 1, [diff~(rx_base(t),t),diff~(ry_base(t),t),diff~(rz_base(t),t),omegax_base(t), omegay_base(t), omegaz_base(t)]):
V_base_s:=Matrix(6, 1, [vxs_base,vys_base,vzs_base,omegaxs_base, omegays_base, omegazs_base]):
# Zeitableitung davon
VD_base_t:=diff~(V_base_t,t):
VD_base_s:=Matrix(6, 1, [vDxs_base,vDys_base,vDzs_base,omegaDxs_base, omegaDys_base, omegaDzs_base]):
# Name der Methode für die Orientierungsdarstellung der Basis
base_method_name := "twist":
# Umrechnung von der Ableitung der Basis-Orientierung zu Winkelgeschwindigkeiten:
# Wird hier nicht weiter beachtet (physikalisch falsch für eine Betrachtung der Dynamik der Basis).
# Ermöglicht eine einfachere Berechnung für Mechanismen, wo die Rückwirkung auf die Basis nicht in der Vorwärtsdynamik simuliert werden muss.
T_basevel := IdentityMatrix(3, 3):
# Verallgemeinerte Koordinaten, gem [2], S. 4, [3], S.1
NQ:=NQJ+NQB:
q_t := Matrix(NQ,1, <X_base_t, qJ_t>):
q_s := Matrix(NQ,1, <X_base_s, qJ_s>):
qD_t:= Matrix(NQ,1, <V_base_t, qJD_t>):
qD_s:= Matrix(NQ,1, <V_base_s, qJD_s>):
qDD_t:= Matrix(NQ,1, <VD_base_t, qJDD_t>):
qDD_s:= Matrix(NQ,1, <VD_base_s, qJDD_s>):
# MDH-Gelenkwinkel neu speichern (Definition der verallg. Koordinaten war dort noch nicht bekannt
theta := value(theta):
# Standard-Werte festlegen
# Anzahl der Körper (Number of Links):
if not assigned(NL) then
  NL := NJ + 1:
  printf("Variable NL ist nicht gegeben. Insgesamt %d Gelenke. Nehme an, dass jedes Gelenk einem Körper zugeordnet ist (keine Schleifen)\n", NJ):
else
  NVJ := NJ - (NL - 1):
  printf("Variable NL=%d ist gegeben. Insgesamt %d Gelenke. Davon sind die ersten %d einem Körper zugeordnet und die letzten %d virtuell.\n", NL, NJ, NJ-NVJ, NVJ):
end if:
# Gelenktyp (0=Revolute, 1=Prismatic, 2=Static). Sollte in der Definition festgelegt sein. Falls nicht, wird alles auf Revolute gesetzt
if not assigned(sigma) then
  sigma := Matrix(NJ,1):
  printf("Variable sigma ist nicht gegeben. Setze alle Gelenke als Drehgelenk (sigma=1)\n"):
end if:
# Aktuierung (1=aktiv, 0=passiv). Sollte in der Definition festgelegt sein. Falls nicht, wird alles auf Aktiv gesetzt
if not assigned(mu) then
  mu := Matrix(NJ,1,1):
  printf("Variable mu ist nicht gegeben. Setze alle Gelenke als aktiv (mu=1)\n"):
end if:
# Segmentreihenfolge der Roboterstruktur. Sollte in der Definition festgelegt sein. Falls nicht, wird eine serielle Kette angenommen
if not assigned(v) then
  v := Matrix(NL-1,1,i->i-1):
  printf("Variable v ist nicht gegeben. Setze serielle Kette (0,1,2,...)\n"):
end if:
# Parameter für Baumstruktur. Sollte in der Definition festgelegt sein. Falls nicht, wird eine serielle Kette angenommen
if not assigned(b) then
  b := Matrix(NJ,1):
  printf("Variable b ist nicht gegeben. zusätzliche Verschiebung auf Null\n"):
end if:
if not assigned(beta) then
  beta := Matrix(NJ,1):
  printf("Variable beta ist nicht gegeben. zusätzliche Verschiebung auf Null\n"):
end if:
# Eingabe prüfen
# Validität der Vorgänger-Indizes. Konvention: Vorgänger hat kleineren Index. Siehe Khalil-Publikationen (dort evtl nur implizit so verwendet).
for i from 1 to NJ do
  if v(i) >= i then
    printf("v(%d)=%d. Verletzt die Regel, dass die Vorgänger-Nummer kleiner sein muss als die aktuelle Nummer.\n", i, v(i)):
    quit: # Funktioniert in GUI nicht richtig...
    robot_name := "": # ...Daher auch Löschung des Roboternamens.
  end if:
end do:

# Dynamics Parameters
# Mass of each link
M_generic := Matrix(NL, 1):
for i from 1 to NL do
  M_generic[i,1]:=parse(sprintf("M%d", i-1)):
end do:
M := copy(M_generic):
if assigned(user_M) then
  if not RowDimension(user_M) = NL then
    printf("Input user_M is not of size %dx1. Error.\n", NL):
  end if:
  # Mass parameter given by user. Some entries can be zero or equal
  for i from 1 to NL do
    if user_M[i,1]<>M_generic[i,1] and user_M[i,1]<>0 then # check input
      printf("User input for mass %d should be either \"%s\" or zero. Not \"%s\". Continue anyway. \n", i-1, M_generic[i,1], user_M[i,1]):
    end if:
    M[i,1]:=user_M[i,1]:
  end do:
  printf("Masseparameter von Benutzer gesetzt. Eingabe für M:\n"):
  print(M);
end if:
# Center of Mass of each link (in link frame)
r_i_i_Si := Matrix(3, NL):
r_i_i_Si_generic := Matrix(3, NL):
for i from 1 to NL do
  r_i_i_Si_generic[1,i]:=parse(sprintf("SX%d", i-1)):
  r_i_i_Si_generic[2,i]:=parse(sprintf("SY%d", i-1)):
  r_i_i_Si_generic[3,i]:=parse(sprintf("SZ%d", i-1)):
end do:
if not assigned(user_CoM) then
  r_i_i_Si := r_i_i_Si_generic:
else
  if RowDimension(user_CoM) <> 3 or ColumnDimension(user_CoM) <> NL then
    printf("Input user_CoM is not of size 3x%d. Error.\n", NL):
  end if:
  for i from 1 to NL do
    for j from 1 to 3 do
      if user_CoM[j,i]<>r_i_i_Si_generic[j,i] and user_CoM[j,i]<>0 then # check input
        printf("User input for CoM %d,%d should be either \"%s\" or zero. Not \"%s\". Continue anyway. \n", j, i-1, r_i_i_Si_generic[j,i], convert(user_CoM[j,i],string)):
      end if:
    end do:
    r_i_i_Si[1..3,i]:=user_CoM[1..3,i]:
  end do:
  printf("Schwerpunktsparameter von Benutzer gesetzt. Eingabe für r_i_i_Si:\n"):
  print(r_i_i_Si);
end if:

# First Moment (mass and center of mass)
mr_i_i_Si := Matrix(3, NL):
mr_i_i_Si_generic := Matrix(3, NL):
xyzstrings := ["X", "Y", "Z"]:
for i from 1 to NL do
  for j from 1 to 3 do # loop over x-, y-, z-coordinates of the CoM
    mr_i_i_Si_generic[j,i]:= parse(sprintf("M%s%d", xyzstrings[j], i-1)):
  end do:
  if M[i,1] = 0 then next; end if: # Zero Mass can be set by the user. Then the first moment stays zero
  for j from 1 to 3 do # loop over x-, y-, z-coordinates of the CoM
    if r_i_i_Si[j,i] = parse(sprintf("S%s%d", xyzstrings[j], i-1)) then
    	 # default value for CoM is set (not overwritten by user input).
    	 # Set default value for first moment
      mr_i_i_Si[j,i]:= parse(sprintf("M%s%d", xyzstrings[j], i-1)):
    else
      # CoM is written by the user. Put this assumption also in the first moment to reduce dynamics parameters
    	 mr_i_i_Si[j,i] := M[i,1]*r_i_i_Si[j,i]:
    end if:
  end do:
end do:

# Inertia of each link (about the center of mass, in link frame)
I_i_Si := Matrix(6, NL):
I_i_Si_generic := Matrix(6, NL):
for i from 1 to NL do
  if M[i,1] = 0 then next; end if: # Zero Mass can be set by the user. Then the inertia stays zero
  I_i_Si_generic[1,i]:=parse(sprintf("XXC%d", i-1)):
  I_i_Si_generic[2,i]:=parse(sprintf("XYC%d", i-1)):
  I_i_Si_generic[3,i]:=parse(sprintf("XZC%d", i-1)):
  I_i_Si_generic[4,i]:=parse(sprintf("YYC%d", i-1)):
  I_i_Si_generic[5,i]:=parse(sprintf("YZC%d", i-1)):
  I_i_Si_generic[6,i]:=parse(sprintf("ZZC%d", i-1)):
end do:
if not assigned(user_inertia) then
  I_i_Si := I_i_Si_generic:
else
  if RowDimension(user_inertia) <> 6 or ColumnDimension(user_inertia) <> NL then
    printf("Input user_inertia is not of size 6x%d. Error.\n", NL):
  end if:
  for i from 1 to NL do
    for j from 1 to 6 do
      if user_inertia[j,i]<>I_i_Si_generic[j,i] and user_inertia[j,i]<>0 then # check input
        printf("User input for inertia %d,%d should be either \"%s\" or zero. Not \"%s\". Continue anyway. \n", j, i-1, convert(I_i_Si_generic[j,i],string), convert(user_inertia[j,i],string)):
      end if:
    end do:
    I_i_Si[1..6,i]:=user_inertia[1..6,i]:
  end do:
  printf("Trägheitsparameter von Benutzer gesetzt. Eingabe für I_i_Si:\n"):
  print(I_i_Si);
end if:

# Inertia of each link (about the origin of body frame, in link frame)
# Berechne die Inertial-Trägheitsmomente mit dem Steinerschen Verschiebungssatz
I_i_i_calc := Matrix(6, NL):
for i from 1 to NL do
  # Trägheitstensor (3x3) um den Körperschwerpunkt in Körper-KS
  I_i_Si_Tensor := Matrix([[I_i_Si[1, i], I_i_Si[2, i], I_i_Si[3, i]], [I_i_Si[2, i], I_i_Si[4, i], I_i_Si[5, i]], [I_i_Si[3, i], I_i_Si[5, i], I_i_Si[6, i]]]):
  # Steinerschen Satz anwenden
  I_i_i_Tensor := I_i_Si_Tensor + M[i,1]*Transpose(vec2skew(r_i_i_Si[1..3,i])) . vec2skew(r_i_i_Si[1..3,i]):
  # Wieder als Matrix abspeichern
  I_i_i_calc[..,i] := <I_i_i_Tensor[1,1]; I_i_i_Tensor[1,2]; I_i_i_Tensor[1,3]; I_i_i_Tensor[2,2]; I_i_i_Tensor[2,3]; I_i_i_Tensor[3,3]>:
end do:

# Allgemeine Form des Trägheitstensors (Einträge sind unabhängige Parameter)
I_i_i_generic := Matrix(6, NL):
for i from 1 to NL do
  I_i_i_generic[1,i]:=parse(sprintf("XX%d", i-1)):
  I_i_i_generic[2,i]:=parse(sprintf("XY%d", i-1)):
  I_i_i_generic[3,i]:=parse(sprintf("XZ%d", i-1)):
  I_i_i_generic[4,i]:=parse(sprintf("YY%d", i-1)):
  I_i_i_generic[5,i]:=parse(sprintf("YZ%d", i-1)):
  I_i_i_generic[6,i]:=parse(sprintf("ZZ%d", i-1)):
end do:

# Prüfe, welche Einträge des Trägheitstensors noch Schwerpunkts-Parameter enthalten. Diese müssen wieder auf die Standard-Werte mit unabhängigen Parametern gesetzt werden
I_i_i := Matrix(6, NL):
compstrings := ["XX", "XY", "XZ", "YY", "YZ", "ZZ"]:
for i from 1 to NL do
  if M[i,1] = 0 then next; end if: # Zero Mass can be set by the user. Then the inertia stays zero
  for j from 1 to 6 do # Alle Komponenten des Tensors
    # Nicht Null. Setze erstmal allgemeinen Eintrag:
    I_i_i[j,i] := I_i_i_generic[j,i]:
    # Prüfe, ob Schwerpunktsparameter vorkommt (durch Verschiebungssatz s.o.)
    IhasCoM := false:
    for k from 1 to 3 do
      if has(I_i_i_calc[j,i], r_i_i_Si_generic[k,i]) then
    	   IhasCoM := true: # Der Eintrag enthält die Schwerpunktskoordinaten. Nicht zulässig für weitere Rechnung.
    	   break:
      end if:
    end do:
    if has(I_i_i_calc[j,i], I_i_Si_generic[j,i]) then
    	 IhasCoM := true: # Es steht der Schwerpunktsbezogene Trägheitsterm drin. Weitere Rechnung damit nicht möglich.
    end if:
    if IhasCoM = true then
      next: # Es kommt ein Schwerpunktsparameter vor. Belasse Parameter auf allgemeinem Wert.
    end if:
    I_i_i[j,i] := I_i_i_calc[j,i]:
  end do:
end do:

# Vector of stacked dynamics parameters for regressor form
# Matrix of link inertial parameters, stacked link parameter vectors.
# Diese Parameter-Matrix wird nur benutzt, um für die Generierung der Regressorform danach abzuleiten.
# Hier stehen also die allgemeinen ("generic") Parameter drin. Es ist egal, ob diese bereits durch den Benutzer zu Null gesetzt sind.
# Dann würden die allgemeinen Parameter nicht in der Dynamik vorkommen und die Ableitung wäre Null.
PV2_mat := Matrix(NL, 10):
for i to NL do
  PV2_mat[i, 1 .. 6] := I_i_i_generic[1 .. 6, i]:
  PV2_mat[i, 7 .. 9] := mr_i_i_Si_generic[1 .. 3, i]:
  PV2_mat[i, 10] := M_generic[i, 1]:
end do:

# Parameter-Vektor Erstellen: vector of link inertial parameters (delta in [1]).
# Gleiche Überlegung wie für Parameter-Matrix
PV2_vec := Matrix(10*(NL), 1):
for i to NL do
  for j to 10 do
    PV2_vec[10*(i-1)+j] := PV2_mat[i, j]:
  end do:
end do:
# Kinematische Zwangsbedingungen
# Prüfe, ob kinematische Zwangsbedingungen in der Roboterkonfiguration genannt sind durch Prüfung der Existenz der entsprechenden Variablen.
if type( kintmp_t, 'Matrix') = false then
  kintmp_t := Matrix(1,1): # Dummy-Werte damit später alles funktioniert
  kintmp_s := Matrix(1,1):
end if:
# Export
save q_t, q_s, qD_t, qD_s, qDD_t, qDD_s, qJ_t, qJ_s, qJD_t, qJD_s, qJDD_t, qJDD_s, g_world, X_base_t, X_base_s, V_base_t, V_base_s, VD_base_t, VD_base_s, qoffset, theta, alpha, d, a,v,b,beta, sigma,mu,M, r_i_i_Si, mr_i_i_Si, I_i_i, I_i_Si, PV2_vec, PV2_mat, robot_name, NQ,NQB,NQJ,NJ,NL, base_method_name, T_basevel, kintmp_t, kintmp_s, sprintf("../codeexport/%s/tmp/tree_floatb_twist_definitions", robot_name):
save q_t, q_s, qD_t, qD_s, qDD_t, qDD_s, qJ_t, qJ_s, qJD_t, qJD_s, qJDD_t, qJDD_s, g_world, X_base_t, X_base_s, V_base_t, V_base_s, VD_base_t, VD_base_s, qoffset, theta, alpha, d, a,v,b,beta, sigma,mu,M, r_i_i_Si, mr_i_i_Si, I_i_i, I_i_Si, PV2_vec, PV2_mat, robot_name, NQ,NQB,NQJ,NJ,NL, base_method_name, T_basevel, kintmp_t, kintmp_s, sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Einzelne DH-Parameter als Matlab-Code exportieren. Damit lässt sich in Matlab ein passender Parametersatz generieren.
# Die folgenden Dateien müssen immer generiert werden, wenn die Bash-Skripte zur Funktionszusammensetzung fehlerfrei laufen sollen.
# Benutze die Funktion convert_t_s, um eventuelle Substitutionsvariablen für konstante Gelenkwinkel zu verwenden, da die Matlab-Terme auch mit substituierten Ausdrücken generiert werden.
# Zur Kennzeichnung von zeitabhängigen und konstanten Ausdrücken kann "delta1(t)", "delta1" und "delta1s" verwendet werden.
MatlabExport(v, sprintf("../codeexport/%s/tmp/parameters_mdh_v_matlab.m", robot_name), 2):
MatlabExport(convert_t_s(a), sprintf("../codeexport/%s/tmp/parameters_mdh_a_matlab.m", robot_name), 2):
d_export := d *~ (1-~sigma):
MatlabExport(convert_t_s(d_export), sprintf("../codeexport/%s/tmp/parameters_mdh_d_matlab.m", robot_name), 2):
theta_export := theta *~ sigma:
MatlabExport(convert_t_s(theta_export), sprintf("../codeexport/%s/tmp/parameters_mdh_theta_matlab.m", robot_name), 2):
MatlabExport(convert_t_s(b), sprintf("../codeexport/%s/tmp/parameters_mdh_b_matlab.m", robot_name), 2):
MatlabExport(convert_t_s(alpha), sprintf("../codeexport/%s/tmp/parameters_mdh_alpha_matlab.m", robot_name), 2):
MatlabExport(convert_t_s(beta), sprintf("../codeexport/%s/tmp/parameters_mdh_beta_matlab.m", robot_name), 2):
MatlabExport(convert_t_s(qoffset), sprintf("../codeexport/%s/tmp/parameters_mdh_qoffset_matlab.m", robot_name), 2):
MatlabExport(sigma, sprintf("../codeexport/%s/tmp/parameters_mdh_sigma_matlab.m", robot_name), 2):
MatlabExport(mu, sprintf("../codeexport/%s/tmp/parameters_mdh_mu_matlab.m", robot_name), 2):

# Einzelne Dynamikparameter als Matlab-Code exportieren. Wenn die Parameter durch Benutzereingaben verändert wurden, lässt sich diese Information so weiter benutzen.
# (z.B. in der Definition von Eingabeparametern in den Testskripten).
# In den Maple-Variablen ist die Spalte der Index der Körper und die Zeilen sind die Indizes für die xyz-Komponenten (einfacherer Aufruf der Variablen)
# In Matlab ist es umgekehrt (führt zu konsistenterem Code in Matlab): Die Zeilen entsprechen dem Körper-Index. Daher hier Transponierung.
# Zusätzlich ist die Reihenfolge der Komponenten der Trägheitstensoren in Matlab und Maple unterschiedlich (daher die Indizierung).
# Matlab: xx, yy, zz, xy, xz, yz (erst Hauptmomente, dann Deviationsmomente)
# Maple: xx, xy, xz, yy, yz, zz (Dreiecksform)
MatlabExport(M, sprintf("../codeexport/%s/tmp/parameters_dyn_mges_matlab.m", robot_name), 2):
MatlabExport(Transpose(r_i_i_Si), sprintf("../codeexport/%s/tmp/parameters_dyn_rSges_matlab.m", robot_name), 2):
MatlabExport(Transpose(I_i_Si([1,4,6,2,3,5],..)), sprintf("../codeexport/%s/tmp/parameters_dyn_Icges_matlab.m", robot_name), 2):
MatlabExport(Transpose(mr_i_i_Si), sprintf("../codeexport/%s/tmp/parameters_dyn_mrSges_matlab.m", robot_name), 2):
MatlabExport(Transpose(I_i_i([1,4,6,2,3,5],..)), sprintf("../codeexport/%s/tmp/parameters_dyn_Ifges_matlab.m", robot_name), 2):

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

# MarControls = 12
# /Series/playbacknet/bool/debug = 0
# /Series/playbacknet/bool/mute = 0
# /Series/playbacknet/bool/probe = 0
# /Series/playbacknet/natural/inObservations = 1
# /Series/playbacknet/natural/inSamples = 512
# /Series/playbacknet/natural/onObservations = 1
# /Series/playbacknet/natural/onSamples = 512
# /Series/playbacknet/real/israte = 22050
# /Series/playbacknet/real/osrate = 22050
# /Series/playbacknet/realvec/input = # MARSYAS realvec
# Size = 0

# type: matrix
# rows: 1
# columns: 0

# Size = 0
# MARSYAS realvec

# /Series/playbacknet/string/inObsNames = audio,
# /Series/playbacknet/string/onObsNames = audio,

# Number of links = 6
# Synonyms of /Series/playbacknet/bool/mute =
# Number of synonyms = 1
# /Series/playbacknet/Gain/gt/bool/mute
# Synonyms of /Series/playbacknet/bool/hasData =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/bool/hasData
# Synonyms of /Series/playbacknet/natural/nChannels =
# Number of synonyms = 2
# /Series/playbacknet/SoundFileSource/src/natural/nChannels
# /Series/playbacknet/AudioSink/dest/natural/nChannels
# Synonyms of /Series/playbacknet/natural/pos =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/natural/pos
# Synonyms of /Series/playbacknet/real/israte =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/real/israte
# Synonyms of /Series/playbacknet/string/filename =
# Number of synonyms = 1
# /Series/playbacknet/SoundFileSource/src/string/filename

# nComponents = 5

# MarSystem
# Type = SoundFileSource
# Name = src

# MarControls = 23
# /SoundFileSource/src/bool/advance = 0
# /SoundFileSource/src/bool/debug = 0
# /SoundFileSource/src/bool/mute = 0
# /SoundFileSource/src/bool/hasData = 0
# /SoundFileSource/src/bool/noteon = 0
# /SoundFileSource/src/natural/inObservations = 1
# /SoundFileSource/src/natural/inSamples = 512
# /SoundFileSource/src/natural/loopPos = 0
# /SoundFileSource/src/natural/nChannels = 1
# /SoundFileSource/src/natural/numFiles = 0
# /SoundFileSource/src/natural/onObservations = 1
# /SoundFileSource/src/natural/onSamples = 512
# /SoundFileSource/src/natural/pos = 0
# /SoundFileSource/src/natural/size = 661500
# /SoundFileSource/src/real/duration = -1
# /SoundFileSource/src/real/frequency = 0
# /SoundFileSource/src/real/israte = 22050
# /SoundFileSource/src/real/osrate = 22050
# /SoundFileSource/src/real/repetitions = 1
# /SoundFileSource/src/string/allfilenames = ,
# /SoundFileSource/src/string/filename = defaultfile
# /SoundFileSource/src/string/inObsNames = audio,
# /SoundFileSource/src/string/onObsNames = audio,

# Number of links = 0

# MarSystem
# Type = Gain
# Name = gt

# MarControls = 11
# /Gain/gt/bool/debug = 0
# /Gain/gt/bool/mute = 0
# /Gain/gt/natural/inObservations = 1
# /Gain/gt/natural/inSamples = 512
# /Gain/gt/natural/onObservations = 1
# /Gain/gt/natural/onSamples = 512
# /Gain/gt/real/gain = 1
# /Gain/gt/real/israte = 22050
# /Gain/gt/real/osrate = 22050
# /Gain/gt/string/inObsNames = audio,
# /Gain/gt/string/onObsNames = audio,

# Number of links = 0

# MarSystem
# Type = AudioSink
# Name = dest

# MarControls = 13
# /AudioSink/dest/bool/debug = 0
# /AudioSink/dest/bool/init = 0
# /AudioSink/dest/bool/mute = 0
# /AudioSink/dest/natural/bufferSize = 512
# /AudioSink/dest/natural/inObservations = 1
# /AudioSink/dest/natural/inSamples = 512
# /AudioSink/dest/natural/nChannels = 1
# /AudioSink/dest/natural/onObservations = 1
# /AudioSink/dest/natural/onSamples = 512
# /AudioSink/dest/real/israte = 22050
# /AudioSink/dest/real/osrate = 22050
# /AudioSink/dest/string/inObsNames = audio,
# /AudioSink/dest/string/onObsNames = audio,

# Number of links = 0

# MarSystem
# Type = Gain
# Name = Default

# MarControls = 11
# /Gain/Default/bool/debug = 0
# /Gain/Default/bool/mute = 0
# /Gain/Default/natural/inObservations = 1
# /Gain/Default/natural/inSamples = 512
# /Gain/Default/natural/onObservations = 1
# /Gain/Default/natural/onSamples = 512
# /Gain/Default/real/gain = 1
# /Gain/Default/real/israte = 22050
# /Gain/Default/real/osrate = 22050
# /Gain/Default/string/inObsNames = audio,
# /Gain/Default/string/onObsNames = audio,

# Number of links = 0

# MarSystem
# Type = Gain
# Name = Default

# MarControls = 11
# /Gain/Default/bool/debug = 0
# /Gain/Default/bool/mute = 0
# /Gain/Default/natural/inObservations = 1
# /Gain/Default/natural/inSamples = 512
# /Gain/Default/natural/onObservations = 1
# /Gain/Default/natural/onSamples = 512
# /Gain/Default/real/gain = 1
# /Gain/Default/real/israte = 22050
# /Gain/Default/real/osrate = 22050
# /Gain/Default/string/inObsNames = audio,
# /Gain/Default/string/onObsNames = audio,

# Number of links = 0

ext: mpl

max_stars_repo_name: jaouahbi/VampPlugins

max_issues_repo_name: jaouahbi/VampPlugins

max_forks_repo_name: jaouahbi/VampPlugins


text:
with(DifferentialThomas):
with(Tools):
Ranking([t], [[x1, x2], [y1] ]):
sys := [
diff(x1(t), t) - (a*x1(t) + b*x1(t) - x2(t)*c*x1(t)),
diff(x2(t), t) - (-a*b*x2(t) + d*x2(t)*x1(t)),
y1(t) - (x1(t))
];
res := CodeTools[CPUTime](ThomasDecomposition(sys));


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

$include "gga_c_scan_e0.mpl"
$include "mgga_x_scan.mpl"

scan_b1c := 0.0285764:
scan_b2c := 0.0889:
scan_b3c := 0.125541:
scan_eclda0 := rs -> -scan_b1c/(1 + scan_b2c*sqrt(rs) + scan_b3c*rs):

scan_chi_infty := 0.12802585262625815:
scan_g_infty := s -> 1/(1 + 4*scan_chi_infty*s^2)^(1/4):

(* in the paper it is 2.3631 *)
scan_G_cnst := 2.363:
scan_Gc := z -> (1 - scan_G_cnst*(2^(1/3) - 1)*f_zeta(z))*(1 - z^12):

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

scan_alpha := (z, xt, ts0, ts1) ->
  (t_total(z, ts0, ts1) - xt^2/8)/(K_FACTOR_C*t_total(z, 1, 1)):

(* set parameters of f_alpha *)
params_a_c1 := 0.64:
params_a_c2 := 1.5:
params_a_d  := 0.7:

scan_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  f_pbe(rs, z, xt, xs0, xs1) + scan_f_alpha(scan_alpha(z, xt, ts0, ts1))*(
    + scan_e0(rs, z, X2S*2^(1/3)*xt)
    - f_pbe(rs, z, xt, xs0, xs1)
  ):

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


text:
######################################################################
# W1(A) is Map(A,RR) modulo translation, except that we take
# W1(emptyset) to be empty.
#
# Elements are represented as tables indexed by A, whose values
# are real numbers.  The equivalence relation
# is built into the definition of the function `is_equal/W1` rather
# than the representation of the elements.

`is_element/W1` := (A::set) -> proc(x)
 local a;
 global reason;

 if A = {} then
  reason := [convert(procname,string),"W1(emptyset) = emptyset"];
  return false;
 fi;

 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/RR`(x[a]) then
   reason := [convert(procname,string),"x[a] is not in R",a,x[a],N];
   return false;
  fi;
 od;

 return true;
end;

`is_equal/W1` := (A::set) -> proc(x,y)
 local a0,u0,a;
 global reason;

 if A = {} then
  reason := [convert(procname,string),"W1(emptyset) = emptyset"];
  return false;
 fi;

 a0 := A[1];
 u0 := x[a0] - y[a0];

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

 return true;
end;

`is_leq/W1` := NULL;

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

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

 return eval(x);
end;

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


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

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

$define xc_dimensions_2d
$include "mgga_x_2d_prhg07.mpl"

# we add here a threshold of 1e-10 for either tau or the Fermi curvature
prp10_f := (rs, z, x, u, t) -> -(X_FACTOR_2D_C*prhg07_v(prhg07_y(prhg07_C(x, u, t)))
  - (2*sqrt(2)/(3*Pi)) * sqrt(2*m_max(t - x^2/8, 1e-10)))*n_spin(rs, z)^(1/2):

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


text:
Q := [b*c+a, a*c-b, b*c-2*c^2+a-2, b*c-c^2+a-1, b*c+c^2+a+1, b*c+2*c^2+a+2, 2*
b*c-c^2+2*a-1, 2*b*c+c^2+2*a+1, a*c-2*c^2-b-2, a*c-c^2-b-1, a*c+c^2-b+1, a*c+2
*c^2-b+2, 2*a*c-c^2-2*b-1, 2*a*c+c^2-2*b+1, a*c-2*b*c-2*a-b, a*c-b*c-a-b, a*c+
b*c+a-b, a*c+2*b*c+2*a-b, 2*a*c-b*c-a-2*b, 2*a*c+b*c+a-2*b, a^2+b^2-3*c^2-3, a
^2+b^2-c^2-1, 3*a^2+3*b^2-c^2-1, a^2+b^2-2*b*c-2*a, a^2+b^2-b*c-a, a^2+b^2+b*c
+a, a^2+b^2+2*b*c+2*a, 2*a^2+2*b^2-b*c-a, 2*a^2+2*b^2+b*c+a, a^2-2*a*c+b^2+2*b
, a^2-a*c+b^2+b, a^2+a*c+b^2-b, a^2+2*a*c+b^2-2*b, 2*a^2-a*c+2*b^2+b, 2*a^2+a*
c+2*b^2-b, a*c-2*b*c-2*c^2-2*a-b-2, a*c-2*b*c-c^2-2*a-b-1, a*c-2*b*c+c^2-2*a-b
+1, a*c-2*b*c+2*c^2-2*a-b+2, a*c-b*c-2*c^2-a-b-2, a*c-b*c-c^2-a-b-1, a*c-b*c+c
^2-a-b+1, a*c-b*c+2*c^2-a-b+2, a*c+b*c-2*c^2+a-b-2, a*c+b*c-c^2+a-b-1, a*c+b*c
+c^2+a-b+1, a*c+b*c+2*c^2+a-b+2, a*c+2*b*c-2*c^2+2*a-b-2, a*c+2*b*c-c^2+2*a-b-\
1, a*c+2*b*c+c^2+2*a-b+1, a*c+2*b*c+2*c^2+2*a-b+2, 2*a*c-2*b*c-c^2-2*a-2*b-1,
2*a*c-2*b*c+c^2-2*a-2*b+1, 2*a*c-b*c-2*c^2-a-2*b-2, 2*a*c-b*c-c^2-a-2*b-1, 2*a
*c-b*c+c^2-a-2*b+1, 2*a*c-b*c+2*c^2-a-2*b+2, 2*a*c+b*c-2*c^2+a-2*b-2, 2*a*c+b*
c-c^2+a-2*b-1, 2*a*c+b*c+c^2+a-2*b+1, 2*a*c+b*c+2*c^2+a-2*b+2, 2*a*c+2*b*c-c^2
+2*a-2*b-1, 2*a*c+2*b*c+c^2+2*a-2*b+1, a^2+b^2-4*b*c-3*c^2-4*a-3, a^2+b^2-4*b*
c-c^2-4*a-1, a^2+b^2-4*b*c+c^2-4*a+1, a^2+b^2-4*b*c+3*c^2-4*a+3, a^2+b^2-2*b*c
-3*c^2-2*a-3, a^2+b^2-2*b*c-c^2-2*a-1, a^2+b^2-b*c-c^2-a-1, a^2+b^2+b*c-c^2+a-\
1, a^2+b^2+2*b*c-3*c^2+2*a-3, a^2+b^2+2*b*c-c^2+2*a-1, a^2+b^2+4*b*c-3*c^2+4*a
-3, a^2+b^2+4*b*c-c^2+4*a-1, a^2+b^2+4*b*c+c^2+4*a+1, a^2+b^2+4*b*c+3*c^2+4*a+
3, 3*a^2+3*b^2-4*b*c-c^2-4*a-1, 3*a^2+3*b^2-4*b*c+c^2-4*a+1, 3*a^2+3*b^2-2*b*c
-c^2-2*a-1, 3*a^2+3*b^2+2*b*c-c^2+2*a-1, 3*a^2+3*b^2+4*b*c-c^2+4*a-1, 3*a^2+3*
b^2+4*b*c+c^2+4*a+1, a^2-4*a*c+b^2-3*c^2+4*b-3, a^2-4*a*c+b^2-c^2+4*b-1, a^2-4
*a*c+b^2+c^2+4*b+1, a^2-4*a*c+b^2+3*c^2+4*b+3, a^2-2*a*c+b^2-3*c^2+2*b-3, a^2-\
2*a*c+b^2-c^2+2*b-1, a^2-a*c+b^2-c^2+b-1, a^2+a*c+b^2-c^2-b-1, a^2+2*a*c+b^2-3
*c^2-2*b-3, a^2+2*a*c+b^2-c^2-2*b-1, a^2+4*a*c+b^2-3*c^2-4*b-3, a^2+4*a*c+b^2-
c^2-4*b-1, a^2+4*a*c+b^2+c^2-4*b+1, a^2+4*a*c+b^2+3*c^2-4*b+3, 3*a^2-4*a*c+3*b
^2-c^2+4*b-1, 3*a^2-4*a*c+3*b^2+c^2+4*b+1, 3*a^2-2*a*c+3*b^2-c^2+2*b-1, 3*a^2+
2*a*c+3*b^2-c^2-2*b-1, 3*a^2+4*a*c+3*b^2-c^2-4*b-1, 3*a^2+4*a*c+3*b^2+c^2-4*b+
1, a^2-2*a*c+b^2-2*b*c-2*a+2*b, a^2-2*a*c+b^2-b*c-a+2*b, a^2-2*a*c+b^2+b*c+a+2
*b, a^2-2*a*c+b^2+2*b*c+2*a+2*b, a^2-a*c+b^2-2*b*c-2*a+b, a^2-a*c+b^2-b*c-a+b,
a^2-a*c+b^2+b*c+a+b, a^2-a*c+b^2+2*b*c+2*a+b, a^2+a*c+b^2-2*b*c-2*a-b, a^2+a*c
+b^2-b*c-a-b, a^2+a*c+b^2+b*c+a-b, a^2+a*c+b^2+2*b*c+2*a-b, a^2+2*a*c+b^2-2*b*
c-2*a-2*b, a^2+2*a*c+b^2-b*c-a-2*b, a^2+2*a*c+b^2+b*c+a-2*b, a^2+2*a*c+b^2+2*b
*c+2*a-2*b, 2*a^2-2*a*c+2*b^2-b*c-a+2*b, 2*a^2-2*a*c+2*b^2+b*c+a+2*b, 2*a^2-a*
c+2*b^2-2*b*c-2*a+b, 2*a^2-a*c+2*b^2-b*c-a+b, 2*a^2-a*c+2*b^2+b*c+a+b, 2*a^2-a
*c+2*b^2+2*b*c+2*a+b, 2*a^2+a*c+2*b^2-2*b*c-2*a-b, 2*a^2+a*c+2*b^2-b*c-a-b, 2*
a^2+a*c+2*b^2+b*c+a-b, 2*a^2+a*c+2*b^2+2*b*c+2*a-b, 2*a^2+2*a*c+2*b^2-b*c-a-2*
b, 2*a^2+2*a*c+2*b^2+b*c+a-2*b, a^2-4*a*c+b^2-4*b*c-3*c^2-4*a+4*b-3, a^2-4*a*c
+b^2-4*b*c-c^2-4*a+4*b-1, a^2-4*a*c+b^2-4*b*c+c^2-4*a+4*b+1, a^2-4*a*c+b^2-4*b
*c+3*c^2-4*a+4*b+3, a^2-4*a*c+b^2-4*b*c+5*c^2-4*a+4*b+5, a^2-4*a*c+b^2-2*b*c-3
*c^2-2*a+4*b-3, a^2-4*a*c+b^2-2*b*c-c^2-2*a+4*b-1, a^2-4*a*c+b^2-2*b*c+c^2-2*a
+4*b+1, a^2-4*a*c+b^2-2*b*c+3*c^2-2*a+4*b+3, a^2-4*a*c+b^2+2*b*c-3*c^2+2*a+4*b
-3, a^2-4*a*c+b^2+2*b*c-c^2+2*a+4*b-1, a^2-4*a*c+b^2+2*b*c+c^2+2*a+4*b+1, a^2-\
4*a*c+b^2+2*b*c+3*c^2+2*a+4*b+3, a^2-4*a*c+b^2+4*b*c-3*c^2+4*a+4*b-3, a^2-4*a*
c+b^2+4*b*c-c^2+4*a+4*b-1, a^2-4*a*c+b^2+4*b*c+c^2+4*a+4*b+1, a^2-4*a*c+b^2+4*
b*c+3*c^2+4*a+4*b+3, a^2-4*a*c+b^2+4*b*c+5*c^2+4*a+4*b+5, a^2-2*a*c+b^2-4*b*c-\
3*c^2-4*a+2*b-3, a^2-2*a*c+b^2-4*b*c-c^2-4*a+2*b-1, a^2-2*a*c+b^2-4*b*c+c^2-4*
a+2*b+1, a^2-2*a*c+b^2-4*b*c+3*c^2-4*a+2*b+3, a^2-2*a*c+b^2-2*b*c-3*c^2-2*a+2*
b-3, a^2-2*a*c+b^2-2*b*c-c^2-2*a+2*b-1, a^2-2*a*c+b^2-2*b*c+c^2-2*a+2*b+1, a^2
-2*a*c+b^2-b*c-c^2-a+2*b-1, a^2-2*a*c+b^2-b*c+c^2-a+2*b+1, a^2-2*a*c+b^2+b*c-c
^2+a+2*b-1, a^2-2*a*c+b^2+b*c+c^2+a+2*b+1, a^2-2*a*c+b^2+2*b*c-3*c^2+2*a+2*b-3
, a^2-2*a*c+b^2+2*b*c-c^2+2*a+2*b-1, a^2-2*a*c+b^2+2*b*c+c^2+2*a+2*b+1, a^2-2*
a*c+b^2+4*b*c-3*c^2+4*a+2*b-3, a^2-2*a*c+b^2+4*b*c-c^2+4*a+2*b-1, a^2-2*a*c+b^
2+4*b*c+c^2+4*a+2*b+1, a^2-2*a*c+b^2+4*b*c+3*c^2+4*a+2*b+3, a^2-a*c+b^2-2*b*c-
c^2-2*a+b-1, a^2-a*c+b^2-2*b*c+c^2-2*a+b+1, a^2-a*c+b^2-b*c-c^2-a+b-1, a^2-a*c
+b^2+b*c-c^2+a+b-1, a^2-a*c+b^2+2*b*c-c^2+2*a+b-1, a^2-a*c+b^2+2*b*c+c^2+2*a+b
+1, a^2+a*c+b^2-2*b*c-c^2-2*a-b-1, a^2+a*c+b^2-2*b*c+c^2-2*a-b+1, a^2+a*c+b^2-
b*c-c^2-a-b-1, a^2+a*c+b^2+b*c-c^2+a-b-1, a^2+a*c+b^2+2*b*c-c^2+2*a-b-1, a^2+a
*c+b^2+2*b*c+c^2+2*a-b+1, a^2+2*a*c+b^2-4*b*c-3*c^2-4*a-2*b-3, a^2+2*a*c+b^2-4
*b*c-c^2-4*a-2*b-1, a^2+2*a*c+b^2-4*b*c+c^2-4*a-2*b+1, a^2+2*a*c+b^2-4*b*c+3*c
^2-4*a-2*b+3, a^2+2*a*c+b^2-2*b*c-3*c^2-2*a-2*b-3, a^2+2*a*c+b^2-2*b*c-c^2-2*a
-2*b-1, a^2+2*a*c+b^2-2*b*c+c^2-2*a-2*b+1, a^2+2*a*c+b^2-b*c-c^2-a-2*b-1, a^2+
2*a*c+b^2-b*c+c^2-a-2*b+1, a^2+2*a*c+b^2+b*c-c^2+a-2*b-1, a^2+2*a*c+b^2+b*c+c^
2+a-2*b+1, a^2+2*a*c+b^2+2*b*c-3*c^2+2*a-2*b-3, a^2+2*a*c+b^2+2*b*c-c^2+2*a-2*
b-1, a^2+2*a*c+b^2+2*b*c+c^2+2*a-2*b+1, a^2+2*a*c+b^2+4*b*c-3*c^2+4*a-2*b-3, a
^2+2*a*c+b^2+4*b*c-c^2+4*a-2*b-1, a^2+2*a*c+b^2+4*b*c+c^2+4*a-2*b+1, a^2+2*a*c
+b^2+4*b*c+3*c^2+4*a-2*b+3, a^2+4*a*c+b^2-4*b*c-3*c^2-4*a-4*b-3, a^2+4*a*c+b^2
-4*b*c-c^2-4*a-4*b-1, a^2+4*a*c+b^2-4*b*c+c^2-4*a-4*b+1, a^2+4*a*c+b^2-4*b*c+3
*c^2-4*a-4*b+3, a^2+4*a*c+b^2-4*b*c+5*c^2-4*a-4*b+5, a^2+4*a*c+b^2-2*b*c-3*c^2
-2*a-4*b-3, a^2+4*a*c+b^2-2*b*c-c^2-2*a-4*b-1, a^2+4*a*c+b^2-2*b*c+c^2-2*a-4*b
+1, a^2+4*a*c+b^2-2*b*c+3*c^2-2*a-4*b+3, a^2+4*a*c+b^2+2*b*c-3*c^2+2*a-4*b-3,
a^2+4*a*c+b^2+2*b*c-c^2+2*a-4*b-1, a^2+4*a*c+b^2+2*b*c+c^2+2*a-4*b+1, a^2+4*a*
c+b^2+2*b*c+3*c^2+2*a-4*b+3, a^2+4*a*c+b^2+4*b*c-3*c^2+4*a-4*b-3, a^2+4*a*c+b^
2+4*b*c-c^2+4*a-4*b-1, a^2+4*a*c+b^2+4*b*c+c^2+4*a-4*b+1, a^2+4*a*c+b^2+4*b*c+
3*c^2+4*a-4*b+3, a^2+4*a*c+b^2+4*b*c+5*c^2+4*a-4*b+5, 3*a^2-4*a*c+3*b^2-4*b*c-
c^2-4*a+4*b-1, 3*a^2-4*a*c+3*b^2-4*b*c+c^2-4*a+4*b+1, 3*a^2-4*a*c+3*b^2-2*b*c-
c^2-2*a+4*b-1, 3*a^2-4*a*c+3*b^2-2*b*c+c^2-2*a+4*b+1, 3*a^2-4*a*c+3*b^2+2*b*c-
c^2+2*a+4*b-1, 3*a^2-4*a*c+3*b^2+2*b*c+c^2+2*a+4*b+1, 3*a^2-4*a*c+3*b^2+4*b*c-
c^2+4*a+4*b-1, 3*a^2-4*a*c+3*b^2+4*b*c+c^2+4*a+4*b+1, 3*a^2-2*a*c+3*b^2-4*b*c-
c^2-4*a+2*b-1, 3*a^2-2*a*c+3*b^2-4*b*c+c^2-4*a+2*b+1, 3*a^2-2*a*c+3*b^2-2*b*c-
c^2-2*a+2*b-1, 3*a^2-2*a*c+3*b^2+2*b*c-c^2+2*a+2*b-1, 3*a^2-2*a*c+3*b^2+4*b*c-
c^2+4*a+2*b-1, 3*a^2-2*a*c+3*b^2+4*b*c+c^2+4*a+2*b+1, 3*a^2+2*a*c+3*b^2-4*b*c-
c^2-4*a-2*b-1, 3*a^2+2*a*c+3*b^2-4*b*c+c^2-4*a-2*b+1, 3*a^2+2*a*c+3*b^2-2*b*c-
c^2-2*a-2*b-1, 3*a^2+2*a*c+3*b^2+2*b*c-c^2+2*a-2*b-1, 3*a^2+2*a*c+3*b^2+4*b*c-
c^2+4*a-2*b-1, 3*a^2+2*a*c+3*b^2+4*b*c+c^2+4*a-2*b+1, 3*a^2+4*a*c+3*b^2-4*b*c-
c^2-4*a-4*b-1, 3*a^2+4*a*c+3*b^2-4*b*c+c^2-4*a-4*b+1, 3*a^2+4*a*c+3*b^2-2*b*c-
c^2-2*a-4*b-1, 3*a^2+4*a*c+3*b^2-2*b*c+c^2-2*a-4*b+1, 3*a^2+4*a*c+3*b^2+2*b*c-
c^2+2*a-4*b-1, 3*a^2+4*a*c+3*b^2+2*b*c+c^2+2*a-4*b+1, 3*a^2+4*a*c+3*b^2+4*b*c-
c^2+4*a-4*b-1, 3*a^2+4*a*c+3*b^2+4*b*c+c^2+4*a-4*b+1, 5*a^2-4*a*c+5*b^2-4*b*c+
c^2-4*a+4*b+1, 5*a^2-4*a*c+5*b^2+4*b*c+c^2+4*a+4*b+1, 5*a^2+4*a*c+5*b^2-4*b*c+
c^2-4*a-4*b+1, 5*a^2+4*a*c+5*b^2+4*b*c+c^2+4*a-4*b+1, -b*c+a, a*b+c, -2*b^2-b*
c+a-2, -b^2-2*b*c+2*a-1, -b^2-b*c+a-1, b^2-2*b*c+2*a+1, b^2-b*c+a+1, 2*b^2-b*c
+a+2, a*b-2*b*c+2*a+c, a*b-b*c+a+c, a*b+b*c-a+c, a*b+2*b*c-2*a+c, 2*a*b-b*c+a+
2*c, 2*a*b+b*c-a+2*c, a*b-2*b^2+c-2, a*b-b^2+c-1, a*b+b^2+c+1, a*b+2*b^2+c+2,
2*a*b-b^2+2*c-1, 2*a*b+b^2+2*c+1, a^2-2*b*c+c^2+2*a, a^2-b*c+c^2+a, a^2+b*c+c^
2-a, a^2+2*b*c+c^2-2*a, 2*a^2-b*c+2*c^2+a, 2*a^2+b*c+2*c^2-a, a^2-3*b^2+c^2-3,
a^2-b^2+c^2-1, 3*a^2-b^2+3*c^2-1, a^2-2*a*b+c^2-2*c, a^2-a*b+c^2-c, a^2+a*b+c^
2+c, a^2+2*a*b+c^2+2*c, 2*a^2-a*b+2*c^2-c, 2*a^2+a*b+2*c^2+c, a*b-2*b^2-2*b*c+
2*a+c-2, a*b-2*b^2-b*c+a+c-2, a*b-2*b^2+b*c-a+c-2, a*b-2*b^2+2*b*c-2*a+c-2, a*
b-b^2-2*b*c+2*a+c-1, a*b-b^2-b*c+a+c-1, a*b-b^2+b*c-a+c-1, a*b-b^2+2*b*c-2*a+c
-1, a*b+b^2-2*b*c+2*a+c+1, a*b+b^2-b*c+a+c+1, a*b+b^2+b*c-a+c+1, a*b+b^2+2*b*c
-2*a+c+1, a*b+2*b^2-2*b*c+2*a+c+2, a*b+2*b^2-b*c+a+c+2, a*b+2*b^2+b*c-a+c+2, a
*b+2*b^2+2*b*c-2*a+c+2, 2*a*b-2*b^2-b*c+a+2*c-2, 2*a*b-2*b^2+b*c-a+2*c-2, 2*a*
b-b^2-2*b*c+2*a+2*c-1, 2*a*b-b^2-b*c+a+2*c-1, 2*a*b-b^2+b*c-a+2*c-1, 2*a*b-b^2
+2*b*c-2*a+2*c-1, 2*a*b+b^2-2*b*c+2*a+2*c+1, 2*a*b+b^2-b*c+a+2*c+1, 2*a*b+b^2+
b*c-a+2*c+1, 2*a*b+b^2+2*b*c-2*a+2*c+1, 2*a*b+2*b^2-b*c+a+2*c+2, 2*a*b+2*b^2+b
*c-a+2*c+2, a^2-3*b^2-4*b*c+c^2+4*a-3, a^2-3*b^2-2*b*c+c^2+2*a-3, a^2-3*b^2+2*
b*c+c^2-2*a-3, a^2-3*b^2+4*b*c+c^2-4*a-3, a^2-b^2-4*b*c+c^2+4*a-1, a^2-b^2-2*b
*c+c^2+2*a-1, a^2-b^2-b*c+c^2+a-1, a^2-b^2+b*c+c^2-a-1, a^2-b^2+2*b*c+c^2-2*a-\
1, a^2-b^2+4*b*c+c^2-4*a-1, a^2+b^2-4*b*c+c^2+4*a+1, a^2+b^2+4*b*c+c^2-4*a+1,
a^2+3*b^2-4*b*c+c^2+4*a+3, a^2+3*b^2+4*b*c+c^2-4*a+3, 3*a^2-b^2-4*b*c+3*c^2+4*
a-1, 3*a^2-b^2-2*b*c+3*c^2+2*a-1, 3*a^2-b^2+2*b*c+3*c^2-2*a-1, 3*a^2-b^2+4*b*c
+3*c^2-4*a-1, 3*a^2+b^2-4*b*c+3*c^2+4*a+1, 3*a^2+b^2+4*b*c+3*c^2-4*a+1, a^2-2*
a*b-2*b*c+c^2+2*a-2*c, a^2-2*a*b-b*c+c^2+a-2*c, a^2-2*a*b+b*c+c^2-a-2*c, a^2-2
*a*b+2*b*c+c^2-2*a-2*c, a^2-a*b-2*b*c+c^2+2*a-c, a^2-a*b-b*c+c^2+a-c, a^2-a*b+
b*c+c^2-a-c, a^2-a*b+2*b*c+c^2-2*a-c, a^2+a*b-2*b*c+c^2+2*a+c, a^2+a*b-b*c+c^2
+a+c, a^2+a*b+b*c+c^2-a+c, a^2+a*b+2*b*c+c^2-2*a+c, a^2+2*a*b-2*b*c+c^2+2*a+2*
c, a^2+2*a*b-b*c+c^2+a+2*c, a^2+2*a*b+b*c+c^2-a+2*c, a^2+2*a*b+2*b*c+c^2-2*a+2
*c, 2*a^2-2*a*b-b*c+2*c^2+a-2*c, 2*a^2-2*a*b+b*c+2*c^2-a-2*c, 2*a^2-a*b-2*b*c+
2*c^2+2*a-c, 2*a^2-a*b-b*c+2*c^2+a-c, 2*a^2-a*b+b*c+2*c^2-a-c, 2*a^2-a*b+2*b*c
+2*c^2-2*a-c, 2*a^2+a*b-2*b*c+2*c^2+2*a+c, 2*a^2+a*b-b*c+2*c^2+a+c, 2*a^2+a*b+
b*c+2*c^2-a+c, 2*a^2+a*b+2*b*c+2*c^2-2*a+c, 2*a^2+2*a*b-b*c+2*c^2+a+2*c, 2*a^2
+2*a*b+b*c+2*c^2-a+2*c, a^2-4*a*b-3*b^2+c^2-4*c-3, a^2-4*a*b-b^2+c^2-4*c-1, a^
2-4*a*b+b^2+c^2-4*c+1, a^2-4*a*b+3*b^2+c^2-4*c+3, a^2-2*a*b-3*b^2+c^2-2*c-3, a
^2-2*a*b-b^2+c^2-2*c-1, a^2-a*b-b^2+c^2-c-1, a^2+a*b-b^2+c^2+c-1, a^2+2*a*b-3*
b^2+c^2+2*c-3, a^2+2*a*b-b^2+c^2+2*c-1, a^2+4*a*b-3*b^2+c^2+4*c-3, a^2+4*a*b-b
^2+c^2+4*c-1, a^2+4*a*b+b^2+c^2+4*c+1, a^2+4*a*b+3*b^2+c^2+4*c+3, 3*a^2-4*a*b-
b^2+3*c^2-4*c-1, 3*a^2-4*a*b+b^2+3*c^2-4*c+1, 3*a^2-2*a*b-b^2+3*c^2-2*c-1, 3*a
^2+2*a*b-b^2+3*c^2+2*c-1, 3*a^2+4*a*b-b^2+3*c^2+4*c-1, 3*a^2+4*a*b+b^2+3*c^2+4
*c+1, a^2-4*a*b-3*b^2-4*b*c+c^2+4*a-4*c-3, a^2-4*a*b-3*b^2-2*b*c+c^2+2*a-4*c-3
, a^2-4*a*b-3*b^2+2*b*c+c^2-2*a-4*c-3, a^2-4*a*b-3*b^2+4*b*c+c^2-4*a-4*c-3, a^
2-4*a*b-b^2-4*b*c+c^2+4*a-4*c-1, a^2-4*a*b-b^2-2*b*c+c^2+2*a-4*c-1, a^2-4*a*b-
b^2+2*b*c+c^2-2*a-4*c-1, a^2-4*a*b-b^2+4*b*c+c^2-4*a-4*c-1, a^2-4*a*b+b^2-4*b*
c+c^2+4*a-4*c+1, a^2-4*a*b+b^2-2*b*c+c^2+2*a-4*c+1, a^2-4*a*b+b^2+2*b*c+c^2-2*
a-4*c+1, a^2-4*a*b+b^2+4*b*c+c^2-4*a-4*c+1, a^2-4*a*b+3*b^2-4*b*c+c^2+4*a-4*c+
3, a^2-4*a*b+3*b^2-2*b*c+c^2+2*a-4*c+3, a^2-4*a*b+3*b^2+2*b*c+c^2-2*a-4*c+3, a
^2-4*a*b+3*b^2+4*b*c+c^2-4*a-4*c+3, a^2-4*a*b+5*b^2-4*b*c+c^2+4*a-4*c+5, a^2-4
*a*b+5*b^2+4*b*c+c^2-4*a-4*c+5, a^2-2*a*b-3*b^2-4*b*c+c^2+4*a-2*c-3, a^2-2*a*b
-3*b^2-2*b*c+c^2+2*a-2*c-3, a^2-2*a*b-3*b^2+2*b*c+c^2-2*a-2*c-3, a^2-2*a*b-3*b
^2+4*b*c+c^2-4*a-2*c-3, a^2-2*a*b-b^2-4*b*c+c^2+4*a-2*c-1, a^2-2*a*b-b^2-2*b*c
+c^2+2*a-2*c-1, a^2-2*a*b-b^2-b*c+c^2+a-2*c-1, a^2-2*a*b-b^2+b*c+c^2-a-2*c-1,
a^2-2*a*b-b^2+2*b*c+c^2-2*a-2*c-1, a^2-2*a*b-b^2+4*b*c+c^2-4*a-2*c-1, a^2-2*a*
b+b^2-4*b*c+c^2+4*a-2*c+1, a^2-2*a*b+b^2-2*b*c+c^2+2*a-2*c+1, a^2-2*a*b+b^2-b*
c+c^2+a-2*c+1, a^2-2*a*b+b^2+b*c+c^2-a-2*c+1, a^2-2*a*b+b^2+2*b*c+c^2-2*a-2*c+
1, a^2-2*a*b+b^2+4*b*c+c^2-4*a-2*c+1, a^2-2*a*b+3*b^2-4*b*c+c^2+4*a-2*c+3, a^2
-2*a*b+3*b^2+4*b*c+c^2-4*a-2*c+3, a^2-a*b-b^2-2*b*c+c^2+2*a-c-1, a^2-a*b-b^2-b
*c+c^2+a-c-1, a^2-a*b-b^2+b*c+c^2-a-c-1, a^2-a*b-b^2+2*b*c+c^2-2*a-c-1, a^2-a*
b+b^2-2*b*c+c^2+2*a-c+1, a^2-a*b+b^2+2*b*c+c^2-2*a-c+1, a^2+a*b-b^2-2*b*c+c^2+
2*a+c-1, a^2+a*b-b^2-b*c+c^2+a+c-1, a^2+a*b-b^2+b*c+c^2-a+c-1, a^2+a*b-b^2+2*b
*c+c^2-2*a+c-1, a^2+a*b+b^2-2*b*c+c^2+2*a+c+1, a^2+a*b+b^2+2*b*c+c^2-2*a+c+1,
a^2+2*a*b-3*b^2-4*b*c+c^2+4*a+2*c-3, a^2+2*a*b-3*b^2-2*b*c+c^2+2*a+2*c-3, a^2+
2*a*b-3*b^2+2*b*c+c^2-2*a+2*c-3, a^2+2*a*b-3*b^2+4*b*c+c^2-4*a+2*c-3, a^2+2*a*
b-b^2-4*b*c+c^2+4*a+2*c-1, a^2+2*a*b-b^2-2*b*c+c^2+2*a+2*c-1, a^2+2*a*b-b^2-b*
c+c^2+a+2*c-1, a^2+2*a*b-b^2+b*c+c^2-a+2*c-1, a^2+2*a*b-b^2+2*b*c+c^2-2*a+2*c-\
1, a^2+2*a*b-b^2+4*b*c+c^2-4*a+2*c-1, a^2+2*a*b+b^2-4*b*c+c^2+4*a+2*c+1, a^2+2
*a*b+b^2-2*b*c+c^2+2*a+2*c+1, a^2+2*a*b+b^2-b*c+c^2+a+2*c+1, a^2+2*a*b+b^2+b*c
+c^2-a+2*c+1, a^2+2*a*b+b^2+2*b*c+c^2-2*a+2*c+1, a^2+2*a*b+b^2+4*b*c+c^2-4*a+2
*c+1, a^2+2*a*b+3*b^2-4*b*c+c^2+4*a+2*c+3, a^2+2*a*b+3*b^2+4*b*c+c^2-4*a+2*c+3
, a^2+4*a*b-3*b^2-4*b*c+c^2+4*a+4*c-3, a^2+4*a*b-3*b^2-2*b*c+c^2+2*a+4*c-3, a^
2+4*a*b-3*b^2+2*b*c+c^2-2*a+4*c-3, a^2+4*a*b-3*b^2+4*b*c+c^2-4*a+4*c-3, a^2+4*
a*b-b^2-4*b*c+c^2+4*a+4*c-1, a^2+4*a*b-b^2-2*b*c+c^2+2*a+4*c-1, a^2+4*a*b-b^2+
2*b*c+c^2-2*a+4*c-1, a^2+4*a*b-b^2+4*b*c+c^2-4*a+4*c-1, a^2+4*a*b+b^2-4*b*c+c^
2+4*a+4*c+1, a^2+4*a*b+b^2-2*b*c+c^2+2*a+4*c+1, a^2+4*a*b+b^2+2*b*c+c^2-2*a+4*
c+1, a^2+4*a*b+b^2+4*b*c+c^2-4*a+4*c+1, a^2+4*a*b+3*b^2-4*b*c+c^2+4*a+4*c+3, a
^2+4*a*b+3*b^2-2*b*c+c^2+2*a+4*c+3, a^2+4*a*b+3*b^2+2*b*c+c^2-2*a+4*c+3, a^2+4
*a*b+3*b^2+4*b*c+c^2-4*a+4*c+3, a^2+4*a*b+5*b^2-4*b*c+c^2+4*a+4*c+5, a^2+4*a*b
+5*b^2+4*b*c+c^2-4*a+4*c+5, 3*a^2-4*a*b-b^2-4*b*c+3*c^2+4*a-4*c-1, 3*a^2-4*a*b
-b^2-2*b*c+3*c^2+2*a-4*c-1, 3*a^2-4*a*b-b^2+2*b*c+3*c^2-2*a-4*c-1, 3*a^2-4*a*b
-b^2+4*b*c+3*c^2-4*a-4*c-1, 3*a^2-4*a*b+b^2-4*b*c+3*c^2+4*a-4*c+1, 3*a^2-4*a*b
+b^2-2*b*c+3*c^2+2*a-4*c+1, 3*a^2-4*a*b+b^2+2*b*c+3*c^2-2*a-4*c+1, 3*a^2-4*a*b
+b^2+4*b*c+3*c^2-4*a-4*c+1, 3*a^2-2*a*b-b^2-4*b*c+3*c^2+4*a-2*c-1, 3*a^2-2*a*b
-b^2-2*b*c+3*c^2+2*a-2*c-1, 3*a^2-2*a*b-b^2+2*b*c+3*c^2-2*a-2*c-1, 3*a^2-2*a*b
-b^2+4*b*c+3*c^2-4*a-2*c-1, 3*a^2-2*a*b+b^2-4*b*c+3*c^2+4*a-2*c+1, 3*a^2-2*a*b
+b^2+4*b*c+3*c^2-4*a-2*c+1, 3*a^2+2*a*b-b^2-4*b*c+3*c^2+4*a+2*c-1, 3*a^2+2*a*b
-b^2-2*b*c+3*c^2+2*a+2*c-1, 3*a^2+2*a*b-b^2+2*b*c+3*c^2-2*a+2*c-1, 3*a^2+2*a*b
-b^2+4*b*c+3*c^2-4*a+2*c-1, 3*a^2+2*a*b+b^2-4*b*c+3*c^2+4*a+2*c+1, 3*a^2+2*a*b
+b^2+4*b*c+3*c^2-4*a+2*c+1, 3*a^2+4*a*b-b^2-4*b*c+3*c^2+4*a+4*c-1, 3*a^2+4*a*b
-b^2-2*b*c+3*c^2+2*a+4*c-1, 3*a^2+4*a*b-b^2+2*b*c+3*c^2-2*a+4*c-1, 3*a^2+4*a*b
-b^2+4*b*c+3*c^2-4*a+4*c-1, 3*a^2+4*a*b+b^2-4*b*c+3*c^2+4*a+4*c+1, 3*a^2+4*a*b
+b^2-2*b*c+3*c^2+2*a+4*c+1, 3*a^2+4*a*b+b^2+2*b*c+3*c^2-2*a+4*c+1, 3*a^2+4*a*b
+b^2+4*b*c+3*c^2-4*a+4*c+1, 5*a^2-4*a*b+b^2-4*b*c+5*c^2+4*a-4*c+1, 5*a^2-4*a*b
+b^2+4*b*c+5*c^2-4*a-4*c+1, 5*a^2+4*a*b+b^2-4*b*c+5*c^2+4*a+4*c+1, 5*a^2+4*a*b
+b^2+4*b*c+5*c^2-4*a+4*c+1, a*c+b, a*b-c, a*c-2*b^2-2*c^2+b, a*c-b^2-c^2+b, a*
c+b^2+c^2+b, a*c+2*b^2+2*c^2+b, 2*a*c-b^2-c^2+2*b, 2*a*c+b^2+c^2+2*b, a*b-2*b^
2-2*c^2-c, a*b-b^2-c^2-c, a*b+b^2+c^2-c, a*b+2*b^2+2*c^2-c, 2*a*b-b^2-c^2-2*c,
2*a*b+b^2+c^2-2*c, a*b-2*a*c-2*b-c, a*b-a*c-b-c, a*b+a*c+b-c, a*b+2*a*c+2*b-c,
2*a*b-a*c-b-2*c, 2*a*b+a*c+b-2*c, a^2-3*b^2-3*c^2+1, a^2-b^2-c^2+1, 3*a^2-b^2-
c^2+3, a^2-2*a*c-2*b+1, a^2-a*c-b+1, a^2+a*c+b+1, a^2+2*a*c+2*b+1, 2*a^2-a*c-b
+2, 2*a^2+a*c+b+2, a^2-2*a*b+2*c+1, a^2-a*b+c+1, a^2+a*b-c+1, a^2+2*a*b-2*c+1,
2*a^2-a*b+c+2, 2*a^2+a*b-c+2, a*b-2*a*c-2*b^2-2*c^2-2*b-c, a*b-2*a*c-b^2-c^2-2
*b-c, a*b-2*a*c+b^2+c^2-2*b-c, a*b-2*a*c+2*b^2+2*c^2-2*b-c, a*b-a*c-2*b^2-2*c^
2-b-c, a*b-a*c-b^2-c^2-b-c, a*b-a*c+b^2+c^2-b-c, a*b-a*c+2*b^2+2*c^2-b-c, a*b+
a*c-2*b^2-2*c^2+b-c, a*b+a*c-b^2-c^2+b-c, a*b+a*c+b^2+c^2+b-c, a*b+a*c+2*b^2+2
*c^2+b-c, a*b+2*a*c-2*b^2-2*c^2+2*b-c, a*b+2*a*c-b^2-c^2+2*b-c, a*b+2*a*c+b^2+
c^2+2*b-c, a*b+2*a*c+2*b^2+2*c^2+2*b-c, 2*a*b-2*a*c-b^2-c^2-2*b-2*c, 2*a*b-2*a
*c+b^2+c^2-2*b-2*c, 2*a*b-a*c-2*b^2-2*c^2-b-2*c, 2*a*b-a*c-b^2-c^2-b-2*c, 2*a*
b-a*c+b^2+c^2-b-2*c, 2*a*b-a*c+2*b^2+2*c^2-b-2*c, 2*a*b+a*c-2*b^2-2*c^2+b-2*c,
2*a*b+a*c-b^2-c^2+b-2*c, 2*a*b+a*c+b^2+c^2+b-2*c, 2*a*b+a*c+2*b^2+2*c^2+b-2*c,
2*a*b+2*a*c-b^2-c^2+2*b-2*c, 2*a*b+2*a*c+b^2+c^2+2*b-2*c, a^2-4*a*c-3*b^2-3*c^
2-4*b+1, a^2-4*a*c-b^2-c^2-4*b+1, a^2-4*a*c+b^2+c^2-4*b+1, a^2-4*a*c+3*b^2+3*c
^2-4*b+1, a^2-2*a*c-3*b^2-3*c^2-2*b+1, a^2-2*a*c-b^2-c^2-2*b+1, a^2-a*c-b^2-c^
2-b+1, a^2+a*c-b^2-c^2+b+1, a^2+2*a*c-3*b^2-3*c^2+2*b+1, a^2+2*a*c-b^2-c^2+2*b
+1, a^2+4*a*c-3*b^2-3*c^2+4*b+1, a^2+4*a*c-b^2-c^2+4*b+1, a^2+4*a*c+b^2+c^2+4*
b+1, a^2+4*a*c+3*b^2+3*c^2+4*b+1, 3*a^2-4*a*c-b^2-c^2-4*b+3, 3*a^2-4*a*c+b^2+c
^2-4*b+3, 3*a^2-2*a*c-b^2-c^2-2*b+3, 3*a^2+2*a*c-b^2-c^2+2*b+3, 3*a^2+4*a*c-b^
2-c^2+4*b+3, 3*a^2+4*a*c+b^2+c^2+4*b+3, a^2-4*a*b-3*b^2-3*c^2+4*c+1, a^2-4*a*b
-b^2-c^2+4*c+1, a^2-4*a*b+b^2+c^2+4*c+1, a^2-4*a*b+3*b^2+3*c^2+4*c+1, a^2-2*a*
b-3*b^2-3*c^2+2*c+1, a^2-2*a*b-b^2-c^2+2*c+1, a^2-a*b-b^2-c^2+c+1, a^2+a*b-b^2
-c^2-c+1, a^2+2*a*b-3*b^2-3*c^2-2*c+1, a^2+2*a*b-b^2-c^2-2*c+1, a^2+4*a*b-3*b^
2-3*c^2-4*c+1, a^2+4*a*b-b^2-c^2-4*c+1, a^2+4*a*b+b^2+c^2-4*c+1, a^2+4*a*b+3*b
^2+3*c^2-4*c+1, 3*a^2-4*a*b-b^2-c^2+4*c+3, 3*a^2-4*a*b+b^2+c^2+4*c+3, 3*a^2-2*
a*b-b^2-c^2+2*c+3, 3*a^2+2*a*b-b^2-c^2-2*c+3, 3*a^2+4*a*b-b^2-c^2-4*c+3, 3*a^2
+4*a*b+b^2+c^2-4*c+3, a^2-2*a*b-2*a*c-2*b+2*c+1, a^2-2*a*b-a*c-b+2*c+1, a^2-2*
a*b+a*c+b+2*c+1, a^2-2*a*b+2*a*c+2*b+2*c+1, a^2-a*b-2*a*c-2*b+c+1, a^2-a*b-a*c
-b+c+1, a^2-a*b+a*c+b+c+1, a^2-a*b+2*a*c+2*b+c+1, a^2+a*b-2*a*c-2*b-c+1, a^2+a
*b-a*c-b-c+1, a^2+a*b+a*c+b-c+1, a^2+a*b+2*a*c+2*b-c+1, a^2+2*a*b-2*a*c-2*b-2*
c+1, a^2+2*a*b-a*c-b-2*c+1, a^2+2*a*b+a*c+b-2*c+1, a^2+2*a*b+2*a*c+2*b-2*c+1,
2*a^2-2*a*b-a*c-b+2*c+2, 2*a^2-2*a*b+a*c+b+2*c+2, 2*a^2-a*b-2*a*c-2*b+c+2, 2*a
^2-a*b-a*c-b+c+2, 2*a^2-a*b+a*c+b+c+2, 2*a^2-a*b+2*a*c+2*b+c+2, 2*a^2+a*b-2*a*
c-2*b-c+2, 2*a^2+a*b-a*c-b-c+2, 2*a^2+a*b+a*c+b-c+2, 2*a^2+a*b+2*a*c+2*b-c+2,
2*a^2+2*a*b-a*c-b-2*c+2, 2*a^2+2*a*b+a*c+b-2*c+2, a^2-4*a*b-4*a*c-3*b^2-3*c^2-\
4*b+4*c+1, a^2-4*a*b-4*a*c-b^2-c^2-4*b+4*c+1, a^2-4*a*b-4*a*c+b^2+c^2-4*b+4*c+
1, a^2-4*a*b-4*a*c+3*b^2+3*c^2-4*b+4*c+1, a^2-4*a*b-4*a*c+5*b^2+5*c^2-4*b+4*c+
1, a^2-4*a*b-2*a*c-3*b^2-3*c^2-2*b+4*c+1, a^2-4*a*b-2*a*c-b^2-c^2-2*b+4*c+1, a
^2-4*a*b-2*a*c+b^2+c^2-2*b+4*c+1, a^2-4*a*b-2*a*c+3*b^2+3*c^2-2*b+4*c+1, a^2-4
*a*b+2*a*c-3*b^2-3*c^2+2*b+4*c+1, a^2-4*a*b+2*a*c-b^2-c^2+2*b+4*c+1, a^2-4*a*b
+2*a*c+b^2+c^2+2*b+4*c+1, a^2-4*a*b+2*a*c+3*b^2+3*c^2+2*b+4*c+1, a^2-4*a*b+4*a
*c-3*b^2-3*c^2+4*b+4*c+1, a^2-4*a*b+4*a*c-b^2-c^2+4*b+4*c+1, a^2-4*a*b+4*a*c+b
^2+c^2+4*b+4*c+1, a^2-4*a*b+4*a*c+3*b^2+3*c^2+4*b+4*c+1, a^2-4*a*b+4*a*c+5*b^2
+5*c^2+4*b+4*c+1, a^2-2*a*b-4*a*c-3*b^2-3*c^2-4*b+2*c+1, a^2-2*a*b-4*a*c-b^2-c
^2-4*b+2*c+1, a^2-2*a*b-4*a*c+b^2+c^2-4*b+2*c+1, a^2-2*a*b-4*a*c+3*b^2+3*c^2-4
*b+2*c+1, a^2-2*a*b-2*a*c-3*b^2-3*c^2-2*b+2*c+1, a^2-2*a*b-2*a*c-b^2-c^2-2*b+2
*c+1, a^2-2*a*b-2*a*c+b^2+c^2-2*b+2*c+1, a^2-2*a*b-a*c-b^2-c^2-b+2*c+1, a^2-2*
a*b-a*c+b^2+c^2-b+2*c+1, a^2-2*a*b+a*c-b^2-c^2+b+2*c+1, a^2-2*a*b+a*c+b^2+c^2+
b+2*c+1, a^2-2*a*b+2*a*c-3*b^2-3*c^2+2*b+2*c+1, a^2-2*a*b+2*a*c-b^2-c^2+2*b+2*
c+1, a^2-2*a*b+2*a*c+b^2+c^2+2*b+2*c+1, a^2-2*a*b+4*a*c-3*b^2-3*c^2+4*b+2*c+1,
a^2-2*a*b+4*a*c-b^2-c^2+4*b+2*c+1, a^2-2*a*b+4*a*c+b^2+c^2+4*b+2*c+1, a^2-2*a*
b+4*a*c+3*b^2+3*c^2+4*b+2*c+1, a^2-a*b-2*a*c-b^2-c^2-2*b+c+1, a^2-a*b-2*a*c+b^
2+c^2-2*b+c+1, a^2-a*b-a*c-b^2-c^2-b+c+1, a^2-a*b+a*c-b^2-c^2+b+c+1, a^2-a*b+2
*a*c-b^2-c^2+2*b+c+1, a^2-a*b+2*a*c+b^2+c^2+2*b+c+1, a^2+a*b-2*a*c-b^2-c^2-2*b
-c+1, a^2+a*b-2*a*c+b^2+c^2-2*b-c+1, a^2+a*b-a*c-b^2-c^2-b-c+1, a^2+a*b+a*c-b^
2-c^2+b-c+1, a^2+a*b+2*a*c-b^2-c^2+2*b-c+1, a^2+a*b+2*a*c+b^2+c^2+2*b-c+1, a^2
+2*a*b-4*a*c-3*b^2-3*c^2-4*b-2*c+1, a^2+2*a*b-4*a*c-b^2-c^2-4*b-2*c+1, a^2+2*a
*b-4*a*c+b^2+c^2-4*b-2*c+1, a^2+2*a*b-4*a*c+3*b^2+3*c^2-4*b-2*c+1, a^2+2*a*b-2
*a*c-3*b^2-3*c^2-2*b-2*c+1, a^2+2*a*b-2*a*c-b^2-c^2-2*b-2*c+1, a^2+2*a*b-2*a*c
+b^2+c^2-2*b-2*c+1, a^2+2*a*b-a*c-b^2-c^2-b-2*c+1, a^2+2*a*b-a*c+b^2+c^2-b-2*c
+1, a^2+2*a*b+a*c-b^2-c^2+b-2*c+1, a^2+2*a*b+a*c+b^2+c^2+b-2*c+1, a^2+2*a*b+2*
a*c-3*b^2-3*c^2+2*b-2*c+1, a^2+2*a*b+2*a*c-b^2-c^2+2*b-2*c+1, a^2+2*a*b+2*a*c+
b^2+c^2+2*b-2*c+1, a^2+2*a*b+4*a*c-3*b^2-3*c^2+4*b-2*c+1, a^2+2*a*b+4*a*c-b^2-
c^2+4*b-2*c+1, a^2+2*a*b+4*a*c+b^2+c^2+4*b-2*c+1, a^2+2*a*b+4*a*c+3*b^2+3*c^2+
4*b-2*c+1, a^2+4*a*b-4*a*c-3*b^2-3*c^2-4*b-4*c+1, a^2+4*a*b-4*a*c-b^2-c^2-4*b-\
4*c+1, a^2+4*a*b-4*a*c+b^2+c^2-4*b-4*c+1, a^2+4*a*b-4*a*c+3*b^2+3*c^2-4*b-4*c+
1, a^2+4*a*b-4*a*c+5*b^2+5*c^2-4*b-4*c+1, a^2+4*a*b-2*a*c-3*b^2-3*c^2-2*b-4*c+
1, a^2+4*a*b-2*a*c-b^2-c^2-2*b-4*c+1, a^2+4*a*b-2*a*c+b^2+c^2-2*b-4*c+1, a^2+4
*a*b-2*a*c+3*b^2+3*c^2-2*b-4*c+1, a^2+4*a*b+2*a*c-3*b^2-3*c^2+2*b-4*c+1, a^2+4
*a*b+2*a*c-b^2-c^2+2*b-4*c+1, a^2+4*a*b+2*a*c+b^2+c^2+2*b-4*c+1, a^2+4*a*b+2*a
*c+3*b^2+3*c^2+2*b-4*c+1, a^2+4*a*b+4*a*c-3*b^2-3*c^2+4*b-4*c+1, a^2+4*a*b+4*a
*c-b^2-c^2+4*b-4*c+1, a^2+4*a*b+4*a*c+b^2+c^2+4*b-4*c+1, a^2+4*a*b+4*a*c+3*b^2
+3*c^2+4*b-4*c+1, a^2+4*a*b+4*a*c+5*b^2+5*c^2+4*b-4*c+1, 3*a^2-4*a*b-4*a*c-b^2
-c^2-4*b+4*c+3, 3*a^2-4*a*b-4*a*c+b^2+c^2-4*b+4*c+3, 3*a^2-4*a*b-2*a*c-b^2-c^2
-2*b+4*c+3, 3*a^2-4*a*b-2*a*c+b^2+c^2-2*b+4*c+3, 3*a^2-4*a*b+2*a*c-b^2-c^2+2*b
+4*c+3, 3*a^2-4*a*b+2*a*c+b^2+c^2+2*b+4*c+3, 3*a^2-4*a*b+4*a*c-b^2-c^2+4*b+4*c
+3, 3*a^2-4*a*b+4*a*c+b^2+c^2+4*b+4*c+3, 3*a^2-2*a*b-4*a*c-b^2-c^2-4*b+2*c+3,
3*a^2-2*a*b-4*a*c+b^2+c^2-4*b+2*c+3, 3*a^2-2*a*b-2*a*c-b^2-c^2-2*b+2*c+3, 3*a^
2-2*a*b+2*a*c-b^2-c^2+2*b+2*c+3, 3*a^2-2*a*b+4*a*c-b^2-c^2+4*b+2*c+3, 3*a^2-2*
a*b+4*a*c+b^2+c^2+4*b+2*c+3, 3*a^2+2*a*b-4*a*c-b^2-c^2-4*b-2*c+3, 3*a^2+2*a*b-\
4*a*c+b^2+c^2-4*b-2*c+3, 3*a^2+2*a*b-2*a*c-b^2-c^2-2*b-2*c+3, 3*a^2+2*a*b+2*a*
c-b^2-c^2+2*b-2*c+3, 3*a^2+2*a*b+4*a*c-b^2-c^2+4*b-2*c+3, 3*a^2+2*a*b+4*a*c+b^
2+c^2+4*b-2*c+3, 3*a^2+4*a*b-4*a*c-b^2-c^2-4*b-4*c+3, 3*a^2+4*a*b-4*a*c+b^2+c^
2-4*b-4*c+3, 3*a^2+4*a*b-2*a*c-b^2-c^2-2*b-4*c+3, 3*a^2+4*a*b-2*a*c+b^2+c^2-2*
b-4*c+3, 3*a^2+4*a*b+2*a*c-b^2-c^2+2*b-4*c+3, 3*a^2+4*a*b+2*a*c+b^2+c^2+2*b-4*
c+3, 3*a^2+4*a*b+4*a*c-b^2-c^2+4*b-4*c+3, 3*a^2+4*a*b+4*a*c+b^2+c^2+4*b-4*c+3,
5*a^2-4*a*b-4*a*c+b^2+c^2-4*b+4*c+5, 5*a^2-4*a*b+4*a*c+b^2+c^2+4*b+4*c+5, 5*a^
2+4*a*b-4*a*c+b^2+c^2-4*b-4*c+5, 5*a^2+4*a*b+4*a*c+b^2+c^2+4*b-4*c+5];


text:
# Proof by PDF that
#	if X <~ standardChiSq(n) then cX <~ gamma(k = n/2, t = 2*c)

PDF_chisq := (n, x) -> x^(n/2-1)*exp(-x/2)/(2^(n/2)*GAMMA(n/2));
PDF_cX := (n, c, x) -> PDF_chisq(n, x/c)/c;
PDF_sub := (k, t, x) -> PDF_cX(2*k, t/2, x);
simplify(PDF_sub(k, t, x));

PDF_gamma := (k, t, x) -> x^(k-1)*exp(-x/t)/(GAMMA(k)*t^k);
simplify(PDF_gamma(k, t, x));

# NOTE: Doesn't evaluate true, but inspection of forms given by
#		the 2 calls to simplify prove it so
evalb(PDF_sub(k, t, x) = PDF_gamma(k, t, x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
#
# Functions related to the CMS model or Cherednik operators
#
###################
#list of functions#
###################
#		Hpoly:=proc(poly)
#		Dunkl_i_n:= proc(i, n, poly, Nin:=0, parameter:= beta)
#		H_n:= proc(n, poly, giveinm:=1)
#		I_k_n:= proc(k, n, poly, giveinm:=1)
#		EigenCMSn:= proc(spart,n, hn:=0, IorNot:=NULL, minN:=NULL)
#		lambdabar_i:= proc(composition, i, N)
#		FermDunkl_i_n:= proc(i,n, poly)
#		Itilde_k_n:= proc(k, n, poly, giveinm:=1)

HcmsPolyM:= proc(expr_in_m)
	return z_to_m(sort(expand(Hpoly(old_LinBaseConvert(expr_in_m, `z`)))));
end proc:
part_Hpoly:= proc(poly, i, nbvars)
	local j, term2, term3, terms, term23;
	terms:= [seq(k, k=(i+1)..nbvars)];
	#for j from i+1 to nbvars do
	#if j > i  then
	term2:=0; term3:=0;
		term2:=map(j->(((z||i+z||j)/(z||i-z||j)*(z||i*diff(
			poly,z||i)-z||j*diff(poly,z||j)))), terms);
		term3:=map(j->(((z||i*z||j)/((z||i-z||j)^2)*(
				poly-subs({theta||i=theta||j,
				theta||j=theta||i,phi||i=phi||j,
				phi||j=phi||i},poly)))), terms) ;
	#end if;
	#end do:
	term23:= term2 -2*term3;
	return sort(expand(add(k, k in term23)));
end proc:
Hpoly:=proc(poly)
	local H_on_poly, term1, term23, term3, i,j, Vars, nbvars, varsindices;
	#Hamiltonian as a direct generalization of K_ij
	#The code is pretty much self-explainatory
    	Vars:=giveVars(poly);
    	nbvars:= max(nops(Vars[1]), nops(Vars[2]), nops(Vars[3]));

	term1:=0; term23:=0; term3:=0;
	varsindices:= [seq(k, k=1..nbvars)];
	#term1:=term1+sort(expand(z||i*(diff(z||i*diff(poly,z||i),z||i))));
	term1:= map(x-> z||x*(diff(z||x*diff(poly, z||x), z||x)), varsindices);
	term1:= sort(expand(add(k, k in term1)));
	term23:= map(x-> part_Hpoly(poly,x,nbvars), varsindices);
	term23:= sort(expand(beta*add(k, k in term23)));
	H_on_poly:= sort(expand(term1+term23));
	#H_on_poly:= z_to_m(H_on_poly);
	return factor(H_on_poly);
end proc:

mH_n:= proc(n, poly_in_m)
	return H_n(n,old_LinBaseConvert(poly_in_m, 'z'), 1);
end proc:

H_n:= proc(n, poly, giveinm:=1)
    local Hpoly, N, Vars;
	# The hamiltonian tower in terms of the sum of the dunkl on every variables
    Vars:= giveVars(poly);
	N:= max(nops(Vars[1]), nops(Vars[2]), nops(Vars[3]));
	if N = 1 then N:= 2; end if;
    Hpoly:= simplify(add(Dunkl_i_n(i, n, poly, N), i=1..N));
    if giveinm = 1 then Hpoly:= z_to_m(Hpoly); end if;
    return Hpoly;
end proc:
#J_n:= proc(n, poly, giveninm:=1)
    #local Hpoly, N, Vars;
	# The hamiltonian tower in terms of the sum of the dunkl on every variables
    #Vars:= giveVars(poly);
	#N:= max(nops(Vars[1]), nops(Vars[2]), nops(Vars[3]));
	#if N = 1 then N:= 2; end if;
    #Hpoly:= add(Dunkl_i_n(i, n, phi||i*theta||i*diff(diff(poly, phi||i), theta||i), N), i=1..N);
    #if giveinm = 1 then Hpoly:= z_to_m(Hpoly); end if;
    #return Hpoly;
#end proc:
#H_n:= proc(n, poly, giveinm:=1)
	#local Hpoly, N, Vars;
	#Vars:= giveVars(poly);
	#N:= max(nops(Vars[1]), nops(Vars[2]), nops(Vars[3]));
	#if N = 1 then N:= 2; end if;
    	#Hpoly:= add(FermDunkl_i_n(i, n, poly), i=1..N);
    	#if giveinm = 1 then Hpoly:= z_to_m(Hpoly); end if;
    	#return Hpoly;
#end proc:
project_var_sector:= proc(poly_z,var, sector)
	local mbarbar, mbar, Ppoly, k,j, mubar, ind, thevar;
	mbarbar:= sector[4];
	mubar:= sector[3];
	#print(mubar);
	mbar := sector[2];
	Ppoly:= poly_z;
	if var = "pt" then
		for k from 1 to mbarbar do
			Ppoly:= sort(expand(phi||k*theta||k*diff(diff(Ppoly, phi||k), theta||k)));
		end do:
	elif var = "p" then
		for j from 1 to mbar do
			k:= j+ mbarbar;
			Ppoly:= sort(expand(phi||k*diff(Ppoly, phi||k)));
		end do:
	elif var = "t" then
		for j from 1 to mubar do
			k:= j + mbarbar+mbar;
			Ppoly:= sort(expand(theta||k*diff(Ppoly, theta||k)));
		end do:
	else return NULL;
	end if;
	return Ppoly;
end proc:

expI_n:= proc(smalli,n, poly_in_m)
	local sector, mbarbar, mbar, nbVars, Ipoly, theIndex, mubar, bigM;
	print(poly_in_m);
	sector:= givePolySector(poly_in_m);
	mbarbar:= sector[4];
	mbar:= sector[2];
	mubar:= sector[3];
	bigM:= [0,sector[4], sector[4]+sector[2], sector[4]+sector[3]+sector[2]];
    nbVars:= max(map(x-> nops(Flatten(x)), genN2spart(op(sector))));
	Ipoly:= old_LinBaseConvert(poly_in_m, 'z');
	Ipoly:= project_var_sector(Ipoly, "pt", sector);
	Ipoly:= project_var_sector(Ipoly, "p", sector);
	Ipoly:= project_var_sector(Ipoly, "t", sector);
	Ipoly:= Dunkl_i_n(bigM[smalli]+1,n,Ipoly, nbVars);
	Ipoly:= factor(Kw(Ipoly, nbVars));
	Ipoly:= z_to_m(Ipoly);
	return Ipoly;
end proc:
Ibarbar_n:= proc(n,poly_in_m)
	local sector, mbarbar, nbVars, Ipoly;
	sector:= givePolySector(poly_in_m);
	mbarbar:= sector[4];
	nbVars:= max(map(x-> nops(Flatten(x)), genN2spart(op(sector))));
	Ipoly:= old_LinBaseConvert(poly_in_m, 'z');
	if mbarbar = 0 then return 0; end if;
	Ipoly:= phi1*theta1*diff(diff(Ipoly, phi1), theta1);
	Ipoly:= Dunkl_i_n(1,n,Ipoly, nbVars);
	Ipoly:= factor(Kw(Ipoly, nbVars));
	Ipoly:= z_to_m(Ipoly);
	return Ipoly;
end proc:

Ibar_n:= proc(n, poly_in_m)
	local sector, mbarbar, mbar, nbVars, Ipoly, theIndex;
	sector:= givePolySector(poly_in_m);
	mbarbar:= sector[4];
	mbar := sector[2];
        nbVars:= max(map(x-> nops(Flatten(x)), genN2spart(op(sector))));
	Ipoly:= old_LinBaseConvert(poly_in_m, 'z');
	Ipoly:= project_var_sector(Ipoly, "pt", sector);
	if mbar = 0 then return 0; end if;
	theIndex:= mbarbar+1;
	Ipoly:= phi||theIndex*diff(Ipoly, phi||theIndex);
	Ipoly:= Dunkl_i_n(theIndex,n,Ipoly, nbVars);
	Ipoly:= factor(Kw(Ipoly, nbVars));
	Ipoly:= z_to_m(Ipoly);
	return Ipoly;
end proc:

Iubar_n:= proc(n, poly_in_m)
	local sector, mbarbar, mbar,mubar, nbVars, Ipoly, theIndex;
	sector:= givePolySector(poly_in_m);
	mbarbar:= sector[4];
	mbar := sector[2];
	mubar:= sector[3];
	if mubar = 0 then return 0; end if;
        nbVars:= max(map(x-> nops(Flatten(x)), genN2spart(op(sector))));
	Ipoly:= old_LinBaseConvert(poly_in_m, 'z');
	Ipoly:= project_var_sector(Ipoly, "pt", sector);
	Ipoly:= project_var_sector(Ipoly, "p", sector);
	theIndex:= mbarbar+mbar+1;
	Ipoly:= theta||theIndex*diff(Ipoly, theta||theIndex);
	Ipoly:= Dunkl_i_n(theIndex,n,Ipoly, nbVars);
	Ipoly:= factor(Kw(Ipoly, nbVars));
	Ipoly:= z_to_m(Ipoly);
	return Ipoly;
end proc:
testJn:= proc(n, poly_in_m)
	local poly_z, nbVars, Jpoly;
	poly_z:= old_LinBaseConvert(poly_in_m, 'z');
	nbVars:= giveNbVars(poly_z);
	Jpoly:= map(x-> theta||x*diff(poly_z, theta||x), [seq(k, k=1..nbVars)]);
	Jpoly:= map(x-> factor(Dunkl_i_n(x, n, Jpoly[x])), [seq(k,k=1..nbVars)]);
	Jpoly:= factor(add(k, k in Jpoly));
	return (Jpoly);
end proc:

VAP_I:= proc(spart, i, N)
	local spart_sector, M_i, pts, phis, thetas, zeds, Rspart, ellComp, deltaNell, term1, term2, term3, j,m,thesign;
	spart_sector:= giveSpartSector(spart);
	spart_sector:= [spart_sector[4], spart_sector[2], spart_sector[3]];
	M_i:= [0,spart_sector[1],spart_sector[1] + spart_sector[2],spart_sector[1] + spart_sector[2]+ spart_sector[3]];
	pts:= spart[1];
	phis:= spart[2];
	thetas:= spart[3];
	zeds:= spart[4];
	Rspart:= [op(Reverse(pts)),op(Reverse(phis)), op(Reverse(thetas)), op(Reverse(zeds))];
	ellComp:= nops(Rspart);
	deltaNell:= N-ellComp;
	if deltaNell<>0 then Rspart:= [op(Rspart), seq(0, k=1..deltaNell)]; end if;
	print(Rspart, "rspart");
	ellComp:= N;
	term1:= factorial(N- (spart_sector[1] + spart_sector[2] + spart_sector[3]));
	print(term1, "globalpref");
	term2:= factorial( spart_sector[1]) * factorial(spart_sector[2]) * factorial(spart_sector[3]);
	term2:= term2/factorial(spart_sector[i]);
	print(term2, "specificpref");
	term3:=
		add( factorial(M_i[i+1]-j)*( (M_i[i+1]-j)*lambdabar_i(Rspart,j, ellComp) - beta), j=(M_i[i]+1)..(M_i[i+1]-1))
		+ lambdabar_i(Rspart, M_i[i+1], ellComp);
	for j from M_i[i]+1 to M_i[i+1] - 1 do
		print("itt",j,"pref", factorial(M_i[i+1]-j));
		print("lambdabar", lambdabar_i(Rspart, j, ellComp));
		print("addedterm",factorial(M_i[i+1]-j)*( (M_i[i+1]-j)*lambdabar_i(Rspart,j, ellComp) - beta)) ;
		print((M_i[i+1]-j), hey);
	end do;
	print(M_i);
	print(lambdabar_i(Rspart, M_i[i+1], ellComp));
	return term1*term2*term3;
end proc:

EigenCMSn:= proc(spart,n, hn:=0, IorNot:=NULL, minN:=NULL)
    local DomMonomial, OpDomMonomial, I_1_1EigVal, I_2_1EigVal, H_2EigVal, EigenValues;
    # The function that returns the Eigen value for a given monomial. The code is probably broken.
    if minN <> NULL then DomMonomial:= m2n(spart, minN); else DomMonomial:= m2n(spart)  end if;

    I_1_1EigVal:=0;
    I_2_1EigVal:=0:
    H_2EigVal:= 0;

    OpDomMonomial:= I_k_n(1,n, DomMonomial, 2);
    if OpDomMonomial <> 0 then
        I_1_1EigVal:=  coeff(OpDomMonomial, m[op(spart)]);
    end if;

    OpDomMonomial:= I_k_n(6,n, DomMonomial, 2);
    if OpDomMonomial <> 0 then
        I_2_1EigVal:= coeff(OpDomMonomial, m[op(spart)]);
    end if;

    OpDomMonomial:= H_n(n+hn, DomMonomial);
    if OpDomMonomial <> 0 then
        H_2EigVal:= coeff(OpDomMonomial, m[op(spart)]);
    end if;
    if IorNot=NULL then
        EigenValues:= [I_1_1EigVal, I_2_1EigVal, H_2EigVal];
    else
        #EigenValues:= [I_1_1EigVal+I_2_1EigVal,I_1_1EigVal-I_2_1EigVal, H_2EigVal];
        EigenValues:= [I_1_1EigVal,I_2_1EigVal, H_2EigVal];
    end if;
    return EigenValues;
end proc:


text:
apply_linear := proc(f,T := rational)
 return
   proc(u)
    local c,v;

    if u = 0 then
     return 0;
    elif type(u,`+`) or type(u,list) or type(u,set) then
     return map(apply_linear(f,T),u);
    elif type(u,`*`) then
     c,v := selectremove(type,u,T);
     return c * f(v);
    elif type(u,T) then
     return u * f(1);
    else
     return f(u);
    fi;
   end:
end:

apply_linear_mod := (f,m::posint) -> proc(u)
 local c,v;

 if u = 0 then
  return 0;
 elif type(u,`+`) or type(u,list) or type(u,set) then
  return modp(map(apply_linear(f),u),m);
 elif type(u,`*`) then
  c,v := selectremove(type,u,integer);
  c := modp(c,m);
  if c = 0 then
   return 0;
  else
   return modp(expand(c * f(v)),m);
  fi;
 elif type(u,integer) then
  return modp(u * f(1),m);
 else
  return f(u);
 fi;
end:

apply_list_linear := (f) -> proc(u)
 local c,v;

 if u = 0 then
  return 0;
 elif type(u,`+`) or type(u,list) or type(u,set) then
  return map(apply_list_linear(f),u);
 elif type(u,`*`) then
  c,v := selectremove(type,u,rational);
  return c * f(v);
 elif type(u,rational) then
  return u * f(1);
 else
  return f(u);
 fi;
end:

apply_list_linear_mod := (f,m::posint) -> proc(u)
 local c,v;

 if u = 0 then
  return 0;
 elif type(u,`+`) or type(u,list) or type(u,set) then
  return modp(map(apply_list_linear(f),u),m);
 elif type(u,`*`) then
  c,v := selectremove(type,u,integer);
  c := modp(c,m);
  if c = 0 then
   return 0;
  else
   return modp(expand(c * f(v)),m);
  fi;
 elif type(u,rational) then
  return modp(u * f(1),m);
 else
  return f(u);
 fi;
end:

apply_bilinear := (f) -> proc(u,v)
 local cu,xu,cv,xv;

 if u = 0 or v = 0 then
  return 0;
 elif type(u,`+`) then
  return map(apply_bilinear(f),u,v);
 elif type(v,`+`) then
  return map2(apply_bilinear(f),u,v);
 fi;

 if type(u,`*`) then
  cu,xu := selectremove(type,u,rational);
 else
  cu,xu := 1,u;
 fi;

 if type(v,`*`) then
  cv,xv := selectremove(type,v,rational);
 else
  cv,xv := 1,v;
 fi;

 return cu * cv * f(xu,xv);
end:

apply_bilinear_mod := (f,m::posint) -> proc(u,v)
 local cu,xu,cv,xv,c;

 if u = 0 or v = 0 then
  return 0;
 elif type(u,`+`) then
  return modp(map(apply_bilinear(f),u,v),m);
 elif type(v,`+`) then
  return modp(map2(apply_bilinear(f),u,v),m);
 fi;

 if type(u,`*`) then
  cu,xu := selectremove(type,u,integer);
 else
  cu,xu := 1,u;
 fi;

 if type(v,`*`) then
  cv,xv := selectremove(type,v,integer);
 else
  cv,xv := 1,v;
 fi;

 c := modp(cu * cv,m);
 if c = 0 then
  return 0;
 else
  return modp(expand(c * f(xu,xv)),m);
 fi;
end:

apply_assoc := (f,e) -> proc()
 if nargs = 0 then
  return e;
 elif nargs = 1 then
  return args[1];
 elif args = 2 then
  return f(args[1],args[2]);
 else
  return f(args[1],apply_assoc(f,e)(args[2..-1]));
 fi;
end:

apply_linear_assoc := (f,e) -> proc()
 if nargs = 0 then
  return e;
 elif nargs = 1 then
  return args[1];
 elif args = 2 then
  return apply_bilinear(f)(args[1],args[2]);
 else
  return apply_bilinear(f)(args[1],apply_linear_assoc(f,e)(args[2..-1]));
 fi;
end:

apply_linear_assoc_mod := (f,e,m::posint) -> proc()
 if nargs = 0 then
  return e;
 elif nargs = 1 then
  return args[1];
 elif args = 2 then
  return apply_bilinear_mod(f,m)(args[1],args[2]);
 else
  return apply_bilinear_mod(f,m)(args[1],apply_linear_assoc_mod(f,e,m)(args[2..-1]));
 fi;
end:

apply_deg := proc(f,deg_type := rational)
 return
  proc(u)
   local c,v,d,e,i,n;

   if u = 0 then
    return 0;
   elif type(u,`+`) then
    d := map(apply_deg(f),{op(u)});
    if nops(d) = 1 then
     return d[1];
    else
     return FAIL;
    fi;
   elif type(u,`*`) then
    d := map(apply_deg(f),[op(u)]);
    if type(d[1],deg_type) then
     return `+`(op(d));
    elif type(d[1],list) then
     e := d[1];
     for i from 2 to nops(d) do e := e +~ d[i]; od;
     return e;
    else
     return FAIL;
    fi;
   elif type(u,`^`) and type(op(2,u),deg_type) then
    d := apply_deg(f)(op(1,u));
    n := op(2,u);
    if type(d,deg_type) then
     return n*d;
    elif type(d,list) then
     return n *~ d;
    else
     return FAIL;
    fi;
   else
    return f(u);
   fi;
  end:
end:

apply_max_deg := (f) -> proc(u)
 local c,v,d,e,i,n;

 if u = 0 then
  return 0;
 elif type(u,`+`) then
  return max(op(map(apply_max_deg(f),{op(u)})));
 elif type(u,`*`) then
  d := map(apply_max_deg(f),[op(u)]);
  if type(d[1],rational) then
   return `+`(op(d));
  elif type(d[1],list) then
   e := d[1];
   for i from 2 to nops(d) do e := e +~ d[i]; od;
   return e;
  else
   return FAIL;
  fi;
 elif type(u,`^`) and type(op(2,u),rational) then
  d := apply_max_deg(f)(op(1,u));
  n := op(2,u);
  if type(d,rational) then
   return n*d;
  elif type(d,list) then
   return n *~ d;
  else
   return FAIL;
  fi;
 else
  return f(u);
 fi;
end:

apply_leibniz := (f) -> proc(u)
 local n,v,uu,i;

 if type(u,`+`) or type(u,list) or type(u,set) then
  return map(apply_leibniz(f),u);
 elif type(u,`*`) then
  uu := [op(u)];
  return expand(add(`*`(op(subsop(i = apply_leibniz(f)(uu[i]),uu))),i=1..nops(uu)));
 elif type(u,`^`) then
  v,n := op(u);
  return n * v^(n-1) * apply_leibniz(f)(v);
 elif type(u,rational) then
  return 0;
 else
  return f(u);
 fi;
end:

apply_leibniz_mod := (f,m::posint) -> proc(u)
 local n,v,uu,i;

 if type(u,`+`) or type(u,list) or type(u,set) then
  return map(apply_leibniz(f),u);
 elif type(u,`*`) then
  uu := [op(u)];
  return modp(expand(add(`*`(op(subsop(i = apply_leibniz(f)(uu[i]),uu))),i=1..nops(uu))),m);
 elif type(u,`^`) then
  v,n := op(u);
  if modp(n,m) = 0 then
   return 0;
  else
   return modp(expand(n * v^(n-1) * apply_leibniz(f)(v)),m);
  fi;
 elif type(u,rational) then
  return 0;
 else
  return f(u);
 fi;
end:


text:
[tracks]
1=Dopekick, Dopekick_soundslikeabreeze_139.ogg
2=bigMovieScore, BigMovieScore_104.ogg

[options]
loop=false
shuffle=false

max_stars_repo_name: vardis/pano

max_issues_repo_name: vardis/pano

max_forks_repo_name: vardis/pano


text:
LogText:= proc(text, target)
	# must pass text as sprintf(text)!!!
	UpdateLog(text, target);
end proc:

LogExpression:= proc(text, target)
	# must pass text as sprintf("%q\n", text)!!!
	UpdateLog(text, target);
end proc:

#LogTextSIAN := ()->UpdateLog(sprintf(_passed), "LogAreaSIAN");
#LogExpressionSIAN := ()->UpdateLog(sprintf("%q\n", _passed), "SIAN");

#LogTextME := ()->UpdateLog(sprintf(_passed), "LogAreaME");
#LogExpressionME := ()->UpdateLog(sprintf("%q\n", _passed), "LogAreaME");

#LogTextSE := ()->UpdateLog(sprintf(_passed), "LogAreaSE");
#LogExpressionSE := ()->UpdateLog(sprintf("%q\n", _passed), "LogAreaSE");

UpdateLog := proc(s, target)
local logsofar;

	logsofar := DocumentTools:-GetProperty(target, value);
	if logsofar <> "" then
		logsofar := logsofar, "\n";
	end if;

	DocumentTools:-SetProperty(target, value, cat(logsofar,s), 'refresh');

end proc:

examples := table([
  	"Biohydrogenation" = [ "Taken from R. Munoz-Tamayo, L. Puillet, J.B. Daniel, D. Sauvant, O. Martin, M. Taghipoor, P. Blavy\n Review: To be or not to be an identifiable model. Is this a relevant question in animal science modelling?\ndoi.org/10.1017/S1751731117002774\nSystem (3) in Supplementary Material 2, initial conditions are assumed to be unknown",
  	[
  "dx4/dt = - k5 * x4 / (k6 + x4);\n",
  "dx5/dt = k5 * x4 / (k6 + x4) - k7 * x5/(k8 + x5 + x6);\n",
  "dx6/dt = k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10;\n",
  "dx7/dt = k9 * x6 * (k10 - x6) / k10;\n",
  "y1 = x4;\n",
  "y2 = x5"]],

  	"Chemical Reaction Network" = ["Taken from  Conradi, C., Shiu, A., Dynamics of post-translational modification systems: recent progress and future directions Eq. 3.4",
	[
  "dx1/dt = -k1 * x1 * x2 + k2 * x4 + k4 * x6;\n",
  "dx2/dt = k1 * x1 * x2 + k2 * x4 + k3 * x4;\n",
  "dx3/dt = k3 * x4 + k5 * x6 - k6 * x3 * x5;\n",
  "dx4/dt = k1 * x1 * x2 - k2 * x4 - k3 * x4;\n",
  "dx5/dt = k4 * x6 + k5 * x6 - k6 * x3 * x5;\n",
  "dx6/dt = -k4 * x6 - k5 * x6 + k6 * x3 * x5;\n",
  "y1 = x3;\n"
  "y2 = x2;\n" ]],

	"DAISY Ex. 3" = ["DAISY Example 3", [
  "dx1/dt = -1 * p1 * x1 + p2 * x2 + u(t);\n",
  "dx2/dt = p3 * x1 - p4 * x2 + p5 * x3;\n",
  "dx3/dt = p6 * x1 - p7 * x3;\n",
  "y1 = x1;\n"]],

	"DAISY_mamil3" = ["DAISY mamil 3",
	[
  "dx1/dt = -(a21 + a31 + a01) * x1 + a12 * x2 + a13 * x3 + u(t);\n",
  "dx2/dt = a21 * x1 - a12 * x2;\n",
  "dx3/dt = a31 * x1 - a13 * x3;\n",
  "y = x1"]],

	"DAISY_mamil4" = ["DAISY mamil 4", [
  "dx1/dt = -k01 * x1 + k12 * x2 + k13 * x3 + k14 * x4 - k21 * x1 - k31 * x1 - k41 * x1 + u(t);\n",
  "dx2/dt = -k12 * x2 + k21 * x1;\n",
  "dx3/dt = -k13 * x3 + k31 * x1;\n",
  "dx4/dt = -k14 * x4 + k41 * x1;\n",
  "y = x1"]],

	"HIV" = ["Example (with initial conditions assumed being unknown) from Section IV of 'DAISY: an Efficient Tool to Test Global Identifiability. Some Case Studies' by G. Bellu, M.P. Saccomani",
	[
  "dx1/dt = -b * x1 * x4 - d * x1 + s;\n",
  "dx2/dt = b * q1 * x1 * x4 - k1 * x2 - mu1 * x2;\n",
  "dx3/dt = b * q2 * x1 * x4 + k1 * x2 - mu2 * x3;\n",
  "dx4/dt = -c * x4 + k2 * x3;\n",
  "y1 = x1;\n",
  "y2 = x4"]],

	"HIV2" = ["The system is taken from Wodarz, D., Nowak, M.\nSpecific therapy regimes could lead to long-term immunological control of HIV\nhttps://doi.org/10.1073/pnas.96.25.14464\nPage 1",
	[
  "dx/dt = lm - d * x - beta * x * v;\n",
  "dy/dt = beta * x * v - a * y;\n",
  "dv/dt = k * y - u * v;\n",
  "dw/dt = c * x * y * w - c * q * y * w - b * w;\n",
  "dz/dt = c * q * y * w - h * z;\n",
  "y1 = w;\n",
  "y2 = z"]],

	"Lipolysis" = ["Taken from R. Munoz-Tamayo, L. Puillet, J.B. Daniel, D. Sauvant, O. Martin, M. Taghipoor, P. Blavy\nReview: To be or not to be an identifiable model. Is this a relevant question in animal science modelling?\ndoi.org/10.1017/S1751731117002774\nSystem (1) in Supplementary Material 2, initial conditions are assumed to be unknown\nbrought to the rational function form by introducing new state variable x5 = k1 e^(-k3 t)",
	[
  "dx1/dt = -x1 * x5 / (k2 + x1);\n",
  "dx2/dt = 2 * x1 * x5 / ((k2 + x1) * 3) - k4 * x2;\n",
  "dx3/dt = k4*(x2)/2 - k4*x3;\n",
  "dx4/dt = x1 * x5 / (3 * (k2 + x1)) + k4 * (x2)/2 + k4 * x3;\n",
  "dx5/dt = -k3 * x5;\n",
  "y1 = x1;\n",
  "y2 = x2 + x3;\n",
  "y3 = x4"]],

  	"LV" = ["Lotka-Volterra System",[
  	"dx1/dt = a*x1 - b*x1*x2;\n",
  	"dx2/dt = -c*x2 + d*x1*x2;\n",
  	"y = x1;\n"]],
	"OralGlucose" = ["Example (with initial conditions assumed being unknown) from Section III of 'DAISY: an Efficient Tool to Test Global Identifiability. Some Case Studies'\nby G. Bellu, M.P. Saccomani",
	[
  "dG/dt = -(p1 + X) * G + p1 * Gb + v * R;\n",
  "dX/dt = -p2 * X + p3 * (u(t) - Ib);\n",
  "dR/dt = k;\n",
  "dIb/dt = 0;\n",
  "dGb/dt = 0;\n",
  "y1 = G;\n",
  "y2 = Ib;\n",
  "y3 = Gb;\n"]],

	"SEIR" = ["Taken from N. Tuncer, T. Le\n'Structural and practical identifiability analysis of outbreak models'\nhttps://doi.org/10.1016/j.mbs.2018.02.004\nEquation (2.2) with prevalence observations",
[
  "dS/dt = -b * S * In / N;\n",
  "dE/dt = b * S * In / N - nu * E;\n",
  "dIn/dt = nu * E - a * In;\n",
  "dN/dt = 0;\n",
  "y1 = In;\n",
  "y2 = N;\n"]],

	"SEIR2" = ["Taken from N. Tuncer, T. Le\n'Structural and practical identifiability analysis of outbreak models'\nhttps://doi.org/10.1016/j.mbs.2018.02.004\nEquation (2.2) with cumulative incidence observations",
	[
  "dS/dt = -b * S * In / N;\n",
  "dE/dt = b * S * In / N - nu * E;\n",
  "dIn/dt = nu * E - a * In;\n",
  "dN/dt = 0;\n",
  "dCu/dt = nu * E;\n",
  "y1 = Cu;\n",
  "y2 = N"]],

	"SIR_R0" = ["SIR R0",[
  "dS/dt = -b * In * S;\n",
  "dIn/dt = b * In * S - g * In;\n",
  "dR/dt = g * In;\n",
  "daux/dt = 0;\n",
  "y1 = In;\n",
  "y2 = b / g + aux;"]],

	"SIRSForced" = ["Taken from Capistran M., Moreles M., Lara B.\n'Parameter Estimation of Some Epidemic Models.\n The Case of Recurrent Epidemics Caused by Respiratory Syncytial Virus'\ndoi.org/10.1007/s11538-009-9429-3\nEquations (7)-(11)",
[
  "ds/dt = mu - mu * s - b0 * (1 + b1 * x1) * i * s + g * r;\n",
  "di/dt = b0 * (1 + b1 * x1) * i * s - (nu + mu) * i;\n",
  "dr/dt = nu * i - (mu + g) * r;\n",
  "dx1/dt = -M * x2;\n",
  "dx2/dt = M * x1;\n",
  "y1 = i;\n",
  "y2 = r;\n"]],

	"SlowFast" = ["Taken from Vajda S., Rabitz H.\n'Identifiability and Distinguishability of First-Order Reaction Systems', p. 701\nWe added an extra output x_C",
	[
  "dxA/dt = -k1 * xA;\n",
  "dxB/dt = k1 * xA - k2 * xB;\n",
  "dxC/dt = k2 * xB;\n",
  "deA/dt = 0;\n",
  "deC/dt = 0;\n",
  "y1 = eA * xA + eB * xB + eC * xC;\n",
  "y2 = xC;\n",
  "y3 = eA;\n",
  "y4 = eC"]],

	"Treatment" = ["Taken from N. Tuncer, T. Le\nStructural and practical identifiability analysis of outbreak models'\nhttps://doi.org/10.1016/j.mbs.2018.02.004\nEquation (2.3) with observed treatment",
	[
  "dS/dt = -b * S * In / N - d * b * S * Tr / N;\n",
  "dIn/dt = b * S * In / N + d * b * S * Tr / N - (a + g) * In;\n",
  "dTr/dt = g * In - nu * Tr;\n",
  "dN/dt = 0;\n",
  "y1 = Tr;\n",
  "y2 = N"]],

	"Tumor" = ["Example (with initial conditions assumed being unknown) from Section 3 of\n'Examples of testing global identifiability of biological and biomedical models with the DAISY software'\nby M.P. Saccomani, S. Audoly, G. Bellu, L. D'Angio",
[ "dx1/dt = -(k3 + k7) * x1 + k4 * x2;\n",
  "dx2/dt = k3 * x1 - (k4 + a * k5 + b * d1 * k5) * x2 + k6 * x3 + k6 * x4 + k5 * x2 * x3 + k5 * x2 * x4;\n",
  "dx3/dt = a * k5 * x2 - k6 * x3 - k5 * x2 * x3;\n",
  "dx4/dt = b * d1 * k5 * x2 - k6 * x4 - k5 * x2 * x4;\n",
  "dx5/dt = k7 * x1;\n",
  "da/dt = 0;\n",
  "db/dt = 0;\n",
  "dd1/dt = 0;\n",
  "y1 = x5;\n",
  "y2 = a;\n",
  "y3 = b;\n",
  "y4 = d1;\n"]]
]):

# Setup

timed_SIAN:=proc(sigma, params_to_assess, p, output_targets_sian, count_solutions, char)
	local output, data, start, finish:
	start:= time():
	output := IdentifiabilityODE(sigma, params_to_assess, p, output_targets_sian, count_solutions, char):
	finish:= time():
	DocumentTools:-SetProperty(output_targets_sian[runningtime], value, convert(finish-start, string), 'refresh'): # time
	return  output:
end proc:

timed_Multi:=proc(model, simplified_generators, no_bound, simplify_bound, max_perms, output_targets_multi)
	local start, output, finish, data, bound, generators:
	start:=time():
	bound, generators := op(MultiExperimentIdentifiableFunctions(model, simplified_generators, no_bound, simplify_bound, max_perms, output_targets_multi)):
	finish:=time():
	DocumentTools:-SetProperty(output_targets_multi[runningtime], value, convert(finish-start, string), 'refresh'):
	return [bound, finish-start, generators]: #table([output=bound, runtime=finish-start]):
end proc:

timed_Single:=proc(model, output_targets_single)
	local start, output, finish, data:
	start:=time():
	output := SingleExperimentIdentifiableFunctions(model, output_targets_single):
	finish:=time():
	DocumentTools:-SetProperty(output_targets_single[runningtime], value, convert(finish-start, string), 'refresh'):
	return output:
end proc:

with(StringTools):

sigmaParser := proc(sigma)
	local states, state_eqs, outputs, output_eqs;
	if SearchText("diff", sigma)>0 then
		sigma := [map(x->parse(x), Split(sigma, ";"))]:
	else
		LogExpression(sprintf("%q \n", Split(sigma, ";")), "LogAreaSIAN"):
		sigma := Split(sigma, ";"):

		states := map(x->Trim(RegSubs("d([a-zA-Z0-9]+)/dt(.*)" = "\\1", x)), select(x->SearchText(x, "/dt")>0, sigma)):
		state_eqs := select(x->Has(x, "/dt"), sigma):

		outputs :=  map(x->Trim(Split(x, "=")[1]), select(x->not Has(x, "/dt"), sigma)):
		output_eqs := select(x->not SearchText(x, "/dt")>0, sigma):

		state_eqs := map(x->convert(subs({seq(parse(states[i])=parse(cat(states[i],"(t)")), i=1..nops(states))}, parse(x)), string),  state_eqs):
		state_eqs := map(x->parse(RegSubs("d([a-zA-Z0-9]+)/dt(.*)" = "diff(\\1(t), t)\\2", x)), state_eqs):

		output_eqs := map(x->parse(Subs({seq(outputs[i]=cat(outputs[i],"(t)"), i=1..nops(outputs))}, x)), output_eqs):
		output_eqs := map(x->subs({seq(parse(states[i])=parse(cat(states[i],"(t)")), i=1..nops(states))}, x),  output_eqs):

		sigma := [op(state_eqs), op(output_eqs)]:
	end if:
	return sigma;
end proc:

DocumentTools:-SetProperty("RunningTimeSingle", value, ""):
DocumentTools:-SetProperty("RunningTimeMulti", value, ""):
DocumentTools:-SetProperty("RunningTimeSIAN", value, ""):

DocumentTools:-SetProperty("run_system", enabled, true):
DocumentTools:-SetProperty("Meter_sian", visible, true):
DocumentTools:-SetProperty("Meter_sian", value, 0):
DocumentTools:-SetProperty("sigma", enabled, true):
DocumentTools:-SetProperty("p", value, "0.99"):
DocumentTools:-SetProperty("params", value, ""):
DocumentTools:-SetProperty("replicas", value, "1"):

DocumentTools:-SetProperty("p", enabled, true):
DocumentTools:-SetProperty("params", enabled, true):
DocumentTools:-SetProperty("replicas", enabled, true):
DocumentTools:-SetProperty("LogAreaSE", value, ""):
DocumentTools:-SetProperty("LogAreaSIAN", value, ""):
DocumentTools:-SetProperty("LogAreaME", value, ""):

DocumentTools:-SetProperty("printSolutions", enabled, true):
DocumentTools:-SetProperty("printSolutions", value, true):

DocumentTools:-SetProperty("GlobalParams1", expression, NULL):
DocumentTools:-SetProperty("LocalParams1", expression, NULL):
DocumentTools:-SetProperty("NoIDParams1", expression, NULL):

DocumentTools:-SetProperty("Parameters", value, ""):

DocumentTools:-SetProperty("Bound", expression, NULL):
DocumentTools:-SetProperty("MultiFunctions", expression, NULL):
DocumentTools:-SetProperty("SingleFunctions", expression, NULL):

DocumentTools:-SetProperty("use_char", enabled, true):
DocumentTools:-SetProperty("use_char", value, false):
DocumentTools:-SetProperty("p_label", enabled, true):

DocumentTools:-SetProperty("ComputeId", enabled, true):
DocumentTools:-SetProperty("ComputeId", value, false):
DocumentTools:-SetProperty("bypass", enabled, false):

DocumentTools:-SetProperty("bypass", value, true):
DocumentTools:-SetProperty("SimplifiedGen", enabled, false):
DocumentTools:-SetProperty("SimplifiedGen", value, true):
DocumentTools:-SetProperty("SkipSingle", enabled, false):
DocumentTools:-SetProperty("SkipSingle", value, false):
DocumentTools:-SetProperty("Refine", enabled, false):
DocumentTools:-SetProperty("NoBound", enabled, false):
DocumentTools:-SetProperty("NoBound", value, false):
DocumentTools:-SetProperty("UsingUpTo", enabled, false):
DocumentTools:-SetProperty("MaxPermutations", enabled, false):
DocumentTools:-SetProperty("Permutations", enabled, false):
DocumentTools:-SetProperty("RunSIAN", enabled, true):
DocumentTools:-SetProperty("RunSIAN", value, true):
DocumentTools:-SetProperty("being_refined", caption, "");
DocumentTools:-SetProperty("sigma", value, "dx1/dt = a*x1 + x2*b + u(t);\ndx2/dt = x2*c + x1;\ny=x2"):
DocumentTools:-SetProperty("example_box", value, "Custom"):
DocumentTools:-SetProperty(reference, value, ""):
DocumentTools:-SetProperty("LocalLabel" , caption, "Locally Identifiable Paramters");
DocumentTools:-SetProperty(TxtOutput, visible, false);
DocumentTools:-SetProperty(TxtOutput, value, "");
DocumentTools:-SetProperty(SaveOutputLabel, visible, false);

readyToSave:=false:
counter:=0:
exname:="Custom":

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
`is_element/E/Fbar` := (N::posint) -> (A::set,B::set) -> proc(tpqx)
 global reason;
 local t,p,q,x,nt,nx;

 if not(type(tpqx,list) and nops(tpqx) = 4) then
  reason := [convert(procname,string),"tpqx is not a list of length 4",tpqx];
  return false;
 fi;

 t,p,q,x := op(tpqx);

 if not (`is_element/real_functions`(A)(t) and
         `is_nonnegative/real_functions`(A)(t)) then
  reason := [convert(procname,string),"t is not a nonnegative real function on A",t,A];
  return false;
 fi;

 nt := `norm/real_functions`(A)(t)^2;

 if simplify(nt - 1) <> 0 then
  reason := [convert(procname,string),"t is not in normalised",t];
  return false;
 fi;

 if not(`is_element/RR`(p) and
        `is_element/RR`(q) and
        p >= 0 and q >= 0 and simplify(p^2+q^2 - 1) = 0) then
  reason := [convert(procname,string),"(p,q) is not in normalised",p,q];
  return false;
 fi;

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

 nx := `norm/prime_W`(N)(B)(x)^2;

 if simplify(nx - 1) <> 0 then
  reason := [convert(procname,string),"x is not in normalised",x];
  return false;
 fi;

 return true;
end:

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

`random_element/E/Fbar` := (N::posint) -> (A::set,B::set) -> proc()
 local t,p,q,x,n,t0,r,i,b;

 n := nops(A);
 t0 := map(abs,`random_element/sphere`(n-1)());
 t := table([seq(A[i]=t0[i],i=1..n)]):
 p,q := op(map(abs,`random_element/sphere`(1)()));
 r := 0;
 while r = 0 do
  x := `random_element/prime_W`(N)(B)();
  r := `norm/prime_W`(N)(B)(x);
 od;
 x := table([seq(b = simplify(x[b] /~ r),b in B)]);
 return [eval(t),p,q,eval(x)];
end:

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

`is_equal/E/Fbar` := (N::posint) -> (A::set,B::set) -> proc(tpqx1,tpqx2)
 local t1,p1,q1,x1,t2,p2,q2,x2;

 t1,p1,q1,x1 := op(tpqx1);
 t2,p2,q2,x2 := op(tpqx2);

 return `is_equal/real_functions`(A)(t1,t2) and
        `is_equal/RR`(p1,p2) and
        `is_equal/RR`(q1,q2) and
        `is_equal/prime_W`(N)(B)(x1,x2);
end:

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

`theta/E/D/Fbar` := (N::posint) -> (A::set,B::set) -> proc(tpqx)
 local t,p,q,x,u,v,a,b;

 t,p,q,x := op(tpqx);
 u := table();
 v := table();
 for a in A do u[a] := p * t[a]; od;
 for b in B do
  v[b] := q *~ x[b];
 od;

 return [u,v];
end;

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

`is_leq/E/Fbar` := NULL:
`list_elements/E/Fbar` := NULL:
`count_elements/E/Fbar` := NULL:


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

pos := parse~(StringTools:-Split(input,",")):
answer1 := min(seq(add(abs~(pos-~i)),i=min(pos)..max(pos)));
answer2 := min(seq(
    add(((pos[i]-j)^2+abs(pos[i]-j))/2, i=1..nops(pos)),
        j=min(pos)..max(pos)));


text:
 # a function returning unary(i + j) * -998
#type $SS <struct {@f1 f32}>
 func $foo (
  var %k <struct {@f1 f32}>,
  var %f f32,
  var %i i32,
   var %j i32)i64  {
   return (
       cvt i64 i32 (constval i32 12345))}

 # 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:
# Lagrange Formalis for Robot based on MDH frames
# Einleitung
# Berechnung des Lagrange-Formalismus
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# floatb -> floating base wird durch base twist (Geschwindigkeit der Basis) oder vollständige Orientierung (Euler-Winkel) berücksichtigt
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# lagrange -> Berechnung des Lagrange-Formalismus
# worldframe -> Berechnung basierend auf Energien aus Welt-KS (KS W)
# par12 -> Parametersatz 1 (Schwerpunkt als Parameter: SX,SY,SZ) oder Parametersatz 2 (1. und 2. Moment MX,MY,MZ,...)
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Sources
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# [KhalilDombre2002] Modeling, Identification and Control of Robots
# [Ortmaier2014] Vorlesungsskript Robotik I
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
# Einstellungen für Code-Export: Optimierungsgrad (2=höchster) und Aktivierung jedes Terms.
codegen_opt := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_Lagrange1":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
# Einstellungen für Term-Vereinfachungen: 0=keine, 1=dTdq_s, 2=auch dTdq_s, 3=auch dTdqDdt_s
if not assigned(simplify_options) or simplify_options(8)=-1 then # Standard-Einstellungen:
  # Die Vereinfachung ist ein Kompromiss aus Laufzeit des Programms und Güte der Ergebnisse
  # Bei zu komplizierten Termen stürzt Maple bei der Vereinfachung unvorhersehbar ab.
  # Annahme: Bei kinematischen Zwangsbedingungen werden die Terme zu kompliziert.
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # Annahme, dass Optimierung der Gravitations-Terme noch möglich ist
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(8): # achter Eintrag ist für Lagrange
end if:

# Kennung des Parametersatzes, für den die Dynamikfunktionen erstellt werden sollen. Muss im Repo und in der mpl-Datei auf 1 gelassen werden, da die folgende Zeile mit einem Skript verarbeitet wird.
codegen_dynpar := 1:
# Ergebnisse der Energie laden
if codegen_dynpar = 1 then
  read sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_worldframe_par1_maple.m", robot_name, base_method_name):
  read sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_worldframe_par1_maple.m", robot_name, base_method_name):
elif codegen_dynpar = 2 then
  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):
else
  printf("Energiefunktionen nur für Parametersatz 1 oder 2 definiert\n"):
end:
T := T:
U_grav := U_grav:
printf("%s. Generiere Terme des Lagrange-Formalismus für %s mit Parametersatz %d und %s\n", \
  FormatTime("%Y-%m-%d %H:%M:%S"), robot_name, codegen_dynpar, base_method_name):

# Lagrange Formalismus (mit Funktion)
tmp_t1:=time():
OutputLagrange := Lagrange1(T, U_grav, NQ):
dTdqDdt_s := OutputLagrange[1]:
dTdq_s := OutputLagrange[2]:
dUdq_s := OutputLagrange[3]:
printf("%s. Lagrange-Formalismus für %s mit Parametersatz %d und %s berechnet. Rechenzeit %1.1fs.\n", \
  FormatTime("%Y-%m-%d %H:%M:%S"), robot_name, codegen_dynpar, base_method_name, time()-tmp_t1):
# Terme vereinfachen: Einzeln und komponentenweise
if use_simplify>=1 then
  # Term dUdq_s
  tmp_t1:=time(): tmp_l1 := length(dUdq_s):
  printf("%s. Beginne Vereinfachung: Lagrange-Term dU/dq (Param. %d). Länge: %d.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l1):
  for i from 1 to NQ do
    dUdq_s[i,1] := simplify2(dUdq_s[i,1]):
  end do:
  tmp_t2:=time(): tmp_l2 := length(dUdq_s):
  printf("%s. Lagrange-Term dU/dq (Param. %d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
end if:
if use_simplify>=2 then # Die Terme können sehr kompliziert sein. Daher nicht immer vereinfachen.
  # Term dTdq_s
  for i from 1 to NQ do
    tmp_t1:=time(): tmp_l1 := length(dTdq_s[i,1]):
    printf("%s. Beginne mit Vereinfachung von Lagrange-Term dT/dq (Param. %d; Komponente %d/%d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, NQ, tmp_l1):
    dTdq_s[i,1] := simplify2(dTdq_s[i,1]):
    tmp_t2:=time(): tmp_l2 := length(dTdq_s[i,1]):
    printf("%s. Lagrange-Term dT/dq (Param. %d; Komponente %d/%d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, NQ, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end do:
end if:
if use_simplify>=3 then
  # Term dTdqDdt_s
  for i from 1 to NQ do
    tmp_t1:=time(): tmp_l1 := length(dTdqDdt_s[i,1]):
    printf("%s. Beginne mit Vereinfachung von Lagrange-Term d(dT/dqD)/dt (Param. %d; Komponente %d/%d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, NQ, tmp_l1):
    dTdqDdt_s[i,1] := simplify2(dTdqDdt_s[i,1]):
    tmp_t2:=time(): tmp_l2 := length(dTdqDdt_s[i,1]):
    printf("%s. Lagrange-Term d(dT/dqD)/dt (Param. %d; Komponente %d/%d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, NQ, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end do:
end if:
save dUdq_s, sprintf("../codeexport/%s/tmp/floatb_%s_lagrange_dUdq_s_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
save dTdq_s, sprintf("../codeexport/%s/tmp/floatb_%s_lagrange_dTdq_s_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
save dTdqDdt_s, sprintf("../codeexport/%s/tmp/floatb_%s_lagrange_dTdqDdt_s_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
printf("%s. Lagrange-Formalismus durchgeführt und Ergebnisse gespeichert (Dynamik-Parametersatz %d)\n", \
  FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar):


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

sa_a := 2.413:
sa_b := 0.348:

params_a_b      := 0.40:
params_a_c      := 1.59096:
params_a_e      := 1.537:
params_a_mu     := 0.21951:

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

(* Equation (8) *)

tpss_ff    := z -> 2:
tpss_kappa := (x, t) -> 2*Pi/(3*sqrt(5)) * \
          sqrt(sa_alpha(x, t) + 1)/sqrt(sa_a + log(sa_alpha(x, t) + sa_b)):

$include "tpss_x.mpl"

sa_a1  := (x, t) -> tpss_kappa(x, t)/(tpss_kappa(x, t) + tpss_fx(x, t)):
sa_f   := (x, u, t) -> 1 + tpss_kappa(x, t)*(1 - sa_a1(x, t)):

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


text:
# Proof by PDF that X <~ standardChiSq(2) == Y <~ exponential(2)

PDF_chisq := (n, x) -> x^(1/2*n-1)*exp(-1/2*x)/(2^(1/2*n))/GAMMA(1/2*n);
simplify(PDF_chisq(2,x));

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

evalb(PDF_chisq(2,x) = PDF_exponential(2,x));


text:
repeat process ACV_gels {

// Water
/scoreboard players set @e[type=item,name=ACV_CleansingGel] OnGround 1 {OnGround:1b}
/execute @e[type=item,name=ACV_CleansingGel,score_OnGround_min=1] ~ ~ ~ fill ~-1 ~ ~-1 ~1 ~ ~1 air 0 replace carpet 1
/execute @e[type=item,name=ACV_CleansingGel,score_OnGround_min=1] ~ ~ ~ fill ~-1 ~ ~-1 ~1 ~ ~1 air 0 replace carpet 11
/kill @e[type=item,name=ACV_CleansingGel,score_OnGround_min=1]

// Repulsion Gel
/execute @a ~ ~ ~ detect ~ ~ ~ carpet 11 effect @p 8 1 3 true
/execute @a ~ ~ ~ detect ~ ~ ~ air 0 effect @p 8 0 3 true

/scoreboard players set @e[type=item,name=ACV_RepulsionGel] OnGround 1 {OnGround:1b}
/execute @e[type=item,name=ACV_RepulsionGel,score_OnGround_min=1] ~ ~ ~ fill ~-1 ~ ~-1 ~1 ~ ~1 carpet 11 keep
/execute @e[type=item,name=ACV_RepulsionGel,score_OnGround_min=1] ~ ~ ~ fill ~-1 ~ ~-1 ~1 ~ ~1 carpet 11 replace carpet 1
/kill @e[type=item,name=ACV_RepulsionGel,score_OnGround_min=1]

// Propulsion Gel
/execute @a ~ ~ ~ detect ~ ~ ~ carpet 1 effect @p 1 1 6 true
/execute @a ~ ~ ~ detect ~ ~ ~ air 0 effect @p 1 0 6 true

/scoreboard players set @e[type=item,name=ACV_PropulsionGel] OnGround 1 {OnGround:1b}
/execute @e[type=item,name=ACV_PropulsionGel,score_OnGround_min=1] ~ ~ ~ fill ~-1 ~ ~-1 ~1 ~ ~1 carpet 1 keep
/execute @e[type=item,name=ACV_PropulsionGel,score_OnGround_min=1] ~ ~ ~ fill ~-1 ~ ~-1 ~1 ~ ~1 carpet 1 replace carpet 11
/kill @e[type=item,name=ACV_PropulsionGel,score_OnGround_min=1]
}


text:
#This test case covers the following items:

#MapleIR Direct read/write
#Direct acfgrecesses are associated with scalar variables and fixed fields inside
#structures. Direct accesses can be mapped to pseudo-register if the variable
#or field has no alias.
#In MAPLE IR, dassign and dread are the opcodes for direct assignments and
#direct references;

#Syntax
#dread<prim-type><var-name><field-id>
#dassign<var-name><field-id>(<rhs-expression>)
#Type syntax
#type <id-name><derived-type>
func $strcpy ( var %dest <*i8>, var %src <*i8> ) <*i8>

type $People < struct{
     @name <[10] i8>,
     @age i8 }>

type $Info < struct{
     @number i32,
     @employee <$People> }>

var $p <$People>
var $temp<[6] i8>= [ 104, 117, 97, 119, 101, 105 ]

func $main ( )i32 {
     var %a i8
     dassign %a ( constval i8 120 )
     var %b i8
     dassign %b ( dread i32 %a )
     var %i i32
     dassign %i ( constval i32 5 )
     var %ii i32
     dassign %ii ( dread i32 %i )
     var %j f32
     dassign %j ( constval f32 12.9f )
     var %jj f32
     dassign %jj ( dread f32 %j )
     var %k f64
     dassign %k ( constval f64 20.5 )
     var %kk f64
     dassign %kk ( dread f64 %k )
     var %p < $People >
     var %pp < $People >
     dassign $p 2 ( constval i32 25)
     # dassign $p 1 ( constval i32 0 )
     #call <pu-name>(<opnd0>,...,<opndn>)
     call &strcpy ( addrof a64 %p, addrof a64 $temp)
     call &strcpy ( addrof a64 %pp, addrof a64 %p)
     var %info < $Info >
     var %info2 < $Info >
     dassign %info 2( dread agg %p )
     dassign %info2 2( dread agg %info 2)
     return ( constval i32 444 )
}

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
$ifndef _LOGOUT_
$define _LOGOUT_

LogLevelHolder:=module()
    option object;
    export logLevel:=1;
end module:

# 设置输出级别
setLogLevel:=proc(v::{0,1,2,3})
    LogLevelHolder:-logLevel:=v;
    return;
end proc:

flog:=proc()
    local v;
    if type(procname,indexed) then
        v:=op(procname);
    else
        v:=1;
    end if;
    if (v>=LogLevelHolder:-logLevel) then
        print(_passed);
    end if;
    return;
end proc:

flogf:=proc()
    local v;
    if type(procname,indexed) then
        v:=op(procname);
    else
        v:=1;
    end if;
    if (v>=LogLevelHolder:-logLevel) then
        printf(_passed);
        printf("\n");
    end if;
    return;
end proc:

$endif


text:
//impulse: /give @p spawn_egg 1 0 {EntityTag:{id:"minecraft:squid"},display:{Name:"spawnMeg",Lore:["Use this to spawn","a Material Emancipation Grid"]}}
//impulse: /give @p spawn_egg 1 0 {EntityTag:{id:"minecraft:squid"},display:{Name:"deleteMeg",Lore:["Use this to delete","a Material Emancipation Grid"]}}

#MMH
repeat process MMH_megs {
  /execute @e[type=squid,name=spawnMeg] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:ACV_MEG,Duration:2147483647,Passengers:[{id:creeper,CustomName:ACV_MEG,powered:1,ActiveEffects:[{Amplifier:1,Duration:1000000,Id:14,ShowParticles:0b}],PersistenceRequired:1,Invulnerable:1,NoAI:1}]}
  /tp @e[type=squid,name=spawnMeg] ~ -100 ~
  /execute @e[type=squid,name=deleteMeg] ~ ~ ~ tp @e[name=ACV_MEG,dy=-1] ~ -100 ~
  /tp @e[type=squid,name=deleteMeg] ~ -100 ~
}


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

sg4_mu1     := 0.042:
sg4_mu2     := 0.26 - sg4_mu1:
sg4_nu_MGE4 := -0.195:
sg4_k2      := -sg4_mu2^2/sg4_nu_MGE4:
sg4_k1      := 0.804 - sg4_k2:

sg4_f0 := s -> 1 + sg4_k1 + sg4_k2
   - sg4_k1*(1 - sg4_mu1*s^2/sg4_k1)/(1 - (sg4_mu1*s^2/sg4_k1)^5)
   - sg4_k2/(1 + sg4_mu2*s^2/sg4_k2):
sg4_f:= x -> sg4_f0(X2S*x):

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