text:
# Inverse Dynamics for Robot based on MDH frames
# Einleitung
# Berechnung der inversen Dynamik
#
# 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
# dynamics -> Berechnung der Dynamik
# 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_act := true:
codegen_debug := false:
codegen_opt := 2:
codeexport_grav := true:
codeexport_corvec := true:
codeexport_cormat := true:
codeexport_inertia := true:
codeexport_inertiaD := true:
codeexport_invdyn := true:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_index_symmat2vector":
read "../helper/proc_symmat2vector":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotx":
read "../transformation/proc_troty":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../transformation/proc_trafo_mdh":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
# Einstellungen für Term-Vereinfachungen: 0=keine, 1=Massenmatrix, 2=Coriolis-Matrix, 3=Massenmatrix-Zeitableitung, Coriolis-Vektor, Invers-Dynamik
# Motivation: Massenmatrix ist Voraussetzung für andere symbolische Berechnungen (lohnt sich damit mehr)
# Coriolis-Matrix lässt sich besser vereinfachen, da viele Einträge getrennt voneinander betrachtet werden.
# Coriolis-Vektor bzw. Invers-Dynamik am kompliziertesten.
if not assigned(simplify_options) or simplify_options(9)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    # pauschale Annahme für kinematische Zwangsbedingungen.
    use_simplify := 1: # nur Massenmatrix
  end if:
