text:
# Implizite kinematische Zwangsbedingungen für 2D-Picker

#
# Erklärung des Dateinamens
# picker2D: 2D-Roboter für Pick&Place
# m2: Modellierung der geschlossenen Ketten als Parallelogramm
# IC: Modellierung des Roboters mit impliziten Zwangsbedingungen
# kinematic_contraints_implizit: Berechnung der kinematischen Zwangsbedingungen (implizite Form)
#
# Quelle
# Aufzeichnungen Schappler, 6.5.2019
# Ursprüngliche Notizen von Abderahman Bejaoui (Hiwi bei Moritz Schappler)
#
# Autor
# Abderahman Bejaoui
# Hiwi bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2019-03
# (C) Institut fuer mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../transformation/proc_trotz":
read "../transformation/proc_transl":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env_IC":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name_OL):
# Ergebnisse der Kinematik laden
read sprintf("../codeexport/%s/tmp/kinematics_floatb_%s_rotmat_maple.m", robot_name_OL, base_method_name);
read "../robot_codegen_definitions/robot_env_IC":
Trf := Trf:
Trf_c := Trf_c:
# VGK Gelb 0-7-13-10-4-2-1-0
# Schleife (0-7)-(7-13)
T_1_13 := combine( Matrix(Trf(1..4,1..4, 7)) . Matrix(Trf(1..4,1..4,13))):
# Schleife (0-1)(1-2)-(2-4)-(4-10)
T_1_10:= combine( Matrix(Trf(1..4,1..4, 1)) . Matrix(Trf(1..4,1..4,2)) . Matrix(Trf(1..4,1..4,4)) . Matrix(Trf(1..4,1..4,10))):
h1t := T_1_13(1..3,4) - T_1_10(1..3,4):

#tmp := Transpose( Matrix(T_1_13(1..3,1..3)) ) . Matrix(T_1_10(1..3,1..3));  nur anzeigen lassen für h1r
;
#combine(tmp);# nur anzeigen lassen für h1r
;
h1r := (-qJ7(t)+qJ1(t)+qJ2(t)+qJ4(t)+qJ10(t))+Pi/2:
# VGK PINK 0-5-14-11-8-1-0
# Schleife (0-5)-(5-14)
T_2_14 := combine(Matrix(Trf(1..4,1..4, 5)) . Matrix(Trf(1..4,1..4,14))):
# Schleife (0-1)-(1-8)-(8-11)
T_2_11:= combine(Matrix(Trf(1..4,1..4, 1)) . Matrix(Trf(1..4,1..4,8)).Matrix(Trf(1..4,1..4, 11))):
h2t := T_2_11(1..3,4) - T_2_14(1..3,4):
#tmp := Transpose( Matrix(T_2_11(1..3,1..3)) ) . Matrix(T_2_14(1..3,1..3));  # nur anzeigen lassen für h1r
;
 #combine(tmp); # nur anzeigen lassen für h1r
;
h2r:=(-qJ1(t)-qJ8(t)-qJ11(t)+phi05+qJ5(t)):
# VGK GRÜN 2-6-12-15-9-3
# Schleife (2-6)-(6-12)
T_2_12 := combine( Matrix(Trf(1..4,1..4, 6)) . Matrix(Trf(1..4,1..4,12))):
# Schleife (2-3)-(3-9)-(9-15)
T_2_15:= combine( Matrix(Trf(1..4,1..4, 3)) . Matrix(Trf(1..4,1..4,9)).Matrix(Trf(1..4,1..4,15)) ):
h3t := T_2_12(1..3,4) - T_2_15(1..3,4):
#tmp := Transpose( Matrix(T_2_12(1..3,1..3)) ) . Matrix(T_2_15(1..3,1..3));  nur anzeigen lassen für h2r
;
#combine(tmp);  nur anzeigen lassen für h3r
;
h3r:=(qJ6(t)+qJ12(t)-qJ3(t)-qJ9(t)):
# Zusätzliche Zwangsbedingung (6 und 8 sind gleiche Körper)
# TODO: In einer neuen Modellierung sollten beide Körper zusammengefasst werden.
h4:=Pi/2 - qJ8(t) - qJ9(t) + qJ2(t):
# Zusammenstellen aller Zwangsbedingungen
implconstr_t := <h1t([1, 2],1);h2t([1, 2],1);h3t([1, 2],1); h1r; h2r; h3r; h4>; # TODO: In h1r, h2r, h3r muss das richtige drinstehen.
implconstr_s := convert_t_s(implconstr_t):
# Exportiere Code für folgende Skripte
kin_constraints_exist:=true:
save implconstr_t, implconstr_s, kin_constraints_exist, sprintf("../codeexport/%s/tmp/kinematic_constraints_implicit_maple.m", robot_name):
# Exportieren des vollständigen Ausdruckes
if codegen_act then
  MatlabExport(implconstr_s, sprintf("../codeexport/%s/tmp/kinconstr_impl_matlab.m", robot_name), 2):
end if:
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( implconstr_s )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_implicit_constraints_symbols_list_matlab.m", robot_name), 2):


text:
# Convert Domain to a KB. Very partial, mainly for backwards compatibility
# for parts of the code which still work with KB.
export Bound := module ()
    export toKB := proc(dom :: DomBound, kb0 :: t_kb, $)::[t_kb_mb, list(name=name)];
        local kb := kb0, vs := op(1, dom), rn := []
            , vn, vt, make, lo, hi, vn_rn, rn_t, v ;
        for v in vs do
            vn, vt, make := op(v);
            lo, hi := ExtBound[make]:-SplitRange(vt);
            vn_rn, kb := ExtBound[make]:-MakeKB( vn, lo, hi, kb );
            rn := [ vn=vn_rn, op(rn) ];
        end do;
        [ kb, rn ]
    end proc;

    export varsOf := proc(dom :: DomBound, output::identical("list","set","seq"):="list")
      local v := map(x->op(1,x), op(1, dom));
      op(table(["list"=[v],"set"=[{op(v)}],"seq"=v])[output]);
    end proc;

    export withVarsIxs := proc(dom :: DomBound, $) :: DomBound;
        DBound( op(1..2,dom), table([ seq(op([1,i,1],dom)=i,i=1..nops(op(1,dom))) ]) );
    end proc;

    export varIx := proc(dom0::DomBound, v::DomBoundVar, $)::nonnegint;
      local i := varIx_mb(dom0,v);
      if i = -1 then
        error "cannot find var %1 in %2", v, dom0;
      else i end if;
    end proc;

    export varIx_mb := proc(dom0::DomBound, v::DomBoundVar, $)::{nonnegint,identical(-1)};
      local dom := `if`(nops(dom0)>=3, (_->_), withVarsIxs)(dom0);
      if assigned(op(3,dom)[v]) then
        op(3,dom)[v]
      else
        -1
      end if;
    end proc;

    export isEmpty := proc(dom :: DomBound, $)::truefalse;
        evalb(nops(op(1,dom))=0);
    end proc;

    export get := proc(dom :: DomBound, var :: DomBoundVar, $)
        local th;
        th := select(x->op(1,x)=var, op(1,dom));
        if nops(th) = 1 then
            op([1,2], th), op([1,3], th)
        elif nops(th) = 0 then
            error "cannot find var %1 in %2", var, dom;
        else
            error "multiple references in %1", dom;
        end if;
    end proc;

    export upd := proc(dom0 :: DomBound, vr :: DomBoundVar, ty :: DomBoundRange, $)
        local dom := `if`(nops(dom0)>=3, (_->_), withVarsIxs)(dom0), i
            , old_lo,old_hi, lo,hi, old_ty, new_ty, _, vr_k;
        i := varIx(dom, vr);

        _, old_ty, vr_k := op(op([1,i], dom));
        old_lo,old_hi := Domain:-ExtBound[vr_k]:-SplitRange(old_ty);
        lo, hi        := Domain:-ExtBound[vr_k]:-SplitRange(ty);

        if lo :: realcons then
          lo := Domain:-ExtBound[vr_k]:-Max(old_lo, lo);
        end if;
        if hi :: realcons then
          hi := Domain:-ExtBound[vr_k]:-Min(old_hi, hi);
        end if;

        new_ty := Domain:-ExtBound[vr_k]:-MakeRange(lo,hi);

        subsop([1,i,2]=new_ty, dom);
    end proc;

    local constrain := proc( opts, vn::DomBoundVar, ty::DomBoundRange, mk :: DomBoundKind, $ )
                    :: set({relation, `::`});
        local lo, hi, noi, bt
            , noinf := evalb('no_infinity' in opts)
            , btys  := evalb('bound_types' in opts);
        lo,hi := Domain:-ExtBound[mk]:-SplitRange(ty);
        bt  := `if`(btys, { vn :: Domain:-ExtBound[mk]:-BoundType }, {});
        bt := { Domain:-ExtBound[mk]:-Constrain( lo, vn )
        , Domain:-ExtBound[mk]:-Constrain( vn, hi )
        , bt[]};
        `if`(noinf,x->remove(c->c::relation and ormap(s->s(c)::identical(infinity,-infinity),[lhs,rhs]), x),_->_)
            (bt);
    end proc;

    local toConstraints_opts := { 'no_infinity', 'bound_types' };
    export toConstraints := proc( bnd :: DomBound )
        local cs, opts, bad_opts, noinf;
        opts := { args[2..-1] } ;
        bad_opts := opts minus toConstraints_opts;
        if bad_opts <> {} then
            error "invalid arguments: %1", bad_opts;
        end if;
        {op(map(b->constrain(opts, op(b))[], op(1,bnd)))
        ,op( `if`(nops(bnd)>1,op(2,bnd),{}) ) };
    end proc;
end module;#Bound

ext: mpl

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


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

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

(* type: work_lda *)

$define lda_x_params
$include "lda_x.mpl"

(* attenuation function *)
aux1_erf := a -> sqrt(Pi)*erf(1/(2*a)):
aux2_erf := a -> exp(-1/(4*a^2)) - 1:
aux3_erf := a -> 2*a^2*aux2_erf(a) + 1/2:
attenuation_erf := a -> 1 - 8/3*a*(aux1_erf(a) + 2*a*(aux2_erf(a) - aux3_erf(a))):

a_cnst := (4/(9*Pi))^(1/3)*p_a_cam_omega/2:

lda_x_erf_s := (rs, z) -> convert(piecewise(z = -1, 0,
  (1 + z)^(4/3)*attenuation_erf(a_cnst*rs/(1 + z)^(1/3))
), 'Heaviside'):

f_lda_x_erf := (rs, z) -> lda_x_ax*(
   lda_x_erf_s(rs, z) + lda_x_erf_s(rs, -z)
)/rs:

f := (rs, z) -> f_lda_x_erf(rs, z):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


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

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

rational_p_f0 := s -> (1 + params_a_C2/params_a_p * s^2)^(-params_a_p):
rational_p_f := x -> rational_p_f0(X2S*x):

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


