1 #ifndef __OPTIMIZATION_H 2 #define __OPTIMIZATION_H 3 #ifdef HAVE_CONFIG_H 4 #include "config.h" 5 #endif 6 #include "first.h" 7 #include "gen.h" 8 #include "unary.h" 9 10 #ifndef NO_NAMESPACE_GIAC 11 namespace giac { 12 #endif // ndef NO_NAMESPACE_GIAC 13 14 enum critical_point_classification { 15 _CPCLASS_UNDECIDED=0, 16 _CPCLASS_MIN=1, 17 _CPCLASS_MAX=2, 18 _CPCLASS_POSSIBLE_MIN=3, 19 _CPCLASS_POSSIBLE_MAX=4, 20 _CPCLASS_SADDLE=5 21 }; 22 23 enum kernel_density_estimation_method { 24 _KDE_METHOD_EXACT, 25 _KDE_METHOD_PIECEWISE, 26 _KDE_METHOD_LIST 27 }; 28 29 enum bandwidth_selection_method { 30 _KDE_BW_METHOD_DPI, 31 _KDE_BW_METHOD_ROT 32 }; 33 34 enum bvp_output_type { 35 _BVP_LIST, 36 _BVP_DIFF, 37 _BVP_PIECEWISE, 38 _BVP_SPLINE 39 }; 40 41 class ipdiff { 42 /* IPDIFF CLASS (Implicit Partial DIFFerentiation) 43 * This class is used for implicit differentiation of f with respect to g=0. 44 * Here, f:R^(n+m)->R, g:R^(n+m)->R^m with rank(g')=m and J(x)=f(x,h(x)) where x in R^n and h is 45 * the implicit function guaranteed to exist (in a neighborhood of some point x0) by the 46 * Implicit Function Theorem if rank(g')=m. Function J may now be differentiated wrt x1,x2,...,xn. */ 47 public: 48 typedef std::vector<int> ivector; 49 typedef ivector::const_iterator ivector_iter; 50 typedef std::vector<ivector> ivectors; 51 typedef std::map<ivector,int> ivector_map; 52 typedef std::pair<ivector,ivector_map> diffterm; 53 typedef std::map<diffterm,int> diffterms; 54 typedef std::map<ivector,gen> pd_map; 55 private: 56 vecteur vars; // list of variables [x1,x2,..,xn,y1,y2,..,ym] where x are free and y are dependent 57 gen f; // the expression f(x1,x2,..,xn,y1,y2,..,ym) 58 vecteur g; // list [g1,g2,..,gm], it is assumed that gj(x1,x2,..,xn,y1,y2,..,ym)=0 for all j=1,2,..,m 59 const context *ctx; // giac context 60 pd_map pdv; // partial derivatives of J 61 pd_map pdf; // partial derivatives of f 62 pd_map pdg; // partial derivatives of g 63 pd_map pdh; // partial derivatives of h 64 std::map<ivector,diffterms> cterms; 65 int ord; // the maximal order of the currently computed partial derivatives of h 66 int nvars; // the number of free variables 67 int nconstr; // the number of constraints = the number of dependent variables 68 diffterms derive_diffterms(const diffterms &terms,ivector &sig); 69 void find_nearest_terms(const ivector &sig,diffterms &match,ivector &excess); 70 const gen &differentiate(const gen &e,pd_map &pds,const ivector &sig); 71 void compute_h(const std::vector<diffterms> &grv,int order); 72 void compute_pd(int order,const ivector &sig=ivector(0)); 73 void raise_order(int order); // compute the missing partial derivatives of h 74 const gen &get_pd(const pd_map &pds,const ivector &sig) const; 75 public: 76 /* construct the ipdiff instance for function f_orig(vars_orig) subject to 77 * the constraints in g_orig, the Jacobian of g_orig must have maximal rank */ 78 ipdiff(const gen &f_orig,const vecteur &g_orig,const vecteur &vars_orig,GIAC_CONTEXT); 79 /* return the partial derivative of J represented by sig */ 80 const gen &derivative(const ivector &sig); 81 /* return the partial derivative of J with respect to the variables in dvars */ 82 const gen &derivative(const vecteur &dvars); 83 /* store all partial derivatives of the specified order to pd_map */ 84 void partial_derivatives(int order,pd_map &pdmap); 85 /* return 1/k!*(sum(i=1)^n (xi-ai)*d/dxi)^k f(x)|a */ 86 gen taylor_term(const vecteur &a,int k,bool shift=true); 87 /* return the Taylor polynomial, with the specified order, of J in the vicinity of a */ 88 gen taylor(const vecteur &a,int order); 89 /* store grad(J) as vector res */ 90 void gradient(vecteur &res); 91 /* store Hessian(J) as matrix res */ 92 void hessian(matrice &res); 93 /* generate all partitions of m into n terms and store them to vector c */ 94 static void ipartition(int m,int n,ivectors &c,const ivector &p=ivector(0)); 95 /* return the sum the elements of a vector v of integers, ignore the last element if drop_last=true */ 96 static int sum_ivector(const ivector &v,bool drop_last=false); 97 }; 98 99 class tprob { 100 /* TPROB CLASS (Transportation PROBlem) 101 * The class implementing the MODI method for balanced TP with degeneracy handling */ 102 public: 103 typedef std::pair<int,int> ipair; 104 typedef std::vector<ipair> ipairs; 105 private: 106 const context *ctx; 107 vecteur demand; 108 vecteur supply; 109 gen eps; // epsilon 110 gen M; // symbol for marking forbidden routes in cost matrix 111 void north_west_corner(matrice &feas); 112 ipairs stepping_stone_path(ipairs &path_orig,const matrice &X); 113 void modi(const matrice &P_orig,matrice &X); 114 public: 115 /* construct the TP with supply s and demand d, m marks forbidden routes */ 116 tprob(const vecteur &s,const vecteur &d,const gen &m,GIAC_CONTEXT); 117 /* solve the transportation problem with the given cost matrix, output in sol */ 118 void solve(const matrice &cost_matrix,matrice &sol); 119 }; 120 121 gen _implicitdiff(const gen &g,GIAC_CONTEXT); 122 gen _minimize(const gen &g,GIAC_CONTEXT); 123 gen _maximize(const gen &g,GIAC_CONTEXT); 124 gen _extrema(const gen &g,GIAC_CONTEXT); 125 gen _minimax(const gen &g,GIAC_CONTEXT); 126 gen _tpsolve(const gen &g,GIAC_CONTEXT); 127 gen _nlpsolve(const gen &g,GIAC_CONTEXT); 128 gen _thiele(const gen &g,GIAC_CONTEXT); 129 gen _triginterp(const gen &g,GIAC_CONTEXT); 130 gen _kernel_density(const gen &g,GIAC_CONTEXT); 131 gen _fitdistr(const gen &g,GIAC_CONTEXT); 132 gen _bvpsolve(const gen &g,GIAC_CONTEXT); 133 gen _euler_lagrange(const gen &g,GIAC_CONTEXT); 134 gen _jacobi_equation(const gen &g,GIAC_CONTEXT); 135 gen _conjugate_equation(const gen &g,GIAC_CONTEXT); 136 gen _convex(const gen &g,GIAC_CONTEXT); 137 138 extern const unary_function_ptr * const at_implicitdiff; 139 extern const unary_function_ptr * const at_minimize; 140 extern const unary_function_ptr * const at_maximize; 141 extern const unary_function_ptr * const at_extrema; 142 extern const unary_function_ptr * const at_minimax; 143 extern const unary_function_ptr * const at_tpsolve; 144 extern const unary_function_ptr * const at_nlpsolve; 145 extern const unary_function_ptr * const at_thiele; 146 extern const unary_function_ptr * const at_triginterp; 147 extern const unary_function_ptr * const at_kernel_density; 148 extern const unary_function_ptr * const at_fitdistr; 149 extern const unary_function_ptr * const at_bvpsolve; 150 extern const unary_function_ptr * const at_euler_lagrange; 151 extern const unary_function_ptr * const at_convex; 152 153 #ifndef NO_NAMESPACE_GIAC 154 } // namespace giac 155 #endif // ndef NO_NAMESPACE_GIAC 156 #endif // __OPTIMIZATION_H 157