else # Benutzer-Einstellungen:
 use_simplify := simplify_options(9): # neunter Eintrag ist für Dynamik
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 des Lagrange-Formalismus laden
read sprintf("../codeexport/%s/tmp/floatb_%s_lagrange_dUdq_s_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
dUdq_s := dUdq_s:
read sprintf("../codeexport/%s/tmp/floatb_%s_lagrange_dTdq_s_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
dTdq_s := dTdq_s:
read sprintf("../codeexport/%s/tmp/floatb_%s_lagrange_dTdqDdt_s_par%d_maple.m", robot_name, base_method_name, codegen_dynpar):
dTdqDdt_s := dTdqDdt_s:

DynString := "Term:":
if codeexport_grav then
  DynString := sprintf("%s g",DynString):
end if:
if codeexport_corvec then
  DynString := sprintf("%s c",DynString):
end if:
if codeexport_cormat then
  DynString := sprintf("%s C",DynString):
end if:
if codeexport_inertia then
  DynString := sprintf("%s M",DynString):
end if:
if codeexport_inertiaD then
  DynString := sprintf("%s MD",DynString):
end if:
if codeexport_invdyn then
  DynString := sprintf("%s tau",DynString):
end if:
printf("%s. Generiere Dynamik (%s) für %s mit Parametersatz %d und %s\n", \
  FormatTime("%Y-%m-%d %H:%M:%S"), DynString, robot_name, codegen_dynpar, base_method_name):

# Gravitational Load
# Generate
taug_s := dUdq_s:
save taug_s, sprintf("../codeexport/%s/tmp/gravload_par%d_maple.m", robot_name, codegen_dynpar):
# Matlab Export
# Belastung der Basis (ist falsch für floatb_twist, da das Moment durch diese Wahl der verallgemeinerten Koordinaten (floatb_twist) nicht berechnet werden kann!)
# Ist korrekt für floatb_eulxyz.
# Die Berechnungen werden deshalb nicht für floatb_twist durchgeführt.
if codegen_act and codeexport_grav and not(base_method_name="twist") then
  MatlabExport(taug_s(1..6), sprintf("../codeexport/%s/tmp/gravload_base_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Belastung der Gelenke
if codegen_act and codeexport_grav then
  MatlabExport(taug_s(7..NQ), sprintf("../codeexport/%s/tmp/gravload_joint_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Kompletter Vektor
if codegen_act and codeexport_grav and not(base_method_name="twist") then
  MatlabExport(taug_s(1..NQ), sprintf("../codeexport/%s/tmp/gravload_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Mass Matrix
# Generate
tauMM_s := dTdqDdt_s:
MM_s := Matrix(NQ, NQ):

if codeexport_inertia then
  for i to NQ do
    for j to NQ do
    	 if i > j then # linke untere Dreiecksmatrix
    	   MM_s[i, j] := MM_s[j,i]: # nehme symmetrischen Teil (wurde schon berechnet)
    	   next:
    	 end if:
      MM_s[i, j] := diff(tauMM_s[i, 1], qDD_s[j, 1]):
    end do:
  end do:
  # Terme vereinfachen
  if use_simplify>=1 then
    tmp_t0:=time():
    tmp_l0 := length(MM_s):
    printf("%s: Beginne Elementweise Vereinfachung Massenmatrix (Param. %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l0):
    for i to NQ do
      for j to NQ do
    	   if i > j then # linke untere Dreiecksmatrix
    	     MM_s[i, j] := MM_s[j,i]: # nehme symmetrischen Teil (wurde schon vereinfacht)
    	     next:
    	   end if:
        tmp_t1:=time():
        tmp_l1 := length(MM_s[i,j]):
        printf("%s: Beginne Vereinfachung Massenmatrix (Param. %d; Elem. %d,%d). Länge: %d.\n", \
          FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, j, tmp_l1):
        MM_s[i,j] := simplify2(MM_s[i,j]):
        tmp_t2:=time():
        tmp_l2 := length(MM_s[i,j]):
        printf("%s: Massenmatrix (Param. %d; Elem. %d,%d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
          FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, j, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
      end do:
    end do:
    tmp_l2 := length(MM_s):
    printf("%s: Massenmatrix (Param. %d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l0, tmp_l2, tmp_t2-tmp_t0):
  end if:
  save MM_s, sprintf("../codeexport/%s/tmp/inertia_par%d_maple.m", robot_name, codegen_dynpar):
end if:
# Matlab Export (nur linke untere Dreiecksmatrix bei vollständigen Matrizen)
if codeexport_inertia and not(base_method_name="twist") then
  MM_s_vek := symmat2vec(MM_s):
  filename_tmp := sprintf("../codeexport/%s/tmp/inertia_floatb_%s_par%d_maple_symvec.m", robot_name, base_method_name, codegen_dynpar):
  save MM_s_vek, filename_tmp:
  read filename_tmp:
  if codegen_act then
    MatlabExport(MM_s_vek, sprintf("../codeexport/%s/tmp/inertia_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  end if:
end if:
if codegen_act and codeexport_inertia and not(base_method_name="twist") then
  MatlabExport(MM_s(7..NQ,1..6), sprintf("../codeexport/%s/tmp/inertia_joint_base_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Gelenk-Massenmatrix.
# Eliminiere die Basis-Orientierung. Hat keinen Einfluss auf die Matrix, wird aber noch Maple teilweise nicht automatisch entfernt, da die Ausdrücke zu kompliziert zum Optimieren sind.
MMjj_s := MM_s(7..NQ,7..NQ):
for i from 1 to NQB do
  MMjj_s := subs({X_base_s[i,1]=0},MMjj_s):
end do:
if codeexport_inertia then
  MMjj_s_vek := symmat2vec(MMjj_s):
  filename_tmp := sprintf("../codeexport/%s/tmp/inertia_joint_joint_floatb_%s_par%d_maple_symvec.m", robot_name, base_method_name, codegen_dynpar):
  save MMjj_s_vek, filename_tmp:
  read filename_tmp:
  if codegen_act then
    MatlabExport(MMjj_s_vek, sprintf("../codeexport/%s/tmp/inertia_joint_joint_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  end if:
end if:
# Basis-Massenmatrix
MMbb_s := MM_s(1..6,1..6):
if codeexport_inertia and not(base_method_name="twist") then
  MMbb_s_vek := symmat2vec(MMbb_s):
  filename_tmp := sprintf("../codeexport/%s/tmp/inertia_base_base_floatb_%s_par%d_maple_symvec.m", robot_name, base_method_name, codegen_dynpar):
  save MMbb_s_vek, filename_tmp:
  read filename_tmp:
  if codegen_act then
    MatlabExport(MMbb_s_vek, sprintf("../codeexport/%s/tmp/inertia_base_base_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  end if:
end if:

# Mass Matrix Time Derivative
# Konvertiere Massenmatrix in zeitabhängige Variablen, um Zeitableitung zu berechnen

if codeexport_inertiaD then
  read sprintf("../codeexport/%s/tmp/inertia_par%d_maple.m", robot_name, codegen_dynpar):
  MM_s := MM_s: # Annahme: Massenmatrix wurde vorher generiert.
  MM_t := convert_s_t(MM_s):
  MMD_t := diff~(MM_t, t):
  MMD_s := convert_t_s(MMD_t):
  # Terme vereinfachen
  # Terme vereinfachen
  if use_simplify>=3 then
    tmp_t0:=time():
    tmp_l0 := length(MMD_s):
    printf("%s: Beginne Elementweise Vereinfachung Massenmatrix-Zeitableitung (Param. %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l0):
    for i to NQ do
      for j to NQ do
    	   if i > j then # linke untere Dreiecksmatrix
    	     MMD_s[i, j] := MMD_s[j,i]: # nehme symmetrischen Teil (wurde schon vereinfacht)
    	     next:
    	   end if:
        tmp_t1:=time():
        tmp_l1 := length(MMD_s[i,j]):
        printf("%s: Beginne Vereinfachung Massenmatrix-Zeitableitung (Param. %d; Elem. %d,%d). Länge: %d.\n", \
          FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, j, tmp_l1):
        MMD_s[i,j] := simplify2(MMD_s[i,j]):
        tmp_t2:=time():
        tmp_l2 := length(MMD_s[i,j]):
        printf("%s: Massenmatrix-Zeitableitung (Param. %d; Elem. %d,%d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
          FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, j, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
      end do:
    end do:
    tmp_l2 := length(MMD_s):
    printf("%s: Massenmatrix-Zeitableitung (Param. %d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l0, tmp_l2, tmp_t2-tmp_t0):
  end if:
end if:
# Matlab Export
if codeexport_inertiaD and not(base_method_name="twist") then
  MMD_s_vek := symmat2vec(MMD_s):
  filename_tmp := sprintf("../codeexport/%s/tmp/inertia_time_derivative_floatb_%s_par%d_maple_symvec.m", \
    robot_name, base_method_name, codegen_dynpar):
  save MMD_s_vek, filename_tmp:
  read filename_tmp:
  if codegen_act then
    MatlabExport(MMD_s_vek, sprintf("../codeexport/%s/tmp/inertia_time_derivative_floatb_%s_par%d_matlab.m", \
      robot_name, base_method_name, codegen_dynpar), codegen_opt):
  end if:
end if:

if codegen_act and codeexport_inertiaD and not(base_method_name="twist") then
  MatlabExport(MMD_s[7..NQ,1..6], sprintf("../codeexport/%s/tmp/inertia_joint_base_time_derivative_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
MMDjj_s := MMD_s(7..NQ,7..NQ):
for i from 1 to NQB do
  MMDjj_s := subs({X_base_s[i,1]=0},MMDjj_s):
end do:
for i from 1 to 6 do
  MMDjj_s := subs({V_base_s[i,1]=0},MMDjj_s):
end do:
if codeexport_inertiaD then
  # Vektor der unteren linken Dreiecksmatrix generieren, speichern und wieder laden. Ohne Neuladen hängt sich maple durch einen Bug auf. TODO: Klären.
  MMDjj_s_vek := symmat2vec(MMDjj_s):
  filename_tmp := sprintf("../codeexport/%s/tmp/inertia_joint_joint_time_derivative_floatb_%s_par%d_maple_symvec.m", robot_name, base_method_name, codegen_dynpar):
  save MMDjj_s_vek, filename_tmp:
  read filename_tmp:
  if codegen_act then
    MatlabExport(MMDjj_s_vek, sprintf("../codeexport/%s/tmp/inertia_joint_joint_time_derivative_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
  end if:
end if:

# Coriolis Vector
# Generiere Coriolis-Vektor aus dem Lagrange-Momenten
tauCC_s := Matrix(NQ,1): # Platzhalter
;

if codeexport_corvec then
  tauCC_s := dTdqDdt_s-dTdq_s:
  for i to NQ do
    tauCC_s := subs({qDD_s(i, 1) = 0}, tauCC_s):
  end do:
  # Terme vereinfachen
  if use_simplify>=3 then
    for i to NQ do
      tmp_t1:=time():
      tmp_l1 := length(tauCC_s[i,1]):
      printf("%s: Beginne Vereinfachung Coriolis-Vektor (Param. %d; Elem. %d). Länge: %d.\n", \
        FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, tmp_l1):
      tauCC_s[i,1] := simplify2(tauCC_s[i,1]):
      tmp_t2:=time():
      tmp_l2 := length(tauCC_s[i,1]):
      printf("%s: Coriolis-Vektor (Param %d; Elem. %d) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
        FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, i, tmp_l1, tmp_l2, tmp_t2-tmp_t1):
   end do:
  end if:
  save tauCC_s, sprintf("../codeexport/%s/tmp/coriolisvec_par%d_maple.m", robot_name, codegen_dynpar):
end if:

# Einzelne Komponenten exportieren (falls der ganze Vektor zu lange dauert).
if codegen_act and codeexport_corvec and not(base_method_name="twist") and codegen_debug then
  for i to NQ do
    MatlabExport(tauCC_s(i), sprintf("../codeexport/%s/tmp/coriolisvec_floatb_%s_%d_par%d_matlab.m", robot_name, base_method_name, i, codegen_dynpar), codegen_opt):
  end do:
end if:
# Matlab Export: Floating base
if codegen_act and codeexport_corvec and not(base_method_name="twist") then
  MatlabExport(tauCC_s[1..6], sprintf("../codeexport/%s/tmp/coriolisvec_base_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
if codegen_act and codeexport_corvec and not(base_method_name="twist") then
  MatlabExport(tauCC_s[7..NQ], sprintf("../codeexport/%s/tmp/coriolisvec_joint_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
if codegen_act and codeexport_corvec and not(base_method_name="twist") then
  MatlabExport(tauCC_s[1..NQ], sprintf("../codeexport/%s/tmp/coriolisvec_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Matlab Export: Fixed base
tauCC_s_fixb:=tauCC_s:
for i from 1 to NQB do
  tauCC_s_fixb := subs({X_base_s[i,1]=0},tauCC_s_fixb):
end do:
for i from 1 to 6 do
  tauCC_s_fixb := subs({V_base_s[i,1]=0},tauCC_s_fixb):
end do:
if codegen_act and codeexport_corvec then
  MatlabExport(tauCC_s_fixb[7..NQ], sprintf("../codeexport/%s/tmp/coriolisvec_joint_fixb_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
end if:

# Coriolis Matrix
# Calculation with Christoffel Symbol approach
# [KhalilDombre2002], equ. (9.7) (p. 195)

cijk := proc (i::integer, j::integer, k::integer, A, qs)
  local c:
  c := (1/2)*(diff(A[i, j], qs(k, 1)))+(1/2)*(diff(A[i, k], qs(j, 1)))-(1/2)*(diff(A[j, k], qs(i, 1))):
  return c:
end proc:

if codeexport_cormat then
  read sprintf("../codeexport/%s/tmp/inertia_par%d_maple.m", robot_name, codegen_dynpar):
  MM_s := MM_s: # Annahme: Massenmatrix wurde vorher generiert.
  Cqs:=Matrix(NQ,NQ):
  for i  from 1 to NQ do
    for j from 1 to NQ do
      Cqs[i,j]:=0:
      for k from 1 to NQ do
        Cqs[i,j]:=Cqs[i,j]+cijk(i,j,k,MM_s,q_s)*qD_s[k,1]:
      end do:
    end do:
  end do:
  # Terme vereinfachen
  if use_simplify>=2 then
    tmp_t1:=time():
    tmp_l1 := length(Cqs):
    printf("%s: Beginne Vereinfachung Coriols-Matrix (Param. %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l1):
    Cqs := simplify2(Cqs):
    tmp_t2:=time():
    tmp_l2 := length(Cqs):
    printf("%s: Coriolis-Matrix vereinfacht (Param. %d). 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:
end if:

# Matlab Export: Floating base
if codegen_act and codeexport_cormat and not(base_method_name="twist") then
  MatlabExport(Cqs[1..NQ,1..NQ], sprintf("../codeexport/%s/tmp/coriolismat_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
if codegen_act and codeexport_cormat and not(base_method_name="twist") then
  MatlabExport(Cqs[7..NQ,1..NQ], sprintf("../codeexport/%s/tmp/coriolismat_joint_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Matlab Export: Fixed base
Cqs_fixb:=Cqs:
if codeexport_cormat then
  for i from 1 to NQB do
    Cqs_fixb := subs({X_base_s[i,1]=0},Cqs_fixb):
  end do:
  for i from 1 to 6 do
    Cqs_fixb := subs({V_base_s[i,1]=0},Cqs_fixb):
  end do:
end if:
if codegen_act and codeexport_cormat then
  MatlabExport(Cqs_fixb[7..NQ,7..NQ], sprintf("../codeexport/%s/tmp/coriolismat_joint_fixb_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
end if:

# Joint Torques

if codeexport_invdyn then
  tau := dTdqDdt_s-dTdq_s+dUdq_s:
  # Terme vereinfachen
  if use_simplify>=3 then
    tmp_t1:=time():
    tmp_l1 := length(tau):
    printf("%s: Beginne Vereinfachung Invers-Dynamik (Param. %d). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), codegen_dynpar, tmp_l1):
    tau := simplify2(tau):
    tmp_t2:=time():
    tmp_l2 := length(tau):
    printf("%s: Invers-Dynamik (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:
end if:
# Matlab Export: Floating base
# Berechnung der Basis-Belastung ist für manche Basis-Darstellungen falsch (siehe oben unter Gravitationslast).
if codegen_act and codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(tau, sprintf("../codeexport/%s/tmp/invdyn_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
if codegen_act and codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(tau[1..6], sprintf("../codeexport/%s/tmp/invdyn_base_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
if codegen_act and codeexport_invdyn and not(base_method_name="twist") then
  MatlabExport(tau[7..NQ], sprintf("../codeexport/%s/tmp/invdyn_joint_floatb_%s_par%d_matlab.m", robot_name, base_method_name, codegen_dynpar), codegen_opt):
end if:
# Matlab Export: Fixed base
if codegen_act and codeexport_invdyn then
  taus_fixb:=tau:
  for i from 1 to NQB do
    taus_fixb := subs({X_base_s[i,1]=0},taus_fixb):
  end do:
  for i from 1 to 6 do
    taus_fixb := subs({V_base_s[i,1]=0},taus_fixb):
    taus_fixb := subs({VD_base_s[i,1]=0},taus_fixb):
  end do:
  save taus_fixb, sprintf("../codeexport/%s/tmp/invdyn_fixb_par%d_maple.m", robot_name, codegen_dynpar):
end if:
if codegen_act and codeexport_invdyn then
  MatlabExport(taus_fixb[7..NQ], sprintf("../codeexport/%s/tmp/invdyn_fixb_par%d_matlab.m", robot_name, codegen_dynpar), codegen_opt):
end if:


text:
#######################################################################
# This file is part of the crlibm library, and is distributed under
# the  LGPL.
# To use:
# restart; read "log-de.mpl";
#TODO output scripts for the Gappa proof out for Paterson/Stockmayer
Digits := 100:

interface(quiet=true):

read "common-procedures.mpl":
read "double-extended.mpl":
mkdir("TEMPLOG"):

log2h,log2l := hiloExt(log(2)):

L := 7: # number of bits used to address the table

MAXINDEX    := round(2^L * (sqrt(2)-1)):

for i from 0 to MAXINDEX-1 do
    center[i] := 1 + i*2^(-L): # center[i] in [1, 2[
    # We want it to fit on 11 bits of mantissa
    r[i] :=  round(evalf(  (1/center[i]) * 2^(11)) ) / 2^(11) ;

od:
for i from MAXINDEX to 2^L do
    # y has been divided by two, center[i] in [0.5, 1[
    center[i]:=(1 + i*2^(-L)) / 2:
    # We want it to fit on 11 bits of mantissa,
    r[i] :=  round(evalf(  (1/center[i]) * 2^(10)) ) / 2^(10) ;
od:

# Note that we go up to 2^L although the case 2^L is wrapped to zero
# in the C code. It could be important for zmax (but it turns out not).

for i from 0 to 2^L do
    (logirh[i], logirl[i]) := hiloExt(-log(r[i])):
od:

#Computation of ZMax.
for i from 0 to MAXINDEX-1 do
    t_x := center[i] + 2^(-L-1) :
    zmax[i] := (t_x*r[i]-1) :
    t_x := center[i] - 2^(-L-1) :
    zmin[i] := (t_x*r[i]-1) :
    zabsmax[i] := max(abs(zmin[i]), abs(zmax[i])):
od:
for i from MAXINDEX to 2^L do
    t_x := center[i] + 2^(-L-2) :
    zmax[i] := (t_x*r[i]-1) :
    t_x := center[i] - 2^(-L-2) :
    zmin[i] := (t_x*r[i]-1) :
    zabsmax[i] := max(abs(zmin[i]), abs(zmax[i])):
od:

zmaxmax:=0:
zminmin:=0:
for i from 0 to 2^L do
    if zmax[i] > zmaxmax then zmaxmax := zmax[i]: fi:
    if zmin[i] < zminmin then zminmin := zmin[i]: fi:
od:
printf("zminmin = -2^(%2f)   zmaxmax = 2^(%2f)\n", log2(-zminmin), log2(zmaxmax) ) :

PolyDegreeQuick:=7:

#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof

if(1+1=3) then # The polynomial used to be computed here,
polyQuick:= polyExact2Ext(x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-1,0], 1 ,  'deltaApprox'), 0):

epsilonApproxQuick := numapprox[infnorm]( 1-polyQuick/log(1+x), x=zminmin..zmaxmax):
deltaApproxQuick := numapprox[infnorm]( polyQuick-log(1+x), x=zminmin..zmaxmax):
else
# magical polynomial given by Sylvain's tool
polyQuick := x * (1 + (x * ((-0.50000000000000000000000000000000000000000000000000000000000000e0) + (x * (0.33333333333333342585191871876304503530263900756835937500000000e0 + (x * ((-0.24999999998423708125194764306797878816723823547363281250000000e0) + (x * (0.19999999996748479835773082413652446120977401733398437500000000e0 + (x * ((-0.16666957260954737285452154083031928166747093200683593750000000e0) + (x * 0.14286056338555042088955815415829420089721679687500000000000000e0)))))))))))):

# validated infinite norms given by Christoph's tool
epsilonApproxQuick := 2^(-64.119667587221):
deltaApproxQuick := 2^(-72.230387592107):
fi:
printf(" rel approximation error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) :
printf(" abs approximation error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) :

# For the Paterson/Stockmayer method
polyQuick0:= x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-1,0], 1 ,  'deltaApprox'):

## Création des nouveaux coefficients
polyQuick0:= x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-1,0], 1 ,  'deltaApprox'):

## Création des nouveaux coefficients
a7 := convert(coeff(polyQuick0,x,7),rational):
a6 := convert(coeff(polyQuick0,x,6)/a7,rational):
a5 := convert(coeff(polyQuick0,x,5)/a7,rational):
a4 := convert(coeff(polyQuick0,x,4)/a7,rational):

a3 := convert(coeff(polyQuick0,x,3),rational):
a2 := convert(coeff(polyQuick0,x,2),rational):
a1 := convert(coeff(polyQuick0,x,1),rational):
a0 := convert(coeff(polyQuick0,x,0),rational):

alpha0 := a5 - 1:
alpha1 := a6:
alpha2 := a4 - alpha0 * a6:

##
a7_hi,a7_lo         := hiloExt(a7):
a3_hi,a3_lo         := hiloExt(a3):
a2_hi,a2_lo         := hiloExt(a2):
a0_hi,a0_lo         := hiloExt(a0):

c_hi,c_lo         := hiloExt(alpha0):
alpha_hi,alpha_lo := hiloExt(alpha1):
beta_hi,beta_lo   := hiloExt(alpha2):

#coef_c := c_hi + c_lo:
#coef_alpha := alpha_hi + alpha_lo:
#coef_beta := beta_hi + beta_lo:
#coef_c7 := a7_hi + a7_lo:
#coef_c3 := a3_hi + a3_lo:
#coef_c2 := a2_hi + a2_lo:
#coef_c0 := a0_hi + a0_lo:

coef_c := c_hi:
coef_alpha := alpha_hi:
coef_beta := beta_hi:
coef_c7 := a7_hi:
coef_c3 := a3_hi:
coef_c2 := a2_hi:
coef_c0 := a0_hi:

polyQuickPat := ((x^2 + coef_c)*(x + coef_alpha) + (x + coef_beta)) * (coef_c7*x^4) + ((coef_c3 * x + coef_c2)*x^2 + (x)):

epsilonApproxQuickPat := numapprox[infnorm]( 1-polyQuickPat/log(1+x), x=zminmin..zmaxmax):
printf("   approximation rel error for Paterson/Stockmayer in the quick phase is 2^(%2f)\n", log2(epsilonApproxQuickPat) ) :
deltaApproxQuickPat := numapprox[infnorm]( polyQuickPat-log(1+x), x=zminmin..zmaxmax):
printf("   approximation abs error for Paterson/Stockmayer in the quick phase is 2^(%2f)\n", log2(deltaApproxQuickPat) ) :

PolyDegreeAccurate:=14:

printf("Computing the polynomial for accurate phase (this may take some time...)\n"):
pe:= x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeAccurate-1,0], 1 ,  'deltaApprox'):

MaxDegreeDDE:=8:  #

polyAccurate := polyExact2Ext(pe, MaxDegreeDDE):
deltaApproxAccurate := numapprox[infnorm](polyAccurate-log(1+x), x=-zmaxmax..zmaxmax):
epsilonApproxAccurate := numapprox[infnorm]( 1-polyAccurate/log(1+x), x=-zmaxmax..zmaxmax):
printf("   approximation error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) :

filename:="TEMPLOG/log-de.h":
fd:=fopen(filename, WRITE, TEXT):

fprintf(fd, "/*File generated by maple/log-de.mpl*/\n\n"):
# some constants
fprintf(fd, "#define L        %d\n", L):
fprintf(fd, "#define MAXINDEX %d\n", MAXINDEX):
fprintf(fd, "#define INDEXMASK %d\n", 2^L-1):
fprintf(fd, "static const double two64 = %1.25e ;\n", evalf(2^64)):

fprintf(fd, "#if 1\n#define ESTRIN\n#else\n#define PATERSON\n#endif\n"):

fprintf(fd, "#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)\n\n"):

fprintf(fd, "#ifdef PATERSON\n"):
  #  polynomial for quick phase
  fprintf(fd,"static const long double c7 = %1.50e;\n",coef_c7):
  fprintf(fd,"static const long double c3 = %1.50e;\n",coef_c3):
  fprintf(fd,"static const long double c2 = %1.50e;\n",coef_c2):
  fprintf(fd,"static const long double ps_alpha = %1.50e;\n",coef_alpha):
  fprintf(fd,"static const long double ps_beta = %1.50e;\n",coef_beta):
  fprintf(fd,"static const long double ps_c = %1.50e;\n",coef_c):
fprintf(fd, "#endif/* PATERSON*/\n\n"):

fprintf(fd, "#ifdef ESTRIN\n"):
 fprintf(fd, "  /* Coefficients are read directly from the array thanks to the following macros */ \n"):
  for i from PolyDegreeQuick to 2 by -1 do
    fprintf(fd, "#define c%d  c[%d]\n", i, PolyDegreeQuick-i):
  od:
   fprintf(fd, "static const double c[%d] =  {\n",PolyDegreeQuick):
   for i from PolyDegreeQuick to 2 by -1 do
     fprintf(fd, "   /* c%d  = */ %1.50eL, \n", i, coeff(polyQuick,x,i)):
   od:
  fprintf(fd, "}; \n \n"):
fprintf(fd, "#endif/* ESTRIN */\n\n"):

  for i from PolyDegreeAccurate to 1 by -1 do
    fprintf(fd, "#define c%dh  ch[%d]\n", i, PolyDegreeAccurate-i):
  od:

  for i from MaxDegreeDDE-1 to 1 by -1 do
    fprintf(fd, "#define c%dl  cl[%d]\n", i, MaxDegreeDDE-1-i):
  od:
  fprintf(fd, "#define PREFETCH_POLY_ACCURATE \n"):
  fprintf(fd, "\n#else /* not(defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)),\n   assuming Itanium, otherwise we shouldn't be there */ \n\n"):
  fprintf(fd, "#define PREFETCH_POLY_QUICK "):
  for i from PolyDegreeQuick to 2 by -1 do
    fprintf(fd, "c%d=c[%d]; ", i, PolyDegreeQuick-i):
  od:
  fprintf(fd, "\n#define PREFETCH_POLY_ACCURATE "):
  for i from PolyDegreeAccurate to 1 by -1 do
    fprintf(fd, "c%dh=ch[%d]; ", i, PolyDegreeAccurate-i):
    if i mod 4 =0 then  fprintf(fd, "\\\n        "): fi:
  od:
  fprintf(fd, "\\\n        "):
  for i from MaxDegreeDDE-1 to 1 by -1 do
    fprintf(fd, "c%dl=cl[%d]; ", i, MaxDegreeDDE-1-i):
  od:

  fprintf(fd, "\n#endif /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */ \n\n"):

  # Various constants
  fprintf(fd, "static const long double log2H  = %1.50eL ;\n", log2h):
  fprintf(fd, "static const long double log2L  = %1.50eL ;\n", log2l):

  # The polynomials
  #  polynomial for quick phase
#  for i from PolyDegreeQuick to 1 by -1 do
#    fprintf(fd, "static const long double c%d =    %1.50eL ;\n", i, coeff(polyQuick,x,i)):
#  od:

  #  polynomial for accurate phase
  fprintf(fd, "static const long double ch[%d] =  {\n",PolyDegreeAccurate):
   for i from PolyDegreeAccurate to 1 by -1 do
     (ch, cl) := hiloExt(coeff(polyAccurate,x,i)):
     fprintf(fd, "   /* ch%d  = */ %1.50eL, \n", i, ch):
   od:
  fprintf(fd, "}; \n \n"):

  fprintf(fd, "static const long double cl[%d] =  {\n", MaxDegreeDDE):
  for i from MaxDegreeDDE-1 to 1 by -1 do
     (ch, cl) := hiloExt(coeff(polyAccurate,x,i)):
     fprintf(fd, "   /* cl%d  = */ %1.50eL, \n", i, cl):
   od:
  fprintf(fd, "}; \n \n"):

  # The tables

  fprintf(fd, "typedef struct rri_tag {float r; long double logirh;  long double logirl; } rri ;  \n"):
  fprintf(fd, "static const rri argredtable[%d] = {\n", 2^L):
  for i from 0 to 2^L-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50eL,   /* r[%d] */ \n", r[i], i):
      fprintf(fd, "    %1.50eL, /* logirh[%d] */ \n", logirh[i], i):
      fprintf(fd, "    %1.50eL, /* logirl[%d] */ \n", logirl[i], i):
      fprintf(fd, "  } "):
      if(i<2^L-1) then  fprintf(fd, ", \n"): fi
  od:
fprintf(fd, "}; \n \n"):

fclose(fd):

printf("\n************ DONE TEMPLOG/log-de.h ************\n");

filename:="TEMPLOG/polynomials.sed":
fd:=fopen(filename, WRITE, TEXT):
for i from PolyDegreeAccurate to 1 by -1 do
    (ch, cl) := hiloExt(coeff(polyAccurate,x,i)):
    fprintf(fd, " s/_c%dh/%1.40e/g\n", i, ch):
    if(i<MaxDegreeDDE) then
        fprintf(fd, " s/_c%dl/%1.40e/g\n", i, cl):
    fi:
od:
for i from PolyDegreeQuick to 1 by -1 do
    fprintf(fd, " s/_c%d/%1.40e/g\n", i, coeff(polyQuick,x,i)):
od:
fprintf(fd, " s/_deltaApproxQuick/%1.40e/g\n", deltaApproxQuick):
fprintf(fd, " s/_epsilonApproxQuick/%1.40e/g\n", epsilonApproxQuick):
fprintf(fd, " s/_deltaApproxAccurate/%1.40e/g\n", deltaApproxAccurate):
fprintf(fd, " s/_epsilonApproxAccurate/%1.40e/g\n", epsilonApproxAccurate):
fprintf(fd, " s/_log2h/%1.40e/g\n", log2h):
fprintf(fd, " s/_log2l/%1.40e/g\n", log2l):
fclose(fd):

printf("\n************ DONE TEMPLOG/polynomials.sed ************\n");

for j from 0 to 2^L-1 do
    filename:=cat("TEMPLOG/log-de_",j,".sed"):
    fd:=fopen(filename, WRITE, TEXT):
    fprintf(fd, " s/_rval/%1.40e/g\n",     r[j]):
    fprintf(fd, " s/_logirh/%1.40e/g\n", logirh[j]):
    fprintf(fd, " s/_logirl/%1.40e/g\n", logirl[j]):
    fprintf(fd, " s/_zabsmax/%1.40e/g\n", zabsmax[j]):
    fclose(fd):
  od:

printf("\n************ DONE TEMPLOG/log-de*.sed ************\n");

# A shell script to use them
filename:="../gappa/log-de/run-log-de-proof.sh":
fd:=fopen(filename, WRITE, TEXT):
fprintf(fd, "#!/bin/sh\n"):
fprintf(fd, "# You probably need to edit the path to the gappa executable\n"):
fprintf(fd, "GAPPA=~/local/src/gappa/src/gappa\n"):
fprintf(fd, "CRLIBMDIR=../..\n"):

fprintf(fd, "echo Accurate phase, case E=0, index=0: 1>&2\n"):
fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index0-E0.gappa | $GAPPA \n"):

fprintf(fd, "echo Accurate phase, case E!=0, index=0 1>&2\n"):
fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index0-E1N.gappa | $GAPPA \n"):

fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
fprintf(fd, "  echo Accurate phase, case E=0, index=$num 1>&2\n"):
fprintf(fd, "  sed -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_$num.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index1N-E0.gappa | $GAPPA \n"):
fprintf(fd, "  echo 1>&2\n"):
fprintf(fd, "done\n"):

fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
fprintf(fd, "  echo Accurate phase, case E!=0, index = $num 1>&2 \n"):
fprintf(fd, "  sed -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_$num.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index1N-E1N.gappa | $GAPPA \n"):
fprintf(fd, "  echo 1>&2\n"):
fprintf(fd, "done\n\n"):

fprintf(fd, "echo Quick phase, case E=0, index=0  1>&2\n"):
fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-index0-E0.gappa | $GAPPA \n"):
fprintf(fd, "  echo 1>&2\n"):

fprintf(fd, "echo Quick phase, case E!=0, index=0 1>&2 \n"):
fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-index0-E1N.gappa | $GAPPA \n"):
fprintf(fd, "  echo 1>&2\n"):

fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
fprintf(fd, "  echo Quick phase, for all E,  index=$num 1>&2\n"):
fprintf(fd, "  sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_$num.sed $CRLIBMDIR/gappa/log-de/log-de-index1N-E0N.gappa | $GAPPA \n"):
fprintf(fd, "  echo 1>&2\n"):
fprintf(fd, "done\n"):

fclose(fd):

printf("\n************ DONE TEMPLOG/run-log-de-proof.sh ************\n"):
printf("Now you should run \n"):
printf(" sh ../../gappa/log-de/run-log-de-proof.sh  2> ../../maple/TEMPLOG/Gappa.out\n"):
printf("  (You probably need to edit the path to the gappa executable within run-log-de-proof.sh)\n"):
printf("Then look at TEMPLOG/Gappa.out. It shouldn't contain 'some enclosures were not satisfied'.\n If it does, report it !\n"):

printf("----DONE---\n") :

ext: mpl

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


text:
repeat process ACV_megs {
/execute @e[type=creeper,name=ACV_MEG] ~ ~ ~ kill @e[type=armor_stand,name=ACV_Cube,dy=0]

/scoreboard players tag @a[tag=ACV_Blue] remove ACV_Blue
/scoreboard players tag @a[tag=ACV_Orange] remove ACV_Orange
/scoreboard players tag @a[tag=ACV_Purple] remove ACV_Purple
/scoreboard players tag @a[tag=ACV_Red] remove ACV_Red

/scoreboard players tag @a add ACV_Blue {Inventory:[{id:"minecraft:bow",tag:{ench:[{id:48,lvl:6}]}}]}
/scoreboard players tag @a add ACV_Orange {Inventory:[{id:"minecraft:bow",tag:{ench:[{id:48,lvl:7}]}}]}
/scoreboard players tag @a add ACV_Purple {Inventory:[{id:"minecraft:bow",tag:{ench:[{id:48,lvl:8}]}}]}
/scoreboard players tag @a add ACV_Red {Inventory:[{id:"minecraft:bow",tag:{ench:[{id:48,lvl:9}]}}]}

/execute @e[type=creeper,name=ACV_MEG] ~ ~ ~ /scoreboard players tag @a[dy=0] add ACV_megDeletePortal_PARAM
conditional: /setblock ${this + 1} redstone_block

skip

impulse: /setblock ${this - 1} stone

/execute @p[tag=ACV_megDeletePortal_PARAM] ~ ~ ~ execute @p[tag=ACV_Blue,r=0] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Blue"],Duration:2147483647}
/execute @p[tag=ACV_megDeletePortal_PARAM] ~ ~ ~ execute @p[tag=ACV_Orange,r=0] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Orange"],Duration:2147483647}
/execute @p[tag=ACV_megDeletePortal_PARAM] ~ ~ ~ execute @p[tag=ACV_Purple,r=0] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Purple"],Duration:2147483647}
/execute @p[tag=ACV_megDeletePortal_PARAM] ~ ~ ~ execute @p[tag=ACV_Red,r=0] ~ ~ ~ summon area_effect_cloud ~ ~ ~ {CustomName:"ACV_remote_deletePortals_PARAM",Tags:["ACV_Red"],Duration:2147483647}
/scoreboard players tag @a remove ACV_megDeletePortal_PARAM
start ACV_remote_deletePortals
}


text:
with(simplex);

Constraints := [
# conservation of units at each node
e13+e14 = 300, # CHI
e23+e24 = 300, # DC
e13+e23 = 300, # HOU
e14+e24 = 300, # BOS
# maximum flow on individual edges
0 <= e13, e13 <= 200,
0 <= e14, e14 <= 200,
0 <= e23, e23 <= 280,
0 <= e24, e24 <= 350
];

Cost := 7*e13 + 6*e14 + 4*e23 + 6*e24;
# Invoke linear programming to solve problem
minimize (Cost, Constraints, NONNEGATIVE);


text:
(*

区间对象
实现运算：
    + 交集
    + 并集
    + 补集
    + 子集
实现了三种初始化方法：
    + 不等式集合
    + RealRange表达式
    + 区间的字符串表达形式
导出了辅助RealRange计算的两个函数

*)
$ifndef _SEG_
$define _SEG_

Seg:=module()
    option object;
    export
            # 运算符，优先级从高到低
            ## 逻辑版本
            `not`::static,          # 区间补集
            `and`::static,          # 区间交集
            `or`::static,           # 区间并集
            `xor`::static,          # 区间差集
            `implies`::static,      # 区间子集
            ## 集合版本
            `&C`::static,           # 区间补集
            `intersect`::static,    # 区间交集
            `union`::static,        # 区间并集
            `minus`::static,        # 区间差集
            `subset`::static,       # 区间属于

            # 属性
            bound:=real,            # 对应的 RealRange 对象

            # 导出工具函数
            formatRange::static,    # 将RealRange转化为区间表示形式
            evalRange::static;      # RealRange彻底求值
    local
            # 初始化
            ModuleApply::static,    # 初始化
            conBuild::static,       # 使用约束集合初始化
            rangeBuild::static,     # 使用 RealRange 初始化
            con2range::static,      # 约束转化为 RealRange 对象
            strBuild::static,       # 使用字符串初始化
            exprOutput::static,     # 字符串解析辅助函数
            # 输出
            ModulePrint::static,    # 输出
            range2str::static,      # RealRange 对象转化为字符串
            sortOps::static,        # 输出时按左值排序
            leftBound::static,      # 获取 RealRange 左端值,用于排序
            # 计算
            expandRange::static,    # 化简 RealRange 对象
            subsRange::static,      # 替换Range表达式中的集合对象和Non表达式
            rangeNot::static;       # RealRange 对象取补集

    # 初始化
    ModuleApply:=proc(x)
        if type(x,set({`=`,`<`,`<=`,`<>`})) then
            return conBuild(x);
        elif type(x,string) then
            return strBuild(x);
        else
            return rangeBuild(x);
        end if;
    end proc:

    # 使用约束集合初始化
    conBuild:=proc(cons::set({`=`,`<`,`<=`,`<>`}))
        if numelems(indets(cons,name))<>1
        or ormap(x->numelems(indets(cons,name))<>1,cons) then
            error "每个不等式（等式）都只能使用同一个变量";
        end if;
        return rangeBuild(AndProp(con2range~(cons)[]));
    end proc:

    # 使用 RealRange 初始化
    rangeBuild:=proc(r)
        local this;
        this:=Object(Seg);
        this:-bound:=evalRange(r);
        return this;
    end proc:

    # 将约束转化为RealRange
    con2range:=proc(_c::{`=`,`<`,`<=`,`<>`})
        local c,v;
        c:=RealDomain:-solve({_c})[];
        if type(c,`<`) then
            if type(lhs(c),numeric) then
                return RealRange(Open(lhs(c)),infinity);
            else
                return RealRange(-infinity,Open(rhs(c)));
            end if;
        elif type(c,`<=`) then
            if type(lhs(c),numeric) then
                return RealRange(lhs(c),infinity);
            else
                return RealRange(-infinity,rhs(c));
            end if;
        else
            if type(c,`=`) then
                return rhs(c);
            else
                v:=Open(rhs(c));
                return OrProp(RealRange(-infinity,v),RealRange(v,infinity));
            end if;
        end if;
    end proc:

    # 交集
    `and`:=proc()
        option overload;
        return rangeBuild(AndProp(map(x->x:-bound,[_passed])[]));
    end proc:
    `intersect`:=`and`;

    # 并集
    `or`:=proc()
        option overload;
        return rangeBuild(OrProp(map(x->x:-bound,[_passed])[]));
    end proc:
    `union`:=`or`;

    # 补集
    `not`:=proc(x::Seg,$)
        option overload;
        return rangeBuild(rangeNot(x:-bound));
    end proc:
    `&C`:=`not`;

    # 差集
    `xor`:=proc(x::Seg,y::Seg,$)
        return x and (not y);
    end proc:
    `minus`:=`xor`;

    # 子集
    `subset`:=proc(x::Seg,y::Seg,$)
        local z;
        z:=x and y;
        return evalb(z:-bound=x:-bound);
    end proc:
    `implies`:=`subset`;

    # RealRange彻底求值
    evalRange:=proc(x)
        return expandRange(subsRange(x));
    end proc:

    # RealRange补集展开
    rangeNot:=proc(b)
        local lv,rv,lb,rb;
        if b=BottomProp then
            return RealRange(-infinity,infinity);
        elif type(b,infinity) then
            return real;
        elif type(evalf(b),numeric) then
            return OrProp(RealRange(-infinity,Open(b)),RealRange(Open(b),infinity));
        elif b=real then
            return BottomProp;
        elif op(0,b)=RealRange then
            if op([1,0],b)=Open then
                lv:=op([1,1],b);
            else
                lv:=Open(op(1,b));
            end if;
            if op([2,0],b)=Open then
                rv:=op([2,1],b);
            else
                rv:=Open(op(2,b));
            end if;
            if has(lv,infinity) then
                lb:=NULL;
            else
                lb:=RealRange(-infinity,lv);
            end if;
            if has(rv,infinity) then
                rb:=NULL;
            else
                rb:=RealRange(rv,infinity);
            end if;
            return OrProp(lb,rb);
        elif op(0,b)=AndProp then
            return OrProp(thisproc~([op(b)])[]);
        elif op(0,b)=OrProp then
            return AndProp(thisproc~([op(b)])[]);
        end if;
    end proc:

    # 把 property 表达式中的集合都替换为 OrProp
    # 把 Non 展开
    subsRange:=proc(x)
        if type(x,set) then
            return OrProp(op(x));
        elif op(0,x)=Non then
            return rangeNot(op(x));
        elif op(0,x)=AndProp or op(0,x)=OrProp then
            return op(0,x)(thisproc~([op(x)])[]);
        else
            return x;
        end if;
    end proc:

    # 展开计算RealRange
    expandRange:=proc(x)
        local orv,nxt,rst,tmp;
        if op(0,x)=AndProp then
            # 所有的交集操作的返回结果都不应该表示为交集的形式
            if has([op(x)],OrProp) then
                # A⋂(B⋃C)=(A⋂B)⋃(A⋂C)
                orv:=op(1,x);
                nxt:=op(2,x);
                rst:=op(3..-1,x);
                tmp:=OrProp(map(x->AndProp(x,nxt),[op(orv)])[]);
                tmp:=AndProp(tmp,rst);
            elif type(op(1,x),extended_numeric) then
                # 计算 AndProp(x,property)
                # 因为集合已经被事先展开，所以这里只可能有一个数值对象
                # 而且经测试，顺序为：numeric < OrProp < RealRange
                if is(op(1,x),op(2,x)) then
                    tmp:=op(1,x);
                else
                    tmp:=BottomProp;
                end if;
            else
                tmp:=x;
            end if;
            return thisproc(tmp);
        elif op(0,x)=OrProp and has([op(x)],AndProp) then
        	return thisproc(OrProp(thisproc~([op(x)])[]));
        else
            return x;
        end if;
    end proc:

    # 将RealRange转化为区间表示形式
    range2str:=proc(b)
        local lc,lv,rc,rv,bs,r;
        if b=BottomProp then
            return "∅";
        elif type(evalf(b),extended_numeric) then
            return sprintf("{%a}",b);
        elif b=real then
            return "( -∞ , +∞ )";
        elif op(0,b)=RealRange then
            if op([1,0],b)=Open then
                lc:="(";
                lv:=op([1,1],b);
            else
                lc:="[";
                lv:=op(1,b);
            end if;
            if op([2,0],b)=Open then
                rc:=")";
                rv:=op([2,1],b);
            else
                rc:="]";
                rv:=op(2,b);
            end if;
            if type(lv,infinity) then
                lv:="-∞";
                lc:="(";
            else
                lv:=convert(lv,string);
            end if;
            if type(rv,infinity) then
                rv:="+∞";
                rc:=")";
            else
                rv:=convert(rv,string);
            end if;
            return sprintf("%s %s , %s %s",lc,lv,rv,rc);
        elif op(0,b)=OrProp then
            bs:=thisproc~(sortOps([op(b)]));
            r:=cat(map(x->cat(x," ⋃ "),bs)[]);
            return sprintf("< %s >",r[1..(-length(" ⋃ ")-1)]);
        elif op(0,b)=AndProp then
            bs:=thisproc~(sortOps([op(b)]));
            r:=cat(map(x->cat(x," ⋂ "),bs)[]);
            return sprintf("< %s >",r[1..(-length(" ⋂ ")-1)]);
        end if;
    end proc:

    # 用于排序
    sortOps:=proc(x)
        return sort(x,key=leftBound);
    end proc:

    # 区间左值
    leftBound:=proc(b)
        if b=BottomProp then
            return evalf(-infinity);
        elif type(evalf(b),extended_numeric) then
            return evalf(b);
        else
            return thisproc(op(1,b));
        end if;
    end proc:

    # 默认输出
    ModulePrint:=proc(x::Seg)
        return formatRange(x:-bound);
    end proc:

    # 将RealRange转化为区间表示形式
    formatRange:=proc(x)
        local r:=range2str(x);
        if r[1]="<" then
            r:=r[3..-3];
        end if;
        return r;
    end proc:

    # 根据区间的字符串表达形式初始化
    strBuild:=proc(_s::string)
        local r,s,cb,cx,cz,i,n,os,st,md,ed;
        s:=_s;
        # 处理无穷和交并符号
        s:=StringTools:-SubstituteAll(s,"\x26\x69\x6E\x66\x69\x6E\x3B","infinity");
        s:=StringTools:-SubstituteAll(s,"\x26\x62\x69\x67\x63\x61\x70\x3B","and");
        s:=StringTools:-SubstituteAll(s,"\x26\x62\x69\x67\x63\x75\x70\x3B","or");
        # 内部默认使用逻辑版本的运算符
        s:=StringTools:-SubstituteAll(s,"&C","not");
        s:=StringTools:-SubstituteAll(s,"intersect","and");
        s:=StringTools:-SubstituteAll(s,"union","or");
        s:=StringTools:-SubstituteAll(s,"minus","xor");
        s:=StringTools:-SubstituteAll(s,"subset","implies");
        # 解析字符串
        r:=StringTools:-StringBuffer();
        n:=length(s);
        cb:=0;# 区间括号标记
        cx:=0;# 小括号标记
        cz:=0;# 中括号标记
        os:=1;# 未处理的字符串的开始下标
        st:=0;# 区间开始下标
        md:=0;# 区间分隔符下标
        ed:=0;# 区间结束下标
        for i from 1 to n do
            if   s[i]="(" then
                if cx=0 and cb=0 then
                    cb:=1;
                    st:=i;
                else
                    cx:=cx+1;
                end if;
            elif s[i]="[" then
                if cz=0 and cb=0 then
                    cb:=1;
                    st:=i;
                else
                    cz:=cz+1;
                end if;
            elif s[i]="," and cb=1 and cx=0 and cz=0 then
            # 这么写保证了只有在其它表达式的()[]匹配的情况下才会将读取的`,`作为区间的分隔符，
            # 考虑到了其它表达式中出现`,`的可能性。
                md:=i;
            elif s[i]=")" then
                if cx=0 and cb=1 then
                    ed:=i;
                    cb:=0;
                    r:-append(exprOutput(s,os,st,md,ed));
                    os:=i+1;
                else
                    cx:=cx-1;
                end if;
            elif s[i]="]" then
                if cz=0 and cb=1 then
                    ed:=i;
                    cb:=0;
                    r:-append(exprOutput(s,os,st,md,ed));
                    os:=i+1;
                else
                    cz:=cz-1;
                end if;
            end if;
        end do;
        r:-append(s[os..-1]);
        s:=r:-value();
        r:-clear();
        # 处理单点区间
        s:=StringTools:-RegSubs("{(.*)}"="Seg({\\1})",s);
        # 处理括号
        s:=StringTools:-SubstituteAll(s,"<","(");
        s:=StringTools:-SubstituteAll(s,">",")");
        # 返回结果
        return eval(subs(Open(infinity)=infinity,Open(-infinity)=-infinity,parse(s)));
    end proc:

    # 辅助转化区间表达式为Maple表达式
    exprOutput:=proc(s,os,st,md,ed)
        local lv,rv,res;
        if s[st]="(" then
            lv:=sprintf("Open(%s)",s[(st+1)..(md-1)]);
        else
            lv:=s[(st+1)..(md-1)];
        end if;
        if s[ed]=")" then
            rv:=sprintf("Open(%s)",s[(md+1)..(ed-1)]);
        else
            rv:=s[(md+1)..(ed-1)];
        end if;
        res:=sprintf("%sSeg(RealRange(%s,%s))",s[os..(st-1)],lv,rv);
        return res;
    end proc:

end module:

$endif


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

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

(* type: work_mgga_x *)

params_a_kappa := 4.8827323:
params_a_mu    := 0.3511128:
$include "gga_x_pbe.mpl"

dldf_a := [1, -0.1637571, -0.1880028, -0.4490609, -0.0082359]:
csi_HF := 1 - 0.6144129:

f := (rs, x, t, u) ->
  + csi_HF*f_pbe(x)*mgga_series_w(dldf_a, 5, t):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
DV := NumberTheory[Divisors]:

ghost_coeff := (n) -> proc(x)
 local d;
 return add(d * x[d]^(n/d), d in DV(n));
end:

ghost_vector := (n::posint) -> proc(x)
 local m;
 return [seq(ghost_coeff(m)(x),m=1..n)];
end:

unghost_vector := (n::posint) -> proc(w)
 local x,m,d;
 x := []:
 for m from 1 to n do
  x := [op(x),expand((w[m] - add(d*x[d]^(m/d),d in DV(m) minus {m}))/m)];
 od:
 return x;
end:

witt_add_term := proc(n::posint)
 option remember;
 local x,y,u,d,rel;

 for d in DV(n) minus {n} do
  u[d] := witt_add_term(d)(x,y);
 od:

 rel := ghost_coeff(n)(u) - ghost_coeff(n)(x) - ghost_coeff(n)(y);
 u[n] := expand(solve(rel,u[n]));
 return unapply(u[n],x,y);
end:

witt_add := (n::posint) -> proc(x,y)
 local m;
 return [seq(witt_add_term(m)(x,y),m=1..n)];
end:

witt_mult_term := proc(n::posint)
 option remember;
 local x,y,u,d,rel;

 for d in DV(n) minus {n} do
  u[d] := witt_mult_term(d)(x,y);
 od:

 rel := ghost_coeff(n)(u) - ghost_coeff(n)(x) * ghost_coeff(n)(y);
 u[n] := expand(solve(rel,u[n]));
 return unapply(u[n],x,y);
end:

witt_mult := (n::posint) -> proc(x,y)
 local m;
 return [seq(witt_mult_term(m)(x,y),m=1..n)];
end:

witt_series := (n::posint) -> proc(x)
 local m;
 return expand(rem(mul(1 - x[m]*t^m,m=1..n),t^(n+1),t));
end:

witt_unseries := (n::posint) -> proc(f)
 local c,g,d;
 c := []:
 g := f;
 for d from 1 to n do
  c := [op(c),-coeff(g,t,d)];
  g := expand(convert(series(g/(1-c[d]*t^d),t=0,n+1),polynom,t));
 od:
 return c;
end:

ghost_p_coeff := (n) -> proc(x)
 local i;
 return add(p^i * x[i]^(p^(n-i)), i=0..n);
end:

ghost_p_vector := (n::nonnegint) -> proc(x)
 local i;
 return table([seq(i = ghost_p_coeff(i)(x),i=0..n-1)]);
end:

ghost_p_add := (n::nonnegint) -> proc(a,b)
 local i;
 return table([seq(i=expand(a[i]+b[i]),i=0..n-1)]);
end:

ghost_p_sub := (n::nonnegint) -> proc(a,b)
 local i;
 return table([seq(i=expand(a[i]-b[i]),i=0..n-1)]);
end:

ghost_p_mult := (n::nonnegint) -> proc(a,b)
 local i;
 return table([seq(i=expand(a[i]*b[i]),i=0..n-1)]);
end:

ghost_p_is_equal := (n::nonnegint) -> proc(a,b)
 local i;
 return `and`(seq(evalb(simplify(a[i]-b[i])=0),i=0..n-1));
end:

unghost_p_vector := (n::nonnegint) -> proc(w)
 local x,m,k;
 x := table:
 for m from 0 to n-1 do
  x[m] := expand((w[m] - add(p^k*x[k]^(p^(m-k)),k=0..m-1))/p^m);
 od:
 return eval(x);
end:

witt_p_add_term := proc(n::nonnegint)
 option remember;
 local x,y,u,m,rel;
 u := table():
 for m from 0 to n-1 do
  u[m] := witt_p_add_term(m)(x,y);
 od:

 rel := ghost_p_coeff(n)(u) - ghost_p_coeff(n)(x) - ghost_p_coeff(n)(y);
 u[n] := expand(solve(rel,u[n]));
 return unapply(u[n],x,y);
end:

witt_p_add := (n::posint) -> proc(x,y)
 local m;
 return table([seq(m=witt_p_add_term(m)(x,y),m=0..n-1)]);
end:

witt_p_mult_term := proc(n::nonnegint)
 option remember;
 local x,y,u,m,rel;
 u := table():
 for m from 0 to n-1 do
  u[m] := witt_p_mult_term(m)(x,y);
 od:

 rel := ghost_p_coeff(n)(u) - ghost_p_coeff(n)(x) * ghost_p_coeff(n)(y);
 u[n] := expand(solve(rel,u[n]));
 return unapply(u[n],x,y);
end:

witt_p_mult := (n::posint) -> proc(x,y)
 local m;
 return table([seq(m=witt_p_mult_term(m)(x,y),m=0..n-1)]);
end:


text:
export Improve := module ()
    export Simplifiers := table();
    export ModuleApply := proc(dom :: HDomain, $)::HDomain_mb;
        local es := map(si->Simplifiers[si]
                       , sort( [indices(Domain:-Improve:-Simplifiers,nolist) ]
                             , key=(si->Simplifiers[si]:-SimplOrder)))
             , bnd, sh;
        bnd, sh := op([1..2], dom);
        sh := foldr(((f,q)->proc() cmp_simp_sh(f,q,args) end proc)
                   ,(proc(_,b,$) b end proc), op(es))(bnd, sh);
        DOMAIN(bnd, sh);
    end proc;

    local ModuleLoad := proc($)
      LoadSimplifiers();
      LoadSimplifiers := (proc() end proc);
    end proc;#ModuleLoad

    local LoadSimplifiers := proc($)
      local lib_path, drs, improve_path, simpls_paths, simpl_path, names0, names1, new_names, nm;
      unprotect(Simplifiers);
      lib_path := LibraryTools:-FindLibrary(Hakaru);
      ASSERT(lib_path<>NULL);
      lib_path := FileTools:-ParentDirectory(lib_path);
      drs := kernelopts(dirsep);
      improve_path := `cat`(lib_path,drs,"Domain",drs,"Improve");

      simpls_paths := FileTools:-ListDirectory(improve_path);
      for simpl_path in simpls_paths do
        names0 := {anames(user)};
        read(`cat`(improve_path,drs,simpl_path));
        names1 := {anames(user)};
        new_names := names1 minus names0;
        for nm in new_names minus {entries(Simplifiers,nolist)} do
          if nm :: `module`(SimplName,SimplOrder) then
            if assigned(Simplifiers[nm:-SimplName])
            and not Simplifiers[nm:-SimplName] = nm then
              WARNING("overwriting old simplifier! (%1)", Simplifiers[nm:-SimplName]):
            end if;
            Simplifiers[nm:-SimplName] := eval(nm);
          else
            WARNING("not recognized as a simplifier: %1", nm);
          end if;
        end do;
      end do;
      protect(Simplifiers);
    end proc;#LoadSimplifiers

    local combine_errs := proc(err::DomNoSol, mb, $)
        if mb :: DomNoSol then
            DNoSol(op(err),op(mb));
        else
            mb
        end if;
    end proc;

    # compose two simplifiers, combining errors if both fail
    local cmp_simp_sh := proc(simp0, simp1, bnd, sh :: {DomShape,DomNoSol}, $)::{DomShape,DomNoSol};
      local res, sh1; sh1 := simp0(bnd, sh); simp0:-SimplName;
      if not sh1 :: DomNoSol then
          simp1(bnd, sh1);
      else
          combine_errs( sh1, simp1(bnd, sh) );
      end if;
    end proc;

    export classify_relation := proc(r0::relation, vars0, $)
      ::{identical(FAIL), [identical(B_LO,B_HI,B_EQ,B_NEQ), satisfies(q->q in {indices(flip_relation,nolist)}), name, algebraic]};
      local r_k, r_s, in_vars, vars := vars0, r := r0;
      if vars :: ({set,list})({name,list(name)}) then
        vars := map(x->`if`(x::list,op(1,x),x),vars);
        vars := {op(vars)}; in_vars := x->x in vars;
      elif vars :: DomBound then
        vars := DomBound:-Bound:-varsOf(vars,"set"); in_vars := x->x in vars;
      elif vars :: procedure then
        in_vars := vars;
      elif vars :: type then
        in_vars := x->type(x,vars);
      else error "unknown argument: %1 (expecting variables, variable membership check, "
                 "or domain bounds)", vars;
      end if;

      if r :: {`<=`,`<`, `=`, `<>`} then
        if in_vars(rhs(r)) then
          r_k := flip_relation[op(0,r)]; r_s := rhs(r), lhs(r);
        elif in_vars(lhs(r)) then
          r_k := op(0,r); r_s := op(r);
        else return FAIL end if;
        [ classify_relop[r_k], r_k, r_s ];
      else FAIL
      end if;
    end proc;
    local flip_relation := table([`=`=`=`,`<>`=`<>`,`<=`=`>=`,`<`=`>`,`>=`=`<=`,`>`=`<`]);
    local classify_relop := table([`=`=B_EQ,`<>`=B_NEQ,`<=`=B_HI,`<`=B_HI,`>=`=B_LO,`>`=B_LO]);
end module;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
read("bssn_spher.mpl");

read("built_eqns.mpl");

read("compactify.mpl");

read("perform_fd.mpl");

read("built_direct_ire_test.mpl");


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

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

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

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

params_a_gamma := (1 - log(2))/Pi^2:
params_a_BB    := 1:
$include "gga_c_pbe.mpl"

zvpbeint_nu := (rs, z, t) ->
  t*mphi(z)*(3/rs)^(1/6):

(* we write (z^2)^(omega/2) instead of z^omega in order to
   avoid the use of abs(z). Max is required not to get float
   exceptions for z->0 *)
zvpbeint_ff := (rs, z, t) ->
  exp(-params_a_alpha*zvpbeint_nu(rs, z, t)^3*m_max(z^2, 1e-20)^(params_a_omega/2)):

f  := (rs, z, xt, xs0, xs1) ->
  f_pw(rs, z) + zvpbeint_ff(rs, z, tp(rs, z, xt))*fH(rs, z, tp(rs, z, xt)):


text:
`is_element/hasse_diagrams` := (A::set) -> proc(G)
 global reason;
 local D,R;

 if not(`is_element/autorel`(A)(G)) then
  reason := ["is_element/hasse_diagrams","G is not a relation on A",G,reason];
  return false;
 fi;

 if not(`is_antisymmetric/autorel`(A)(G)) then
  reason := ["is_element/hasse_diagrams","G is not antisymmetric",G];
  return false;
 fi;

 D := `id/autorel`(A);
 R := `o/autorel`(A)(G,G);

 while R <> {} do
  if R intersect G <> {} then
   reason := ["is_element/hasse_diagrams","G has an arrow parallel to a longer path",G,R intersect G];
   return false;
  fi;

  if R intersect D <> {} then
   reason := ["is_element/hasse_diagrams","G has a loop",G,R intersect D];
   return false;
  fi;

  R := `o/autorel`(A)(R,G);
 od;

 return true;
end:

`is_element/ranked_hasse_diagrams` := (A::set) -> proc(Gr)
 local G,r,AA,rA,d,A1,r1,a,b,AA1,G1,M1,M2,Aa,A0,Ma,rMa,i;
 global reason;

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

 G,r := op(Gr);

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

 if not(type(G,set) and G minus AA = {}) then
  reason := [convert(procname,string),"G is not a subset of A^2 minus Delta",G,A];
  return false;
 fi;

 if not(type(r,table) and [indices(r)] = map(a -> [a],[op(A)])) then
  reason := [convert(procname,string),"r is not a table indexed by A",r,A];
  return false;
 fi;

 if A = {} then return true; fi;

 rA := map(a -> r[a],A);
 if not(type(rA,list(integer)) and min(op(rA)) = 0) then
  reason := [convert(procname,string),"The values of r are not integers starting from 0",r,A];
  return false;
 fi;

 d := max(rA);
 if rA <> {seq(i,i=0..d)} then
  reason := [convert(procname,string),"The values of r do not form an interval starting from 0",r,A];
  return false;
 fi;

 if d = 0 then
  if G = {} then
   return true;
  else
   reason := [convert(procname,string),"All ranks are zero but G is not empty",G];
   return false;
  fi;
 fi;

 A1 := select(a -> (r[a] < d),A);
 r1 := table();
 for a in A1 do r1[a] := r[a]; od;
 AA1 := {seq(seq([a,b],b in A1),a in A1)} minus {seq([a,a],a in A1)};
 G1 := G intersect AA1;

 if not(`is_element/ranked_hasse_diagrams`(A1)([G1,r1])) then
  return false;
 fi;

 M1 := NULL;
 for a in A1 do
  Aa := select(b -> member([a,b],G1),A1);
  if Aa = {} then M1 := M1,a; fi;
 od:
 M1 := {M1};
 M2 := select(a -> (r[a] = d-1),M1);

 A0 := A minus A1;

 for a in A0 do
  Ma := select(b -> member([b,a],G),A);
  rMa := map(b -> r[b],Ma);
  if not(member(d-1,rMa)) then
   reason := [convert(procname,string),"a has rank d but does not cover anything of rank d-1",a,d,Ma];
   return false;
  fi;
 od:
end:


text:
(*
 Copyright (C) 2021 Susi Lehtola

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

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

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

$include "gga_x_optx.mpl"
$include "gga_x_ityh.mpl"

ityh_enhancement := xs -> optx_f(xs):


text:
program TeaBreak;
begin
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    writeln('Say "BREAK!!" ');
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    ;
    ;

max_stars_repo_name: Hiroya-W/lang-processing

max_issues_repo_name: Hiroya-W/lang-processing

max_forks_repo_name: Hiroya-W/lang-processing


text:
#include "utils.mpl" // you can include mpl files
$creepAmt:10 // just like normal CRPL
$addPerTime:1
$Delay:15
$maxAmt:1000
$incrementInterval:10

amt()inline, CRPL{<-creepAmt 1 add} // also you can create CRPL functions

main()inline{ // the function that calls every tick
	/*
	recommended to be inline.
	all inline function will be implaced to their calls.
	*/
	once{
		RemoveImages(Self());
		SetImage(Self(), "main", "CustomEmitter");
		SetTextSize(0.6);
		SetTimer1(incrementInterval + 1);
	}
	if(GetTimer1() == 0){
		SetText("Amt " $ creepAmt $ endl() $ "Interval: " $ round(Delay/30.0, 2));
		SetCreeper(CurrentCoords(), amt());
		SetTimer1(Delay);
	}
	if (GetTimer0()==1 && creepAmt < maxAmt){
		creepAmt = creepAmt + addPerTime;
		SetTimer0(incrementInterval + 1);
	}
	if (GetTimer2() == 0){
		f(RandCoords());
		SetTimer2(60);
	}
	else if(creepAmt > maxAmt) once creepAmt = maxAmt;
}


text:
!Properties from TestCase_0 that are present in this test case:
! - no matrix A
! - only bin vars on the first stage
! - standard bounds for 1st stage cont vars (if any)
! - only vector rhs is stochastic

!Properties specific for this TestCase:
! all scenarios have the same RHS
! integer vars on the 2nd stage

!These properties of are present in all TestCases:
! - MIN obj - always the case for smps files
! - 2 stages - only such problems are in scope of the research
! - only cont vars on the 2nd stage - we will solve the relaxation otherwise
! - >= constraints in subproblems - we do not care about subproblems
! - bin and cont (if any) vars on the first - int vars are prohibited

! Solution value: -17600
! Solution value 2nd stage relaxation: -17600.4
! Solution value all stages relaxation: -17900.1

TITLE
	SNDP_TestCase_14;

STOCHASTIC

SCENARIO
	SCEN := 1..3;

PROBABILITIES
	Prob[SCEN] := (0.25, 0.5, 0.25);

RANDOM DATA
	Demand[SCEN] := (5000, 5000, 5000);

STAGE1 INTEGER VARIABLES
	yc2;
	yc3;

STAGE2 VARIABLES
	x_c24;
	x_c34;
	x_a12;
	x_b12;
	x_a13;

STAGE2 INTEGER VARIABLES
	x_b13;
	x_a32;
	x_a33;

MIN	! should be min for solver
	- 13*(x_c24 + x_c34) + 2*x_a12 + 2*x_b12 + 2*x_a13 + 2*x_b13 + 1*x_a32 + 2*x_c24 + 3*x_c34 + 2000*yc2 + 2000*yc3;

SUBJECT TO
	x_a12 + x_a32 >= 2*x_c24;
	-x_a12 + -x_a32 >= -2*x_c24;

	x_a13 + x_a33 >= 2*x_c34;
	-x_a13 + -x_a33 >= -2*x_c34;

	x_b12 >= 3*x_c24;
	-x_b12 >= -3*x_c24;

	x_b13 >= 3*x_c34;
	-x_b13 >= -3*x_c34;

	-x_c24 >= -5000*yc2;
	-x_c34 >= -5000*yc3;
	-x_c24 + -x_c34 >= -Demand;

	-x_b13 + -x_a33 >= -24500.5; ! in the original problem these vars get 15K and 10K values

BOUNDS
	!x_b13 <= 100000;  ! integer should have UB to be treated as not binary in MPS
	!x_a32 <= 100000;
	!x_a33 <= 100000;
	yc2 <= 1;
	yc3 <= 1;
END

max_stars_repo_name: pashtetGP/optconvert

max_issues_repo_name: pashtetGP/optconvert

max_forks_repo_name: pashtetGP/optconvert


text:
# # Comuting BGG operator and characteristic homomorphism
# # by Shizuo Kaji
# # version: 2007/12/10 12:30
# # See the included demo maple worksheet for usage
# # LIMITATION: due to computational complexity,
# #       for type-E this script deal with modulo
# #         maximal parabolic P2 by default

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

# Data files
# E8data.m: load reduced words of E8 up to length 15, into a variable "redu"
# E8data-P2.m: load reduced words of E8/P2, into a variable "redu"
# E8gen.mpl: load definition for ring generators of E8/T

# global variables:
# R: Lie type
# B: simple roots
# redu[deg]: reduced words
# reduw: flattened list of reduced words
# w[i]: i-th weight
# a[i]: i-th root
# s[[]]: weyl group element
# X[weyl]: Schubert class corresponding to weyl
# S[weyl]: Schubert polynomial corresponding to weyl
# gam[deg],x[deg]: special generator
# rho[deg]: relation

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

# convert poly of weights into that of e
w2e := proc (wpol::polynom) global R,B;
    return eval(wpol, [seq(w[j]=weights(R)[j], j=1..nops(B))]);
end proc;

# convert poly of special generators to that of weights "w[1..r]"
gen2w := proc (f::polynom) global R,B;
    return eval(f, [seq(G[i] = gam[i], i = [3,4,5,6,9,10,15]),
		seq(c[i] = c(i), i = 1 .. nops(B)),
		seq(T[i] = t[i], i = 0 .. 8)]);
end proc;

# listup monomials of degree "deg" consisting of "gens", whose degrees are specified by "degs"
list_monom := proc (gens::list, degs::list, deg::integer) local i, L;
 L := {};
 for i to nops(gens) do
   if 0 < deg-degs[i] then
     L := `union`(L, map2(`*`, gens[i], list_monom(gens, degs, deg-degs[i])))
   elif deg-degs[i] = 0 then
     L := `union`(L, {gens[i]});
   end if;
 end do;
 return L;
end proc:

# BGG operator for w-poly wrt any root
BGG_root := proc (vec::linear, f::polynom) global R;
    return -simplify((f - subs([seq(w[i]=e2w(reflect(vec, weights(R)[i])),i=1..rank(R))], f))
	      / e2w(vec));
end proc;

# BGG for w-poly by recursive computation using Leibniz-rule
BGG_w := proc (rt::integer, f::polynom) local i, j, term;
	if degree(f) < 1 then return 0
    elif type(f, monomial) then
		return -sum(ref_weight(rt)^i*w[rt]^(-i-1)*f, i=0..degree(f, w[rt])-1);
    elif op(0, f) = `+` then
		return expand(map2(procname, rt, f));
    elif op(0, f) = `*` then
		if type(op(1,f), 'rational') then
			return op(1,f)*procname(rt, f/op(1,f));
		end if;
		j := 1;
		for i to nops(f) do
		    if degree(op(i, f)) <= degree(op(j, f)) and length(op(i, f)) < length(op(j, f)) then
			j := i
		    end if
		end do;
		term[1] := op(j, f);
		term[2] := mul(op(i, f), i = 1 .. j-1)*mul(op(i, f), i = j+1 .. nops(f));
		return procname(rt, term[1])*term[2]+ref_rt(rt, term[1])*procname(rt, term[2]);
    elif op(0, f) = `^` then
	    term[1] := op(1, f);
	    return procname(rt, term[1])*add(term[1]^(i-1)*ref_rt(rt, term[1])^(op(2,f)-i),i=1..op(2,f));
    else
        print("cannot compute BGG:",f);
        return infinity;
    end if
end proc;

# BGG operator for e-poly
BGG_e := proc (rt::integer, f::polynom) global R;
    return -simplify((f - subs([seq(e(i)=reflect(B[rt], e(i)),i=1..dimr(R))], f))
	      / B[rt]);
end proc;

# a simple reflection on fundamental weight
ref_weight := proc (rt::integer) option remember; global R,B;
	return w[rt]-e2w(B[rt]);
end proc;

# simple reflection on Schubert cycle
ref_weyl := proc(rt::integer,weyl::list) local a,j,v,u,w,temp; global B; option remember;
	temp := X[weyl];
	w := weyl.[rt];
    if nops(w) > nops(weyl) then
		return X[weyl];
	end if;
	for a in pos_roots(B) do
       	 v := reflect(a,interior_pt(B));
     	 v := reflect(seq(B[j],j=w),v);
	  	  vec2fc(v,R,'u');    # u = w.s_a
		  if nops(u)=nops(weyl) then
			   temp := temp - iprod(2/iprod(a,a)*a,B[rt])*X[u];
          end if;
     end do;
     return temp;
end proc:

# simple reflection on polynomial
ref_rt := proc (rt::integer, f::polynom) local ff;
    if type(f, 'rational') then
        return f;
    else
        ff:=subs(w[rt]=ref_weight(rt), f);
        if has(f,'X') then
            if type(ff, 'indexed') and op(0,ff)='X' then return ref_weyl(rt,op(1,ff));
            elif op(0,ff) =`+` or op(0,ff)=`*` then return map2(procname, rt, ff);
            elif op(0,ff) =`^` then return procname(rt,op(1,ff))^op(2,ff);
            end if;
        end if;
        return ff;
    end if;
end proc;

# weyl group action on polynomial
ref := proc(weyl, f::polynom);
    return WeylFunc(weyl,f,ref_rt);
end proc:

# Schubert basis expansion
Char := proc (f::polynom, func::procedure, roots::set, R) local weyl;
    return add( WeylFunc(weyl, f, func)*X[weyl],
	       weyl=PReducedWd(R, roots, degree(f))[degree(f)]);
end proc;

# Schubert basis expansion for poly in special generators
Char_w := proc (f::polynom) local temp, deg; global redu; global R;
	temp:=collect(gen2w(f), {seq(w[i], i=1..rank(R))});
    if not type(temp, polynom) then
        print('not polynomial', f); return -infinity;
    end if;
    deg := degree(temp);
    if deg < 1 then return temp;
	elif deg <> ldegree(temp) then
		print("not homogeneous:", temp); return -infinity;
    elif deg>nops(redu) then return 0;
    end if;
	return add(WeylFunc(weyl, temp, BGG_w)*X[weyl], weyl=redu[deg]);
# for multi thread
#	return Threads:-Add( WeylFunc(weyl, temp)*X[weyl], weyl=redu[deg] );
end proc;

# Schubert basis expansion for poly in "X[weyl]"
Char_X := proc(f::polynom) local weyl,deg; global redu;
    deg := deg_X(f);
    if deg < 1 then return f;
    elif deg>nops(redu) then return 0;
    end if;
  return add(WeylFunc(weyl, f, BGG_X)*X[weyl], weyl=redu[deg]);
end proc:

# Schubert basis expansion for poly in "e"
Char_e := proc(f::polynom) local weyl,deg; global redu;
    deg := degree(f);
    if deg <1 then return f;
	elif deg <> ldegree(f) then
		print("not homogeneous:", f); return -infinity;
    elif deg>nops(redu) then return 0;
    end if;
  return add(WeylFunc(weyl, f, BGG_e)*X[weyl], weyl=redu[degree(f)]);
end proc:

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

# top (point) cohomology class
top := proc (roots::set := {}) local k, j, i, RP, temp, lie_type; global R,P,reduw;
	if not roots={} then
       RP := relative_roots(R,roots);
	   return (-1)^(nops(RP))*mul(i,i=RP)/nops(reduw);
    end if;
    temp := 1;
	lie_type := name_of(R);
    if lie_type = F4 then
        return (-255887*gam[3]^4*gam[4]^3+178123*gam[4]^6)/173408256000;
    elif lie_type = E7 then
		temp := mul(t[i]-t[0], i = 1 .. 7);
		for i to 6 do temp := temp*mul(t[i]-t[j], j = i+1 .. 7) end do;
	    for i to 5 do
		 for j from i+1 to 6 do
			temp := temp*mul(t[i]+t[j]+t[k]-t[0], k = j+1 .. 7);
	     end do;
		end do;
		return temp;
    elif lie_type = E8 then
		temp := mul(t[i], i = 1 .. 8);
		for i to 7 do temp := temp*mul(t[i]-t[j], j = i+1 .. 8) end do;
		for i to 7 do temp := temp*mul(t[i]+t[j]-t[0], j = i+1 .. 8) end do;
		for i to 6 do
			for j from i+1 to 7 do
				temp := temp*mul(t[i]+t[j]+t[k]-t[0], k = j+1 .. 8)
			end do;
		end do;
		return temp;
    elif lie_type = G2 then
		return (-9*w[1]^6+6*w[1]^5*w[2])/12;
    elif lie_type = cat('A',rank(R)) then
	     return mul(e(i)^(rank(R)-i+1), i = 1 .. rank(R));
	elif lie_type = cat('B',rank(R)) then
		 return (-1)^rank(R)*mul(e(i)^(2*i-1), i = 1 .. rank(R))/2^rank(R);
    elif lie_type = cat('C',rank(R)) then
		 return (-1)^rank(R)*mul(e(i)^(2*i-1), i = 1 .. rank(R));
    else
		return (-1)^(nops(P))*mul(i, i = P)/mul(i,i=degrees(R));
    end if
end proc;

# enumerate Schubert polynomials by applying BGG to the top degree class
schubertpolys := proc() global S, R, reduw; local w,func,top_class,top_weyl;
  top_class := top();
  top_weyl  := longest_elt(R);
  S[[]] := 1;
  if evalb(name_of(R) in {'G2','F4','E6','E7','E8'}) then
  	 func:='BGG_w';
  else
	func:='BGG_e';
  end if;
  for w in reduw do
  	   S[w] := WeylFunc(inv(w).top_weyl, top_class, func);
  end do;
  print("The top class is", top_weyl, top_class);
end proc:

# polynomial representatives in given generators for Schubert classes in a given degree
Schubert_Polynomials := proc (gens::list, deg::integer) local i,j,v,temp,M,L,A,V,degs;
 L:=[]; A:=[];
 degs := map(degree, map(gen2w, gens));
 M := convert(list_monom(gens, degs, deg), `list`);
 for i from 1 to nops(M) do
  temp := Char_w(M[i]);
  L:=[op(L), M[i] = temp];
  v := vector(nops(redu[deg])+nops(M), 0);
  v[nops(redu[deg])+i] := 1;
  for j to nops(redu[deg]) do
     v[j] := coeff(temp, X[redu[deg][j]]);
  end do;
  A := [op(A), convert(v, list)];
 end do;
 A := LinearAlgebra[HermiteForm](Matrix(A));
# for printing
 V := [];
 for v in convert(A, `listlist`) do
	V := [op(V),
	    add(v[i]*X[redu[deg][i]], i=1..nops(redu[deg]))
        = add(v[nops(redu[deg])+i]*M[i], i=1..nops(M))];
  end do;
  return V;
end proc:

## initialization for Lie types (using Toda-Watanabe basis for exceptional types)
setup := proc (lie_type, roots::set := {}, upto::integer := 0) local i,S;
		global R, P, B, x, u, z, t, udim, gam, rho, redu, reduw;
  forget(ref_weight);
  forget(ref_weyl);
  unassign('gam','t','x','S','rho');

# type G2
  if lie_type = G2 then
    R := lie_type;
    t[1] := -w[1];
    t[2] := -w[1]+w[2];
    t[3] := 2*w[1]-w[2];
    rho[2] := w[2]^2+3*w[1]^2-3*w[1]*w[2];
    gam[3] := (1/2)*t[1]*t[2]*t[3];
	x[3] := 1/2*(w[1]*w[2]^2-w[1]^2*w[2]); # [2,1,2]
# type F4
  elif lie_type = F4 then
    R := linalg[transpose](cartan_matrix(lie_type));
    t[1] := -w[4];
    t[2] := w[1]-w[4];
    t[3] := -w[1]+w[2]-w[4];
    t[4] := -w[2]+2*w[3]-w[4];
    t[0] := w[3]-2*w[4];
    gam[3] := (1/2)*c(3);
    gam[4] := (1/3)*c(4)-(2/3)*t[0]*c(3)+(8/3)*t[0]^4;
	rho[2] := c(2)-2*t[0]^2;
	rho[6] := gam[3]^2-3*t[0]^2*gam[4]-4*t[0]^3*gam[3]+8*t[0]^6;
	rho[8] := 3*gam[4]^2+6*t[0]*gam[3]*gam[4]-3*t[0]^4*gam[4]-13*t[0]^8;
# type E6 / P2
  elif lie_type = E6 then
    R := linalg[transpose](cartan_matrix(lie_type));
    t[1] := -w[1]+w[2];
    t[2] := w[1]+w[2]-w[3];
    t[3] := w[2]+w[3]-w[4];
    t[4] := w[4]-w[5];
    t[5] := w[5]-w[6];
    t[6] := w[6]; t[0] := w[2];
    gam[3] := c(3);
    gam[4] := c(4)+2*t[0]^4;
# type E7 / P2
  elif lie_type = E7 then
    R := linalg[transpose](cartan_matrix(lie_type));
    t[1] := -w[1]+w[2];
    t[2] := w[1]+w[2]-w[3];
    t[3] := w[2]+w[3]-w[4];
    t[4] := w[4]-w[5];
    t[5] := w[5]-w[6];
    t[6] := w[6]-w[7];
    t[7] := w[7];
    t[0] := w[2];
    gam[3] := 1/2*c(3);
    gam[4] := 1/3*(c(4)+2*t[0]^4);
    gam[5] := 1/2*(c(5)-t[0]*c(4)+t[0]^2*c(3)-2*t[0]^5);
    gam[9] := 1/2*(c(6)*c(3)+t[0]^2*c(7)-3*t[0]^3*c(6));
# type E8 / P2
  elif lie_type = E8 then
    R := linalg[transpose](cartan_matrix(lie_type));
	read "E8data-p2.m";
	read "E8gen.mpl";
  else
    R := lie_type;
  end if;

  B := base(R);
  P := pos_roots(lie_type);
  udim := `if`(upto>0,upto,nops(P));
  if not lie_type = E8 then
      redu := PReducedWd(R, roots, udim);
  end if;
  reduw := [[],seq(op(i),i=redu)];
  print(lie_type, "mod the parabolic w.r.t", roots);
  return diagram(R);
end proc;

max_stars_repo_name: shizuo-kaji/Maple-flag-cohomology

max_issues_repo_name: shizuo-kaji/Maple-flag-cohomology

max_forks_repo_name: shizuo-kaji/Maple-flag-cohomology


text:
(*
 Copyright (C) 2020 Susi Lehtola

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

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

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

$include "mgga_x_rscan.mpl"
$include "mgga_x_scan.mpl"

(* eqn S6 *)
r2scan_alpha := (x, t) -> (t - x^2/8)/(K_FACTOR_C + params_a_eta*x^2/8):

(* f(alpha) replaced with a polynomial for alpha in [0, 2.5], eqn S7 *)
r2scan_f_alpha_neg := a -> exp(-params_a_c1*a/(1 - a)):
r2scan_f_alpha := (a, ff) -> my_piecewise5(a <= 0, r2scan_f_alpha_neg(m_min(a, 0)), a <= 2.5, rscan_f_alpha_small(m_min(a, 2.5), ff), rscan_f_alpha_large(m_max(a, 2.5))):

(* eqn S11 *)
Cn := 20/27 + params_a_eta*5/3:
(* eqn S12 *)
C2 := ff -> -add(i*ff[9-i], i=1..8) * (1-scan_h0x):

(* eqn S10; this is analogous to scan_y *)
r2scan_x := (p, ff) -> (Cn*C2(ff)*exp(-p^2/params_a_dp2^4)+MU_GE)*p:

r2scan_f := (x, u, t) -> (scan_h1x(r2scan_x(scan_p(x), rscan_fx)) + r2scan_f_alpha(r2scan_alpha(x, t), rscan_fx) * (scan_h0x - scan_h1x(r2scan_x(scan_p(x), rscan_fx))))*scan_gx(x):

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


text:
############################################################
# Generic conformal 3 metric for BSSN in spherical symmetry
# See: BSSN_Spher.mw or BSSN_Spher.mpl
############################################################

Ndim_ := 3:
x1_ := x:
x2_ := theta:
x3_ := phi:
complex_ := {}:
g11_ := (A(t,x)):
g12_ := 0:
g13_ := 0:
g22_ := x^2*(B(t,x)):
g23_ := 0:
g33_ := x^2*(B(t,x))*(sin(theta))^2:


text:
# Vierachsroboter implizit
# Initialisierung
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env_IC":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name_OL):
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name_OL, base_method_name);
read "../robot_codegen_definitions/robot_env_IC":
Trf := Trf:
Trf_c := Trf_c:
Trf
;
# VGK Gelb 2-7-11-14-6
# Schleife (1-6)-(6-14)
T_1_14 := combine( Matrix(Trf(1..4,1..4, 6)) . Matrix(Trf(1..4,1..4,14))):
# Schleife (1-2)-(2-7)-(7-11)
T_1_11:= combine( Matrix(Trf(1..4,1..4, 2)) . Matrix(Trf(1..4,1..4,7)) . Matrix(Trf(1..4,1..4,11))):
h1t := T_1_14(1..3,4) - T_1_11(1..3,4);

#tmp := Transpose( Matrix(T_1_14(1..3,1..3)) ) . Matrix(T_1_11(1..3,1..3));  nur anzeigen lassen für h1r
;
#combine(tmp);# nur anzeigen lassen für h1r
;
h1r := -(-qJ6(t)+qJ2(t)+qJ7(t)+phi711+qJ11(t))+Pi/2:
# VGK PINK 2-3-12-15-9-8-2
# Schleife (2-3)-(3-12)
T_2_12 := combine(Matrix(Trf(1..4,1..4, 3)) . Matrix(Trf(1..4,1..4,12))):
# Schleife (2-8)-(8-9)-(9-15)
T_2_15:= combine(Matrix(Trf(1..4,1..4, 8)) . Matrix(Trf(1..4,1..4,9)).Matrix(Trf(1..4,1..4, 15))):
h2t := T_2_12(1..3,4) - T_2_15(1..3,4):
tmp := Transpose( Matrix(T_2_12(1..3,1..3)) ) . Matrix(T_2_15(1..3,1..3));  # nur anzeigen lassen für h1r
;
 combine(tmp); # nur anzeigen lassen für h1r
;
h2r:=(qJ3(t)+phi312+qJ12(t)-qJ8(t)-qJ9(t))+3*Pi/2:
# VGK GRÜN 2-3-4-13-16-10-7-2
# Schleife (2-3)-(3-4)-(4-13)
T_2_13 := combine( Matrix(Trf(1..4,1..4, 3)) . Matrix(Trf(1..4,1..4,4)) . Matrix(Trf(1..4,1..4,13))):
# Schleife (2-7)-(7-10)-(10-16)
T_2_16:= combine( Matrix(Trf(1..4,1..4, 7)) . Matrix(Trf(1..4,1..4,10)).Matrix(Trf(1..4,1..4,16)) ):
h3t := T_2_13(1..3,4) - T_2_16(1..3,4):
#tmp := Transpose( Matrix(T_2_13(1..3,1..3)) ) . Matrix(T_2_16(1..3,1..3));  nur anzeigen lassen für h2r
;
#combine(tmp);  nur anzeigen lassen für h3r
;
h3r:=-(qJ3(t)+qJ4(t)+phi413+qJ13(t)-qJ7(t)+phi710-qJ10(t))+5*Pi/2:
# Zusammenstellen aller Zwangsbedingungen
implconstr_t := <h1t([1, 3],1);h2t([1, 2],1);h3t([1, 2],1); h1r; h2r; h3r>; # TODO: In h1r, h2r, h3r muss das richtige drinstehen.
implconstr_s := convert_t_s(implconstr_t):
# Exportiere Code für folgende Skripte
kin_constraints_exist:=true:
save implconstr_t, implconstr_s, kin_constraints_exist, sprintf("../codeexport/%s/tmp/kinematic_constraints_implicit_maple.m", robot_name):
# Exportieren des vollständigen Ausdruckes
if codegen_act then
  MatlabExport(implconstr_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_matlab.m", robot_name), 2):
end if:
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( implconstr_s )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_matlab.m", robot_name), 2):

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


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

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

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

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

AA :=  0.757211:
BB := -0.106364:
CC := -0.118649:
DD :=  0.609650:

fH := s -> add(params_a_a[i]*s^(1+i), i=1..6)/(1 + add(params_a_b[i]*s^i, i=1..9)):

mzeta   := s -> s^2*fH(s):
meta    := s -> AA + mzeta(s):
mlambda := s -> DD + mzeta(s):
mchi    := (rs, z, s) -> nu(rs, z)/sqrt(mlambda(s) + nu(rs, z)^2):

fF := (rs, z, s) ->
  1 - s^2/(27*CC*(1 + s^2/4)) - mzeta(s)/(2*CC):

fG := (rs, z, s) ->
  - 2/5  * CC*fF(rs, z, s)*mlambda(s)
  - 4/15 * BB*mlambda(s)^2
  - 6/5  * AA*mlambda(s)^3
  - mlambda(s)^(7/2)*(4/5*sqrt(Pi) + 12/5*(sqrt(mzeta(s)) - sqrt(meta(s)))):

f1 := (rs, z, s) ->
   + AA
   - 4/9 * BB*(1 - mchi(rs, z, s))/mlambda(s)
   - 2/9 * CC*fF(rs, z, s)*(2 - 3*mchi(rs, z, s) + mchi(rs, z, s)^3)/mlambda(s)^2
   - 1/9 * fG(rs, z, s)*(8 - 15*mchi(rs, z, s) + 10*mchi(rs, z, s)^3 - 3*mchi(rs, z, s)^5)/mlambda(s)^3
   + 2*nu(rs, z)*(sqrt(mzeta(s) + nu(rs, z)^2) -  sqrt(meta(s) + nu(rs, z)^2))
   + 2*mzeta(s)*log((nu(rs, z) + sqrt(mzeta(s) + nu(rs, z)^2))/(nu(rs, z) + sqrt(mlambda(s) + nu(rs, z)^2)))
   - 2*meta(s)*log((nu(rs, z) + sqrt(meta(s) + nu(rs, z)^2))/(nu(rs, z) + sqrt(mlambda(s) + nu(rs, z)^2))):

f_lda := (rs, z) -> -X_FACTOR_C*RS_FACTOR*((1 + z)/2)^(4/3)/rs:

f := (rs, z, xt, xs0, xs1) ->
  f_lda(rs, z)*f1(rs, z, X2S*xs0) + f_lda(rs, -z)*f1(rs, -z, X2S*xs1):


text:
var $sconst1 <struct {
        @f11 i32,
        @f22:15 i32,
        @f33 i32,
        @f44 i32,
        @f55 f32,
        @f66 f32}> = [ 6=10.07f, 2=699, 3=-273, 4=75, 1=31425926, 5=6.02e+23f ]
func &printf (var %p1 <* i8>)void
func &main ( ) i32 {
   call &printf (addrof a32 $sconst1)
   return (constval i32 0) }

 # 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:
# Proof by PDF that X <~ rayleigh(1) == X^2 <~ weibull(1,1)

PDF_rayleigh := (a, x) -> 2*x*exp(-x^2/a)/a;

PDF_transform := (x) -> PDF_rayleigh(1,sqrt(x))*diff(sqrt(x),x);
simplify(PDF_transform(x));

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

evalb(PDF_weibull(1, 1, x) = PDF_transform(x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


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

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

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

 return true;
end;

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

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

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

 return true;
end;

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

`random_element/real_functions` := (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/real_functions` := NULL;
`count_elements/real_functions` := NULL;

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

`is_nonnegative/real_functions` := (A::set) -> proc(x)
 local a;

 for a in A do
  if is(x[a] < 0) then
   return false;
  fi;
 od;

 return true;
end:

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

`is_leq/real_functions` := (A::set) -> proc(x,y)
 local a;

 for a in A do
  if is(x[a] > y[a]) then
   return false;
  fi;
 od;

 return true;
end:

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

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

 n := 0;
 for a in A do n := n + x[a]^2; od;

 return sqrt(n);
end:

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

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

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

 return u;
end:

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

`max/real_functions` := (A::set) -> proc(x)
 local a,u;

 u := -infinity;
 for a in A do u := max(u,x[a]); od;

 return u;
end:

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

`min/real_functions` := (A::set) -> proc(x)
 local a,u;

 u := infinity;
 for a in A do u := min(u,x[a]); od;

 return u;
end:

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

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

 n := 0;
 for a in A do n := n + (x[a]-y[a])^2; od;

 return sqrt(n);
end:

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

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

 d := 0;
 for a in A do d := d + x[a]*y[a]; od;

 return d;
end:

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

`width/real_functions` := (A::set) -> proc(x)
 local V;
 V := map(a -> x[a],A);
 return max(op(V)) - min(op(V));
end:

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

`gap/real_functions` := (A::set) -> proc(x)
 local V,n,i;
 V := sort([op(map(a -> x[a],A))]);
 n := nops(V);
 return max(0,seq(V[i+1]-V[i],i=1..n-1));
end:

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

`relative_gap/real_functions` := (A::set) -> proc(x)
 return `gap/real_functions`(A)(x)/`width/real_functions`(A)(x);
end:


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

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

(* type: gga_exc *)

$include "gga_x_pw86.mpl"

params_a_aa := 2.208:
params_a_bb := 9.27:
params_a_cc := 0.2:

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


text:
var s : string;
print "Write something: ";
read s;
print "You wrote: ";
print s;
print "\n";

max_stars_repo_name: nicohi/mini-pl-interpreter

max_issues_repo_name: nicohi/mini-pl-interpreter

max_forks_repo_name: nicohi/mini-pl-interpreter


text:
var $nan_f32 <f32> = nanf
var $nan_f64 <f64> = nan
var $inf_f32 <f32> = -inff
var $inf_f64 <f64> = inf
var $g_f32 <f32> = -1.2f
var $g_f64 <f64> = 1.2
func $add_f32 () f32 {
  return (
    add f32(dread f32 $g_f32, constval f32 inff))}
func $add_f64 () f64 {
  return (
    add f64(dread f64 $g_f64, constval f64 -inf))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
# Kinematik-Berechnung 3. Arm KAS5
# Beschreibung
#
# Modell m5:
# * Kurbel N7 ist starr mit Körper K2 gekoppelt (über die Achse) (wie m3),
# * Zahnradkopplung zwischen Körpern Z2 und Z3 (Unterschied zu m3)
#
# Setze die Zwangsbedingung (Zahnradkopplung) in die bereits berechneten Zwangsbedingungen für das KAS5m3 ein
#
# Berechnung der Kopplung der Zahnräder am Ellenbogen
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-02
# Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(LinearAlgebra):
with(StringTools): # Für Zeitausgabe
;
with(codegen):
with(CodeGeneration):
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
codegen_act := true:
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Lese Ergebnisse der Kinematik von KAS5m3 aus.
# Siehe KAS5m3_sym_codegen_kinematic_constraints.mw
# kintmp_qs, kintmp_qt, lpar_qs, lpar_qt, kintmp_subsexp
read  sprintf("../codeexport/KAS5m5/tmp/kinematic_constraints_maple_inert.m", robot_name):
kintmp_qs_KAS5m5 := kintmp_qs:
kintmp_qt_KAS5m5 := kintmp_qt:
lpar_qs_KAS5m5 := lpar_qs:
lpar_qt_KAS5m5 := lpar_qt:
kintmp_subsexp_KAS5m5 := kintmp_subsexp:
read  sprintf("../codeexport/KAS5m5/tmp/tree_floatb_definitions", robot_name):
kintmp_s_KAS5m5 := kintmp_s:
kintmp_t_KAS5m5 := kintmp_t:
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
kintmp_s := kintmp_s:
kintmp_t := kintmp_t:
kintmp_qs := kintmp_s:
kintmp_qt := kintmp_t:
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:
printf("Beginn der Berechnungen. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):
st := time():
# Übernehme die Ersetzungsausdrücke von KAS5m5
#
for i from 1 to RowDimension(kintmp_s_KAS5m5) do
  for j from 1 to RowDimension(kintmp_s) do
    if kintmp_s_KAS5m5(i,1) = kintmp_s(j,1) then
      # printf("KAS5m5 %d = KAS5m7 %d\n", i, j):
      kintmp_qs(j,1) := kintmp_qs_KAS5m5(i,1):
      kintmp_qt(j,1) := kintmp_qt_KAS5m5(i,1):
    end if:
  end do:
end do:
for i from 1 to RowDimension(kintmp_subsexp_KAS5m5) do
  for j from 1 to RowDimension(kintmp_subsexp) do
    if kintmp_subsexp_KAS5m5(i,1) = kintmp_subsexp(j,1) then
      # printf("KAS5m5 %d = KAS5m7 %d\n", i, j):
      kintmp_subsexp(j,2) := kintmp_subsexp_KAS5m5(i,2):
    end if:
  end do:
end do:
# Hierbei werden die Definitionen für kintmp_s und kintmp_t überschrieben, da die Elemente in kintmp_qs und kintmp_qt drin stehen.
# Umgehung des Problems: Definitionen neu laden.
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
kintmp_s := kintmp_s:
kintmp_t := kintmp_t:
interface(rtablesize=100):
# Korrigiere die Bedingungen. Dort wurden die konstanten Winkel auch mit substitutierten Variablen "s" bezeichnet, was eigentlich keinen Sinn ergibt.
# Die mit Suffix "s" gekennzeichneten substituierten Winkel beziehen sich nur auf Zeitableitungen, die mit einem konstanten Ausdruck ersetzt wurden.
# Damit nicht delta8 und delta8s als unterschiedliche Variablen auftauchen, werden diese ersetzt.
subs_kp1 := <delta8s; delta9s; delta10s; delta12s; delta17s; delta18s>:
subs_kp2 := <delta8;  delta9;  delta10;  delta12;  delta17;  delta18>:
for i from 1 to RowDimension(subs_kp1) do
  for j from 1 to RowDimension(kintmp_subsexp) do
    kintmp_subsexp(j,2) := subs({subs_kp1(i,1)=subs_kp2(i,1)}, kintmp_subsexp(j,2));
  end do:
  for j from 1 to RowDimension(kintmp_qs) do
    kintmp_qs(j,1) := subs({subs_kp1(i,1)=subs_kp2(i,1)}, kintmp_qs(j,1));
    kintmp_qt(j,1) := subs({subs_kp1(i,1)=subs_kp2(i,1)}, kintmp_qt(j,1));
  end do:
  lpar_qs_KAS5m5 := subs({subs_kp1(i,1)=subs_kp2(i,1)}, lpar_qs_KAS5m5);
end do:
# Zahnrad Abwälzbedingung (delta19)
# Siehe Aufzeichnungen vom 3.6.2016
kintmp_s(24);
kintmp_subsexp(47,2);
kintmp_subsexp(48,2);
delta19_qs := convert_t_s(theta(4)):
kintmp_qs(24) := delta19_qs:
kintmp_subsexp(47,2) := sin(delta19_qs):
kintmp_subsexp(48,2) := cos(delta19_qs):
# Zahnrad Abwälzbedingung (rho6)
# Siehe Aufzeichnungen vom 3.6.2016
kintmp_s(28);
kintmp_subsexp(55,2);
kintmp_subsexp(56,2);
rho6_qs := convert_t_s(theta(5)) - convert_t_s(theta(4)) + delta18:
kintmp_qs(28) := rho6_qs:
kintmp_subsexp(55,2) := sin(rho6_qs):
kintmp_subsexp(56,2) := cos(rho6_qs):
# Berechne delta7
# Gl. (30)
kintmp_s(12);
kintmp_subsexp(23,2);
kintmp_subsexp(24,2);
delta7_qs := delta17 + delta19_qs: # Geändert ggü. KAS5
;
kintmp_qs(12) := delta7_qs:
kintmp_subsexp(23,2) := sin(delta7_qs):
kintmp_subsexp(24,2) := cos(delta7_qs):

# Dummy-Einträge für Schnitt-Gelenke
# Die Winkel der Schnittgelenke sind egal, daher keine Berechnung.

# delta 14
#
kintmp_s(19);
kintmp_subsexp(37,2);
kintmp_subsexp(38,2);
kintmp_qs(19) := 0:
# delta 21
#
kintmp_s(26);
kintmp_subsexp(51,2);
kintmp_subsexp(52,2);
kintmp_qs(26) := 0:

# Setze Federlänge als Variable
kintmp_s(29);
kintmp_qs(29) := lpar_qs_KAS5m5:
# Nachverarbeiten
#
kintmp_qt := convert_s_t(kintmp_qs):
# Speichern der Ergebnisse
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
# Exportieren der Zwangsbedingungen
# Matlab-Funktionen für ZB
if codegen_act then
  MatlabExport(kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematik_parallel_kintmp_subsexp_matlab_opt1.m", robot_name), 1):
end if:
printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste der Parameter speichern
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp(..,2) )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m", robot_name), 2):
printf("Zwangsbedingungen der Parallelstruktur von KAS5m5 nach KAS5m7 angepasst. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):


text:
program Write;
begin
    writeln('It''s OK?')
end.


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

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

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

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

sogga11_alpha := params_a_mu*X2S*X2S/params_a_kappa:

sogga11_f0 := x -> 1 - 1/(1 + sogga11_alpha*x^2):
sogga11_f1 := x -> 1 - exp(-sogga11_alpha*x^2):

sogga11_f  := x -> add(params_a_a[i]*sogga11_f0(x)^(i-1), i=1..6) + add(params_a_b[i]*sogga11_f1(x)^(i-1), i=1..6):

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


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

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

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

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

mpbe_f0 := s -> s^2/(1 + params_a_a*s^2):
mpbe_f := x -> 1
  + params_a_c1*mpbe_f0(X2S*x)
  + params_a_c2*mpbe_f0(X2S*x)^2
  + params_a_c3*mpbe_f0(X2S*x)^3:

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


text:
input := StringTools:-Trim(FileTools:-Text:-ReadFile("AoC-2021-15-input.txt"));

grid := Matrix(map(s->parse~(StringTools:-Explode(s)),StringTools:-Split(input,"\n"))):

(gridw,gridl) := upperbound(grid);

nbs := proc(i,j,gridsize)
local out := NULL;
    if i < gridsize then
       out := out, [i+1, j];
    end if;
    if j < gridsize then
       out := out, [i, j+1];
    end if;

    if i > 1 then
      out := out, [i-1,j];
    end if;
    if j > 1 then
       out := out, [i, j-1];
    end if;

    return [out];
end proc:

Dijkstra := proc(start, theend, grid)
local pq, costsofar, curr, nb, newcost, pt, gridsize;
uses priqueue;

gridsize := upperbound(grid)[1];
costsofar := table(sparse=infinity):

initialize(pq):
insert([0,start], pq):
costsofar[start] := 0:

while pq[0] <> 0 do
   curr := extract(pq)[2];
   if curr = theend then
       return costsofar[theend];
   end if;
   nb := nbs(curr[], gridsize);
   for pt in nb do
       newcost := costsofar[curr] + grid[pt[]];
       if newcost < costsofar[pt] then
           costsofar[pt] := newcost;
           insert([-newcost,pt], pq);
       end if;
   end do;
end do:

end proc:

answer1 := Dijkstra([1,1], [gridl,gridw], grid);
# should really do this lazily
inc := a -> (a mod 9) +~ 1;
biggrid:= Matrix([seq([seq((inc@@(i+j-2))(grid),j=1..5)],i=1..5)] );
(biggridw,biggridl) := upperbound(biggrid);
answer2 := Dijkstra([1,1], [biggridl,biggridw], biggrid);