text:
function [x,f,exitflag,output] = minFunc(funObj,x0,options,varargin)
% minFunc(funObj,x0,options,varargin)
%
% Unconstrained optimizer using a line search strategy
%
% Uses an interface very similar to fminunc
%   (it doesn't support all of the optimization toolbox options,
%       but supports many other options).
%
% It computes descent directions using one of ('Method'):
%   - 'sd': Steepest Descent
%       (no previous information used, not recommended)
%   - 'csd': Cyclic Steepest Descent
%       (uses previous step length for a fixed length cycle)
%   - 'bb': Barzilai and Borwein Gradient
%       (uses only previous step)
%   - 'cg': Non-Linear Conjugate Gradient
%       (uses only previous step and a vector beta)
%   - 'scg': Scaled Non-Linear Conjugate Gradient
%       (uses previous step and a vector beta,
%           and Hessian-vector products to initialize line search)
%   - 'pcg': Preconditionined Non-Linear Conjugate Gradient
%       (uses only previous step and a vector beta, preconditioned version)
%   - 'lbfgs': Quasi-Newton with Limited-Memory BFGS Updating
%       (default: uses a predetermined nunber of previous steps to form a
%           low-rank Hessian approximation)
%   - 'newton0': Hessian-Free Newton
%       (numerically computes Hessian-Vector products)
%   - 'pnewton0': Preconditioned Hessian-Free Newton
%       (numerically computes Hessian-Vector products, preconditioned
%       version)
%   - 'qnewton': Quasi-Newton Hessian approximation
%       (uses dense Hessian approximation)
%   - 'mnewton': Newton's method with Hessian calculation after every
%   user-specified number of iterations
%       (needs user-supplied Hessian matrix)
%   - 'newton': Newton's method with Hessian calculation every iteration
%       (needs user-supplied Hessian matrix)
%   - 'tensor': Tensor
%       (needs user-supplied Hessian matrix and Tensor of 3rd partial derivatives)
%
% Several line search strategies are available for finding a step length satisfying
%   the termination criteria ('LS'):
%   - 0: Backtrack w/ Step Size Halving
%   - 1: Backtrack w/ Quadratic/Cubic Interpolation from new function values
%   - 2: Backtrack w/ Cubic Interpolation from new function + gradient
%   values (default for 'bb' and 'sd')
%   - 3: Bracketing w/ Step Size Doubling and Bisection
%   - 4: Bracketing w/ Cubic Interpolation/Extrapolation with function +
%   gradient values (default for all except 'bb' and 'sd')
%   - 5: Bracketing w/ Mixed Quadratic/Cubic Interpolation/Extrapolation
%   - 6: Use Matlab Optimization Toolbox's line search
%           (requires Matlab's linesearch.m to be added to the path)
%
%   Above, the first three find a point satisfying the Armijo conditions,
%   while the last four search for find a point satisfying the Wolfe
%   conditions.  If the objective function overflows, it is recommended
%   to use one of the first 3.
%   The first three can be used to perform a non-monotone
%   linesearch by changing the option 'Fref'.
%
% Several strategies for choosing the initial step size are avaiable ('LS_init'):
%   - 0: Always try an initial step length of 1 (default for all except 'cg' and 'sd')
%       (t = 1)
%   - 1: Use a step similar to the previous step (default for 'cg' and 'sd')
%       (t = t_old*min(2,g'd/g_old'd_old))
%   - 2: Quadratic Initialization using previous function value and new
%   function value/gradient (use this if steps tend to be very long)
%       (t = min(1,2*(f-f_old)/g))
%   - 3: The minimum between 1 and twice the previous step length
%       (t = min(1,2*t)
%   - 4: The scaled conjugate gradient step length (may accelerate
%   conjugate gradient methods, but requires a Hessian-vector product)
%       (t = g'd/d'Hd)
%
% Inputs:
%   funObj is a function handle
%   x0 is a starting vector;
%   options is a struct containing parameters
%  (defaults are used for non-existent or blank fields)
%   all other arguments are passed to funObj
%
% Outputs:
%   x is the minimum value found
%   f is the function value at the minimum found
%   exitflag returns an exit condition
%   output returns a structure with other information
%
% Supported Input Options
%   Display - Level of display [ off | final | (iter) | full | excessive ]
%   MaxFunEvals - Maximum number of function evaluations allowed (1000)
%   MaxIter - Maximum number of iterations allowed (500)
%   TolFun - Termination tolerance on the first-order optimality (1e-5)
%   TolX - Termination tolerance on progress in terms of function/parameter changes (1e-9)
%   Method - [ sd | csd | bb | cg | scg | pcg | {lbfgs} | newton0 | pnewton0 |
%       qnewton | mnewton | newton | tensor ]
%   c1 - Sufficient Decrease for Armijo condition (1e-4)
%   c2 - Curvature Decrease for Wolfe conditions (.2 for cg methods, .9 otherwise)
%   LS_init - Line Search Initialization -see above (2 for cg/sd, 4 for scg, 0 otherwise)
%   LS - Line Search type -see above (2 for bb, 4 otherwise)
%   Fref - Setting this to a positive integer greater than 1
%       will use non-monotone Armijo objective in the line search.
%       (20 for bb, 10 for csd, 1 for all others)
%   numDiff - compute derivative numerically
%       (default: 0) (this option has a different effect for 'newton', see below)
%   useComplex - if 1, use complex differentials when computing numerical derivatives
%       to get very accurate values (default: 0)
%   DerivativeCheck - if 'on', computes derivatives numerically at initial
%       point and compares to user-supplied derivative (default: 'off')
%   outputFcn - function to run after each iteration (default: []).  It
%       should have the following interface:
%       outputFcn(x,infoStruct,state,varargin{:})
%   useMex - where applicable, use mex files to speed things up (default: 1)
%
% Method-specific input options:
%   newton:
%       HessianModify - type of Hessian modification for direct solvers to
%       use if the Hessian is not positive definite (default: 0)
%           0: Minimum Euclidean norm s.t. eigenvalues sufficiently large
%           (requires eigenvalues on iterations where matrix is not pd)
%           1: Start with (1/2)*||A||_F and increment until Cholesky succeeds
%           (an approximation to method 0, does not require eigenvalues)
%           2: Modified LDL factorization
%           (only 1 generalized Cholesky factorization done and no eigenvalues required)
%           3: Modified Spectral Decomposition
%           (requires eigenvalues)
%           4: Modified Symmetric Indefinite Factorization
%           5: Uses the eigenvector of the smallest eigenvalue as negative
%           curvature direction
%       cgSolve - use conjugate gradient instead of direct solver (default: 0)
%           0: Direct Solver
%           1: Conjugate Gradient
%           2: Conjugate Gradient with Diagonal Preconditioner
%           3: Conjugate Gradient with LBFGS Preconditioner
%           x: Conjugate Graident with Symmetric Successive Over Relaxation
%           Preconditioner with parameter x
%               (where x is a real number in the range [0,2])
%           x: Conjugate Gradient with Incomplete Cholesky Preconditioner
%           with drop tolerance -x
%               (where x is a real negative number)
%       numDiff - compute Hessian numerically
%                 (default: 0, done with complex differentials if useComplex = 1)
%       LS_saveHessiancomp - when on, only computes the Hessian at the
%       first and last iteration of the line search (default: 1)
%   mnewton:
%       HessianIter - number of iterations to use same Hessian (default: 5)
%   qnewton:
%       initialHessType - scale initial Hessian approximation (default: 1)
%       qnUpdate - type of quasi-Newton update (default: 3):
%           0: BFGS
%           1: SR1 (when it is positive-definite, otherwise BFGS)
%           2: Hoshino
%           3: Self-Scaling BFGS
%           4: Oren's Self-Scaling Variable Metric method
%           5: McCormick-Huang asymmetric update
%       Damped - use damped BFGS update (default: 1)
%   newton0/pnewton0:
%       HvFunc - user-supplied function that returns Hessian-vector products
%           (by default, these are computed numerically using autoHv)
%           HvFunc should have the following interface: HvFunc(v,x,varargin{:})
%       useComplex - use a complex perturbation to get high accuracy
%           Hessian-vector products (default: 0)
%           (the increased accuracy can make the method much more efficient,
%               but gradient code must properly support complex inputs)
%       useNegCurv - a negative curvature direction is used as the descent
%           direction if one is encountered during the cg iterations
%           (default: 1)
%       precFunc (for pnewton0 only) - user-supplied preconditioner
%           (by default, an L-BFGS preconditioner is used)
%           precFunc should have the following interfact:
%           precFunc(v,x,varargin{:})
%   lbfgs:
%       Corr - number of corrections to store in memory (default: 100)
%           (higher numbers converge faster but use more memory)
%       Damped - use damped update (default: 0)
%   pcg:
%       cgUpdate - type of update (default: 2)
%   cg/scg/pcg:
%       cgUpdate - type of update (default for cg/scg: 2, default for pcg: 1)
%           0: Fletcher Reeves
%           1: Polak-Ribiere
%           2: Hestenes-Stiefel (not supported for pcg)
%           3: Gilbert-Nocedal
%       HvFunc (for scg only)- user-supplied function that returns Hessian-vector
%           products
%           (by default, these are computed numerically using autoHv)
%           HvFunc should have the following interface:
%           HvFunc(v,x,varargin{:})
%       precFunc (for pcg only) - user-supplied preconditioner
%           (by default, an L-BFGS preconditioner is used)
%           precFunc should have the following interfact:
%           precFunc(v,x,varargin{:})
%   bb:
%       bbType - type of bb step (default: 1)
%           0: min_alpha ||delta_x - alpha delta_g||_2
%           1: min_alpha ||alpha delta_x - delta_g||_2
%           2: Conic BB
%           3: Gradient method with retards
%   csd:
%       cycle - length of cycle (default: 3)
%
% Supported Output Options
%   iterations - number of iterations taken
%   funcCount - number of function evaluations
%   algorithm - algorithm used
%   firstorderopt - first-order optimality
%   message - exit message
%   trace.funccount - function evaluations after each iteration
%   trace.fval - function value after each iteration
%
% Author: Mark Schmidt (2006)
% Web: http://www.cs.ubc.ca/~schmidtm
%
% Sources (in order of how much the source material contributes):
%   J. Nocedal and S.J. Wright.  1999.  "Numerical Optimization".  Springer Verlag.
%   R. Fletcher.  1987.  "Practical Methods of Optimization".  Wiley.
%   J. Demmel.  1997.  "Applied Linear Algebra.  SIAM.
%   R. Barret, M. Berry, T. Chan, J. Demmel, J. Dongarra, V. Eijkhout, R.
%   Pozo, C. Romine, and H. Van der Vost.  1994.  "Templates for the Solution of
%   Linear Systems: Building Blocks for Iterative Methods".  SIAM.
%   J. More and D. Thuente.  "Line search algorithms with guaranteed
%   sufficient decrease".  ACM Trans. Math. Softw. vol 20, 286-307, 1994.
%   M. Raydan.  "The Barzilai and Borwein gradient method for the large
%   scale unconstrained minimization problem".  SIAM J. Optim., 7, 26-33,
%   (1997).
%   "Mathematical Optimization".  The Computational Science Education
%   Project.  1995.
%   C. Kelley.  1999.  "Iterative Methods for Optimization".  Frontiers in
%   Applied Mathematics.  SIAM.

if nargin < 3
    options = [];
end

% Get Parameters
[verbose,verboseI,debug,doPlot,maxFunEvals,maxIter,tolFun,tolX,method,...
    corrections,c1,c2,LS_init,LS,cgSolve,qnUpdate,cgUpdate,initialHessType,...
    HessianModify,Fref,useComplex,numDiff,LS_saveHessianComp,...
    DerivativeCheck,Damped,HvFunc,bbType,cycle,...
    HessianIter,outputFcn,useMex,useNegCurv,precFunc] = ...
    minFunc_processInputOptions(options);

if isfield(options, 'logfile')
    logfile = options.logfile;
else
    logfile = [];
end

% Constants
SD = 0;
CSD = 1;
BB = 2;
CG = 3;
PCG = 4;
LBFGS = 5;
QNEWTON = 6;
NEWTON0 = 7;
NEWTON = 8;
TENSOR = 9;

% Initialize
p = length(x0);
d = zeros(p,1);
x = x0;
t = 1;

% If necessary, form numerical differentiation functions
funEvalMultiplier = 1;
if numDiff && method ~= TENSOR
    varargin(3:end+2) = varargin(1:end);
    varargin{1} = useComplex;
    varargin{2} = funObj;
    if method ~= NEWTON
        if debug
            if useComplex
                fprintf('Using complex differentials for gradient computation\n');
            else
                fprintf('Using finite differences for gradient computation\n');
            end
        end
        funObj = @autoGrad;
    else
        if debug
            if useComplex
                fprintf('Using complex differentials for gradient computation\n');
            else
                fprintf('Using finite differences for gradient computation\n');
            end
        end
        funObj = @autoHess;
    end

    if method == NEWTON0 && useComplex == 1
        if debug
            fprintf('Turning off the use of complex differentials\n');
        end
        useComplex = 0;
    end

    if useComplex
        funEvalMultiplier = p;
    else
        funEvalMultiplier = p+1;
    end
end

% Evaluate Initial Point
if method < NEWTON
    [f,g] = feval(funObj, x, varargin{:});
else
    [f,g,H] = feval(funObj, x, varargin{:});
    computeHessian = 1;
end
funEvals = 1;

if strcmp(DerivativeCheck,'on')
    if numDiff
        fprintf('Can not do derivative checking when numDiff is 1\n');
    end
    % Check provided gradient/hessian function using numerical derivatives
    fprintf('Checking Gradient:\n');
    [f2,g2] = autoGrad(x,useComplex,funObj,varargin{:});

    fprintf('Max difference between user and numerical gradient: %f\n',max(abs(g-g2)));
    if max(abs(g-g2)) > 1e-4
        fprintf('User NumDif:\n');
        [g g2]
        diff = abs(g-g2)
        pause;
    end

    if method >= NEWTON
        fprintf('Check Hessian:\n');
        [f2,g2,H2] = autoHess(x,useComplex,funObj,varargin{:});

        fprintf('Max difference between user and numerical hessian: %f\n',max(abs(H(:)-H2(:))));
        if max(abs(H(:)-H2(:))) > 1e-4
            H
            H2
            diff = abs(H-H2)
            pause;
        end
    end
end

% Output Log
if verboseI
    fprintf('%10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');
end

if logfile
    fid = fopen(logfile, 'a');
    if (fid > 0)
        fprintf(fid, '-- %10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');
        fclose(fid);
    end
end

% Output Function
if ~isempty(outputFcn)
    callOutput(outputFcn,x,'init',0,funEvals,f,[],[],g,[],sum(abs(g)),varargin{:});
end

% Initialize Trace
trace.fval = f;
trace.funcCount = funEvals;

% Check optimality of initial point
if sum(abs(g)) <= tolFun
    exitflag=1;
    msg = 'Optimality Condition below TolFun';
    if verbose
        fprintf('%s\n',msg);
    end
    if nargout > 3
        output = struct('iterations',0,'funcCount',1,...
            'algorithm',method,'firstorderopt',sum(abs(g)),'message',msg,'trace',trace);
    end
    return;
end

% Perform up to a maximum of 'maxIter' descent steps:
for i = 1:maxIter

    % ****************** COMPUTE DESCENT DIRECTION *****************

    switch method
        case SD % Steepest Descent
            d = -g;

        case CSD % Cyclic Steepest Descent

            if mod(i,cycle) == 1 % Use Steepest Descent
                alpha = 1;
                LS_init = 2;
                LS = 4; % Precise Line Search
            elseif mod(i,cycle) == mod(1+1,cycle) % Use Previous Step
                alpha = t;
                LS_init = 0;
                LS = 2; % Non-monotonic line search
            end
            d = -alpha*g;

        case BB % Steepest Descent with Barzilai and Borwein Step Length

            if i == 1
                d = -g;
            else
                y = g-g_old;
                s = t*d;
                if bbType == 0
                    yy = y'*y;
                    alpha = (s'*y)/(yy);
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                elseif bbType == 1
                    sy = s'*y;
                    alpha = (s'*s)/sy;
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                elseif bbType == 2 % Conic Interpolation ('Modified BB')
                    sy = s'*y;
                    ss = s'*s;
                    alpha = ss/sy;
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                    alphaConic = ss/(6*(myF_old - f) + 4*g'*s + 2*g_old'*s);
                    if alphaConic > .001*alpha && alphaConic < 1000*alpha
                        alpha = alphaConic;
                    end
                elseif bbType == 3 % Gradient Method with retards (bb type 1, random selection of previous step)
                    sy = s'*y;
                    alpha = (s'*s)/sy;
                    if alpha <= 1e-10 || alpha > 1e10
                        alpha = 1;
                    end
                    v(1+mod(i-2,5)) = alpha;
                    alpha = v(ceil(rand*length(v)));
                end
                d = -alpha*g;
            end
            g_old = g;
            myF_old = f;

        case CG % Non-Linear Conjugate Gradient

            if i == 1
                d = -g; % Initially use steepest descent direction
            else
                gtgo = g'*g_old;
                gotgo = g_old'*g_old;

                if cgUpdate == 0
                    % Fletcher-Reeves
                    beta = (g'*g)/(gotgo);
                elseif cgUpdate == 1
                    % Polak-Ribiere
                    beta = (g'*(g-g_old)) /(gotgo);
                elseif cgUpdate == 2
                    % Hestenes-Stiefel
                    beta = (g'*(g-g_old))/((g-g_old)'*d);
                else
                    % Gilbert-Nocedal
                    beta_FR = (g'*(g-g_old)) /(gotgo);
                    beta_PR = (g'*g-gtgo)/(gotgo);
                    beta = max(-beta_FR,min(beta_PR,beta_FR));
                end

                d = -g + beta*d;

                % Restart if not a direction of sufficient descent
                if g'*d > -tolX
                    if debug
                        fprintf('Restarting CG\n');
                    end
                    beta = 0;
                    d = -g;
                end

                % Old restart rule:
                %if beta < 0 || abs(gtgo)/(gotgo) >= 0.1 || g'*d >= 0

            end
            g_old = g;

        case PCG % Preconditioned Non-Linear Conjugate Gradient

            % Apply preconditioner to negative gradient
            if isempty(precFunc)
                % Use L-BFGS Preconditioner
                if i == 1
                    old_dirs = zeros(length(g),0);
                    old_stps = zeros(length(g),0);
                    Hdiag = 1;
                    s = -g;
                else
                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);

                    if useMex
                        s = lbfgsC(-g,old_dirs,old_stps,Hdiag);
                    else
                        s = lbfgs(-g,old_dirs,old_stps,Hdiag);
                    end
                end
            else % User-supplied preconditioner
                s = precFunc(-g,x,varargin{:});
            end

            if i == 1
                d = s;
            else

                if cgUpdate == 0
                    % Preconditioned Fletcher-Reeves
                    beta = (g'*s)/(g_old'*s_old);
                elseif cgUpdate < 3
                    % Preconditioned Polak-Ribiere
                    beta = (g'*(s-s_old))/(g_old'*s_old);
                else
                    % Preconditioned Gilbert-Nocedal
                    beta_FR = (g'*s)/(g_old'*s_old);
                    beta_PR = (g'*(s-s_old))/(g_old'*s_old);
                    beta = max(-beta_FR,min(beta_PR,beta_FR));
                end
                d = s + beta*d;

                if g'*d > -tolX
                    if debug
                        fprintf('Restarting CG\n');
                    end
                    beta = 0;
                    d = s;
                end

            end
            g_old = g;
            s_old = s;
        case LBFGS % L-BFGS

            % Update the direction and step sizes

            if i == 1
                d = -g; % Initially use steepest descent direction
                old_dirs = zeros(length(g),0);
                old_stps = zeros(length(d),0);
                Hdiag = 1;
            else
                if Damped
                    [old_dirs,old_stps,Hdiag] = dampedUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                else
                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                end

                if useMex
                    d = lbfgsC(-g,old_dirs,old_stps,Hdiag);
                else
                    d = lbfgs(-g,old_dirs,old_stps,Hdiag);
                end
            end
            g_old = g;

        case QNEWTON % Use quasi-Newton Hessian approximation

            if i == 1
                d = -g;
            else
                % Compute difference vectors
                y = g-g_old;
                s = t*d;

                if i == 2
                    % Make initial Hessian approximation
                    if initialHessType == 0
                        % Identity
                        if qnUpdate <= 1
                            R = eye(length(g));
                        else
                            H = eye(length(g));
                        end
                    else
                        % Scaled Identity
                        if debug
                            fprintf('Scaling Initial Hessian Approximation\n');
                        end
                        if qnUpdate <= 1
                            % Use Cholesky of Hessian approximation
                            R = sqrt((y'*y)/(y'*s))*eye(length(g));
                        else
                            % Use Inverse of Hessian approximation
                            H = eye(length(g))*(y'*s)/(y'*y);
                        end
                    end
                end

                if qnUpdate == 0 % Use BFGS updates
                    Bs = R'*(R*s);
                    if Damped
                        eta = .02;
                        if y'*s < eta*s'*Bs
                            if debug
                                fprintf('Damped Update\n');
                            end
                            theta = min(max(0,((1-eta)*s'*Bs)/(s'*Bs - y'*s)),1);
                            y = theta*y + (1-theta)*Bs;
                        end
                        R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                    else
                        if y'*s > 1e-10
                            R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                        else
                            if debug
                                fprintf('Skipping Update\n');
                            end
                        end
                    end
                elseif qnUpdate == 1 % Perform SR1 Update if it maintains positive-definiteness

                    Bs = R'*(R*s);
                    ymBs = y-Bs;
                    if abs(s'*ymBs) >= norm(s)*norm(ymBs)*1e-8 && (s-((R\(R'\y))))'*y > 1e-10
                        R = cholupdate(R,-ymBs/sqrt(ymBs'*s),'-');
                    else
                        if debug
                            fprintf('SR1 not positive-definite, doing BFGS Update\n');
                        end
                        if Damped
                            eta = .02;
                            if y'*s < eta*s'*Bs
                                if debug
                                    fprintf('Damped Update\n');
                                end
                                theta = min(max(0,((1-eta)*s'*Bs)/(s'*Bs - y'*s)),1);
                                y = theta*y + (1-theta)*Bs;
                            end
                            R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                        else
                            if y'*s > 1e-10
                                R = cholupdate(cholupdate(R,y/sqrt(y'*s)),Bs/sqrt(s'*Bs),'-');
                            else
                                if debug
                                    fprintf('Skipping Update\n');
                                end
                            end
                        end
                    end
                elseif qnUpdate == 2 % Use Hoshino update
                    v = sqrt(y'*H*y)*(s/(s'*y) - (H*y)/(y'*H*y));
                    phi = 1/(1 + (y'*H*y)/(s'*y));
                    H = H + (s*s')/(s'*y) - (H*y*y'*H)/(y'*H*y) + phi*v*v';

                elseif qnUpdate == 3 % Self-Scaling BFGS update
                    ys = y'*s;
                    Hy = H*y;
                    yHy = y'*Hy;
                    gamma = ys/yHy;
                    v = sqrt(yHy)*(s/ys - Hy/yHy);
                    H = gamma*(H - Hy*Hy'/yHy + v*v') + (s*s')/ys;
                elseif qnUpdate == 4 % Oren's Self-Scaling Variable Metric update

                    % Oren's method
                    if (s'*y)/(y'*H*y) > 1
                        phi = 1; % BFGS
                        omega = 0;
                    elseif (s'*(H\s))/(s'*y) < 1
                        phi = 0; % DFP
                        omega = 1;
                    else
                        phi = (s'*y)*(y'*H*y-s'*y)/((s'*(H\s))*(y'*H*y)-(s'*y)^2);
                        omega = phi;
                    end

                    gamma = (1-omega)*(s'*y)/(y'*H*y) + omega*(s'*(H\s))/(s'*y);
                    v = sqrt(y'*H*y)*(s/(s'*y) - (H*y)/(y'*H*y));
                    H = gamma*(H - (H*y*y'*H)/(y'*H*y) + phi*v*v') + (s*s')/(s'*y);

                elseif qnUpdate == 5 % McCormick-Huang asymmetric update
                    theta = 1;
                    phi = 0;
                    psi = 1;
                    omega = 0;
                    t1 = s*(theta*s + phi*H'*y)';
                    t2 = (theta*s + phi*H'*y)'*y;
                    t3 = H*y*(psi*s + omega*H'*y)';
                    t4 = (psi*s + omega*H'*y)'*y;
                    H = H + t1/t2 - t3/t4;
                end

                if qnUpdate <= 1
                    d = -R\(R'\g);
                else
                    d = -H*g;
                end

            end
            g_old = g;

        case NEWTON0 % Hessian-Free Newton

            cgMaxIter = min(p,maxFunEvals-funEvals);
            cgForce = min(0.5,sqrt(norm(g)))*norm(g);

            % Set-up preconditioner
            precondFunc = [];
            precondArgs = [];
            if cgSolve == 1
                if isempty(precFunc) % Apply L-BFGS preconditioner
                    if i == 1
                        old_dirs = zeros(length(g),0);
                        old_stps = zeros(length(g),0);
                        Hdiag = 1;
                    else
                        [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                        if useMex
                            precondFunc = @lbfgsC;
                        else
                            precondFunc = @lbfgs;
                        end
                        precondArgs = {old_dirs,old_stps,Hdiag};
                    end
                    g_old = g;
                else
                    % Apply user-defined preconditioner
                    precondFunc = precFunc;
                    precondArgs = {x,varargin{:}};
                end
            end

            % Solve Newton system using cg and hessian-vector products
            if isempty(HvFunc)
                % No user-supplied Hessian-vector function,
                % use automatic differentiation
                HvFun = @autoHv;
                HvArgs = {x,g,useComplex,funObj,varargin{:}};
            else
                % Use user-supplid Hessian-vector function
                HvFun = HvFunc;
                HvArgs = {x,varargin{:}};
            end

            if useNegCurv
                [d,cgIter,cgRes,negCurv] = conjGrad([],-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFun,HvArgs);
            else
                [d,cgIter,cgRes] = conjGrad([],-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFun,HvArgs);
            end

            funEvals = funEvals+cgIter;
            if debug
                fprintf('newtonCG stopped on iteration %d w/ residual %.5e\n',cgIter,cgRes);

            end

            if useNegCurv
                if ~isempty(negCurv)
                    %if debug
                    fprintf('Using negative curvature direction\n');
                    %end
                    d = negCurv/norm(negCurv);
                    d = d/sum(abs(g));
                end
            end

        case NEWTON % Newton search direction

            if cgSolve == 0
                if HessianModify == 0
                    % Attempt to perform a Cholesky factorization of the Hessian
                    [R,posDef] = chol(H);

                    % If the Cholesky factorization was successful, then the Hessian is
                    % positive definite, solve the system
                    if posDef == 0
                        d = -R\(R'\g);

                    else
                        % otherwise, adjust the Hessian to be positive definite based on the
                        % minimum eigenvalue, and solve with QR
                        % (expensive, we don't want to do this very much)
                        if debug
                            fprintf('Adjusting Hessian\n');
                        end
                        H = H + eye(length(g)) * max(0,1e-12 - min(real(eig(H))));
                        d = -H\g;
                    end
                elseif HessianModify == 1
                    % Modified Incomplete Cholesky
                    R = mcholinc(H,debug);
                    d = -R\(R'\g);
                elseif HessianModify == 2
                    % Modified Generalized Cholesky
                    if useMex
                        [L D perm] = mcholC(H);
                    else
                        [L D perm] = mchol(H);
                    end
                    d(perm) = -L' \ ((D.^-1).*(L \ g(perm)));

                elseif HessianModify == 3
                    % Modified Spectral Decomposition
                    [V,D] = eig((H+H')/2);
                    D = diag(D);
                    D = max(abs(D),max(max(abs(D)),1)*1e-12);
                    d = -V*((V'*g)./D);
                elseif HessianModify == 4
                    % Modified Symmetric Indefinite Factorization
                    [L,D,perm] = ldl(H,'vector');
                    [blockPos junk] = find(triu(D,1));
                    for diagInd = setdiff(setdiff(1:p,blockPos),blockPos+1)
                        if D(diagInd,diagInd) < 1e-12
                            D(diagInd,diagInd) = 1e-12;
                        end
                    end
                    for blockInd = blockPos'
                        block = D(blockInd:blockInd+1,blockInd:blockInd+1);
                        block_a = block(1);
                        block_b = block(2);
                        block_d = block(4);
                        lambda = (block_a+block_d)/2 - sqrt(4*block_b^2 + (block_a - block_d)^2)/2;
                        D(blockInd:blockInd+1,blockInd:blockInd+1) = block+eye(2)*(lambda+1e-12);
                    end
                    d(perm) = -L' \ (D \ (L \ g(perm)));
                else
                    % Take Newton step if Hessian is pd,
                    % otherwise take a step with negative curvature
                    [R,posDef] = chol(H);
                    if posDef == 0
                        d = -R\(R'\g);
                    else
                        if debug
                            fprintf('Taking Direction of Negative Curvature\n');
                        end
                        [V,D] = eig(H);
                        u = V(:,1);
                        d = -sign(u'*g)*u;
                    end
                end

            else
                % Solve with Conjugate Gradient
                cgMaxIter = p;
                cgForce = min(0.5,sqrt(norm(g)))*norm(g);

                % Select Preconditioner
                if cgSolve == 1
                    % No preconditioner
                    precondFunc = [];
                    precondArgs = [];
                elseif cgSolve == 2
                    % Diagonal preconditioner
                    precDiag = diag(H);
                    precDiag(precDiag < 1e-12) = 1e-12 - min(precDiag);
                    precondFunc = @precondDiag;
                    precondArgs = {precDiag.^-1};
                elseif cgSolve == 3
                    % L-BFGS preconditioner
                    if i == 1
                        old_dirs = zeros(length(g),0);
                        old_stps = zeros(length(g),0);
                        Hdiag = 1;
                    else
                        [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,t*d,corrections,debug,old_dirs,old_stps,Hdiag);
                    end
                    g_old = g;
                    if useMex
                        precondFunc = @lbfgsC;
                    else
                        precondFunc = @lbfgs;
                    end
                    precondArgs = {old_dirs,old_stps,Hdiag};
                elseif cgSolve > 0
                    % Symmetric Successive Overelaxation Preconditioner
                    omega = cgSolve;
                    D = diag(H);
                    D(D < 1e-12) = 1e-12 - min(D);
                    precDiag = (omega/(2-omega))*D.^-1;
                    precTriu = diag(D/omega) + triu(H,1);
                    precondFunc = @precondTriuDiag;
                    precondArgs = {precTriu,precDiag.^-1};
                else
                    % Incomplete Cholesky Preconditioner
                    opts.droptol = -cgSolve;
                    opts.rdiag = 1;
                    R = cholinc(sparse(H),opts);
                    if min(diag(R)) < 1e-12
                        R = cholinc(sparse(H + eye*(1e-12 - min(diag(R)))),opts);
                    end
                    precondFunc = @precondTriu;
                    precondArgs = {R};
                end

                % Run cg with the appropriate preconditioner
                if isempty(HvFunc)
                    % No user-supplied Hessian-vector function
                    [d,cgIter,cgRes] = conjGrad(H,-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs);
                else
                    % Use user-supplied Hessian-vector function
                    [d,cgIter,cgRes] = conjGrad(H,-g,cgForce,cgMaxIter,debug,precondFunc,precondArgs,HvFunc,{x,varargin{:}});
                end
                if debug
                    fprintf('CG stopped after %d iterations w/ residual %.5e\n',cgIter,cgRes);
                    %funEvals = funEvals + cgIter;
                end
            end

        case TENSOR % Tensor Method

            if numDiff
                % Compute 3rd-order Tensor Numerically
                [junk1 junk2 junk3 T] = autoTensor(x,useComplex,funObj,varargin{:});
            else
                % Use user-supplied 3rd-derivative Tensor
                [junk1 junk2 junk3 T] = feval(funObj, x, varargin{:});
            end
            options_sub.Method = 'newton';
            options_sub.Display = 'none';
            options_sub.TolX = tolX;
            options_sub.TolFun = tolFun;
            d = minFunc(@taylorModel,zeros(p,1),options_sub,f,g,H,T);

            if any(abs(d) > 1e5) || all(abs(d) < 1e-5) || g'*d > -tolX
                if debug
                    fprintf('Using 2nd-Order Step\n');
                end
                [V,D] = eig((H+H')/2);
                D = diag(D);
                D = max(abs(D),max(max(abs(D)),1)*1e-12);
                d = -V*((V'*g)./D);
            else
                if debug
                    fprintf('Using 3rd-Order Step\n');
                end
            end
    end

    if ~isLegal(d)
        fprintf('Step direction is illegal!\n');
        pause;
        return
    end

    % ****************** COMPUTE STEP LENGTH ************************

    % Directional Derivative
    gtd = g'*d;

    % Check that progress can be made along direction
    if gtd > -tolX
        exitflag=2;
        msg = 'Directional Derivative below TolX';
        break;
    end

    % Select Initial Guess
    if i == 1
        if method < NEWTON0
            t = min(1,1/sum(abs(g)));
        else
            t = 1;
        end
    else
        if LS_init == 0
            % Newton step
            t = 1;
        elseif LS_init == 1
            % Close to previous step length
            t = t*min(2,(gtd_old)/(gtd));
        elseif LS_init == 2
            % Quadratic Initialization based on {f,g} and previous f
            t = min(1,2*(f-f_old)/(gtd));
        elseif LS_init == 3
            % Double previous step length
            t = min(1,t*2);
        elseif LS_init == 4
            % Scaled step length if possible
            if isempty(HvFunc)
                % No user-supplied Hessian-vector function,
                % use automatic differentiation
                dHd = d'*autoHv(d,x,g,0,funObj,varargin{:});
            else
                % Use user-supplid Hessian-vector function
                dHd = d'*HvFunc(d,x,varargin{:});
            end

            funEvals = funEvals + 1;
            if dHd > 0
                t = -gtd/(dHd);
            else
                t = min(1,2*(f-f_old)/(gtd));
            end
        end

        if t <= 0
            t = 1;
        end
    end
    f_old = f;
    gtd_old = gtd;

    % Compute reference fr if using non-monotone objective
    if Fref == 1
        fr = f;
    else
        if i == 1
            old_fvals = repmat(-inf,[Fref 1]);
        end

        if i <= Fref
            old_fvals(i) = f;
        else
            old_fvals = [old_fvals(2:end);f];
        end
        fr = max(old_fvals);
    end

    computeHessian = 0;
    if method >= NEWTON
        if HessianIter == 1
            computeHessian = 1;
        elseif i > 1 && mod(i-1,HessianIter) == 0
            computeHessian = 1;
        end
    end

    % Line Search
    f_old = f;
    if LS < 3 % Use Armijo Bactracking
        % Perform Backtracking line search
        if computeHessian
            [t,x,f,g,LSfunEvals,H] = ArmijoBacktrack(x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,LS_saveHessianComp,funObj,varargin{:});
        else
            [t,x,f,g,LSfunEvals] = ArmijoBacktrack(x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,1,funObj,varargin{:});
        end
        funEvals = funEvals + LSfunEvals;

    elseif LS < 6
        % Find Point satisfying Wolfe

        if computeHessian
            [t,f,g,LSfunEvals,H] = WolfeLineSearch(x,t,d,f,g,gtd,c1,c2,LS,25,tolX,debug,doPlot,LS_saveHessianComp,funObj,varargin{:});
        else
            [t,f,g,LSfunEvals] = WolfeLineSearch(x,t,d,f,g,gtd,c1,c2,LS,25,tolX,debug,doPlot,1,funObj,varargin{:});
        end
        funEvals = funEvals + LSfunEvals;
        x = x + t*d;

    else
        % Use Matlab optim toolbox line search
        [t,f_new,fPrime_new,g_new,LSexitFlag,LSiter]=...
            lineSearch({'fungrad',[],funObj},x,p,1,p,d,f,gtd,t,c1,c2,-inf,maxFunEvals-funEvals,...
            tolX,[],[],[],varargin{:});
        funEvals = funEvals + LSiter;
        if isempty(t)
            exitflag = -2;
            msg = 'Matlab LineSearch failed';
            break;
        end

        if method >= NEWTON
            [f_new,g_new,H] = funObj(x + t*d,varargin{:});
            funEvals = funEvals + 1;
        end
        x = x + t*d;
        f = f_new;
        g = g_new;
    end

    % Output iteration information
    if verboseI
        fprintf('%10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g)));
    end

    if logfile
        fid = fopen(logfile, 'a');
        if (fid > 0)
            fprintf(fid, '-- %10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g)));
            fclose(fid);
        end
    end

    % Output Function
    if ~isempty(outputFcn)
        callOutput(outputFcn,x,'iter',i,funEvals,f,t,gtd,g,d,sum(abs(g)),varargin{:});
    end

    % Update Trace
    trace.fval(end+1,1) = f;
    trace.funcCount(end+1,1) = funEvals;

    % Check Optimality Condition
    if sum(abs(g)) <= tolFun
        exitflag=1;
        msg = 'Optimality Condition below TolFun';
        break;
    end

    % ******************* Check for lack of progress *******************

    if sum(abs(t*d)) <= tolX
        exitflag=2;
        msg = 'Step Size below TolX';
        break;
    end

    if abs(f-f_old) < tolX
        exitflag=2;
        msg = 'Function Value changing by less than TolX';
        break;
    end

    % ******** Check for going over iteration/evaluation limit *******************

    if funEvals*funEvalMultiplier > maxFunEvals
        exitflag = 0;
        msg = 'Exceeded Maximum Number of Function Evaluations';
        break;
    end

    if i == maxIter
        exitflag = 0;
        msg='Exceeded Maximum Number of Iterations';
        break;
    end

end

if verbose
    fprintf('%s\n',msg);
end
if nargout > 3
    output = struct('iterations',i,'funcCount',funEvals*funEvalMultiplier,...
        'algorithm',method,'firstorderopt',sum(abs(g)),'message',msg,'trace',trace);
end

% Output Function
if ~isempty(outputFcn)
    callOutput(outputFcn,x,'done',i,funEvals,f,t,gtd,g,d,sum(abs(g)),varargin{:});
end

end

max_stars_repo_name: ebertolazzi/NLtoolbox

max_issues_repo_name: ebertolazzi/NLtoolbox

max_forks_repo_name: ebertolazzi/NLtoolbox


text:
# Postprocess Definition File for Code Generation
# Init
# Nachverarbeitung der Datei robot_env. Falls dort mit Maple-Programmcode Teile der DH-Parameter definiert werden, wird diese Datei ausgewertet und die Variablen neu gespeichert.
# Benötigt, da in robot_codegen_tmpvar_bash.sh ein definiertes Muster für die robot_env erwartet wird.
#
# Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2020-11
# (C) Institut fuer Mechatronische Systeme, Leibniz Universitaet Hannover
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
# Lese Umgebungsvariable für Codegenerierung.
read "../robot_codegen_definitions/robot_env":
printf("Nachverarbeitung der Eingabedatei für %s\n",robot_name):

# Postprocess
# Speichere die Variablen neu ab
NQJ := NQJ:
NJ := NJ:
robot_name := robot_name:
save NQJ, NJ, robot_name, "../robot_codegen_definitions/robot_env2":

max_stars_repo_name: SchapplM/robsynth-modelgen

max_issues_repo_name: SchapplM/robsynth-modelgen

max_forks_repo_name: SchapplM/robsynth-modelgen


text:
(*
 Copyright (C) 2020 Susi Lehtola

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

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

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

rda_s := x -> X2S*x:
rda_p := u -> X2S^2*u:

(* Equation (61) *)
rda_k4 := (s, p, b) -> sqrt(s^4 + b*p^2):
(* Equation (63) *)
rda_k2 := (s, p, b) -> s^2 + b*p:

(* Equation (71); first term is von Weiszäcker according to equation (13) *)
rda_f0 := (s, p) ->
       5/3*s^2 + params_a_A0
       + params_a_A1 * (rda_k4(s,p,params_a_a) / (1 + params_a_beta1*rda_k4(s,p,params_a_a)))^2
       + params_a_A2 * (rda_k4(s,p,params_a_b) / (1 + params_a_beta2*rda_k4(s,p,params_a_b)))^4
       + params_a_A3 * (rda_k2(s,p,params_a_c) / (1 + params_a_beta3*rda_k2(s,p,params_a_c))):

(* Complete functional *)
rda_f := (xs, us) -> rda_f0(rda_s(xs), rda_p(us)):

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


text:
# Berechne kinematische Zwangsbedingungen für Fünfgelenkkette

# Einleitung
# Die kinematischen Zwangsbedingungen werden als Ersetzungsausdruck für die abhängigen Winkel aufgestellt.
#
# fivebar1TE -> Fünfgelenkkette, Modellierung der Zwangsbedingungen mit trigonometrischer Elimination der Winkel
# kinematic_constraint -> Berechnungen zu kinematischen Zwangsbedingungen

# Quelle
# SA Bejaoui: Bejaoui2018_S749; "Modellierung kinematischer Zwangsbedingungen für hybride serielle Roboter mit planaren Parallelmechanismen"
# (Fünfgelenkkette S. 23 in [Bejaoui2018_S749])
# Autor
# Abderahman Bejaoui, Moritz Schappler
# Studienarbeit bei: Moritz Schappler, moritz.schappler@imes.uni-hannover.de, 2019-01
# (C) Institut fuer Mechatronische Systeme, Leibniz Universitaet Hannover
# Initialisierung
interface(warnlevel=0): # Unterdrücke die folgende Warnung.
restart: # Gibt eine Warnung, wenn über Terminal-Maple mit read gestartet wird.
interface(warnlevel=3):
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
#with(ListTools):
codegen_act := true:
codegen_opt := 2:
read "../helper/proc_MatlabExport":
read "../helper/proc_simplify2":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions",robot_name):
use_simplify := 0: # keine zusätzlichen simplify-Befehle
;
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:
printf("%s. Beginn der Berechnungen.\n", FormatTime("%Y-%m-%d %H:%M:%S")):
st := time():
# Zwangsbedingungen mit Kreisschnittpunkt D zur Bestimmung von rho
r_p1:=<AE*cos(qJ1s), AE*sin(qJ1s)>:
r_p2:=<AB+BC*cos(qJ2s), BC*sin(qJ2s)>:
r1:=ED:
r2:=CD:
A:=intersect_circle(r_p1, r_p2, r1, r2):
# Es gibt zwei Lösungen für den Kreisschnittpunkt. Bei der gewählten Definition der Kreise ist die erste Lösung in der Fünfgelenkkette und die zweite außen (für die gewählten Beispiel-Parameter).
# Die zweite Lösung entspricht der gewählten Lösung in [Bejaoui2018_S749]
r_s1x:=simplify(A(1,2),trig):
r_s1y:=simplify(A(2,2),trig):
# Zwangsbedingungen für eta, zweiter Weg zum Kreisschnittpunkt D über B
radx1:=AB+BC*cos(qJ2s)+CD*(cos(qJ2s)*cos(eta_s)-sin(qJ2s)*sin(eta_s)):
rady1:=BC*sin(qJ2s)+CD*(sin(qJ2s)*cos(eta_s)+cos(qJ2s)*sin(eta_s)):
f11:=r_s1x-radx1=0: #GL3
;
f12:=r_s1y-rady1=0: #GL4
;
# Zwangsbedingungen für xi , zwei Wege Punkt E
racx1:=AE*cos(qJ1s):
racy1:=AE*sin(qJ1s):
racx2:=AB+BC*cos(qJ2s)+CD*(cos(qJ2s)*cos(eta_s)-sin(qJ2s)*sin(eta_s))+ED*(cos(qJ2s)*(cos(eta_s)*cos(xi_s)-sin(eta_s)*sin(xi_s))-sin(qJ2s)*(cos(xi_s)*sin(eta_s)+cos(eta_s)*sin(xi_s))):
racy2:=BC*sin(qJ2s)+CD*(sin(qJ2s)*cos(eta_s)+cos(qJ2s)*sin(eta_s))+ED*(sin(qJ2s)*(cos(eta_s)*cos(xi_s)-sin(eta_s)*sin(xi_s))+cos(qJ2s)*(cos(xi_s)*sin(eta_s)+cos(eta_s)*sin(xi_s))):
f21:=racx1-racx2=0:  # GL5
;
f22:=racy1-racy2=0:  # GL6
;
# Bestimme Winkel
# eta bestimmen
sol1:=solve({f11,f12},{cos(eta_s),sin(eta_s)}):
assign(sol1):
cos_eta_s:=simplify(rhs(sol1[1]),trig):
sin_eta_s:=simplify(rhs(sol1[2]),trig):
# xi bestimmen
sol2:=solve({f21,f22},{cos(xi_s),sin(xi_s)}):
sol2:=op(sol2):
cos_xi_s:=simplify(rhs(sol2[1]),trig):
sin_xi_s:=simplify(rhs(sol2[2]),trig):
# rho bestimmen
cos_psi_phi:=cos(qJ2s)*cos(qJ1s)+sin(qJ2s)*sin(qJ1s): # cos(psi-phi)
;
cos_eta_xi:=cos_eta_s*cos_xi_s-sin_eta_s*sin_xi_s: # cos(eta+xi)
;
sin_psi_phi:=sin(qJ2s)*cos(qJ1s)-cos(qJ2s)*sin(qJ1s): # sin(psi-phi)
;
sin_eta_xi:=sin_eta_s*cos_xi_s+cos_eta_s*sin_xi_s: #sin(eta+xi)
;
cos_rho_s:=cos_psi_phi*cos_eta_xi-sin_psi_phi*sin_eta_xi:
sin_rho_s:=sin_psi_phi*cos_eta_xi+cos_psi_phi*sin_eta_xi:
# Cos-und Sinus-Terme zuweisen

kintmp_subsexp(1,2):=sin_rho_s:
kintmp_subsexp(2,2):=cos_rho_s:
kintmp_subsexp(3,2):=sin_eta_s:
kintmp_subsexp(4,2):=cos_eta_s:
kintmp_subsexp(5,2):=sin_xi_s:
kintmp_subsexp(6,2):=cos_xi_s:
# Terme vereinfachen

if use_simplify>=1 then # bringt leider nichts
  tmp_t0 := time():
  tmp_l0 := length(kintmp_subsexp):
  printf("%s: Vereinfache Ersetzungsausdrücke für Zwangsbedingungen. Länge: %d.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l0):
  for i from 1 to RowDimension(kintmp_subsexp) do
    tmp_t1 := time():
    tmp_l1 := length(kintmp_subsexp(i,2)): # es kann sowieso nur einer der beiden Informationen enthalten
    printf("%s: Vereinfache Ersetzungsausdruck %d (%s). Länge: %d.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i, String(kintmp_subsexp(i,1)), tmp_l1):
    kintmp_subsexp(i,2) := simplify2(kintmp_subsexp(i,2)):
    tmp_t2 := time():
    tmp_l2 := length(kintmp_subsexp(i,2)):
    printf("%s: Ersetzungsausdruck %d (%s) vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
      FormatTime("%Y-%m-%d %H:%M:%S"), i, String(kintmp_subsexp(i,1)), tmp_l1, tmp_l2, tmp_t2-tmp_t1):
  end do:
  tmp_l3 := length(kintmp_subsexp):
  printf("%s: Ersetzungsausdrücke vereinfacht. Länge: %d->%d. Rechenzeit %1.1fs.\n", \
    FormatTime("%Y-%m-%d %H:%M:%S"), tmp_l0, tmp_l3, tmp_t2-tmp_t0):
end if:
# Winkel zuweisen:
kintmp_qs(1,1):=%arctan(sin_rho_s,cos_rho_s):
kintmp_qs(2,1):=%arctan(sin_eta_s,cos_eta_s):
kintmp_qs(3,1):=%arctan(sin_xi_s,cos_xi_s):
# Exportiere Code für folgende Skripte
#
# Die Annahmen sind im Ausdruck bereits in den Variablen gespeichert. Ersetze das "~"-Zeichen.
# Quelle: http://www.mapleprimes.com/questions/207601-Remove-Assumptions
winkel := <<sin_xi|sin(xi_s)>; <cos_xi|cos(xi_s)>; <sin_rho | sin(rho_s)>; <cos_rho | cos(rho_s)>;<sin_eta | sin(eta_s)>; <cos_eta | cos(eta_s)>>:
save winkel, sprintf("../codeexport/fivebar1TE/tmp/kinematic_constraints_maple_inert"):
kintmp_qt := convert_s_t(kintmp_qs):
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
printf("%s. Ausdrücke für kintmp_subsexp gespeichert (Maple). CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d", robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt,kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
printf("%s. Ausdrücke mit Inert-Arctan exportiert (Matlab). CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp ));
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m",robot_name),2);
printf("%s. Zwangsbedingungen berechnet. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):

#

max_stars_repo_name: SchapplM/robsynth-serhybroblib

max_issues_repo_name: SchapplM/robsynth-serhybroblib

max_forks_repo_name: SchapplM/robsynth-serhybroblib


text:
`is_element/cacti` := (A::set) -> proc(C)
 global reason;
 local tag,n,i,j,a,b,U;

 tag := "is_element/cacti";

 if not(is_table_on(A)(C)) then
  reason := [tag,"C is not a table indexed by A",C,A];
  return false;
 fi;

 for a in A do
  if not(`is_element/cactus_lobes`(C[a])) then
   reason := [tag,"C[a] is not a cactus lobe",a,C[a],reason];
   return false;
  fi;
 od;

 n := nops(A);
 for i from 1 to n do
  for j from i+1 to n do
   a := A[i];
   b := A[j];
   if `are_crossed/RZ_set`(C[a][2],C[b][2]) then
    reason := [tag,"lobes C[a] and C[b] are crossed",a,b,C[a],C[b]];
    return false;
   fi;
  od;
 od;

 U := `empty/RZ_set`;
 for a in A do
  U := `union/RZ_set`(U,C[a][2]);
 od;

 if U <> `top/RZ_set` then
  reason := [tag,"The lobes do not fill the whole circle",C,U];
  return false;
 fi;

 return true;
end:

`is_equal/cacti` := (A::set) -> proc(C0,C1)
 local a;
 `and`(seq(evalb(C0[a]=C1[a]),a in A));
end:

`is_leq/cacti` := NULL:

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

`random_element/cacti` := (A) -> proc()
 local J,s,s1,n,l,m,c,T,i,j,a,C;

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

 J,s := op(`random_element/cactus_planar_trees`(A)());

 if nops(A) = 1 then
  return table([op(A) = `id/cactus_lobes`]);
 fi;

 s1 := map(e -> e[1],s);
 n := nops(s1);
 l := [seq(rand(1..6)(),i=1..n)];
 l := l /~ `+`(op(l));
 m := [seq(add(l[i],i=1..j),j=0..n)];
 c := table([seq(a = `random_element/RZ`(),a in A)]);
 T := table([seq(a = [],a in A)]):
 for i from 1 to n do
  T[s1[i]] := [op(T[s1[i]]),m[i],m[i+1]];
 od:
 C := table():
 for a in A do C[a] := [c[a],T[a]]; od:
 C := `source_shift/cacti`(A)(`random_element/RZ`())(C):
 return eval(C);
end:

`list_elements/cacti` := NULL:
`count_elements/cacti` := NULL:

######################################################################
# Here t should be in R/Z

`source_shift/cacti` := (A::set) -> (t) -> proc(C)
 local D,a;

 D := table():
 for a in A do
  D[a] := `source_shift/cactus_lobes`(t)(C[a]);
 od:
 return eval(D);
end:

######################################################################
# Here t should be a table with indices A and entries in R/Z

`target_shift/cacti` := (A::set) -> (t) -> proc(C)
 local D,a;

 D := table():
 for a in A do
  D[a] := `target_shift/cactus_lobes`(t[a])(C[a]);
 od:
 return eval(D);
end:

######################################################################
# The standard basepoint for the space of cacti on {1,..,n}

`basepoint/cacti` := proc(n::posint)
 local i,C;
 C := table():
 for i from 1 to n do
  C[i] := [0,[(i-1)/n,i/n]];
 od:
 return eval(C);
end:

######################################################################
# A path in the space of cacti on {1,..,n}, from the basepoint to
# its image under the transposition (m,m+1).

`omega/cacti` := (n::posint) -> (m::posint) -> proc(t)
 local A,C,i,k;
 if m >= n then error("m is out of range"); fi;
 A := {seq(i,i=1..n)};
 C := table():
 for k from 1 to n do
  if k = m then
   C[k] := [0,[m-1+t,m+t] /~ n]
  elif k = m+1 then
   if t = 0 then
    C[k] := [0,[m,m+1] /~ n];
   elif t = 1 then
    C[k] := [0,[m-1,m] /~ n];
   else
    C[k] := [0,[m-1,m-1+t,m+t,m+1] /~ n];
   fi;
  else
   C[k] := [0,[k-1,k] /~ n];
  fi;
 od:
 return eval(C);
end:

######################################################################
# A family of paths in the space of cacti on {1,2,3}

`alpha/cacti` := proc(i,t)
 local L,C;
 if t < 0 or t > 1 then
  error("t out of range");
 fi;
 if t = 0 then
  L := [[1,2] /~ 3,[0,1] /~ 3,[2,3] /~ 3];
 elif t = 1 then
  L := [[0,1] /~ 3,[1,2] /~ 3,[2,3] /~ 3];
 else
  L := [[0,t,1+t,2] /~ 3,[t,1+t] /~ 3,[2,3] /~ 3];
 fi;
 L := map(u -> [0,u],L);
 C := table():
 if i = 1 then
  C[1] := L[1]; C[2] := L[2]; C[3] := L[3];
 elif i = 2 then
  C[1] := L[3]; C[2] := L[1]; C[3] := L[2];
 elif i = 3 then
  C[1] := L[2]; C[2] := L[3]; C[3] := L[1];
 else
  error("i out of range");
 fi;
 return eval(C);
end:

######################################################################
# A family of paths in the space of cacti on {1,2,3}

`beta/cacti` := proc(z,i,t)
 local C,w;

 C := `alpha/cacti`(i,t);
 if i = 1 then
  w := z;
 elif i = 2 then
  w := z - (1+t)/3;
 elif i = 3 then
  w := z + (1+t)/3;
 else
  error("i is out of range");
 fi;

 return `source_shift/cacti`({1,2,3})(w)(C);
end:

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

`planar_tree/cacti` := (A::set) -> proc(C)
 local S,SI,D,E,E0,E1,EI,J,t0,t1,t,j,jt,u,v,a,ok,s;

 S := table():
 D := table():
 SI := table():

 for a in A do
  S[a] := `starts/RZ_set`(C[a][2]);
  D[a] := `boundary/RZ_set`(C[a][2]);
  for t in S[a] do SI[t] := a; od;
 od:

 E1 := map(op,{seq(S[a],a in A)});
 E1 := sort([op(E1)]);

 J := {};
 E0 := E1;
 EI := table():
 while E0 <> [] do
  t0 := E0[1];
  jt := {t0};
  j := {SI[t0]};
  for t1 in [op(2..-1,E0)] do
   u := `interval/RZ_set`(t0,t1);
   v := `interval/RZ_set`(t1,t0);
   ok := true;
   for a in A do
    if not (`intersect/RZ_set`(u,C[a][2]) = `empty/RZ_set` or
            `intersect/RZ_set`(v,C[a][2]) = `empty/RZ_set`) then
     ok := false;
     break;
    fi;
   od;
   if ok then
    jt := {op(jt),t1};
    j := {op(j),SI[t1]};
   fi;
  od:
  J := {op(J),j};
  E0 := select(t -> not(member(t,jt)),E0);
  for t in jt do EI[t] := [SI[t],j]; od;
 od:

 s := [seq(EI[t],t in E1)];
 return [J,s];
end:

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

# This computes the matrix x as in lem-cactus-x

`x/cacti` := (A::set) -> proc(C)
 local x,a,a1;

 x := table();
 for a1 in A do
  for a in A do
   if a1 <> a then
    x[a,a1] := `eval/cactus_lobes`(C[a])(C[a1][2][1]);
   fi;
  od:
 od:

 return eval(x);
end:

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

`eta/cacti` := (A::set) ->
 `if`(nops(A) = 1,table([op(A) = `id/cactus_lobes`]),FAIL):

`gamma/cacti` := (A::set,B::set) -> (p) -> proc(D,C)
 local a,E;

 E := table():
 for a in A do
  E[a] := `o/cactus_lobes`(C[p[a]][a],D[p[a]]);
 od:
 return eval(E);
end:

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

`phi/ord_simplex_interior/cacti` := (A::set) -> proc(Rx)
 local R,x,n,C,y,i,j;
 R,x := op(Rx);
 n := nops(R);
 C := table():
 y := [seq(add(x[R[i]],i=1..j),j=0..n)];
 for i from 1 to n do
  C[R[i]] := [0,[y[i],y[i+1]]];
 od:
 return eval(C);
end:

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

`phi/cacti/simplex_interior` := (A::set) -> proc(C)
 local x,a;
 x := table();
 for a in A do
  x[a] := `length/RZ_set`(C[a][2]);
 od:
 return eval(x):
end:

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

`plot/cacti` := (A::set) -> proc(C,r := 1,c := [0,0])
 local n,z,i;
 n := nops(A);
 z := table([seq(i = `centroid_R2/RZ_set`(C[A[i]][2]),i=1..n)]);
 display(
  seq(`plot/RZ_set`(C[A[i]][2],r,c,colour=standard_colour(i)),i=1..n),
  seq(point(z[i],colour=standard_colour(i)),i=1..n)
 );
end:


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

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

(* type: mgga_exc *)

(* prefix:
  mgga_x_mvs_params *params;

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

(* equation 10 *)
mvs_fa := a -> (1 - a) / ((1 + params_a_e1*a^2)^2 + params_a_c1*a^4)^(1/4):

(* alpha *)
mvs_alpha := (t,x) -> (t - x^2/8)/K_FACTOR_C:

(* eq 7 *)
mvs_f := (x, u, t) -> (1 + params_a_k0*mvs_fa(mvs_alpha(t,x)))
      / (1 + params_a_b*(X2S*x)^4)^(1/8):

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


text:
# Andrews-Paule identity. See §7.1 in “Multiple binomial sums”.

S := Sum(z^n*Sum(Sum(Multinomial([i,j])^2*Multinomial([2*n-2*j,2*n-2*i]), i=0..infinity), j=0..infinity), n=0..infinity);

R, ord := BinomSums[sumtores](S, u);

# Checking
BinomSums[rser](R, ord, 6);
BinomSums[computesum](S, 5);

# We obtain the following differential operator annihilating the generating function S:
dop := 1048576*D^6*z^8 + 2883584*D^6*z^7-40960*D^6*z^6 +
22020096*D^5*z^7-29696*D^6*z^5 + 64225280*D^5*z^6 + 1296*D^6*z^4 +
1269760*D^5*z^5 + 146407424*D^4*z^6-605184*D^5*z^4 + 455041024*D^4*z^5 +
16848*D^5*z^3 + 24412672*D^4*z^4 + 363069440*D^3*z^5-3632352*D^4*z^3 +
1211465728*D^3*z^4 + 59292*D^4*z^2 + 106845184*D^3*z^3 +
305827840*D^2*z^4-7352832*D^3*z^2 + 1109626112*D^2*z^3 + 58320*D^3*z +
139138736*D^2*z^2 + 60272640*D*z^3-4247073*D^2*z + 244005120*D*z^2 + 9720*D^2 +
42117840*D*z + 691200*z^2-374625*D + 3369600*z + 996300;

# Conversion to a differential equation with initial condition
init := BinomSums[computesum](S, 6):
deq := {
  DETools[diffop2de](dop, y(z), [D,z]),
  seq( (D@@i)(y)(0) = i!*coeff(init, z, i), i=0..4 )
  };

# Conversion to a recurrence.
rec := gfun[diffeqtorec](deq, y(z), u(n));

# Resolution of the recurrence.
rsolve(rec, u(n));

max_stars_repo_name: lairez/binomsum

max_issues_repo_name: lairez/binomsum

max_forks_repo_name: lairez/binomsum


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

$include "mgga_c_scan.mpl"

(* Coefficients of the rSCAN switching function from SI, in reversed(!) order: 7, 6, ..., 0 *)
rscan_fc := [-0.051848879792, 0.516884468372, -1.915710236206, 3.061560252175, -1.535685604549, -0.4352, -0.64, 1]:

np53 := rs -> n_total(rs)^(5/3):
(* First regularization: tau^u -> tau^u + tau^r *)

rscan_alpha0 := (rs, z, xt, ts0, ts1) ->
  (np53(rs)*m_max(t_total(z, ts0, ts1) - xt^2/8, 0))/((K_FACTOR_C*np53(rs) + 2^(2/3)*params_a_taur)*t_total(z, 1, 1)):

(* Second regularization: alpha -> alpha^3/(alpha^2 + alpha_r) *)
rscan_alpha := (rs, z, xt, ts0, ts1) -> rscan_alpha0(rs, z, xt, ts0, ts1)^3/(rscan_alpha0(rs, z, xt, ts0, ts1)^2 + params_a_alphar):

(* f(alpha) replaced with a polynomial for alpha in [0, 2.5] *)
rscan_f_alpha_small := (a,ff) -> add(ff[8-i]*a^i, i=0..7):
rscan_f_alpha_large := a -> -params_a_d*exp(params_a_c2/(1 - a)):
rscan_f_alpha := (a, ff) -> my_piecewise3( a <= 2.5, rscan_f_alpha_small(m_min(a, 2.5),ff), rscan_f_alpha_large(m_max(a, 2.5)) ):

(* set parameters of f_alpha *)
params_a_alphar := 1e-3:
params_a_taur := 1e-4:

rscan_f := (rs, z, xt, xs0, xs1, ts0, ts1) ->
  f_pbe(rs, z, xt, xs0, xs1) + rscan_f_alpha(rscan_alpha(rs, z, xt, ts0, ts1), rscan_fc)*(
    + scan_e0(rs, z, X2S*2^(1/3)*xt)
    - f_pbe(rs, z, xt, xs0, xs1)
  ):

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


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

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

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

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

Cf := 3/10 * (3*Pi^2)^(2/3):

omega := rs -> params_a_B*exp(-params_a_c*rs)/(1 + params_a_d*rs):
delta := rs -> (params_a_c + params_a_d/(1 + params_a_d*rs))*rs:

aux6 := 1/2^(8/3):
aux4 := aux6/4:
aux5 := aux4/(9*2):

t1 := (rs, z) ->
  -(1 - z^2)/(1 + params_a_d*rs):
t2 := (rs, z, xt) ->
  -xt^2*((1 - z^2)*(47 - 7*delta(rs))/(4*18) - 2/3):
t3 := (z) ->
  -Cf/2*(1 - z^2)*((1 + z)^(8/3) + (1 - z)^(8/3)):
t4 := (rs, z, xs0, xs1) ->
  aux4*(1 - z^2)*(5/2 - delta(rs)/18)*(xs0^2*(1 + z)^(8/3) + xs1^2*(1 - z)^(8/3)):
t5 := (rs, z, xs0, xs1) ->
  aux5*(1 - z^2)*(delta(rs) - 11)*(xs0^2*(1 + z)^(11/3) + xs1^2*(1 - z)^(11/3)):
t6 := (z, xs0, xs1) ->
  -aux6*(2/3*(xs0^2*(1 + z)^(8/3) + xs1^2*(1 - z)^(8/3))
  -(1 + z)^2*xs1^2*(1 - z)^(8/3)/4 - (1 - z)^2*xs0^2*(1 + z)^(8/3)/4):

f_lyp := (rs, z, xt, xs0, xs1) -> params_a_A*(t1(rs/RS_FACTOR, z) + omega(rs/RS_FACTOR)*(
  + t2(rs/RS_FACTOR, z, xt) + t3(z) + t4(rs/RS_FACTOR, z, xs0, xs1)
  + t5(rs/RS_FACTOR, z, xs0, xs1) + t6(z, xs0, xs1)
)):

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


text:
#!/bin/bash maple
# Use `maple -q expand7.mpl -D n=20` to run

e := (1 + sqrt(3) * x + sqrt(5) * y) ^ n:
f := e * (e + sqrt(7)):

st := time[real]():
f := expand(f):
1000*(time[real]() - st);

done

max_issues_repo_name: HQSquantumsimulations/symengine

max_forks_repo_name: HQSquantumsimulations/symengine


text:
func &foo ( var %i i32 ) void {
   var %p ptr
   dassign %p ( malloc ptr ( constval i32 8 ))
   dassign %p ( alloca ptr ( constval i32 16 ))
   dassign %p ( gcmalloc ref i32 )
   dassign %p ( gcmallocjarray ref <[] i32> ( constval i32 5 ))
   free ( dread ref %p)
}

 # 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:
SemiAlgebraic := module()
    uses Domain, Domain_Type, SolveTools, Utilities;
    export SimplName  := "SemiAlgebraic";
    export SimplOrder := 6+(1/2);

    export ModuleApply := proc(dbnds :: DomBound, dshape :: DomShape, $)
        local vs;
        vs := ListTools[Reverse](Domain:-Bound:-varsOf(dbnds));
        vs := map(v->if v::list then op(1,v) else v end if, vs);
        do_sh( dshape , dbnds, vs );
    end proc;

    local extra_sol := (c -> c::`=` and (lhs(c)::DomBoundVar and evalb(lhs(c)=rhs(c))));
    local postproc := proc(sol, $)
      foldl((x,k)->subsindets(x,op(k)), sol
           ,[specfunc('piecewise')
            ,x->DSplit(Partition:-PWToPartition(x,'check_valid'))]
           ,[list(list),DSum@op]
           ,[list(relation),x->DConstrain(remove(extra_sol,x)[])]);
    end proc;

    # Check for trivial solutions, which are rejected because they are not an
    # improvement in terms of Hakaru semantics (but SemiAlgebraic doesn't know
    # anything about that)
    local trivial_sol := proc(sh :: DomShape, ctx :: DomBound, $)
      local ps, vs, vsc;
      if sh::DomConstrain then
        nops(sh)=0
      elif sh::DomSplit then
        # we assume that the conditions of a solution of SemiAlgebraic cover the
        # entire domain, so it is left to check only the values
        ps := KB:-kb_Partition( op(1,sh), Domain:-Bound:-contextOf(ctx),
                               ((x_,y_)->true), #unused dummy to pass typechecks
                               ((v,kb)->trivial_sol
                                  (v,Domain:-Bound:-onContext(_->kb,ctx))) );
        `and`(op(map(evalb@Partition:-valOf,Partition:-piecesOf(ps))));
      elif sh::DomSum then
        if nops(sh)<>2 then return false end if;
        if ormap(s->not(s::DomConstrain) or nops(s)<>1,sh)
        then return false; end if;
        vs := Domain:-Bound:-varsOf(ctx,"set");
        vsc := map((v->(v,-v)), vs);
        ps := map2(op,1,[op(sh)]);
        ps := map(p->classify_relation(p,v->v in vsc), ps);
        if not(has(ps, {B_LO,B_HI})) then return false; end if;
        if not(the(map2(op,3,ps))) then return false; end if;
        if not(the(map(abs,map2(op,4,ps)))) then return; false end if;
        return true;
      elif sh::DomInto then
        error "solution should not have DIntos: %1", sh;
      else
        error "unrecognized DomShape: %1", sh;
      end if;
    end proc;

    local do_sh := proc( sh :: DomShape, ctx :: DomBound, vs, $)
        local sol;
        if sh :: DomConstrain then
            sol := do_Constrain(sh, ctx, vs);
            sol := `if`(sol::DomShape, sol, postproc(sol));
            if not(sol::DomNoSol) and trivial_sol(sol, ctx) then
              userinfo(3, Domain:-Improve,
                       printf("Discarding solution: %a\n", sol));
              sh;
            else sol end if;
        elif sh :: DomSplit then
            DSplit(Partition:-Amap(
              [ (c,_)->c
              , (p,c)-> p->do_sh(p, applyop(cs->{c,op(cs)},2,ctx), vs)
              , c->c ], op(1, sh) ));
        elif sh :: DomSum then
            map(s->do_sh(s, ctx, vs), sh);
        else
            DNoSol(sprintf("Don't know how to solve DOMAIN(%a, %a)", ctx, sh));
        end if;
    end proc;

    local do_Constrain := proc( sh :: DomConstrain , ctx, vs_, $ )
        local ts, vs := vs_, cs;
        cs := { op( Domain:-Bound:-toConstraints(ctx,'no_infinity') ), op(sh) } ;
        ts, cs := selectremove(type, cs, `::`);
        ts := select(t->op(2,t)<>real, ts);
        if ts <> {} and has(cs,ts) then
          WARNING("Variables of unrecognized types being passed to SemiAlgebraic: %1",ts);
        end if;
        try SemiAlgebraic(cs, vs);
        catch: DNoSol(StringTools[FormatMessage](lastexception[2..-1]));
        end try;
    end proc;
end module; #SemiAlgebraic

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
MikMod playlist
"/usr/share/mikmod/mod/spacedeb.mod" "" 0 0
"/usr/share/mikmod/mod/GSLINGER.MOD" "" 0 0
"/usr/share/mikmod/mod/jt_911.xm" "" 0 0
"/usr/share/mikmod/mod/aryx.s3m" "" 0 0
"/usr/share/mikmod/mod/celestial_fantasia.s3m" "" 0 0
"/usr/share/mikmod/mod/aceman_-_my_first_console.xm" "" 0 0
"/usr/share/mikmod/mod/beyond_music.mod" "" 0 0
"/usr/share/mikmod/mod/ELYSIUM.MOD" "" 0 0
"/usr/share/mikmod/mod/stardstm.mod" "" 0 0
"/usr/share/mikmod/mod/eternity-1553.xm" "" 0 0
"/usr/share/mikmod/mod/strshine.s3m" "" 0 0
"/usr/share/mikmod/mod/BUTTERFL.XM" "" 0 0
"/usr/share/mikmod/mod/CTGOBLIN.S3M" "" 0 0
"/usr/share/mikmod/mod/sweetdre.xm" "" 0 0
"/usr/share/mikmod/mod/enigma.mod" "" 0 0
"/usr/share/mikmod/mod/aurora.mod" "" 0 0
"/usr/share/mikmod/mod/external.xm" "" 0 0
"/usr/share/mikmod/mod/SATELL.S3M" "" 0 0
"/usr/share/mikmod/mod/jt_letgo.xm" "" 0 0
"/usr/share/mikmod/mod/DEADLOCK.XM" "" 0 0


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

$define xc_dimensions_2d

_2d_b86_beta  := 0.002105:
_2d_b86_gamma := 0.000119:

_2d_b86_f := x -> (1 + _2d_b86_beta*x^2)/(1 + _2d_b86_gamma*x^2):

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


text:
TITLE
	SNDP_10_5_0_1

STOCHASTIC

DATA
	NrOfScen := DATAFILE("10_5_0_1_ScalarData.dat");
	NrOfLocations := DATAFILE("10_5_0_1_ScalarData.dat"); !last location is a final location (market)
	Market := NrOfLocations; !last location is a final location (market)
	NrOfProducts := DATAFILE("10_5_0_1_ScalarData.dat");
	EndProduct := NrOfProducts;	!last product is the end product
	SalesPrice := DATAFILE("10_5_0_1_ScalarData.dat");
	PlantCost := DATAFILE("10_5_0_1_ScalarData.dat");
	PlantCapacity := DATAFILE("10_5_0_1_ScalarData.dat");

SCENARIO
	SCEN := 1..NrOfScen;

PROBABILITIES
	Prob[SCEN] := SPARSEFILE("10_5_0_1_Prob.dat");

RANDOM DATA
	Demand[SCEN] := SPARSEFILE("10_5_0_1_Demand.dat");

INDEX
	location := 1..NrOfLocations;
	start := location;
	finish := location;
	product := 1..NrOfProducts;
	material[product] WHERE (product <> EndProduct);
	arc[start, finish]:= INDEXFILE("10_5_0_1_arc.dat");

DATA
	ShipCost[start, finish]:= SPARSEFILE("10_5_0_1_ShipCost.dat");
	ArcProduct[product, start, finish] := SPARSEFILE("10_5_0_1_ArcProduct.dat");
	MaterialReq[material] := SPARSEFILE("10_5_0_1_MaterialReq.dat");

INDEX
	plantLocation[location] WHERE (ArcProduct[product:=EndProduct, start:=location, finish:=Market] = 1);

STAGE1 BINARY VARIABLES
	OpenProduction[plantLocation];

STAGE2 VARIABLES
	Ship[product, start, finish] WHERE (ArcProduct[product, start, finish] = 1);

MACROS
	TotalRevenue := SUM(start:  SalesPrice * Ship[product:=EndProduct, start, finish:=Market]);
	FixedCost := SUM(plantLocation: PlantCost * OpenProduction);
	TotalShipCost := SUM(product,start,finish: ShipCost * Ship);

MAX
	Profit = TotalRevenue - FixedCost - TotalShipCost;

SUBJECT TO
	BOMConstr[plantLocation, material]:
		MaterialReq[material] * Ship[product:=EndProduct, start:=plantLocation, finish:=Market] = SUM(start: Ship[product:=material, start, finish:=plantLocation]);

	DemandConstr -> DemConstr: SUM(start: Ship[product:=EndProduct,start,finish:=Market]) <= Demand;

	PlantConstr[plantLocation]:
		Ship[product:=EndProduct, start:=plantLocation, finish:=Market] <= PlantCapacity * OpenProduction[plantLocation];

max_stars_repo_name: pashtetGP/sndpgen

max_issues_repo_name: pashtetGP/sndpgen

max_forks_repo_name: pashtetGP/sndpgen


text:
######################################################################
# simplex_interior(A) is the set of maps x : A -> (0,1] with sum = 1.

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

 if not `is_element/simplex`(A)(x) then
  reason := [convert(procname,string),"x is not in the A-simplex",x,A,reason];
  return false;
 fi;

 for a in A do if x[a] = 0 then
  reason := [convert(procname,string),"x[a] = 0",a];
  return false;
 fi; od;

 return true;
end;

`is_equal/simplex_interior` := (A::set) -> proc(x,y)
 `is_equal/simplex`(A)(x,y);
end;

`is_leq/simplex_interior` := NULL;

`random_element/simplex_interior` := (A::set) -> proc()
 local x,a,r1,r2,u;
 if A = {} then return FAIL; fi;

 r1 := rand(1..5);
 r2 := rand(1..5);
 x := table();
 u := 0;
 for a in A do
  x[a] := r1()/r2();
  u := u + x[a];
 od;
 for a in A do x[a] := x[a]/u; od;

 return eval(x);
end;

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


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

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

(* type: work_gga_x *)

ad  := 1e-8:
b   := 0.024974:
a2  := (ad + 0.1234)/b:
a4  := 29.790:
a6  := 22.417:
a8  := 12.119:
a10 := 1570.1:
a12 := 55.944:

f0 := s-> 1 + a2*s^2 + a4*s^4 + a6*s^6 + a8*s^8 + a10*s^10 + a12*s^12:
f1 := s-> f0(s)^b/(1 + ad*s^2):

f  := x->f1(X2S*x):


text:
read "../IdentifiabilityODE.mpl";

sys := [
diff(S(t), t) = -b*S(t)*Ninv(t)*A(t)*q - b*S(t)*Ninv(t)*II(t) - b*S(t)*Ninv(t)*J(t),
diff(A(t), t) = -E(t)*k*r + E(t)*k - A(t)*g1,
diff(E(t), t) = b*S(t)*Ninv(t)*A(t)*q + b*S(t)*Ninv(t)*II(t) + b*S(t)*Ninv(t)*J(t) - E(t)*k,
diff(C(t), t) = alpha*II(t),
diff(Ninv(t), t) = 0,
diff(J(t), t) = alpha*II(t) - g2*J(t),
diff(II(t), t) = -alpha*II(t) + E(t)*k*r - g1*II(t),
y2(t) = Ninv(t),
y(t) = C(t)
];
CodeTools[CPUTime](IdentifiabilityODE(sys, GetParameters(sys)));


text:
######################################
# Implementation of the lapse choices
# in standard G-BSSN
# Gamma drive and 1 + log
# Also some experimental other lapse
# choices
######################################

read "/home/arman/FD/FD.mpl":

CFD();
MFD();

grid_functions:= {alpha,K,beta,Lamx,DLamx,BB,ctfmp};

# A generalization of 1 + log
eq1 := diff(alpha(t,x),t) = -gma*alpha(t,x)^gmb*K(t,x);

AVGT := f -> (FD(f,[[1],[0,0,0]]) + FD(f,[[0],[0,0,0]]))/2 ;
FD_table[t] := [ [0] , [0,1] ];

spec := [
  { i = [1,1,1] , b=xmin } = Gen_Sten(lhs(eq1)) - AVGT(Gen_Sten(rhs(eq1))) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq1)) - AVGT(Gen_Sten(rhs(eq1))),
  { i=[Nx,Nx,1], b=xmax  } = alpha(n+1,i) + myzero*x(i) - 1
];

A_Gen_Solve_Code(spec,{alpha(n+1,i)},proc_name="evol_alphaopl");
A_Gen_Res_Code(spec,input="d",proc_name="resid_alphaopl");

# Some kinda generalization of 1 + log
eq2:= diff(alpha(t,x),t)  + epsal*diff(K(t,x),t) =  -epsal*ck*K(t,x);

spec2 := [
  { i = [1,1,1] , b=xmin } = Gen_Sten(lhs(eq2)) - AVGT(Gen_Sten(rhs(eq2))) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq2)) - AVGT(Gen_Sten(rhs(eq2))),
  { i=[Nx,Nx,1], b=xmax  } = alpha(n+1,i) + myzero*x(i) - 1
];

A_Gen_Solve_Code(spec2,{alpha(n+1,i)},proc_name="evol_alphakd");
A_Gen_Res_Code(spec2,input="d",proc_name="resid_alphakd");

# This is gamma drive

eq3 := diff(beta(t,x),t) = mus*Lamx(t,x) - eta*beta(t,x);
eq4 := diff(alpha(t,x),t) = -2*alpha(t,x)*K(t,x);

spec := [
  { i = [1,1,1] , b=xmin } = Gen_Sten(lhs(eq4)) - AVGT(Gen_Sten(rhs(eq4))) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq4)) - AVGT(Gen_Sten(rhs(eq4))),
  { i=[Nx,Nx,1], b=xmax  } = alpha(n+1,i) + myzero*x(i) - 1
];

A_Gen_Solve_Code(spec,{alpha(n+1,i)},proc_name="evol_alphadyopl");
A_Gen_Res_Code(spec,input="d",proc_name="resid_alphadyopl");

spec := [
  { i = [1,1,1] , b=xmin } = Gen_Sten(lhs(eq3)) - AVGT(Gen_Sten(rhs(eq3))) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq3)) - AVGT(Gen_Sten(rhs(eq3))),
  { i=[Nx,Nx,1], b=xmax  } = beta(n+1,i) + myzero*x(i)
];

A_Gen_Solve_Code(spec,{beta(n+1,i)},proc_name="evol_betagd");
A_Gen_Res_Code(spec,input="d",proc_name="resid_betagd");

# Some strange driver, for testing
eq5 := diff(beta(t,x),t) = mus*diff(K(t,x),x) - eta*beta(t,x);

spec := [
  { i = [1,1,1] , b=xmin } = beta(n+1,i) + myzero*x(i) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq5)) - AVGT(Gen_Sten(rhs(eq5))),
  { i=[Nx,Nx,1], b=xmax  } = beta(n+1,i) + myzero*x(i)
];

A_Gen_Solve_Code(spec,{beta(n+1,i)},proc_name="evol_betakd");
A_Gen_Res_Code(spec,input="d",proc_name="resid_betakd");

# More playing around with shift condition
eq6:= diff(beta(t,x),t)  + epsal*diff(Lamx(t,x),t) =  -epsal*ck*Lamx(t,x);

spec2 := [
  { i = [1,1,1] , b=xmin } = Gen_Sten(lhs(eq6)) - AVGT(Gen_Sten(rhs(eq6))) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq6)) - AVGT(Gen_Sten(rhs(eq6))),
  { i=[Nx,Nx,1], b=xmax  } = beta(n+1,i) + myzero*x(i)
];

A_Gen_Solve_Code(spec2,{beta(n+1,i)},proc_name="evol_betagd2");
A_Gen_Res_Code(spec2,input="d",proc_name="resid_betagd2");

eq7 := diff(beta(t,x),t) = advc*beta(t,x)*diff(beta(t,x),x)/ctfmp(x) + mus*BB(t,x);
eq8 := diff(BB(t,x),t) = advc*beta(t,x)*diff(BB(t,x),x)/ctfmp(x) - eta*BB(t,x) + DLamx(t,x);
eq9 := diff(alpha(t,x),t) = advc*beta(t,x)*diff(alpha(t,x),x)/ctfmp(x) - 2*alpha(t,x)*K(t,x);

spec := [
  { i = [1,1,1] , b=xmin } = beta(n+1,i)+myzero*x(i) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq7)) - AVGT(Gen_Sten(rhs(eq7))),
  { i=[Nx,Nx,1], b=xmax  } = beta(n+1,i) + myzero*x(i)
];

A_Gen_Solve_Code(spec,{beta(n+1,i)},proc_name="evol_beta_hgd_adv");
A_Gen_Res_Code(spec,input="d",proc_name="resid_beta_hgd_adv");

spec := [
  { i = [1,1,1] , b=xmin } = BB(n+1,i)+myzero*x(i) ,
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq8)) - AVGT(Gen_Sten(rhs(eq8))),
  { i=[Nx,Nx,1], b=xmax  } = BB(n+1,i) + myzero*x(i)
];

A_Gen_Solve_Code(spec,{BB(n+1,i)},proc_name="evol_BB_hgd_adv");
A_Gen_Res_Code(spec,input="d",proc_name="resid_BB_hgd_adv");

spec := [
  { i = [1,1,1] , b=xmin } =  Gen_Sten(lhs(eq9)) - AVGT(Gen_Sten(rhs(eq9))),
  { i = [2,Nx-1,1]       } = Gen_Sten(lhs(eq9)) - AVGT(Gen_Sten(rhs(eq9))),
  { i=[Nx,Nx,1], b=xmax  } = alpha(n+1,i) - 1 + myzero*x(i)
];

A_Gen_Solve_Code(spec,{alpha(n+1,i)},proc_name="evol_alpha_hgd_adv");
A_Gen_Res_Code(spec,input="d",proc_name="resid_alpha_hgd_adv");


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

a := -0.74860:
b :=  0.06001:
c :=  3.60073:
d :=  0.90000:

f_num := (z, xt) -> sqrt(1 - z^2)*(a + b*xt):
f_den := (rs, xs0, xs1) -> c + d*(xs0 + xs1) + rs:

f := (rs, zeta, xt, xs0, xs1) -> f_num(zeta, xt)/f_den(rs, xs0, xs1):


text:
knights_tour_map := [
 [ 0, 5, 2,57,62,15,18,21],
 [ 3,56,63,14, 9,20,61,16],
 [ 6, 1, 4,43,58,17,22,19],
 [55,44,27, 8,13,10,39,60],
 [28, 7,42,45,40,59,12,23],
 [51,54,49,26,11,36,33,38],
 [48,29,52,41,46,31,24,35],
 [53,50,47,30,25,34,37,32]
];

knights_tour_pos := table():
knights_tour_opos := table():

for i from 0 to 7 do
 for j from 0 to 7 do
  knights_tour_pos[knights_tour_map[i+1,j+1]] := [i,j];
  knights_tour_opos[knights_tour_map[i+1,j+1]] := [i,j] +~ [0.5,0.5];
 od:
od:
knights_tour_pos[64] := knights_tour_pos[0]:
knights_tour_opos[64] := knights_tour_opos[0]:

chess_grid := display(
 seq(line([i,0],[i,8],colour=grey),i=0..8),
 seq(line([0,i],[8,i],colour=grey),i=0..8),
 scaling=constrained,axes=none
):

shaded_chess_grid := display(
 seq(seq(rectangle([2*i+1,2*j+1],[2*i+2,2*j],colour=grey),j=0..3),i=0..3),
 seq(seq(rectangle([2*i,2*j+2],[2*i+1,2*j+1],colour=grey),j=0..3),i=0..3),
 chess_grid
);

knights_tour_plot := display(
 shaded_chess_grid,
 seq(line(knights_tour_opos[i],knights_tour_opos[i+1],colour=red),i=0..63),
 scaling=constrained,axes=none
);

is_knight_move := (u,v) -> evalb({op(map(abs,u -~ v))} = {1,2});

check_knights_tour := proc()
 evalb(sort(map(op,knights_tour_map)) = [seq(i,i=0..63)]) and
 `and`(seq(is_knight_move(knights_tour_pos[i],knights_tour_pos[i+1]),i=0..63));
end: