text:
# Calculate Jacobi-Matrix based on MDH-Frames
# Introduction
# Jacobi-Matrix für Roboterkinematik basierend auf MDH-Parametern berechnen.
# Die Berechnung erfolgt durch symbolische Ableitung der Transformationsmatrizen für jeden Körper
#
# robot -> Berechnung für allgemeinen Roboter
# tree -> Berechnung für eine beliebige Baumstruktur (ohne Schleifen)
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# jacobian -> Berechnung der Jacobi-Matrix
# baseframe -> Berechnung der Jacobi-Matrix bezogen auf das Basis-Koordinatensystem des Roboters (nicht: Welt-KS)

# Tests
#
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-10
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
#
# Sources
# [SchapplerTapOrt2019] Schappler, M. et al.: Resolution of Functional Redundancy for 3T2R Robot Tasks using Two Sets of Reciprocal Euler Angles, Proc. of the 15th IFToMM World Congress, 2019
# [Ortmaier2014] Vorlesungsskript Robotik I
# [Ortmaier2014a] Vorlesungsskript Robotik II
# Initialization
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_debug := false:
codegen_opt := 2:
codegen_jacobi := 3: # 1=nur Position, 2=Rotationsmatrix-Jacobi, 3=mit Zeitableitung
;
# Funktionen aus Robotik-Repo
read("../robotics_repo_path"): # Lädt Variable robotics_repo_path
;
read(sprintf("%s/transformation/maple/proc_r2eulxyz_inert", robotics_repo_path)):
# Funktionen aus diesem Repo
read sprintf("../helper/proc_MatlabExport"):
read sprintf("../transformation/proc_transl"):
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
# Roboter-Definitionen
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name, base_method_name):
Trf := Trf:
Trf_c := Trf_c:
# Eliminiere die Basis-Position und Orientierung. Darf keinen Einfluss auf die Jacobi-Matrix haben.
# Bei komplizierten Ausdrücken (insbesondere bei kinematischen Zwangsbedingungen) können die Symbole manchmal nicht eliminiert werden
for i from 1 to NQB do
  Trf_c := subs({X_base_t[i,1]=0},Trf_c):
end do:
# Link-Index, für den die Jacobi-Matrix aufgestellt wird. Hier wird angenommen, dass der Endeffektor das letzte Segment (=Link) ist. Die Jacobi-Matrix kann hier aber für beliebige Segmente aufgestellt werden. (0=Basis)
LIJAC:=NL-1:
printf("Generiere Jacobi-Matrix für %s (Körper %d)\n", robot_name, LIJAC):
# Jacobi-Matrix analytisch (Translatorisch)
# Ortmaier2014a Gl. (1.15), S.14: Geometrische Zwangsbedingungen in impliziter Form
# Gleichung enthält verallgemeinerte Koordinaten und Endeffektorposition und ergibt Null.
x_EE := Matrix(3,1,<r_xEE; r_yEE; r_zEE>):
# Transformationsmatrix von der Basis zum gegebenen Punkt auf dem Körper. Die Koordinaten px,py,pz für den Punkt werden später eingesetzt.
T_p := Trf_c(1 .. 4, 1..4, LIJAC+1) . transl(<px;py;pz>):
h_transl := Matrix(T_p(1..3,4)) - x_EE:
h_transl_s :=convert_t_s(h_transl):
# Jacobi-Matrix der inversen Kinematik
b_transl := Matrix(3, NQJ):
for i from 1 to 3 do
  for j from 1 to NQJ do
    b_transl(i,j) := diff(h_transl_s(i,1), qJ_s(j,1)):
  end do:
end do:
# Export
save b_transl, sprintf("../codeexport/%s/tmp/jacobia_transl_%d_maple.m", robot_name, LIJAC):
# Ausdruck muss nochmal geladen werden, ansonsten hängt sich die Code-Optimerung mit "tryhard" auf.
# TODO: Klären warum das so ist und Problem beheben.
read sprintf("../codeexport/%s/tmp/jacobia_transl_%d_maple.m", robot_name, LIJAC):
b_transl := b_transl:
if codegen_act and codegen_debug then
  for i from 1 to 3 do
    for j from 1 to NQJ do
      MatlabExport(b_transl(i,j), sprintf("../codeexport/%s/tmp/jacobia_transl_%d_floatb_%s_%d_%d_matlab.m", robot_name, LIJAC, base_method_name, i, j), codegen_opt):
    end do:
  end do:
end if:
if codegen_act then
  MatlabExport(b_transl, sprintf("../codeexport/%s/tmp/jacobia_transl_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
end if:
# Jacobi-Matrix analytisch (Rotatorisch)
# Ortmaier2014a Gl. (1.15), S.14: Geometrische Zwangsbedingungen in impliziter Form
# Rotationsdarstellung des Endeffektors in RPY-Winkeln
xo_EE := Matrix(3,1,<phi_xEE; phi_yEE; phi_zEE>):
# r2eulxyz_inert liefert inert-arctan (%). Auswertung vor Export dauert zu lange. In exportiertem Matlab-Code muss "%arctan" händisch gegen "atan2" getauscht werden.
h_rota_rpy := r2eulxyz_inert(T_p) - xo_EE:
h_rota_s :=convert_t_s(h_rota_rpy):
# Jacobi-Matrix der inversen Kinematik
# Falls die RPY-Darstellung für dieses Körper-KS singulär ist, existiert die analytische Jacobi-Matrix nicht. Maple gibt dann einen Fehler aus. Die Kennzeichnung mit NaN sorgt dafür, dass dies in Matlab erkennbar ist.
# Die Singularität kann bei kinematischen Strukturen am Anfang auftreten.
b_rota := Matrix(3, NQJ):
for i from 1 to 3 do
  for j from 1 to NQJ do
    if T_p(3,3) = 0 then
      b_rota(i,j) := NaN: # Singulär
    else
      b_rota(i,j) := diff(h_rota_s(i,1), qJ_s(j,1)):
    end if:
  end do:
end do:
# Export
save b_rota, sprintf("../codeexport/%s/tmp/jacobia_rot_%d_maple.m", robot_name, LIJAC):
# Ausdruck nochmal laden.
read sprintf("../codeexport/%s/tmp/jacobia_rot_%d_maple.m", robot_name, LIJAC):
b_rota := b_rota:
# Einzelne Ausdrücke als Code exportieren (falls Gesamtausdruck zu komplex ist)
if codegen_act and codegen_debug then
  for i from 1 to 3 do
    for j from 1 to NQJ do
      MatlabExport(b_rota(i,j), sprintf("../codeexport/%s/tmp/jacobia_rot_%d_floatb_%s_%d_%d_matlab.m", robot_name, LIJAC, base_method_name, i, j), codegen_opt):
    end do:
  end do:
end if:
if codegen_act then
  interface(warnlevel=0): # Unterdrücke die Warnung für arctan
  MatlabExport(b_rota, sprintf("../codeexport/%s/tmp/jacobia_rot_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
  interface(warnlevel=3):
end if:
# Jacobi-Matrix geometrisch (Rotatorisch)
# Zusammenhang zwischen Geschwindigkeit der verallgemeinerten Koordinaten und Winkelgeschwindigkeit des Endeffektors ausgedrückt im Basis-Koordinatensystems
read sprintf("../codeexport/%s/tmp/velocity_worldframe_floatbase_%s_par1_maple.m", robot_name, base_method_name):
omega_W_i := omega_W_i:
# Basis-Terme entfernen (falls durch Maple später nicht möglich, s.o.)
for i from 1 to NQB do
  omega_W_i := subs({X_base_t[i,1]=0},omega_W_i):
  omega_W_i := subs({V_base_t[i,1]=0},omega_W_i):
end do:
# Winkelgeschwindigkeit als allgemeiner Ausdruck (nicht über Gelenk-Konfiguration)
# (fällt später sowieso weg)
omega_EE := Matrix(3,1,<omega_xEE; omega_yEE; omega_zEE>):
# Differenz als kinematische Zwangsbedingung aufstellen
h_rotg := Matrix(omega_W_i(1..3,LIJAC+1)) - omega_EE:
h_rotg_s :=convert_t_s(h_rotg):
# Jacobi-Matrix der inversen Kinematik
# Leite die Winkelgeschwindigkeiten des Endeffektors nach den Gelenkwinkelgeschwindigkeiten ab.
# Damit wird der Zusammenhang nach Ortmaier2014 Gl. (4.17), S.51 umgeformt.
b_rotg := Matrix(3, NQJ):
for i from 1 to 3 do
  for j from 1 to NQJ do
    b_rotg(i,j) := diff(h_rotg_s(i,1), qJD_s(j,1)):
  end do:
end do:
B_s := b_rotg:
# Export
save b_rotg, sprintf("../codeexport/%s/tmp/jacobig_rot_%d_maple.m", robot_name, LIJAC):
# Ausdruck nochmal laden.
read sprintf("../codeexport/%s/tmp/jacobig_rot_%d_maple.m", robot_name, LIJAC):
b_rotg := b_rotg:
if codegen_act and codegen_debug then
  for i from 1 to 3 do
    for j from 1 to NQJ do
      MatlabExport(b_rotg(i,j), sprintf("../codeexport/%s/tmp/jacobig_rot_%d_floatb_%s_%d_%d_matlab.m", robot_name, LIJAC, base_method_name, i, j), codegen_opt):
    end do:
  end do:
end if:
if codegen_act then
  MatlabExport(b_rotg, sprintf("../codeexport/%s/tmp/jacobig_rot_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
end if:
# Jacobi-Matrix der Rotationsmatrix (Rotatorisch)
# Ableitung der Rotationsmatrix nach den Gelenkwinkel (kann z.B. für die inverse Kinematik benutzt werden)
# Siehe [SchapplerTapOrt2019], Gl. 31, Term III und Gl. 27
R_0_i_s := convert_t_s( Trf_c(1 .. 3, 1..3, LIJAC+1) ):
Rb_0_i_s :=Reshape(R_0_i_s, 9,1):
J_Rb_s := Matrix(9, NQJ):
for i from 1 to 9 do
  for j from 1 to NQJ do
    J_Rb_s(i,j) := diff(Rb_0_i_s(i,1), qJ_s(j,1)):
  end do:
end do:
save J_Rb_s, sprintf("../codeexport/%s/tmp/jacobiR_rot_%d_floatb_%s_maple.m", robot_name, LIJAC, base_method_name):
if codegen_act and codegen_jacobi > 1 then
  MatlabExport(J_Rb_s, sprintf("../codeexport/%s/tmp/jacobiR_rot_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
end if:
# Jacobi-Zeitableitungen
# Die Jacobi-Zeitableitungen werden z.B. zur Umrechnung zwischen EE-Beschleunigungen und Gelenk-Beschleunigungen benötigt.
JD_rotg_t := diff~(convert_s_t(b_rotg), t):
JD_rotg_s := convert_t_s(JD_rotg_t):
JD_rota_t := diff~(convert_s_t(b_rota), t):
JD_rota_s := convert_t_s(JD_rota_t):
JD_transla_t := diff~(convert_s_t(b_transl), t):
JD_transla_s := convert_t_s(JD_transla_t):
# Zeitableitung der Ableitung der Rotationsmatrix nach den Gelenkwinkeln. Dieser Term wird zur Berechnung der Zeitableitung der kinematischen Zwangsbedingungen für PKM benötigt.
# Das wird bei der Berechnung der Coriolis-Terme aus dieser Darstellung benötigt.
JD_Rb_t := diff~(convert_s_t(J_Rb_s), t):
JD_Rb_s := convert_t_s(JD_Rb_t):
# Matlab-Export
if codegen_act and codegen_jacobi > 2 then
  MatlabExport(JD_rotg_s, sprintf("../codeexport/%s/tmp/jacobigD_rot_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
  interface(warnlevel=0): # Unterdrücke die Warnung für arctan
  MatlabExport(JD_rota_s, sprintf("../codeexport/%s/tmp/jacobiaD_rot_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
  interface(warnlevel=3):
  MatlabExport(JD_transla_s, sprintf("../codeexport/%s/tmp/jacobiaD_transl_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
  MatlabExport(JD_Rb_s, sprintf("../codeexport/%s/tmp/jacobiRD_rot_%d_floatb_%s_matlab.m", robot_name, LIJAC, base_method_name), codegen_opt):
end if:

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


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

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

(* type: work_gga_x *)

a1 := 1.647127:
a2 := 0.980118:
a3 := 0.017399:
b1 := 1.523671:
b2 := 0.367229:
b3 := 0.011282:

f0 := s -> (1 + a1*s^2 + a2*s^4 + a3*s^6)/(1 + b1*s^2 + b2*s^4 + b3*s^6):
f  := x -> f0(X2S*x):

ext: mpl

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


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

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

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

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

params_a_n := 19:

params_a_a := [
    7/6,  8/6,  9/6, 10/6, 17/12, 9/6, 10/6, 11/6,
   10/6, 11/6, 12/6, 10/6, 11/6, 12/6,  7/6,  8/6,
    9/6, 10/6, 13/12.0
]:

params_a_b := [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0]:
params_a_c := [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0]:
params_a_d := [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0]:

$include "th.mpl"

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


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

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

(* type: mgga_exc *)

(* Equations are from the Appendix of Kurth1999_889 *)

$include "gga_c_gapc.mpl"

gap_par0[10] = 0.06483*((9*Pi)/4)^(2/3): (* this is approximately equal to 0.23878 *)

(* Eq. (A4) *)
kcis_G := (rs, xt) -> xt^2*n_total(rs)^(2/3)/8:

(* Eq. (A9) *)
kcis_t := (rs, xt) -> 2^(2/3)*xt/(8*sqrt(rs)):

kcis_beta := 0.066725:

(* Eq. (A7) *)
kcis_gga0 := (rs, xt) -> f_pw(rs, 0)/(1 + kcis_beta*log(1 + kcis_t(rs, xt)^2/m_abs(f_pw(rs, 0)))):

(* Eq. (A8) *)
kcis_gga1 := (rs, xt) -> f_pw(rs, 1)/(1 + kcis_beta*log(1 + 2^(-1/3)*kcis_t(rs, xt)^2/m_abs(f_pw(rs, 1)))):

(* Eq. (A5) *)
(* The polarized parameters are the same as the unpolarized ones *)
(* except that c_1, c_2, c_3 are multiplied by 0.7. 1.5, and 2.59 respectively *)
kcis_eps_0 := (rs, xt) ->
  + (kcis_gga0(rs, xt) + gap_c1(rs, 0, gap_par0)*kcis_G(rs, xt))
  / (1 + gap_c2(rs, 0, gap_par0)*kcis_G(rs, xt) + gap_c3(rs, 0, gap_par0)*kcis_G(rs, xt)^2):

(* Eq. (A6) *)
kcis_eps_1 := (rs, xt) ->
  + (kcis_gga1(rs, xt) + 0.7*gap_c1(rs, 0, gap_par0)*kcis_G(rs, xt))
  / (1 + 1.5*gap_c2(rs, 0, gap_par0)*kcis_G(rs, xt) + 2.59*gap_c3(rs, 0, gap_par0)*kcis_G(rs, xt)^2):

(* Eq. (A2) *)
gap_f := (rs, z, xt) ->
  + kcis_eps_0(rs, xt)
  + f_zeta(z)*(kcis_eps_1(rs, xt) - kcis_eps_0(rs, xt)):

(* Eq. (A1) *)
kcis_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  + gap_f(rs, z, xt)
  - xs0^2/(8*ts0) * opz_pow_n( z,1)/2 * gap_f(rs*(2/(1 + z))^(1/3),  1, xs0)
  - xs1^2/(8*ts1) * opz_pow_n(-z,1)/2 * gap_f(rs*(2/(1 - z))^(1/3), -1, xs1):

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


text:
# Proof by PDF that X <~ rayleigh(1) == X^2 <~ gamma(1,2)

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

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

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

evalb(PDF_gamma(1, 2, x) = PDF_transform(x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
! Tsp36.mpl

TITLE
    Tsp36;

DATA
    NumCity := 36;

INDEX
    city := 1..NumCity;
    city2 := city;

DATA
    Cost[city, city2] := DATAFILE("Tsp36.dat");

BINARY VARIABLE
    Route[city, city2] WHERE (city <> city2) -> x;

INTEGER VARIABLE
    SeqNr[city] -> s;

MODEL

    MIN TotalCost = SUM(city, city2: Cost * Route);

SUBJECT TO

    Horz[city]: SUM(city2: Route) = 1;

    Vert[city2]: SUM(city: Route) = 1;

    SetSeq[city, city2 > 1]  WHERE (city <> city2):

        SeqNr[city] - SeqNr[city:=city2] + NumCity * Route <= NumCity - 1;

BOUNDS
    Tour[city]:  1 <= SeqNr <= NumCity;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


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

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

$include "gga_x_b88.mpl"

(* equation 11 *)
y := (x,u) -> x^2 - u:
(* equation 5 *)
gBecke := x -> b88_f(x)-1:
(* equation 24 *)
jk_f := (x,u,t) -> 1 + gBecke(x)/(1 + 2*y(x,u)/x^2):

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


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 "op.mpl"

(* Xalpha *)
op_qab         := 2.5654:
op_enhancement := xs -> 1:


text:
var $g i32
var $h i32
func $multiwayfunc ( var %n i32) i32 {
 multiway (dread i32 %n) @labdft {
   (add i32 (dread i32 $g, dread i32 $h)): goto @lab0
   (dread i32 $h): goto @lab1
   (constval i32 888): goto @lab0
   (constval i32 'Y'): goto @lab0
   (neg i32 (dread i32 $h)) : goto @lab9 }
@lab0
 return (constval i32 -3)
@labdft
 return (constval i32 100)
@lab9
 return (constval i32 9)
@lab1
 return (constval i32 1) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


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

`eta/full_trees` := (A::set) -> `eta/trees`(A);

`gamma/full_trees` := (A::set,B::set) -> (p) -> (TT,UU) ->
  `gamma/trees`(A,B)(p)(TT,UU);


text:
__BEGIN__

def xjoin(tgt, between, after, join)
	tgt.map {|a, b| "#{a}#{between}#{b}#{after}" }.join(join)
end

def sep(i, c, between = ",")
	between if c.length != i+1
end

__END__
%def genproto
%self.each do |msg|
struct [%msg.name%] {
	static const int method = [%msg.id%];  %>if msg.id

	%# constructor
	%unless msg.body.empty?
	[%msg.name%](
		const [%b.type%]& [%b.name%]_[%sep(i, msg.body)%]  %|b,i| msg.body.each_with_index
	) :
		[%b.name%]([%b.name%]_)[%sep(i, msg.body)%]  %|b,i| msg.body.each_with_index
	{ }
	%end

	%# default constructor
	[%msg.name%]() { }

	%# unpack
	void msgpack_unpack(msgpack::object o)
	{
		if(o.type != msgpack::type::ARRAY ||
				o.via.array.size < [%msg.required.size%] ||
				o.via.array.size > [%msg.body.size%]) {
			throw msgpack::type_error();
		}
		[%b.name%] = o.via.array.ptr[[%i%]].as<[%b.type%] >();  %|b,i| msg.required.each_with_index

		%unless msg.optional.empty?
		switch(o.via.array.size) {
		%msg.optional.length.times do |i|
		case [%msg.required.size + i%]:
			[%x.name%] = o.via.array.ptr[[%i%]].as<[%x.type%] >(); %|x| msg.optional[0..i].each
			[%x.name%] = [%x.default%];  %|x| msg.optional[i+1..-1].each
			break;
		%end
		}
		%end
	}

	%# pack
	template <typename Packer>
	void msgpack_pack(Packer& pk) const
	{
		pk.pack_array([%msg.body.size%]);
		pk.pack([%b.name%]); %|b| msg.body.each
	}

	%# members
	[%b.type%] [%b.name%];  %|b| msg.body.each
};

%end
%end  # def genproto
%# vim: syntax=mplex

max_stars_repo_name: frsyuki/ccf

max_issues_repo_name: frsyuki/ccf

max_forks_repo_name: frsyuki/ccf


text:
(*
 Copyright (C) 2019 Susi Lehtola

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

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

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

s12g_f := x -> params_a_bx*(params_a_A + params_a_B*(1 - 1/(1 + params_a_C*x^2 + params_a_D*x^4))*(1 - 1/(1 + params_a_E*x^2))):

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


text:
# Use this to define the canonical_labeling function in maple:

canonical_labeling := define_external('canonical_labeling', \
colored_edges::(ARRAY(1..3*m, integer[4])), \ # colored edges as array
\ # Format: ( va_1, vb_1, w_1, va_2, vb_2, w_2, ... )
\ # For a graph with edges (va_k, vb_k) having colors w_k.
ext_vertices::(ARRAY(1..l,integer[4])), \ # external vertices
\ # Array of vertices with another color ( ve_1, ve_2, ... )
lab::(ARRAY(1..n,integer[4])), \ # canonical labeling returned
\ # Canonical labelling is written into this.
n::(integer[4]), \ # number of vtcs
m::(integer[4]), \ # number of edges
w::(integer[4]), \ # number of colors
l::(integer[4]), \ # number of ext vertices
LIB="nauty_wrapper.so" ):

# Example:

lab:=Array([0,1,2,3,4],datatype=integer[4]):
ext:=Array([4],datatype=integer[4]):

G:=Array([0,3,0, 0,2,0, 2,1,0, 3,1,0, 3,4,0], datatype=integer[4]):

canonical_labeling(G, ext, lab, 5, 5, 1, 1);

lab;

G:=Array(subs([1=2,2=3,3=4,4=1],[0,3,0, 0,2,0, 2,1,0, 3,1,0, 3,4,0]), datatype=integer[4]):

canonical_labeling(G, ext, lab, 5, 5, 1, 1);

lab;

max_stars_repo_name: michibo/maple-graph-canonicalization

max_issues_repo_name: michibo/maple-graph-canonicalization

max_forks_repo_name: michibo/maple-graph-canonicalization


text:
######################################################################
# Binary operations on finite sets

`is_element/binops` := (A::set,B::set,C::set) -> proc(p)
 local AB,dom,cod,a,b;
 global reason;

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

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

 if {indices(p)} <> AB then
  reason := [convert(procname,string),"wrong domain",{indices(p)},AB];
  return false;
 fi;

 cod := map(op,{entries(p)});
 if cod minus C <> {} then
  reason := [convert(procname,string),"wrong codomain",cod,C];
  return false;
 fi;

 return true;
end;

`is_equal/binops` := (A::set,B::set,C::set) -> proc(p,q)
 local a,b;
 global reason;

 for a in A do
  for b in B do
   if p[a,b] <> q[a,b] then
    reason := [convert(procname,string),"different value",a,b,p[a,b],q[a,b]];
    return false;
   fi;
  od;
 od;

 return true;
end;

`is_leq/binops` := NULL;

`random_element/binops` := (A::set,B::set,C::set) -> proc()
 local p,n,r,a,b;

 p := table();
 if B = {} then
  if A = {} then
   return eval(p);
  else
   return FAIL;
  fi;
 fi;

 n := nops(C);

 r := rand(1..n);

 for a in A do
  for b in B do
   p[a,b] := C[r()];
  od;
 od;

 return eval(p);
end;

`list_elements/binops` := proc(A::set,B::set,C::set)
 local nA,nB,nC,J,L,i,j,k,u;

 nA := nops(A);
 nB := nops(B);
 nC := nops(C);
 J := [seq(seq([A[i],B[j]],j=1..nB),i=1..nA)];
 L := [[]];
 for i from 1 to nA*nB do
  L := [seq(seq([op(u),j],j=1..nC),u in L)];
 od:

 L := map(u -> table([seq(op(J[k])=C[u[k]],k=1..nA*nB)]),L);
 return L;
end;

`count_elements/binops` := (A::set,B::set,C::set) -> nops(C) ^ (nops(A)*nops(B));

`is_associative/binops` := (A::set) -> proc(m)
 local a,b,c;
 for a in A do
  for b in A do
   for c in A do
    if m[a,m[b,c]] <> m[m[a,b],c] then
     return false;
    fi;
   od;
  od;
 od;
 return true;
end:

`is_commutative/binops` := (A::set) -> proc(m)
 local a,b;
 for a in A do
  for b in A do
   if m[a,b] <> m[b,a] then
    return false;
   fi;
  od;
 od;
 return true;
end:

`is_left_proj/binops` := (A::set) -> proc(m)
 local a,b;
 for a in A do
  for b in A do
   if m[a,b] <> a then
    return false;
   fi;
  od;
 od;
 return true;
end:

`is_right_proj/binops` := (A::set) -> proc(m)
 local a,b;
 for a in A do
  for b in A do
   if m[a,b] <> b then
    return false;
   fi;
  od;
 od;
 return true;
end:

`is_idempotent/binops` := (A::set) -> proc(m)
 local a;
 for a in A do
  if m[a,a] <> a then
   return false;
  fi;
 od;
 return true;
end:

`is_left_neutral/binops` := (A::set) -> proc(m,e)
 local a;

 if not(member(e,A)) then return false; fi;
 for a in A do
  if m[e,a] <> a then return false; fi;
 od:
 return true;
end:

`is_right_neutral/binops` := (A::set) -> proc(m,e)
 local a;

 if not(member(e,A)) then return false; fi;
 for a in A do
  if m[a,e] <> a then return false; fi;
 od:
 return true;
end:

`is_neutral/binops` := (A::set) -> proc(m,e)
 return `is_left_neutral/binops`(A)(m,e) and
        `is_right_neutral/binops`(A)(m,e);
end:

`has_left_neutral/binops` := (A::set) -> proc(m)
 local e;

 for e in A do
  if `is_left_neutral/binops`(A)(m,e) then
   return true;
  fi;
 od:
 return false;
end:

`has_right_neutral/binops` := (A::set) -> proc(m)
 local e;

 for e in A do
  if `is_right_neutral/binops`(A)(m,e) then
   return true;
  fi;
 od:
 return false;
end:

`has_neutral/binops` := (A::set) -> proc(m)
 local e;

 for e in A do
  if `is_neutral/binops`(A)(m,e) then
   return true;
  fi;
 od:
 return false;
end:

`is_left_absorbing/binops` := (A::set) -> proc(m,e)
 local a;

 if not(member(e,A)) then return false; fi;
 for a in A do
  if m[e,a] <> e then return false; fi;
 od:
 return true;
end:

`is_right_absorbing/binops` := (A::set) -> proc(m,e)
 local a;

 if not(member(e,A)) then return false; fi;
 for a in A do
  if m[a,e] <> e then return false; fi;
 od:
 return true;
end:

`is_absorbing/binops` := (A::set) -> proc(m,e)
 return `is_left_absorbing/binops`(A)(m,e) and
        `is_right_absorbing/binops`(A)(m,e);
end:

`has_left_absorbing/binops` := (A::set) -> proc(m)
 local e;

 for e in A do
  if `is_left_absorbing/binops`(A)(m,e) then
   return true;
  fi;
 od:
 return false;
end:

`has_right_absorbing/binops` := (A::set) -> proc(m)
 local e;

 for e in A do
  if `is_right_absorbing/binops`(A)(m,e) then
   return true;
  fi;
 od:
 return false;
end:

`has_absorbing/binops` := (A::set) -> proc(m)
 local e;

 for e in A do
  if `is_absorbing/binops`(A)(m,e) then
   return true;
  fi;
 od:
 return false;
end:


text:
program nwhilest;
var i, j, k : integer;
begin
    i := 1;
    while i < 10 do begin
        j := 1;
        while j < 10 do begin
            k := 1;
            while k < 10 do begin
                if (k div 2) * 2 = k then
                        k := k + 1
                    end
                else
                        k := k + 1
                    end
            end;
            j := j + 1;
        end;
        i := i + 1
    end;
    writeln('All End')
end.

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:
### Domain abstraction

# A domain is abstractly a pair of a "square" domain (a product of square
# domains, or an interval), the domain "bounds", and a domain "shape". The
# domain "bounds" are essentially the variables, each with the type of domain
# (Int,Sum, etc) ; the domains "shape" is sums of products of constraints
# (relations) on the domain, interspersed with 'makes' (DInto) which parametrize
# a subdomain by one of the domain variables ('picking' a particular
# parametrization of the domain, in which one variable is fixed in subsequent
# "shape" constructors).

# the seperation between the two 'stages' of domain exists mainly because
# certain simplifications need to happen after extracting the domain
# bounds but before extracting the domain shape, as the shape ends up
# becoming part of the bounds of an inner simplification.

# Note that large parts of this are not very nice. It is a very literal
# translation of the code which it replaced, which a very thin amount
# of abstraction on top. But it does the hard work of factoring
# out all of the domain-related code into a single module.
# At some point, the interface should be improved, and the implementation
# should get rid of all sorts of unnecessary conversions.

# Broad TODOs:
#
# Mechanism for 'checking if' and 'making' constraints about variables should
#   be abstracted out into its own module (it is needed in multiple places, and
#   done sort of ad-hoc) (this could be fixed by a broader design - see "merging
#   KB with Domain")
#
# Shape extraction needs to be reworked; we "flatten" constraints multiple
#   times; This should be done once after shape extraction, not at every step
#
# DInto should also optionally omit the bounds if they are identical to the
#   bounds in the a-priori domain bounds (i.e. just `DInto(x)`); DInto sort of
#   means 'we've solved this bound in this subcontext' but the `DInto`s where
#   the bound is identically the a-priori bound are erased (as they should be)
#   - we want a slightly seperate form to mean "we are done improving this", but
#   we still want to be able to omit bounds which would be identical. It should
#   also allow specifying only a subset of the bounds.
#
# A more granular interface for composing/re-composing simplifiers, including
#   call simplifier A {before/after} B
#
# Domain bounds should store a representation of 'dependancy' of variables
#   (a more refined representation; the current one is a linear order of variables)
#   Basically we want a list of graphs, nodes are variables and edges are
#   variable dependencies. If each graph is acyclic, the variable dependency
#   in the bounds is valid. e.g.
#     [x=0..1, y=0..1, z=0..1]  =>  [ {[x, {}]}, {[y, {}]}, {[z, {}]} ]
#     [x=0..1, y=x..1, z=0..y]  =>  [ {[z, {y}], [y, {x}], [x, {}]} ]
#   This generalizes the previous notion of variable order, although in certain
#     cases it is ambiguous what the 'original' order was. This is only when
#     there are variable dependencies, so perhaps this isn't a big deal - we may
#     want to rearrange the dependant bounds even when there is no 'interesting'
#     constraint to solve in the shape, because we know a particular order will
#     be better in some way.
#   The graph is probably not a convenient place to store the variable bounds.
#    (this is one of the big reasons we want to view the variables as a table -
#    we have lots of different information for which the convenient formats are
#    all different)
#
# A better interface for `simpl_relation` - something that allows it to treat
#   arbitrary constructor forms as And, Or, Not. Typically, we do a `subs` right
#   before the call to `simpl_relation` to put things in the right form, then
#   another `subs` immediately inside `simpl_relation` to turn `And,Or` into `&and,&or`.
#
# Long term, KB and Domain should be merged, as their functionality overlaps
#   quite a bit; basically
#     - Introduce becomes a variable in DBound
#     - Constrain becomes DConstrain
#     - DBound becomes DInto
#     - Let becomes (probably?) a one piece partition,
#           or maybe a new variable type
#          (it doesn't seem to appear in present code?)
#  Logically, this is at most an ephemeral change - Domain already generalizes
#   KB. But mechanically it is no small undertaking.
#  assert_deny becomes a "DConstrain" smart constructor (note that in the
#   current AST for Domain, DConstrain appears only at leaves, i.e. products
#   are pushed all the way down; whereas KB is 'flat') but the shape of domains
#   means it is probably much simpler to 'normalize' the constraints all at
#   once, as opposed to when every constrain is inserted, as KB currently does.
#  kb_subtract becomes a similar `domain_combine` function (kb_subtract is really
#   an 'combine' as well as a 'deconstruct into a way to rebuild the KB' together in
#   one - we don't really need the 'deconstruct' part anymore, but 'add' is
#   still useful) - basically this function takes two domains, checks if the
#   variables of the first are a subset of the second (non commutative!) and if
#   so, combines the shapes according to the passed 'shape constructor'; all of
#   the constructors (aside from DInto)
#  gen{whatever} becomes a part of `Domain:-Extract` - basically we will not
#   have `gen` as part of the "interface" to KB+Domain, but rather an entry in
#   the `Ext` function tables.
#  most of the other functions in KB can be basically literally translated to
#    the subset of Domain which corresponds to KB

Domain := module()
    uses Hakaru, Partition, SolveTools[Inequality] ;
    global DSplit; global DInto; global DNoSol;

    local ModuleLoad := proc($)
      local ty_nm, g;
      for ty_nm in [ indices(DomainTypes, nolist) ] do
        TypeTools[AddType]( ty_nm, DomainTypes[ty_nm] );
      end do;

      #op([2,6], ...) of a module is its globals.
      for g in op([2,6], thismodule) do
        if g <> eval(g) then
          unassign(g);
          WARNING("Previous value of global name '%1' erased.", g)
        end if;
        protect(g);
      end do;
    end proc;

    local ModuleUnload := proc($)
      local ty_nm;
      for ty_nm in [ indices(DomainTypes, nolist) ] do
        if TypeTools[Exists](ty_nm) then TypeTools[RemoveType](ty_nm) end if;
      end do;
    end proc;

    # Extending domain extraction and replacement.
    export ExtBound := table();
    export ExtShape := table();

    export Set_ExtBound := proc(nm,val,$)
      unprotect(Domain:-ExtBound);
      Domain:-ExtBound[nm] := val;
      protect(Domain:-ExtBound);
    end proc;

    export Set_ExtShape := proc(nm,val,$)
      unprotect(Domain:-ExtShape);
      Domain:-ExtShape[nm] := val;
      protect(Domain:-ExtShape);
    end proc;

$include "Domain/Has.mpl"
$include "Domain/Bound.mpl"
$include "Domain/Shape.mpl"
$include "Domain/Types.mpl"
$include "Domain/Extract.mpl"
$include "Domain/Apply.mpl"
$include "Domain/Improve.mpl"

    export simpl_relation :=
    proc( expr_ :: {relation, boolean, specfunc({`And`,`Not`,`Or`}), `and`, `not`, `or`}
        , { norty := 'DNF' }
        , $) :: { specfunc(specfunc({relation, specfunc(relation, Not)}, `Or`), `And`)
                , specfunc(specfunc({relation, specfunc(relation, Not)}, `And`), `Or`)
                };
        local expr := expr_, outty, outmk, inty, inmk, ty_ord ;

        expr := foldr( proc(v,e) subsindets(e, op(v)) end proc
                     , expr
                     , [ { specfunc(relation, `Not`), `not`(relation) }
                       , x-> KB:-negate_rel(op(1,x)) ]
                     , [ { specfunc(`Not`), `not` }
                       , x->Logic:-`&not`(op(1,x)) ]
                     , [ { specfunc(`Or`), `or` }
                       , x->Logic:-`&or`(op(x)) ]
                     , [ { specfunc(`And`), `and` }
                       , x->Logic:-`&and`(op(x)) ] );
        expr := Logic:-Normalize(expr, form=norty);
        expr := foldr( proc(v,e) subsindets(e, op(v)) end proc
                     , expr
                     , [ specfunc(Logic:-`&and`), x->`And`(op(x)) ]
                     , [ specfunc(Logic:-`&or`) , x->`Or`(op(x)) ]
                     , [ specfunc(Logic:-`&not`), x->KB:-negate_rel(op(1,x))  ] );

        if expr :: identical(false) then
            return `if`(norty='DNF', `Or`(), `And`(`Or`()));
        elif expr :: identical(true) then
            return `if`(norty='DNF', `Or`(`And`()), `And`());
        end if;

        ty_ord := `if`(norty='DNF', [1,2], [2,1]);
        outty, inty := [ 'specfunc(Or)', 'specfunc(And)' ][ty_ord][];
        outmk, inmk := [ `Or`, `And` ][ty_ord][];

        if not expr :: outty then expr := outmk(expr) end if;
        map(x -> if not x :: inty then inmk(x) else x end if, expr);
    end proc;

    # The main interface to Domain
    export Reduce := proc(e0,kb::t_kb,f_into,f_body,f_nosimp:=(_->FAIL),$)
      local e := e0, dom_specb, dom_specw, dom_ctx, dom_spec;
      # Build the domain
      dom_specb, e := op(Domain:-Extract:-Bound(e));
      if Domain:-Bound:-isEmpty(dom_specb) then return f_nosimp(e0) end if;
      dom_specw, e := op(Domain:-Extract:-Shape(e));
      dom_ctx := {op(KB:-kb_to_assumptions(kb))};
      dom_specb := DBound(op(1,dom_specb), dom_ctx);
      dom_spec := DOMAIN(dom_specb, dom_specw);

      # Improve, if necessary, then apply back to the expression
      if dom_specw <> DConstrain() then dom_spec := Domain:-Improve(dom_spec) end if;
      Domain:-Apply(dom_spec, kb, f_into, f_body)(e);
    end proc;

end module;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
# convert from Gonnet's FPAccuracy data sets to mpfr format
# http://www.inf.ethz.ch/personal/gonnet/FPAccuracy/all.tar.Z

# 1 - cut the lines from (say) C/acos.c, remove the 3rd (eps) field,
#     replace the commata ',' by spaces, and remove the final '};'
#     (hint: cut -d" " -f1,2,4,5 /tmp/acos.c > /tmp/acos2.c)
# 2 - edit the infile and outfile lines below, and run
#     maple -q < gonnet.mpl

infile := "/tmp/acos2.c":
outfile := "acos":

###################### don't edit below this line #############################

foo := proc(arg_m, val_m, arg_e, val_e, fp)
   fprintf (fp, "53 53 n ", 53);
   to_hex(arg_m, arg_e, fp);
   fprintf (fp, " ");
   # warning: Gonnet stores -val_e
   to_hex(val_m, -val_e, fp);
   fprintf (fp, "\n");
end:

to_hex := proc(m, e, fp)
   if m<0 then fprintf (fp, "-") fi;
   fprintf (fp, "0x%sp%d", convert(abs(m),hex), e);
end:

copyright := proc(fp)
   fprintf (fp, "# This file was generated from the FPAccuracy package\n# http://www.inf.ethz.ch/personal/gonnet/FPAccuracy/all.tar.Z:\n# Copyright (C) Gaston H. Gonnet\n# This program is free software; you can redistribute it and/or\n# modify it under the terms of the GNU General Public License\n# as published by the Free Software Foundation; either version 2\n# of the License, or (at your option) any later version.\n# This program is distributed in the hope that it will be useful,\n# but WITHOUT ANY WARRANTY; without even the implied warranty of\n# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n# GNU General Public License for more details.\n# You should have received a copy of the GNU General Public License\n# along with this program; if not, write to the Free Software\n# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.\n")
end:

fp := fopen (outfile, WRITE):

l := readdata(infile, integer, 4):
copyright(fp):
for e in l do foo(op(e), fp) od:

fclose (fp);

quit;


text:
# Maple 2018.2 file to be called with
# read "log_trans_prob.mpl"

# we assume that b > 0!

with(LinearAlgebra):

trans_prob_unequal := proc(a, b, t, lambda, mu)
  local h, theta, psi, phi, alpha, beta, gamma:

  theta := lambda - mu:
  psi := lambda + mu:

  phi := (exp(theta * t) - 1) / (lambda * exp(theta * t) - mu):
  alpha := mu * phi:
  beta := lambda * phi:
  gamma := 1 - psi * phi:

  binomial(a + b - 1, a - 1) * alpha^a * beta^b + sum(binomial(a, h) * binomial(a + b - h - 1, a - 1) * alpha^(a - h) * beta^(b - h) * gamma^h, h = 1 .. min(a, b))
end proc:

trans_prob_equal := proc(a, b, t, lambda)
  local h, theta, alpha, gamma:

  theta := lambda * t:

  alpha := theta / (1 + theta):
  gamma := (1 - theta) / (1 + theta):

  binomial(a + b - 1, a - 1) * alpha^(a + b) + sum(binomial(a, h) * binomial(a + b - h - 1, a - 1) * alpha^(a + b - 2 * h) * gamma^h, h = 1 .. min(a, b))
end proc:

log_trans_prob := proc(a, b, t, lambda, mu)
  local s, lp:

  s := 0:
  lp := 0:

  if (lambda <> mu) then
    s := trans_prob_unequal(a, b, t, lambda, mu):
  else
    s := trans_prob_equal(a, b, t, lambda):
  end if:

  if evalf[20](s) > 0 then
    lp := ln(s):
  else
    lp := ln(normal(s)):
  end if:

  lp
end proc:

log_trans_prob_float := proc(a, b, t, lambda, mu)
  evalf[100](log_trans_prob(a, b, t, lambda, mu))
end proc:

# generate data for Figures 1 and 3
N := 1000:
A := 25:
B := 35:
T := 2:
L := 1:
M := 3:

values := Array([seq(0, i = 1 .. N)]):

for n to N do
  values[n] := log_trans_prob_float(A, B, T, L, n * M / N):
end do:

values := Transpose(Matrix(convert(values, list))):

Export("log_trans_prob_values_univariate.csv", values):

# generate data for Figure 4 (beware: it takes a long time)
N := 100:
A := 200:
B := 100:
T := 1:
L := 10:
M := 10:

values := Array([seq(0, i = 1 .. (N * N))]):

c0 := 1:
for c1 to N do
  for c2 to N do
    values[c0] := log_trans_prob_float(A, B, T, c1 * L / N, c2 * M / N):
    c0 := c0 + 1:
  end do:
end do:

values := Transpose(Matrix(convert(values, list))):

Export("log_trans_prob_values_bivariate.csv", values):


text:
# Base Parameter Energy Regressor for Robot based on MDH frames
# Einleitung
# Erstellung einer parameterlinearen und -minimalen Regressorform
#
# Dateiname:
# robot -> Berechnung für allgemeinen Roboter
# chain -> Berechnung für eine serielle Struktur (nicht: Baumstruktur)
# fixb -> fixed base. Kein Floating base Modell. Dort ist diese Form der Minimalparameterform nicht möglich.
# rotmat -> Kinematik wird mit Rotationsmatrizen berechnet
# energy -> Berechnung der Energie
# regressor -> Regressorform (parameterlinear)
# Authors
# Alexander Tödtheide, toedtheide@irt.uni-hannover.de
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-03
#
# (C) Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Quellen
# [HRL_IDR] Skript Humanoid Robotics Lab - Serial Chain Robot Identification
# [GautierKhalil1988] A Direct Determination of Minimum Inertial Parameters of Robots
# [GautierKhalil1990] Direct Calculation of Minimum Set of Inertial Parameters of Serial Robots
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
interface(rtablesize=100): # Zur Anzeige von größeren Vektoren
;
with(LinearAlgebra):
with(ArrayTools):
with(codegen):
with(CodeGeneration):
with(StringTools):
codegen_act := true:
codegen_opt := 2:
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../robot_codegen_definitions/robot_env":
printf("%s. Generiere Minimalparameterregressor der Energie für %s (geometrischer Ansatz)\n", \
  FormatTime("%Y-%m-%d %H:%M:%S"), robot_name, codegen_dynpar):
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name, base_method_name):
read sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
kin_constraints_exist := kin_constraints_exist: # nur zum Abschätzen der Komplexität
;
# Term-Vereinfachungen einstellen
if not assigned(simplify_options) or simplify_options(7)=-1 then # Standard-Einstellungen:
  if not kin_constraints_exist then # normale serielle Ketten und Baumstrukturen
    use_simplify := 0: # Standardmäßig aus
  else # mit kinematischen Zwangsbedingungen
    use_simplify := 1: # standardmäßig simplify-Befehle anwenden
  end if:
else # Benutzer-Einstellungen:
  use_simplify := simplify_options(7): # siebter Eintrag ist für Energie-Regressor
end if:

# Ergebnisse der Energie laden
read sprintf("../codeexport/%s/tmp/energy_potential_floatb_%s_worldframe_par2_maple.m", robot_name, base_method_name):
read sprintf("../codeexport/%s/tmp/energy_kinetic_floatb_%s_linkframe_par2_maple.m", robot_name, base_method_name):
T_floatb := T:
U_floatb := U_grav:
# Einfluss der Basis-Geschwindigkeit und -Position entfernen
T_fixb:=T_floatb:
U_fixb:=U_floatb:
for i from 1 to 6 do
  T_fixb := subs({V_base_t[i,1]=0},T_fixb):
  U_fixb := subs({X_base_t[i,1]=0},U_fixb):
end do:
# Die kinetischen und potentiellen Energien aus (2) und (3) stehen ab hier durch T_fixb und U_fixb zur Verfügung.
# Der Parametervektor 'PV2_vec' aus (13) wurde in 'robot_tree_floatb_twist_definitions.mw' aufgestellt.

# Parameterlinearisierung
# Parameterlinearisierung auf Basis von [HRL_IDR] (14) und (15)
# Linearisierung
t_ges := Matrix(1, 10*(NL-1)):
u_ges := Matrix(1, 10*(NL-1)):
for i to 10*(NL-1) do
  t_ges[1,i] := diff(T_fixb,PV2_vec[10+i,1]);
  u_ges[1,i] := diff(U_fixb,PV2_vec[10+i,1]);
end do:
# Terme vereinfachen
if use_simplify=1 then
  tmp_t1:=time():
  tmp_l11 := length(t_ges):
  tmp_l12 := length(u_ges):
  t_ges := simplify2(t_ges):
  u_ges := simplify2(u_ges):
  tmp_t2:=time():
  tmp_l21 := length(t_ges):
  tmp_l22 := length(u_ges):
  printf("%s: Terme für Energie-Regressor vereinfacht. Länge: %d->%d / %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l11, tmp_l21, tmp_l12, tmp_l22, tmp_t2-tmp_t1):
end if:
# Export
save t_ges, sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_regressor_maple.m", robot_name):
save u_ges, sprintf("../codeexport/%s/tmp/energy_potential_fixb_regressor_maple.m", robot_name):
if codegen_act then
  MatlabExport(convert_t_s(t_ges), sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_regressor_matlab.m", robot_name), codegen_opt):
end if:
if codegen_act then
  MatlabExport(convert_t_s(u_ges), sprintf("../codeexport/%s/tmp/energy_potential_fixb_regressor_matlab.m", robot_name), codegen_opt):
end if:
# Parameterminimierung
# Prüfe, ob der folgende Algorithmus funktionieren kann
# Bestimme, ob es eine Baumstruktur ist. Wenn ja, funktioniert der Algorithmus (noch) nicht (nicht implementiert).
tree:=false:
for i from 1 to NJ do
  if not v(i) = i-1 then
    tree := true: break:
  end if:
end do:

if assigned(user_CoM) or assigned(user_M) or assigned(user_inertia) \
  or kin_constraints_exist or tree then
  # es gibt einen Sonderfall, diese Berechnung der Minimalparameter funktioniert voraussichtlich nicht.
  printf("%s. Geometrische Berechnung der Minimalparameter nicht möglich.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
  quit: # Funktioniert in GUI nicht richtig...
  robot_name := "": # ...Daher auch Löschung des Roboternamens.
end if:
# Minimalparametervekor
# Definiere Parametermatrix
# Nehme nur die Inertialparameter der bewegten Segmente, nicht die Basis.
XX := Matrix(NJ, 1, PV2_mat[2 .. NL, 1]):
XY := Matrix(NJ, 1, PV2_mat[2 .. NL, 2]):
XZ := Matrix(NJ, 1, PV2_mat[2 .. NL, 3]):
YY := Matrix(NJ, 1, PV2_mat[2 .. NL, 4]):
YZ := Matrix(NJ, 1, PV2_mat[2 .. NL, 5]):
ZZ := Matrix(NJ, 1, PV2_mat[2 .. NL, 6]):
mX := Matrix(NJ, 1, PV2_mat[2 .. NL, 7]):
mY := Matrix(NJ, 1, PV2_mat[2 .. NL, 8]):
mZ := Matrix(NJ, 1, PV2_mat[2 .. NL, 9]):
m :=  Matrix(NJ, 1, PV2_mat[2 .. NL, 10]):

# Rekursive Berechnung der Minimalparameter
# Rekursive Berechnung der Minimalparameter mit [HRL_IDR] (23), siehe Abbildung 1.
# [GautierKhalil1990] (eq. 15,16)
for i from NJ by -1 to 2 do
  # printf("i=%d, sigma=%d\n", i, sigma[i,1]):
  if sigma[i,1] = 0 then # Drehgelenk (eq. 15)
    XX[i,  1] := XX[i,  1]-YY[i,1]:
    XX[i-1,1] := XX[i-1,1] + YY[i,1] + d[i,1]^2*m[i,1]+2*d[i,1]*mZ[i,1]:
    XY[i-1,1] := XY[i-1,1] + a[i,1]*sin(alpha[i,1])*mZ[i,1]+a[i,1]*d[i,1]*sin(alpha[i,1])*m[i,1]:
    XZ[i-1,1] := XZ[i-1,1] - a[i,1]*cos(alpha[i,1])*mZ[i,1]-a[i,1]*d[i,1]*cos(alpha[i,1])*m[i,1]:
    YY[i-1,1] := YY[i-1,1] + cos(alpha[i,1])^2*YY[i,1]+2*d[i,1]*cos(alpha[i,1])^2*mZ[i,1]+(a[i,1]^2+d[i,1]^2*cos(alpha[i,1])^2)*m[i,1]:
    YZ[i-1,1] := YZ[i-1,1] + cos(alpha[i,1])*sin(alpha[i,1])*YY[i,1]+2*d[i,1]*cos(alpha[i,1])*sin(alpha[i,1])*mZ[i,1]+d[i,1]^2*cos(alpha[i,1])*sin(alpha[i,1])*m[i,1]:
    ZZ[i-1,1] := ZZ[i-1,1] + sin(alpha[i,1])^2*YY[i,1]+2*d[i,1]*sin(alpha[i,1])^2*mZ[i,1]+(a[i,1]^2+d[i,1]^2*sin(alpha[i,1])^2)*m[i,1]:
    mX[i-1,1] := mX[i-1,1] + a[i,1]*m[i,1]:
    mY[i-1,1] := mY[i-1,1] - sin(alpha[i,1])*mZ[i,1]-d[i,1]*sin(alpha[i,1])*m[i,1]:
    mZ[i-1,1] := mZ[i-1,1] + cos(alpha[i,1])*mZ[i,1]+d[i,1]*cos(alpha[i,1])*m[i,1]:
    m[i-1 ,1] := m[i-1 ,1] + m[i,1]:
  else: # Schubgelenk (eq. 16)
    # im Normalfall nicht geänderte Parameter initialisieren.
    # Wird durch Sonderregeln weiter unten überschrieben.
    mR_i := m[i-1 ,1]:
    mXR_i := mX[i-1,1]:
    mYR_i :=mY[i-1,1]:
    mZR_i:=mZ[i-1,1]:
    XXR_i := XX[i-1,1] + cos(theta[i,1])^2*XX[i,1]
                       - 2*cos(theta[i,1])*sin(theta[i,1])*XY[i,1]+sin(theta[i,1])^2*YY[i,1]:
    XYR_i := XY[i-1,1] + cos(theta[i,1])*sin(theta[i,1])*cos(alpha[i,1])*XX[i,1]
                       + (cos(theta[i,1])^2-sin(theta[i,1])^2)*cos(alpha[i,1])*XY[i,1]
                       - cos(theta[i,1])*sin(alpha[i,1])*XZ[i,1]
                       - cos(theta[i,1])*sin(theta[i,1])*cos(alpha[i,1])*YY[i,1]
                       + sin(theta[i,1])*sin(alpha[i,1])*YZ[i,1]:
    XZR_i := XZ[i-1,1] + cos(theta[i,1])*sin(theta[i,1])*sin(alpha[i,1])*XX[i,1]
                       + (cos(theta[i,1])^2-sin(theta[i,1])^2)*sin(alpha[i,1])*XY[i,1]
                       + cos(theta[i,1])*cos(alpha[i,1])*XZ[i,1]
                       - cos(theta[i,1])*sin(theta[i,1])*sin(alpha[i,1])*YY[i,1]
                       - sin(theta[i,1])*cos(alpha[i,1])*YZ[i,1]:
    YYR_i := YY[i-1,1] + sin(theta[i,1])^2*cos(alpha[i,1])^2*XX[i,1]
                       + 2*cos(theta[i,1])*sin(theta[i,1])*cos(alpha[i,1])^2*XY[i,1]
                       - 2*sin(theta[i,1])*cos(alpha[i,1])*sin(alpha[i,1])*XZ[i,1]
                       + cos(theta[i,1])^2*cos(alpha[i,1])^2*YY[i,1]
                       - 2*cos(theta[i,1])*cos(alpha[i,1])*sin(alpha[i,1])*YZ[i,1]
                       + sin(alpha[i,1])^2*ZZ[i,1]:
    YZR_i := YZ[i-1,1] + sin(theta[i,1])^2*cos(alpha[i,1])*sin(alpha[i,1])*XX[i,1]
                       + 2*cos(theta[i,1])*sin(theta[i,1])*cos(alpha[i,1])*sin(alpha[i,1])*XY[i,1]
                       + sin(theta[i,1])*(cos(alpha[i,1])^2-sin(alpha[i,1])^2)*XZ[i,1]
                       + cos(theta[i,1])^2*cos(alpha[i,1])*sin(alpha[i,1])*YY[i,1]
                       + cos(theta[i,1])*(cos(alpha[i,1])^2-sin(alpha[i,1])^2)*YZ[i,1]
                       - cos(alpha[i,1])*sin(alpha[i,1])*ZZ[i,1]:
    ZZR_i := ZZ[i-1,1] + sin(theta[i,1])^2*sin(alpha[i,1])^2*XX[i,1]
                       + 2*cos(theta[i,1])*sin(theta[i,1])*sin(alpha[i,1])^2*XY[i,1]
                       + 2*sin(theta[i,1])*cos(alpha[i,1])*sin(alpha[i,1])*XZ[i,1]
                       + cos(theta[i,1])^2*sin(alpha[i,1])^2*YY[i,1]
                       + 2*cos(theta[i,1])*cos(alpha[i,1])*sin(alpha[i,1])*YZ[i,1]
                       + cos(alpha[i,1])^2*ZZ[i,1]:
    # Im Kap. IV.C. sind Sonderfälle nur bezogen auf Schubgelenke definiert.
    # Diese scheinen aber noch überhaupt nicht gut zu funktionieren. Unklar, warum.
    # Sonderfall parallele Achsen (eq. 19)
    if sin(alpha[i,1]) = 0 and false then # deaktiviert. TODO: Klären, ob der Sonderfall sinnvoll ist
      ZZR_i := ZZ[i-1,1] + 2*a[i,1]*cos(theta[i,1])*mX[i,1] - 2*a[i,1]*sin(theta[i,1])*mY[i,1]:
      mXR_i := mX[i-1,1] + cos(theta[i,1])*mX[i,1] - sin(theta[i,1])*mY[i,1]:
      mYR_i := mY[i-1,1] + sin(theta[i,1])*cos(alpha[i,1])*mX[i,1] + cos(theta[i,1])*cos(alpha[i,1])*mY[i,1]:
    end if:
    # Sonderfall senkrechte Achsen (eq. 20)
    if cos(alpha[i,1]) = 0 and false then # deaktiviert.
      YYR_i := YY[i-1,1] + 2*a[i,1]*cos(theta[i,1])*mX[i,1] - 2*a[i,1]*sin(theta[i,1])*mY[i,1]:
      mXR_i := mX[i-1,1] + cos(theta[i,1])*mX[i,1] - sin(theta[i,1])*mY[i,1]:
      mZR_i := mZ[i-1,1] + sin(alpha[i,1])*sin(theta[i,1])*mX[i,1] + sin(alpha[i,1])*cos(theta[i,1])*mY[i,1]:
    end if:
    # Zuweisung der neuen Variable (wegen Fallunterscheidung parallel/senkrecht)
    # Nur für Schubgelenke notwendig
    m[i-1 ,1] := mR_i:
    mX[i-1,1] := mXR_i:
    mY[i-1,1] := mYR_i:
    mZ[i-1,1] := mZR_i:
    XX[i-1,1] := XXR_i:
    XY[i-1,1] := XYR_i:
    XZ[i-1,1] := XZR_i:
    YY[i-1,1] := YYR_i:
    YZ[i-1,1] := YZR_i:
    ZZ[i-1,1] := ZZR_i:
  end if: # Drehgelenk/Schubgelenk
end do:

# Parameter ohne Einfluss "Markieren"
# Markierung von Parametern ohne Einfluss an Hand von [HRL_IDR] (18)
printf("The following parameters don't have any effect:\n"):

for i to 10*NJ do
  effect := 0:
  if not (t_ges[1,i]=0) then
    # printf("t_ges[1,%d] not Null\n", i):
    effect := 1:
  else
    for j to NQJ do
      if has(u_ges[1,i],{qJ_t(j,1)}) then
        # printf("u_ges[1,%d] not constant\n", i):
        effect := 1:
        break:
      end if:
    end do:
  end if:
  if effect = 0 then
       t_ges[1,i] := REMOVE:
       u_ges[1,i] := REMOVE:

       if `mod`(i, 10) = 1 then XX[trunc((1/10)*i)+1, 1] := 0; printf("XX%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 2 then XY[trunc((1/10)*i)+1, 1] := 0; printf("XY%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 3 then XZ[trunc((1/10)*i)+1, 1] := 0; printf("XZ%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 4 then YY[trunc((1/10)*i)+1, 1] := 0; printf("YY%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 5 then YZ[trunc((1/10)*i)+1, 1] := 0; printf("YZ%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 6 then ZZ[trunc((1/10)*i)+1, 1] := 0; printf("ZZ%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 7 then mX[trunc((1/10)*i)+1, 1] := 0; printf("MX%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 8 then mY[trunc((1/10)*i)+1, 1] := 0; printf("MY%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 9 then mZ[trunc((1/10)*i)+1, 1] := 0; printf("MZ%d \n", trunc((1/10)*i)+1): end if:
       if `mod`(i, 10) = 0 then m[trunc((1/10)*i), 1] := 0;    printf(" M%d \n", trunc((1/10)*i)):   end if:
  end if:
end do:
# Speichere markierte Regressor-Vektoren
save t_ges, sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_regressor_marked_maple.m", robot_name):
save u_ges, sprintf("../codeexport/%s/tmp/energy_potential_fixb_regressor_marked_maple.m", robot_name):
# Parametervektor in Richtiger Reihenfolge aufstellen und Entfernen von mZ, YY, m
# Anwendung von [HRL_IDR] Regel 3 (S. 5) auf den Parametervektor
# Bei Drehgelenken werden die Parameter YY, mZ und m mit denen des folgenden Körpers gruppiert. Bei Schubgelenken werden alle Trägheitsmomente gruppiert (da die Trägheitsmomente keinen Einfluss auf die translatorische Bewegung des Schubgelenkes haben).
# [GautierKhalil1990] (eq. 18)
nt:=add(sigma[i,1],i=1..NJ): # Anzahl Translatorische Gelenke
nr:=NJ-nt: # Anzahl rotatorische Gelenke
;
#Anzahl der Einträge im Minimalparametervektor
MPV_n_max := 7*nr+4*nt: # ignoriere Abzüge in der Formel
;
Paramvec := Matrix(MPV_n_max, 1):
ii := 1:
for i to NJ do
  if sigma[i,1] = 0 then # Drehgelenk (eq. 15)
    Paramvec[ii, 1] := XX[i, 1]: ii:= ii+1:
    Paramvec[ii, 1] := XY[i, 1]: ii:= ii+1:
    Paramvec[ii, 1] := XZ[i, 1]: ii:= ii+1:
    Paramvec[ii, 1] := YZ[i, 1]: ii:= ii+1:
    Paramvec[ii, 1] := ZZ[i, 1]: ii:= ii+1:
    Paramvec[ii, 1] := mX[i, 1]: ii:= ii+1:
    Paramvec[ii, 1] := mY[i, 1]: ii:= ii+1:
  else: # Schubgelenk (eq. 16)
    if sin(alpha[i,1]) = 0 and false then # Sonderregel (eq. 19).
      # Parameter wurden mit denen der Basisnäheren Achse zusammengelegt und werden hier nicht betrachtet.
      # Parameter mZ fällt weg, also bleiben mX und mY übrig
      # TODO: Unklar, was das für das Programm heißt. Folgendes noch falsch.
      Paramvec[ii, 1] := mX[i, 1]: ii:= ii+1:
      Paramvec[ii, 1] := mY[i, 1]: ii:= ii+1:
    elif cos(alpha[i,1]) = 0 and false then # Sonderregel (eq. 20)
      # Parameter mY fällt weg, also bleiben mX und mZ übrig
      # TODO: Unklar, was das für das Programm heißt. Folgendes noch falsch.
      Paramvec[ii, 1] := mX[i, 1]: ii:= ii+1:
      Paramvec[ii, 1] := mZ[i, 1]: ii:= ii+1:
    else
      Paramvec[ii, 1] := mX[i, 1]: ii:= ii+1:
      Paramvec[ii, 1] := mY[i, 1]: ii:= ii+1:
      Paramvec[ii, 1] := mZ[i, 1]: ii:= ii+1:
    end if:
    Paramvec[ii, 1] := m[i, 1]:  ii:= ii+1:
  end if:
end do:

# Entfernung von Nullelementen
# Durch Entfernung der Nullelemente erfolgt 'beta_b' aus [HRL_IDR] (22)
p := 0:
for i to MPV_n_max do
  if Paramvec[i, 1] = 0 then
    p := p+1
  end if:
end do:
# Determine size of matrix dependent on number of Zeroes
Paramvec_size := MPV_n_max-p:
Paramvec2 := Matrix(Paramvec_size, 1):
p := 0:
for i to MPV_n_max do
  # count Zeroes
  if Paramvec[i, 1] = 0 then
    p := p+1:
  else
    # Shift by number of currently found Zeroes
    Paramvec2[i-p, 1] := Paramvec[i, 1]:
  end if:
end do:
printf("Dimension des Minimalparametervektors: %dx1\n", Paramvec_size):
Paramvec2;
# Export - Minimalparametervektor
save Paramvec2, sprintf("../codeexport/%s/tmp/minimal_parameter_vector_fixb_maple", robot_name):
save Paramvec2, sprintf("../codeexport/%s/tmp/minimal_parameter_vector_fixb_maple_geom", robot_name): # zum Testen gegen andere Implementierung
if codegen_act then
   MatlabExport(Paramvec2, sprintf("../codeexport/%s/tmp/minimal_parameter_vector_fixb_matlab.m", robot_name), codegen_opt):
end if;
# Minimal geometrievektor t_i und u_i
# Markieren von t_mZ, t_YY, t_m, u_mZ, u_YY, u_m_j
# Anwendung von [HRL_IDR] Regel 3 (S. 5) auf die Geometrievektoren t und u der Energien T und U
# Der Basis-Körper wird nicht gezählt. Daher ignorieren der ersten zehn Einträge in t_ges
# [GautierKhalil1990], p.369 (rechts: "As a conclusion ..."
for i from 0 to NJ-1 do
  if sigma[i+1,1] = 0 then: # Drehgelenk
    t_ges[1,i*10+ 4]:=REMOVE; # YY
    t_ges[1,i*10+ 9]:=REMOVE; # mZ
    t_ges[1,i*10+10]:=REMOVE; # m
    u_ges[1,i*10+ 4]:=REMOVE; # YY
    u_ges[1,i*10+ 9]:=REMOVE; # mZ
    u_ges[1,i*10+10]:=REMOVE; # m
  else: # Schubgelenk
    u_ges[1,i*10+ 1]:=REMOVE; # XX
    u_ges[1,i*10+ 2]:=REMOVE; # XY
    u_ges[1,i*10+ 3]:=REMOVE; # XZ
    u_ges[1,i*10+ 4]:=REMOVE; # YY
    u_ges[1,i*10+ 5]:=REMOVE; # YZ
    u_ges[1,i*10+ 6]:=REMOVE; # ZZ
    t_ges[1,i*10+ 1]:=REMOVE; # XX
    t_ges[1,i*10+ 2]:=REMOVE; # XY
    t_ges[1,i*10+ 3]:=REMOVE; # XZ
    t_ges[1,i*10+ 4]:=REMOVE; # YY
    t_ges[1,i*10+ 5]:=REMOVE; # YZ
    t_ges[1,i*10+ 6]:=REMOVE; # ZZ
    # Sonderregel (eq. 19-20). Parameter zusammengelegt in Parametervektor. Werden daher hier entfernt.
    # Fälle vorerst auskommentiert (siehe oben). TODO: Das hat alles nicht funktioniert.
    if sin(alpha[i+1,1]) = 0 and false then # deaktiviert!
     u_ges[1,i*10+ 7]:=REMOVE; # mX
     t_ges[1,i*10+ 7]:=REMOVE; # mX
     u_ges[1,i*10+ 8]:=REMOVE; # mY
     t_ges[1,i*10+ 8]:=REMOVE; # mY
    end if:
    if cos(alpha[i+1,1]) = 0 and false then # deaktiviert!
     u_ges[1,i*10+ 7]:=REMOVE; # mX
     t_ges[1,i*10+ 7]:=REMOVE; # mX
     u_ges[1,i*10+ 9]:=REMOVE; # mZ
     t_ges[1,i*10+ 9]:=REMOVE; # mZ
    end if:
  end if:
end do:
# Entfernungen von markierten Elementen
# Durch Entfernung der markierten Elemente, äquivalent zum Minimalparametervektor beta_b, haben t und u die gleiche Spaltenanzahl wie 'C_b'.
p:=0:
for i from 1 to 10*NJ do       #Nullen Zählen um die Matrixgröße von t_ges zu bestimmen
  if t_ges[1,i]=REMOVE and u_ges[1,i]=REMOVE then
    p:=p+1;
  end if
end do:
size_Matrix:=10*NJ-p:
t_ges_minpar:=Matrix(1,size_Matrix):
u_ges_minpar:=Matrix(1,size_Matrix):
printf("Dimension der Regressormatrix: 1x%d\n", size_Matrix):

if size_Matrix <> Paramvec_size then
  printf("Achtung: Dimension der Regressormatrix passt nicht zu der des Parametervektors!\n"):
end if:

p:=0:
for i from 1 to 10*NJ do       #Nullen Zählen
   if t_ges[1,i]=REMOVE and u_ges[1,i]=REMOVE then
      p:=p+1;
   else
      t_ges_minpar[1,i-p]:=t_ges[1,i]; #Um die Anzahl der im Iterationsschritt gezählten Nullen verschieben
      u_ges_minpar[1,i-p]:=u_ges[1,i]; #Um die Anzahl der im Iterationsschritt gezählten Nullen verschieben
   end if
end do:
save t_ges_minpar, sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_regressor_minpar_maple.m", robot_name):
save u_ges_minpar, sprintf("../codeexport/%s/tmp/energy_potential_fixb_regressor_minpar_maple.m", robot_name):
# Zum symbolischen Testen des minimierten Regressors gegen den ursprünglichen
(*
read sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_regressor_maple.m", robot_name):
t_ges := t_ges: # Variable wurde im Arbeitsblatt verändert. Neu laden.
read sprintf("../codeexport/%s/tmp/energy_potential_fixb_regressor_maple.m", robot_name):
u_ges := u_ges:
# Muss Null sein.
test_u := simplify(u_ges_minpar.Paramvec2 - u_ges.PV2_vec(11..,..)):
test_t := simplify(t_ges_minpar.Paramvec2 - t_ges.PV2_vec(11..,..)):
*)
;
# Export
if codegen_act then
  MatlabExport(convert_t_s(t_ges_minpar), sprintf("../codeexport/%s/tmp/energy_kinetic_fixb_regressor_minpar_matlab.m", robot_name), codegen_opt):
end if:
if codegen_act then
  MatlabExport(convert_t_s(u_ges_minpar), sprintf("../codeexport/%s/tmp/energy_potential_fixb_regressor_minpar_matlab.m", robot_name), codegen_opt):
end if:


text:
with(DocumentTools):
with(StringTools):

SetProperty("Meter_sian", value, 100):
SetProperty("LogAreaSIAN", value, ""):
SetProperty("LogAreaSE", value, ""):
SetProperty("LogAreaME", value, ""):
SetProperty("GlobalParams1", expression, NULL):
SetProperty("LocalParams1", expression, NULL):
SetProperty("NoIDParams1", expression, NULL):
SetProperty("Bound", expression, NULL):
SetProperty("MultiFunctions", expression, NULL):
SetProperty("SingleFunctions", expression, NULL):
SetProperty("being_refined", caption, ""):

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

LogParameters := proc(s, target)
	local logsofar;
	logsofar := DocumentTools:-GetProperty(target, value);
	if logsofar <> "" then
		logsofar := logsofar, "\n";
	end if;
	DocumentTools:-SetProperty(target, value, cat(logsofar,s), 'refresh');
end proc:

ClearLog := proc(target)

		DocumentTools:-SetProperty(target, value, "");

end proc:
ClearLog("Parameters"):

skip_single:= parse(GetProperty("SkipSingle", value)):
computeId:= parse(GetProperty("ComputeId", value)):

simplified_generators := parse(GetProperty("SimplifiedGen", value)):
simplify_bound := parse(GetProperty("Refine", value)):
no_bound:= parse(GetProperty("NoBound", value)):
use_char := parse(GetProperty("use_char", value)):
err := 0:
#try:
	sigma:="":
	sigma:=GetProperty("sigma", value):
	if SearchText("diff", sigma)>0 then
		LogText(sprintf("%q\n", "Using Maple input format"), "LogAreaSIAN"):
		sigma := map(x->parse(x), Split(sigma, ";")):
	else
		LogText(sprintf("Using text-based input format"), "LogAreaSIAN"):
		sigma := Split(sigma, ";"):
		states := map(x->Trim(RegSubs("d([a-zA-Z0-9]+)/dt(.*)" = "\\1", x)), select(x->Search("/dt", x)>0, sigma)):
		state_eqs := select(x->Search("/dt", x)>0, sigma):

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

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

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

		sigma := [op(state_eqs), op(output_eqs)]:
	end if:
	LogExpression(sprintf("%q \n", sigma), "LogAreaSIAN"):
(*catch:
	SetProperty("LogAreaSIAN", value, "Error parsing system: check syntax"): # LogText(
	err:=1:
	SetProperty("GlobalParams1", expression, SystemParseError("see Logs")):
	SetProperty("LocalParams1", expression,  SystemParseError("see Logs")):
	SetProperty("NoIDParams1", expression, SystemParseError("see Logs")):
	SetProperty("Bound", expression,  SystemParseError("see Logs")):
	SetProperty("MultiFunctions", expression,  SystemParseError("see Logs")):
	SetProperty("SingleFunctions", expression,  SystemParseError("see Logs")):
	SetProperty("Meter_sian", value, 0):
	error "Error parsing system: check syntax":
end try:*)

#try:
	params_to_assess := [seq(parse(expr), expr in StringTools[Split](GetProperty("params", value), ";"))]:
	LogExpression(sprintf("%q \n", params_to_assess), "LogAreaSIAN"):
(*catch:
	SetProperty("LogAreaSIAN", value, "Error Parsing Parameters"): # LogText(
	err:=1:
	SetProperty("GlobalParams1", expression, ParameterParseError("see Logs")):
	SetProperty("LocalParams1", expression,  ParameterParseError("see Logs")):
	SetProperty("NoIDParams1", expression, ParameterParseError("see Logs")):
	SetProperty("Bound", expression,  ParameterParseError("see Logs")):
	SetProperty("MultiFunctions", expression,  ParameterParseError("see Logs")):
	SetProperty("SingleFunctions", expression,  ParameterParseError("see Logs")):
	SetProperty("Meter_sian", value, 0):
	error "Error parsing parameters: check syntax":
end try:*)

#try:
	replicas := parse(GetProperty("replicas", value)):
	if replicas <=1 then
		initial_conditions:= true:
	else
		sigma := GenerateReplica(sigma, replicas):
		initial_conditions:=false:
	end if:

	if nops(params_to_assess)=0 then
		LogText(sprintf("%q\n", "Checking all parameters"), "LogAreaSIAN"):
		params_to_assess := GetParameters(sigma, initial_conditions):
	end if:
	p:=parse(GetProperty("p", value)):

(*catch:
	SetProperty("LogAreaSIAN", value, "SYNTAX ERROR PLEASE CHECK INPUTS"): # LogText(
	err:=1:
	SetProperty("GlobalParams1", expression, Error("see Logs")):
	SetProperty("LocalParams1", expression, Error("see Logs")):
	SetProperty("NoIDParams1", expression, Error("see Logs")):
	SetProperty("Bound", expression, Error("see Logs")):
	SetProperty("MultiFunctions", expression, Error("see Logs")):
	SetProperty("SingleFunctions", expression, Error("see Logs")):
	SetProperty("Meter_sian", value, 0):
end try:*)

if err=0 then
	run_sian := parse(GetProperty("RunSIAN", value)):
	if run_sian then
		SetProperty("GlobalParams1", expression, Running):
		SetProperty("LocalParams1", expression, Running):
		SetProperty("NoIDParams1", expression, Running):
	fi:
	if computeId then
		if not no_bound then
			SetProperty("Bound", expression, Running):
		fi:
		if not skip_single then
			SetProperty("SingleFunctions", expression, Running):
		fi:
		SetProperty("MultiFunctions", expression, Running):
	end if:

	output_targets_multi:=table([log="LogAreaME", multi="MultiFunctions", bound_="Bound", runningtime="RunningTimeMulti", single="SingleFunctions",
					  	    run_single_=run_single, runningtimesingle="RunningTimeSingle"]):
	output_targets_single:=table([log="LogAreaSE", single="SingleFunctions", runningtime="RunningTimeSingle"]):
	output_targets_sian:=table([globalparams="GlobalParams1", localparams= "LocalParams1", noidparams="NoIDParams1",
						   runningtime="RunningTimeSIAN", log="LogAreaSIAN", parameters="Parameters"]):
	bypass_ := parse(GetProperty("bypass", value)):

	count_solutions := parse(GetProperty("printSolutions", value)):

	if run_sian then
		if use_char then
			char := 2^29 -3:
		else
			char := 0:
		end if:
		out_sian := timed_SIAN(sigma,  params_to_assess, 0.99, output_targets_sian, count_solutions, char):
		basis_sian:=out_sian[basis]:
		ordering_sian:=out_sian[varordering]:
		all_other:=out_sian[allvars]:
		identified_params:=select(x->evalb(x in out_sian[parameters]), out_sian[globally]):
		max_perms:= parse(GetProperty("MaxPermutations", value)):
		if bypass_ then
			if nops(identified_params)=nops(out_sian[parameters]) then
				if computeId then
					LogText(sprintf("Bypassing Multi-Experiment check via SIAN"), "LogAreaME"):
				end if:
				computeId:=false;
				se_me_bypass_time := time():
				SetProperty("Bound", value, 1):
				SetProperty("MultiFunctions", value, convert(out_sian[parameters], list)):
				SetProperty("RunningTimeMulti", value, convert(time() - se_me_bypass_time, string), 'refresh'): # time

				if not skip_single then
					LogText(sprintf("Bypassing Single-Experiment check via SIAN"), "LogAreaSE"):
				end if:
				skip_single:=true:
				se_me_bypass_time := time():
				SetProperty("SingleFunctions", value, convert(out_sian[parameters], list)):
				SetProperty("RunningTimeSingle", value, convert(time() - se_me_bypass_time, string), 'refresh'): # time
			fi:
		end if:
	fi:
	if computeId then
		out_multi := timed_Multi(sigma, simplified_generators, no_bound, simplify_bound, max_perms, output_targets_multi):
		generators := convert(out_multi[3],list):
		finish_bypass_multi:=out_multi[2]:
		bound := out_multi[1]: # bound value
		if bound=1  then
			if not skip_single then
				LogText(sprintf("Deterministically bypassing Single-Experiment via Multi-Experiment check, bound from Multi-Experiment = 1"), "LogAreaSE"):
			end if:
			skip_single:=true:
			SetProperty("RunningTimeSingle", value, finish_bypass_multi, 'refresh'): # time
		fi:
		if not skip_single then
			out_single := timed_Single(sigma, output_targets_single):
			if out_single = generators and bound>1 then
				LogText(sprintf("Deterministically overriding Multi-Experiment via Single-Experiment check,\n bound from Multi-Experiment > 1, single- and multi-experiment generators are the same"), "LogAreaSE"):
				SetProperty("MultiFunctions", value, convert(out_single, list)):
				SetProperty("Bound", value, 1):
			fi:
			LogExpression(sprintf("%q %q\n",out_single, generators), "LogAreaSE");
		fi:
	fi:
	SetProperty("Meter_sian", value, 0):
end if:

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
# simplicial_interval(n) is just the set [n] = {0,1,...,n}

`is_element/simplicial_interval` := (n::nonnegint) -> proc(k)
 type(k,nonnegint) and k <= n;
end:

`is_equal/simplicial_interval` := (n::nonnegint) -> (a,b) -> evalb(a = b);

`is_leq/simplicial_interval` := (n::nonnegint) -> (a,b) -> evalb(a <= b);

`list_elements/simplicial_interval` := proc(n::nonnegint)
 local i;
 [seq(i,i=0..n)];
end:

`count_elements/simplicial_interval` := (n::nonnegint) -> n+1;

`random_element/simplicial_interval` := (n::nonnegint) -> proc()
 rand(0..n)();
end:

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

`is_element/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 global reason;
 local N,M,i;

 N := {seq(i,i=0..n)};
 M := {seq(i,i=0..m)};

 if not(`is_element/maps`(N,M)(f)) then
  reason := [convert(procname,string),"not a map from [n] to [m]",reason];
  return false;
 fi;

 for i from 0 to n-1 do
  if not(f[i] <= f[i+1]) then
   reason := [convert(procname,string),"f[i+1] > f[i]",i,f[i+1],f[i]];
   return false;
  fi;
 od;

 return true;
end:

`to_list/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 local i;
 return [seq(f[i],i=0..n)];
end:

`is_equal/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f,g)
 evalb(`to_list/simplicial_maps`(n,m)(f) = `to_list/simplicial_maps`(n,m)(g));
end:

`is_leq/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f,g)
 local u;
 u := `to_list/simplicial_maps`(n,m)(g) -~ `to_list/simplicial_maps`(n,m)(f);
 evalb( min(op(u)) >= 0 );
end:

# This assumes that
# *  p is a subset of {1,...,n} of size k
# *  q is a subset of {0,...,m} of size k+1
# It returns the nondecreasing map f : {0,...,n} -> {0,...,m} such that
# *  The image of f is q
# * { i > 0 : f(i) > f(i-1) } = p.

`build/simplicial_maps` :=  (n::nonnegint,m::nonnegint) -> proc(k,p,q)
 local pp,i,j,f;
 pp := [0,op(p),n+1];
 f := table();
 i := 0;
 for i from 0 to k do
  for j from pp[i+1] to pp[i+2]-1 do
   f[j] := q[i+1];
  od;
 od;

 return eval(f);
end:

`random_element/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc()
 local p,q,i,j,k,l,f;

 k := rand(0..min(n,m))();
 p := combinat[randcomb](n,k);
 q := combinat[randcomb](m+1,k+1) -~ 1;
 return `build/simplicial_maps`(n,m)(k,p,q);
end:

`list_elements/simplicial_maps` := proc(n::nonnegint,m::nonnegint)
 option remember;
 local i;

 if n = 0 then
  return [seq(`C/simplicial_maps`(0,i),i=0..m)];
 else
  if m = 0 then
   return [`C/simplicial_maps`(n,0)];
  else
   return [op(map(`I/simplicial_maps`(n,m-1),`list_elements/simplicial_maps`(n,m-1))),
           op(map(`T/simplicial_maps`(n-1,m),`list_elements/simplicial_maps`(n-1,m)))];
  fi;
 fi;
end:

`count_elements/simplicial_maps` := (n::nonnegint,m::nonnegint) -> binomial(n+m+1,m);

# Constant maps
`C/simplicial_maps` := proc(n::nonnegint,k::nonnegint)
 local i;
 table([seq(i=k,i=0..n)]);
end:

# Inclusion Delta(n,m) -> Delta(n,m+1)
`I/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 local i;
 table([seq(i=f[i],i=0..n)]);
end:

# Map Delta(n,m) -> Delta(n+1,m): extend by sending n+1 to m
`T/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 local i;
 table([seq(i=f[i],i=0..n),n+1=m]);
end:

# `delta/simplicial_maps`(n)(i) : [n] >-> [n+1]; image omits i
`delta/simplicial_maps` := proc(n::nonnegint,i::nonnegint)
 local j;
 if i > n+1 then return FAIL; fi;
 return table([seq(j=j,j=0..i-1),seq(j=j+1,j=i..n)]);
end:

# `sigma/simplicial_maps`(n)(i) : [n] ->> [n-1]; takes the value i twice
`sigma/simplicial_maps` := proc(n::nonnegint,i::nonnegint)
 local j;
 if i > n-1 then return FAIL; fi;
 return table([seq(j=j,j=0..i),seq(j=j-1,j=i+1..n)]);
end:

`id/simplicial_maps` := proc(n::nonnegint)
 local i;
 return table([seq(i=i,i=0..n)]);
end:

# `compose/simplicial_maps`(n,m,p)(f,g) assumes that f:[n] -> [m] and g:[m] -> [p],
# and it returns the composite g o f : [n] -> [p].

`compose/simplicial_maps` := (n::nonnegint,m::nonnegint,p::nonnegint) -> proc(f,g)
 local i;
 table([seq(i = g[f[i]],i=0..n)]);
end:

`is_wide/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 return evalb(f[0] = 0 and f[n] = m);
end:

`is_mono/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 local i;

 if m < n then return false; fi;
 for i from 0 to n-1 do
  if f[i] = f[i+1] then return false; fi;
 od;
 return true;
end:

`is_epi/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 local i,j;

 if n < m then return false; fi;
 if f[0] <> 0 then return false; fi;
 if f[n] <> m then return false; fi;

 j := 0;
 for i from 1 to n do
  if (f[i] <> j and f[i] <> j+1) then
   return false;
  fi;
  j := f[i];
 od;
 return true;
end:

# This assumes that f : [n] -> [m], and it returns [k,g,h] such that
# g : [n] ->> [k] and h : [k] >-> [m] and f = h o g

`factor/simplicial_maps` := (n::nonnegint,m::nonnegint) -> proc(f)
 local K,k,g,h,hi,i;

 K := sort([op({seq(f[i],i=0..n)})]);
 k := nops(K) - 1;
 h := table([seq(i=K[i+1],i=0..k)]);
 hi := table();
 for i from 0 to k do hi[h[i]] := i; od;
 g := table([seq(i = hi[f[i]],i=0..n)]);
 return [k,eval(g),eval(h)];
end:

`describe/simplicial_maps` := (n::nonnegint) -> proc(f)
 local i;
 cat(seq(nat_code[f[i]],i=0..n));
end:

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

`is_element/simplicial_epi` := (n::nonnegint,m::nonnegint) -> proc(f)
 global reason;

 if not(`is_element/simplicial_maps`(n,m)(f)) then
  reason := [convert(procname,string),"not a simplicial map from [n] to [m]",reason];
  return false;
 fi;

 return `is_epi/simplicial_maps`(n,m)(f);
end:

`is_equal/simplicial_epi` :=  (n::nonnegint,m::nonnegint) -> proc(f,g)
 `is_equal/simplicial_maps`(n,m)(f,g);
end:

`is_leq/simplicial_epi` :=  (n::nonnegint,m::nonnegint) -> proc(f,g)
 `is_leq/simplicial_maps`(n,m)(f,g);
end:

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

 if n < m then return FAIL; fi;

 k := m;
 p := combinat[randcomb](n,k);
 q := [seq(i,i=0..m)];
 return `build/simplicial_maps`(n,m)(k,p,q);
end:

`list_elements/simplicial_epi` := proc(n::nonnegint,m::nonnegint)
 local P,L,S,f,i,j;

 P := combinat[choose]([seq(i,i=1..n)],m);
 L := [];
 for S in P do
  f := table();
  for i from 0 to S[1]-1 do f[i] := 0; od;
  for j from 1 to m-1 do
   for i from S[j] to S[j+1]-1 do f[i] := j; od;
  od;
  for i from S[m] to n do f[i] := m; od;
  L := [op(L),eval(f)];
 od:

 return L;
end:

`count_elements/simplicial_epi` := (n::nonnegint,m::nonnegint) -> binomial(n,m);

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

`is_element/simplicial_mono` := (n::nonnegint,m::nonnegint) -> proc(f)
 global reason;

 if not(`is_element/simplicial_maps`(n,m)(f)) then
  reason := [convert(procname,string),"not a simplicial map from [n] to [m]",reason];
  return false;
 fi;

 return `is_mono/simplicial_maps`(n,m)(f);
end:

`is_equal/simplicial_mono` :=  (n::nonnegint,m::nonnegint) -> proc(f,g)
 `is_equal/simplicial_maps`(n,m)(f,g);
end:

`is_leq/simplicial_mono` :=  (n::nonnegint,m::nonnegint) -> proc(f,g)
 `is_leq/simplicial_maps`(n,m)(f,g);
end:

`random_element/simplicial_mono` := (n::nonnegint,m::nonnegint) -> proc()
 local p,q,i,j,k,l,f;

 if n > m then return FAIL; fi;

 k := n;
 p := [seq(i,i=1..n)];
 q := combinat[randcomb](m+1,k+1) -~ 1;
 return `build/simplicial_maps`(n,m)(k,p,q);
end:

`list_elements/simplicial_mono` := proc(n::nonnegint,m::nonnegint)
 local P,L,S,f,i,j;

 P := combinat[choose]([seq(i,i=0..m)],n+1);
 L := [];
 for S in P do
  f := table();
  for i from 0 to n do f[i] := S[i+1]; od;
  L := [op(L),eval(f)];
 od:

 return L;
end:

`count_elements/simplicial_mono` := (n::nonnegint,m::nonnegint) -> binomial(m+1,n+1);


text:
(*
 Copyright (C) 2017 M.A.L. Marques
               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: lda_exc *)

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

bn05_A  := 3.4602:
bn05_C0 := 3.2:
bn05_C1 := -0.9:

f := (rs, z) ->
  f_lda_x_yukawa(rs, z) + f_pw(rs, z)*bn05_A/(bn05_C0 + bn05_C1*rs + rs^2):


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

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

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

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

$ifdef lda_c_pz_params
params_a_gamma := [-0.1423, -0.0843]:
params_a_beta1 := [ 1.0529,  1.3981]:
params_a_beta2 := [ 0.3334,  0.2611]:
params_a_a     := [ 0.0311,  0.01555]:
params_a_b     := [-0.048,  -0.0269]:
params_a_c     := [ 0.0020,  0.0007]:
params_a_d     := [-0.0116, -0.0048]:
$endif

(* Equation C3 *)
ec_low  := (i, rs) -> params_a_gamma[i] / \
        (1 + params_a_beta1[i]*sqrt(rs) + params_a_beta2[i]*rs):

(* Equation [1].C5 *)
ec_high := (i, rs) -> params_a_a[i]*log(rs) + params_a_b[i] \
        + params_a_c[i]*rs*log(rs) + params_a_d[i]*rs:

ec := (i, x) -> my_piecewise3(x >= 1, ec_low(i, x), ec_high(i, x)):

f_pz := (rs, zeta) -> \
 ec(1, rs) + (ec(2, rs) - ec(1, rs))*f_zeta(zeta):

f := (rs, zeta) -> f_pz(rs, zeta):


text:
func $foo ( var %i i32, var %j i32, var %k i32) i32 {
   return (
     intrinsicop i32 JSOP_STRICTEQ (
      dread i32 %i,
      intrinsicop i32 JS_NUMBER (dread i32 %k)))
}

max_stars_repo_name: MapleSystem/OpenArkCompiler

max_issues_repo_name: MapleSystem/OpenArkCompiler

max_forks_repo_name: MapleSystem/OpenArkCompiler


text:
make_octahedron_complex := proc()
 global octahedron_complex,cube_complex;
 local OC,CC,vO,vC,vOi,vCi,i,j,k,dp,CF,f,G,G0,G1,r,g,h,eqs,sol,sc,OS,s,t;

 CC := table():
 OC := table():

 CC["vertices"] := [seq(i,i=1..8)];
 OC["vertices"] := [seq(i,i=1..6)];

 vO := table([
  1 = [1,0,0], 2 = [0,1,0], 3 = [0,0,1], 4 = [0,0,-1], 5 = [0,-1,0], 6 = [-1,0,0]
 ]);

 vC := table([
  1 = [ 1, 1, 1], 2 = [ 1, 1,-1], 3 = [ 1,-1, 1], 4 = [ 1,-1,-1],
  5 = [-1, 1, 1], 6 = [-1, 1,-1], 7 = [-1,-1, 1], 8 = [-1,-1,-1]
 ]);

 vOi := table():
 vCi := table():

 for i from 1 to 6 do vOi[vO[i]] := i; od;
 for i from 1 to 8 do vCi[vC[i]] := i; od;

 OC["embedding_dim"] := 3;
 CC["embedding_dim"] := 3;

 OC["embedding"]            := eval(vO);
 OC["normalised_embedding"] := copy(vO);
 OC["dual_embedding"]       := copy(vO);

 CC["embedding"]            := eval(vC);
 CC["normalised_embedding"] := table([seq(i = vC[i] *~ (sqrt(3)/3), i = 1 .. 8)]);
 CC["dual_embedding"]       := table([seq(i = vC[i] /~ 3, i = 1 .. 8)]);

 dp := (u,v) -> add(u[i] * v[i],i=1..3);

 OC["faces"]   := [seq(select(i -> dp(vO[i],vC[j]) = 1,OC["vertices"]),j = 1..8)];
 CC["squares"] := [seq(select(j -> dp(vO[i],vC[j]) = 1,CC["vertices"]),i = 1..6)];

 OC["edges"] := {op(map(op,map(f -> [[f[1],f[2]],[f[1],f[3]],[f[2],f[3]]],OC["faces"])))};
 OC["edges"] := sort([op(OC["edges"])]);

 CC["short_edges"] := [[1,2],[1,3],[1,5],[1,6],
                      [2,4],[2,6],[3,4],[3,7],[5,6],[5,7],
		      [4,8],[5,8],[6,8],[7,8]];

 CC["long_edges"]  := [[1,4],[1,6],[1,7],[2,8],[3,8],[5,8]];

 CC["edges"] := sort([op(CC["short_edges"]),op(CC["long_edges"])]);
 sc := [1,6,4,2,7,5,3,8];
 CC["six_cycle"] := sc;

 CF := NULL;
 OS := NULL;
 for s in CC["squares"] do
  i := op({op(s)} minus {seq(sc[j],j in s)});
  if member(1,s) then
   CF := CF,[1,s[2],s[4]],[1,s[3],s[4]];
   OS := OS,[1,i,sc[i],sc[sc[i]]];
  else
   CF := CF,[s[1],s[2],8],[s[1],s[3],8];
   OS := OS,[sc[sc[i]],sc[i],i,8];
  fi;
 od:
 OS := [OS];

 CC["faces"] := sort([CF]);
 CC["oriented_squares"] := OS;
 CC["square_plot"] :=
   display(
    seq(polygon(map(j -> vC[j],OS[i]),colour=standard_colour(i)),i=1..6),
    scaling=constrained,axes=none);

 OC["max_simplices"] := OC["faces"];
 CC["max_simplices"] := CC["faces"];

 OC["edge_index"]   := make_index(OC["edges"]):
 OC["face_index"]   := make_index(OC["faces"]):

 CC["edge_index"]   := make_index(CC["edges"]):
 CC["face_index"]   := make_index(CC["faces"]):

 `plot/simplicial_complex`(OC);
 `plot/simplicial_complex`(CC);

 OC["edge_centres"] :=
  map(e -> (vO[e[1]] +~ vO[e[2]]) /~ 2,OC["edges"]);

 OC["face_centres"] :=
  map(f -> (vO[f[1]] +~ vO[f[2]] +~ vO[f[3]]) /~ 3, OC["faces"]);

 CC["short_edge_centres"] :=
  map(e -> (vC[e[1]] +~ vC[e[2]]) /~ 2,CC["short_edges"]);

 CC["square_centres"] :=
  map(f -> (vC[f[1]] +~ vC[f[2]] +~ vC[f[3]] +~ vC[f[4]]) /~ 4, CC["squares"]);

 f := (u) -> simplify(rationalize(u /~ sqrt(add(u[i]^2,i=1..3))));

 OC["poles"] := table([
  2 = map(f,OC["edge_centres"]),
  3 = map(f,OC["face_centres"]),
  4 = map(i -> f(vO[i]),OC["vertices"])
 ]);

 CC["poles"] := copy(OC["poles"]);

 OC["pole_plots"] := table([
  2 = map(u -> line(1.1 *~ u, -1.1 *~ u,color=green,thickness=4),CC["poles"][2]),
  3 = map(u -> line(1.1 *~ u, -1.1 *~ u,color=red  ,thickness=4),CC["poles"][3]),
  4 = map(u -> line(1.2 *~ u, -1.2 *~ u,color=blue ,thickness=4),CC["poles"][4])
 ]);

 OC["all_poles_plot"] := display(
  op(OC["pole_plots"][2]),op(OC["pole_plots"][3]),op(OC["pole_plots"][4]),
  scaling=constrained,axes=none
 );

 CC["pole_plots"] := table([
  2 = map(u -> line(1.6 *~ u, -1.6 *~ u,color=green,thickness=4),CC["poles"][2]),
  3 = map(u -> line(1.9 *~ u, -1.9 *~ u,color=red  ,thickness=4),CC["poles"][3]),
  4 = map(u -> line(1.2 *~ u, -1.2 *~ u,color=blue ,thickness=4),CC["poles"][4])
 ]);

 CC["all_poles_plot"] := display(
  op(CC["pole_plots"][2]),op(CC["pole_plots"][3]),op(CC["pole_plots"][4]),
  scaling=constrained,axes=none
 );

 G0 := [];
 G1 := [[2,4,1,3,6,8,5,7],[1,3,5,7,2,4,6,8]];
 G := [[1,2,3,4,5,6,7,8]];

 while G <> G0 do
  G0 := G;
  G := sort([op({seq(seq(`o/permutations`(8)(g,h),g in G1),h in G)})]);
 od:

 r := table([1 = 1, 2 = 2, 3 = 3, 4 = 4, 5 = 4, 6 = 3, 7 = 2, 8 = 1]):

 CC["vertex_action"]     := table():
 CC["short_edge_action"] := table():
 CC["square_action"]     := table():
 CC["rotation_matrix"]   := table():

 OC["vertex_action"]     := table():
 OC["edge_action"]       := table():
 OC["face_action"]       := table():
 OC["rotation_matrix"]   := table():

 for s in G do
  t := [r[s[1]],r[s[2]],r[s[3]],r[s[4]]];
  CC["vertex_action"][t] := s;
  OC["face_action"][t]   := s;

  CC["short_edge_action"][t] :=
   map(e -> CC["short_edge_index"][sort([s[e[1]],s[e[2]]])],CC["short_edges"]);

  CC["square_action"][t] :=
   map(e -> CC["square_index"][sort([s[e[1]],s[e[2]],s[e[3]],s[e[4]]])],CC["squares"]);

  OC["vertex_action"][t] := CC["square_action"][t];

  OC["edge_action"][t] :=
   map(e -> OC["edge_index"][sort([s[e[1]],s[e[2]]])],OC["edges"]);

  g := Matrix(3,3,[seq(x[i],i=1..9)]);
  eqs := map(op,[seq(convert(g . Vector(vC[i]) - Vector(vC[s[i]]),list),i = 1..8)]);
  sol := solve(eqs);
  CC["rotation_matrix"][t] := subs(sol,convert( g,listlist));
  OC["rotation_matrix"][t] := CC["rotation_matrix"][s];
 od:

 octahedron_complex := eval(OC);
 cube_complex := eval(CC);
 return eval(OC);
end:

make_octahedron_complex():

octosphere_complex := proc(n::nonnegint)
 local T,T0,E,v,x,i,j;

 T := eval(octahedron_complex);

 for i from 1 to n do
  T0 := eval(T);
  T := eval(`triangular_subdivision/simplicial_complex`(T0));
  T0 := eval(T);
  T := eval(`condense/simplicial_complex`(T0));
  E := T["embedding"];
  for v in T["vertices"] do
   x := E[v];
   x := evalf(x /~ sqrt(add(x[j]^2,j=1..3)));
   E[v] := x;
  od:
 od:

 return eval(T):
end:


text:
`eta/BW` := eval(`eta/prime_simplex`);
`gamma/BW` := eval(`gamma/prime_simplex`);
`gamma_min/BW` := eval(`gamma_min/prime_simplex`);


text:
type $Ljava/lang/Object; <class  {
  @DEXINFO_isuserclass 0}>
type $Leh; <class <$Ljava/lang/Object;> {
  @DEXINFO_isuserclass 1,
  &foo public static (<* <$Leh;>>,<* <[0] <* <$Ljava/lang/String;>>>>) void}>
type $Ljava/lang/String; <class  {
  @DEXINFO_isuserclass 0}>
type $Ljava/io/PrintStream; <class  {
  @DEXINFO_isuserclass 0}>
type $Ljava/lang/System; <class  {
  @DEXINFO_isuserclass 0,
  @out <* <$Ljava/io/PrintStream;>>}>
type $Ljava/lang/Exception; <class  {
  @DEXINFO_isuserclass 0}>
  javaclass $javaclass|Leh; <$Leh;>

func &foo public static (var %_this <* <$Leh;>>, var %Reg2_R39 <* <[0] <* <$Ljava/lang/String;>>>>) void
var $const_array_0 <[3] i32> readonly = [1, 2, 3]
var $static|Ljava/lang/System;|out|Ljava/io/PrintStream;| <* <$Ljava/io/PrintStream;>>
func &Ljava/io/PrintStream;|print|_I_V| (var %_this <* <$Ljava/io/PrintStream;>>, var %Reg1_I i32) void
var $const_chars_1 <[11] u8> const readonly = [0, 0, 7, 0, 101, 114, 114, 111, 114, 49, 10]
func &Ljava/io/PrintStream;|print|_Ljava/lang/String;_V| (var %_this <* <$Ljava/io/PrintStream;>>, var %Reg1_R37 <* <$Ljava/lang/String;>>) void
var $const_chars_2 <[11] u8> const readonly = [0, 0, 7, 0, 101, 114, 114, 111, 114, 50, 10]

func &foo public static (var %_this <* <$Leh;>>, var %Reg2_R39 <* <[0] <* <$Ljava/lang/String;>>>>) void {
  var %Reg0_I i32
  var %Reg0_R45 <* <[3] i32>>
  var %Reg1_I i32
  var %Reg1_R47 <* <$Ljava/io/PrintStream;>>
  var %Reg0_R51 <* <$Ljava/lang/Exception;>>
  var %Reg0_R47 <* <$Ljava/io/PrintStream;>>
  var %Reg1_R53 <* u8>
  var %Reg0_R55 <* dynany>
  dassign %Reg0_I 0 (cvt i32 i8 (constval i8 3))
  javatry { @label1 }
  dassign %Reg0_R45 0 (malloc ptr (mul i32 (dread i32 %Reg0_I, sizeoftype u32 i32)))
  endtry
  intrinsiccallassigned DEX_ARRAY_SETELEM (dread ptr %Reg0_R45, constval i32 3, addrof ptr $const_array_0) {}
  dassign %Reg1_I 0 (cvt i32 i8 (constval i8 4))
  javatry { @label1 @label2 }
  dassign %Reg0_I 0 (iread i32 <* <[3] i32>> 0 (array ptr <* <[3] i32>> (dread ptr %Reg0_R45, dread i32 %Reg1_I)))
  dassign %Reg1_R47 0 (dread ptr $static|Ljava/lang/System;|out|Ljava/io/PrintStream;|)
  virtualcallassigned &Ljava/io/PrintStream;|print|_I_V| (dread ptr %Reg1_R47, dread i32 %Reg0_I) {}
  endtry
@label0   return ()
@label2   javacatch <* <$Ljava/lang/Exception;>>
  dassign %Reg0_R51 0 (regread ptr %%thrownval)
  javatry { @label1 }
  dassign %Reg0_R47 0 (dread ptr $static|Ljava/lang/System;|out|Ljava/io/PrintStream;|)
  dassign %Reg1_R53 0 (addrof ptr $const_chars_1)
  virtualcallassigned &Ljava/io/PrintStream;|print|_Ljava/lang/String;_V| (dread ptr %Reg0_R47, dread ptr %Reg1_R53) {}
  endtry
  goto @label0
@label1   javacatch <* <$Ljava/lang/Exception;>>
  dassign %Reg0_R51 0 (regread ptr %%thrownval)
  dassign %Reg0_R47 0 (dread ptr $static|Ljava/lang/System;|out|Ljava/io/PrintStream;|)
  dassign %Reg1_R53 0 (addrof ptr $const_chars_2)
  virtualcallassigned &Ljava/io/PrintStream;|print|_Ljava/lang/String;_V| (dread ptr %Reg0_R47, dread ptr %Reg1_R53) {}
  goto @label0
}


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

C0 := (1 - log(2))/(2*Pi^2): (* Equation (9) *)
C1 := 4*C0/3:                (* Equations (13), (28), (33) *)
C2 := RS_FACTOR:             (* Equation (8) *)
C3 := C2/3:

(* several cutoffs *)
big     := 1.0e4:
cutoff  := 1.0e7:
ei_xmax := 7.0183341467e+02:

(* Equation (39) *)
kssp0_k0 := 1.291551074:
kssp0_k1 := 0.349064173:
kssp0_r1 := 0.08327588:

kssp0 := rs ->
  kssp0_k0 - kssp0_k1*(1 - exp(-kssp0_r1*rs^(4/5))):

(* Equation (45) *)
fssp_A1 := 1.622118767:
fssp_A2 := 0.489958076:
fssp_A3 := 1.379021941:

fssp := (rs, gr) ->
  (1 + fssp_A1*gr + fssp_A2^2*gr^2)*exp(-fssp_A2^2*gr^2)/sqrt(1 + fssp_A3*gr/rs):

(* Equation (34) *)
fa_a1 := 0.939016:
fa_a2 := 1.733170:

f_factor := rs -> exp(-rs^2/(fa_a1*sqrt(rs) + fa_a2*rs)^2):

(* Equation (40) *)
kss0_k0 :=  1.200801774:
kss0_k1 :=  0.859614445:
kss0_k2 := -0.812904345:
kss0_r1 :=  1.089338848:
kss0_r2 :=  0.655638823:

kss0 := (rs, gr) ->
  + kss0_k0
  + kss0_k1*(1 - exp(-kss0_r1*sqrt(rs)))
  + kss0_k2*(1 - exp(-kss0_r2*rs^(2/5))):

fss_A4 := 4.946281353:
fss_A5 := 3.600612059:

fss := (rs, gr) ->
  (1 + fss_A4^2*gr^2)*exp(-fss_A4^2*gr^2)/sqrt(1 + fss_A5*gr/rs):

(* Equation (15) *)
eq15 := mu -> (3 + 2*(sqrt(mu) + mu))/(3 + 6*(sqrt(mu) + mu)):

f_eab := mu ->
  C0*(-xc_E1_scaled(mu)*(1 + 2*mu*eq15(mu)) + 2*eq15(mu)):

(* Equation (13) *)
(*
   This is numerically suboptimal - the max function cuts off f_eab to
   around 10^-6. This is too large. The other possibility that was used
   before in the code was to put term to zero if mu > ei_xmax
*)
mu_ba := (rsa, ga2) -> C1*rsa/m_max(kssp0(rsa)^2*fssp(rsa, ga2)^2, 1e-60):
term1 := (rsa, z, ga2) -> my_piecewise3(mu_ba(rsa, ga2) >= cutoff, 0,
  f_eab(mu_ba(rsa, ga2))*(1 - z)/2):

mu_aa := (rsa, ga2) -> C1*rsa/m_max(kss0(rsa)^2*fss(rsa, ga2)^2, 1e-60):
term2 := (rsa, z, ga2) -> my_piecewise3(mu_aa(rsa, ga2) >= cutoff, 0,
  f_eab(mu_aa(rsa, ga2))*f_factor(rsa)*(1 + z)/2):

f_ft97 := (rs, z, xs) ->
  + term1(rs*(2/(1 + z))^(1/3), z, C3^2*xs^2)
  + term2(rs*(2/(1 + z))^(1/3), z, C3^2*xs^2):

f  := (rs, z, xt, xs0, xs1) -> f_ft97(rs, z, xs0) + f_ft97(rs, -z, xs1):


text:
single_case_Partition := module()
  uses Domain;

  export SimplName  := "Single case partition";
  export SimplOrder := 11;

    export ModuleApply := proc(vs :: DomBound, sh :: DomShape, $)
        subsindets(sh, DomSplit, x->do_simp(op(x)));
    end proc;

    local do_simp := proc(p:: Partition,$)::DomShape;
        local r := Partition:-Simpl:-single_nonzero_piece_cps(
            proc(c,v) if v::DomConstrain then DConstrain(conv_bool(c),op(v)) else p end if
            end proc,p,_testzero=(x->x=DSum()));
        if r :: Partition then DSplit(r) else r end if;
    end proc;

    local conv_bool := proc(r, $)
        if r :: {specfunc(`And`), `and`} then
            op(map(conv_bool,r))
        else
            r
        end if;
    end proc;
end module;

redundant_Partition_Pieces := module()
  uses Domain;

  export SimplName  := "Redundant Partition pieces";
  export SimplOrder := (10+1/2);

  export ModuleApply := proc(vs :: DomBound, sh :: DomShape, $)
    local as := Domain:-Bound:-toConstraints(vs, 'bound_types');
    subsindets(sh, DomSplit, proc(pr)
        local r; r := Partition:-Simpl(op(1,pr)) assuming op(as);
        if not r :: Partition then r else pr end if;
      end proc);
  end proc;
end module;


text:
[tracks]
1=This is on the stairs, Neurowaxx_-_This_is_on_the_stairs.ogg
2=Independence Day, DoKashiteru_-_Independence_Day.ogg
3=Revolve, freestylemix_-_hisboyelroy_-_Revolve.ogg

[options]
loop=false
shuffle=false

max_stars_repo_name: vardis/pano

max_issues_repo_name: vardis/pano

max_forks_repo_name: vardis/pano


text:
######################################################################
# proper_nonempty_subsets(A) is the set of all proper nonempty subsets of A

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

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

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

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

 n := nops(A);

 if n < 2 then return FAIL; fi;

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

 while nops(B) = 0 or nops(B) = n do
  B := select(a -> (r() = 1),A);
 od;

 return B;
end;

`list_elements/proper_nonempty_subsets` := (A::set) ->
 sort(map(sort,[op(combinat[powerset](A) minus {{},A})])):

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


text:
# 计算不变量阶数
$ifndef _INVORDER_
$define _INVORDER_

(*
    * 计算不变量的次数
*)
findInvariantsOrder:=proc(ee)
    if not type(denom(ee),numeric) then
        return findOrder(numer(ee))-findOrder(denom(ee));
    else
        return findOrder(ee);
    end if;
end proc:

(*
    * 计算齐次表达式的的次数（不含分母）
*)
findOrder:=proc(ee)
    local e;
    e:=expand(ee);
    if type(e,`+`) then
        e:=op(1,e);
    end if;
    e:=remove(type,e,numeric);
    if (e=NULL) then
        return 0;
    end if:
    return findItemOrder(e);
end proc:

(*
    * 计算一项的次数
*)
findItemOrder:=proc(ee)
    local s,_e;
    s:=0;
    for _e in ee do
        if type(_e,`^`) then
            s:=s+op(2,_e);
        else
            s:=s+1;
        end if;
    end do;
end proc:

$endif


text:
$ifndef _INVSOL_
$define _INVSOL_

InvSol:=module()
    option object;
    export
            # 公用变量
            As::static,     # 每个生成元的伴随矩阵
            A::static,      # 总的伴随矩阵
            nvars::static,  # 总的变量个数
            # 实例变量
            state,          # 状态代码，主要为人工干预提供入口
                            # 0：还未生成不变量方程
                            # 1：不变量方程求解失败
                            # 2：取特解失败
                            # 3：变换方程求解失败
                            # 4：等待选择最优代表元
                            # 5：求解完成
            oeq::static,    # 偏微分方程组，
                            # 迭代过程中oeq保持不变
                            # 通过变量代换来获取实际的oeq
            Deltas:=[],     # 不变量
            orders:=[],     # 不变量的阶数
            ieq:=[],        # 不变量方程程组，按不变量排序
            ieqCode:=0,     # 不变量方程的编号
            isols:=[],      # 不变量方程组的解
            icons:=[],      # 不变量方程组对应的条件
            isolInd:=1,     # 通解的下标
            useBranch:=false,# 使用分支求解则假设每个特解都是不等价的
            rsols:=[],      # 不变量方程的特解
            tsols:=[],      # TeqSol对象list
            tsolsList,      # TsolsList 对象
            rep,            # 代表元
            # 最优变换方程的解
            teq,
            tInd,
            tsol,
            tcon,

            vars,           # 需要求解的系数
            # 条件相关
            # 不变量方程+附加约束+展示约束，描述了代表元所代表的元素的范围
            # 后期只维护以下变量，而不修改icons,tcons等原始信息
            addcons:={},    # 附加约束，能够参与计算，其中，
                            # + 等式是对不变量方程的扩充，例如a[1]=0
                            # + 不等式式对取值范围的描述，并能够进行处理，例如a[1]>0
            discons:={},    # 展示约束，仅用于展示信息，不能参与计算，包括
                            # + a[1]*a[3]>0 这种不能处理的多变量约束
            # 导出函数
            getSubs::static,     # 在求解新的不变量时代入的条件
            updateVars::static,  # 更新求解变量
            getRealOeq::static,  # 获取代换后的oeq
            getZeroCons::static, # 获取为零约束
            addZeroCons::static, # 添加为零约束
            getIsolCons::static, # 返回通解的约束
            displayIeq::static,  # 显示不变量方程
            printSol::static,
            getDisplayIeq::static,
            getDisplayCons::static,
            uniqueKey::static,
            ModulePrint::static; # 显示函数

    getIsolCons:=proc(s::InvSol)
        return s:-addcons union s:-discons union s:-icons[s:-isolInd];
    end proc:

    ModulePrint:=proc(s::InvSol)
        if s:-state<5 then
            return sprintf("unsolved %d",s:-state);
        else
            return s:-rep;
        end if;
    end proc:

    uniqueKey:=proc(s::InvSol)
        return ModulePrint(s);
    end proc:

    printSol:=proc(s::InvSol)
        if s:-state<5 then
            printf("求解失败 %d",s:-state);
        else
            printf("--------------------------------\n");
            print(s:-rep);
            print(getDisplayCons(s));
            if s:-isols<>[] then
                print(s:-isols[s:-isolInd]);
                print(s:-tsol);
                print(s:-tcon);
            end if;
        end if;
    end proc:

    getDisplayIeq:=proc(s::InvSol)
        return [seq(Delta[i]=rhs(s:-ieq[i]),i=1..numelems(s:-ieq))];
    end proc:

    getDisplayCons:=proc(s::InvSol)
        local res:={};
        if s:-isols<>[] then
            res:=res union s:-icons[s:-isolInd];
        end if;
        res:=res union s:-addcons union s:-discons;
        res:=res minus convert(s:-ieq,set);
        return [getDisplayIeq(s)[],res[]];
    end proc:

    getSubs:=proc(s::InvSol)
        local r;
        r:=getZeroCons(s);
        r:=[r[]];
        if s:-isols<>[] then
            r:=[r[],s:-isols[s:-isolInd][]];
        end if;
        return remove(x->evalb(lhs(x)=rhs(x)),{r[]});
    end proc:

    updateVars:=proc(s::InvSol)
        s:-vars:=s:-vars minus indets(lhs~(getSubs(s)),name);
    end proc:

    getRealOeq:=proc(s::InvSol)
        local oeq:=s:-oeq;
        updateVars(s);
        oeq:=PDETools:-dsubs(phi(seq(a[i],i=1..s:-nvars))=phi(s:-vars[]),oeq);
        oeq:=eval(subs(getSubs(s)[],oeq)) minus {0};
        return oeq;
    end proc:

    getZeroCons:=proc(s::InvSol)
        return select(type,s:-addcons,equation);
    end proc:

    addZeroCons:=proc(s::InvSol,c::set)
        s:-addcons:=s:-addcons union c;
        s:-vars:=s:-vars minus indets(c,name);
    end proc:

    displayIeq:=proc(s::InvSol)
        local n;
        n:=numelems(s:-ieq);
        flog[1]([seq(Delta[i]=rhs(s:-ieq[i]),i=1..n)]);
        flog[1]([seq(Delta[i]=s:-Deltas[i],i=1..n)]);
    end proc:

end module:

$endif


text:
type $SS <class {@f1 i32, @f2 i8, @f3 i8, &fun(i32) i32}>
javaclass $TT <$SS> public static

func &bar (var %y i32) i32 {
  return (dread i32 %y)
}

func &foo (var %x <$SS>) i32 {
  var %ret i32
  polymorphiccallassigned <func(i32) i32> &bar (dread i32 %x 1) {dassign %ret 0}
  return (dread i32 %ret)
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
 # a function returning (i + j) * -998
 func $foo (
# var %i xxx,
  var %i i32, var %j i32, var %k i32,
  var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   assertge (dread i32 %i, constval i32 0)
   assertlt (dread i32 %i, dread i32 %k)
   return (
     ashr i32(
      mul i32 (
       add i32 (
        band i32 (dread i32 %i, dread i32 %j),
        bior i32 (
         bxor i32(dread i32 %i1, dread i32 %k1),
         eq i32 i32 (
          depositbits i32 1 23 (
           land i32 (dread i32 %j1, dread i32 %k),
           lior i32 (
             min i32 (
              shl i32 (dread i32 %k1, dread i32 %j),
              sub i32 (dread i32 %j1, dread i32 %i)),
             rem i32 (dread i32 %i, dread i32 %k))),
          div i32(
           lshr i32 ( dread i32 %k, dread i32 %i),
           max i32 ( dread i32 %i, dread i32 %j))))),
       constval i32 -998),
      dread i32 %i))}

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.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: work_gga_x *)
(* prefix:
  gga_x_pbeint_params *params;

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

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

(* this is the gga_x_pbe expression *)
f0 := s -> 1 + params_a_kappa * (1 - params_a_kappa/(params_a_kappa + mu(s)*s^2)):
f  := x -> f0(X2S * x):


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

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

(* type: work_gga_x *)

mbeta  :=  0.0042:
mgamma :=  6:
c0    :=  1.09878:
c1    := -2.51173:
c2    :=  0.0156233:

f_aux := x -> 1 + mgamma*mbeta*x*arcsinh(x):

f := x -> c0 + mbeta/X_FACTOR_C*x^2*(c1/f_aux(x) + c2/(mbeta*f_aux(x)^2)):


text:
program ifst;
var ch : char;
begin
    readln(ch);
    if ch = 'a' then
        writeln('It is ''a'' ');