ADMB Documentation: adrndeff.h (* $Id$) Source File

Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  *
00007  * ADModelbuilder and associated libraries and documentations are
00008  * provided under the general terms of the "BSD" license.
00009  *
00010  * License:
00011  *
00012  * Redistribution and use in source and binary forms, with or without
00013  * modification, are permitted provided that the following conditions are
00014  * met:
00015  *
00016  * 1. Redistributions of source code must retain the above copyright
00017  * notice, this list of conditions and the following disclaimer.
00018  *
00019  * 2.  Redistributions in binary form must reproduce the above copyright
00020  * notice, this list of conditions and the following disclaimer in the
00021  * documentation and/or other materials provided with the distribution.
00022  *
00023  * 3.  Neither the name of the  University of California, Otter Research,
00024  * nor the ADMB Foundation nor the names of its contributors may be used
00025  * to endorse or promote products derived from this software without
00026  * specific prior written permission.
00027  *
00028  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00029  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00030  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00031  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
00032  * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00033  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
00034  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00035  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
00036  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00037  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00038  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00039  *
00040  */
00046 #if !defined(__RANDEFFECTS__)
00047 #  define __RANDEFFECTS__
00048 
00049 #  include <admodel.h>
00050 
00051 class dcompressed_triplet;
00052 class hs_symbolic;
00053 class dvar_compressed_triplet;
00054 
00058 class random_effects_number : public param_init_number
00059 {
00060   virtual void set_random_effects_active();
00061   virtual void set_random_effects_inactive();
00062   virtual void set_only_random_effects_active();
00063   virtual void set_only_random_effects_inactive();
00064 };
00065 
00069 class random_effects_bounded_number : public param_init_bounded_number
00070 {
00071   virtual void set_random_effects_active();
00072   virtual void set_random_effects_inactive();
00073   virtual void set_only_random_effects_active();
00074   virtual void set_only_random_effects_inactive();
00075 };
00076 
00080 class random_effects_vector : public param_init_vector
00081 {
00082   virtual void set_random_effects_active();
00083   virtual void set_random_effects_inactive();
00084   virtual void set_only_random_effects_active();
00085   virtual void set_only_random_effects_inactive();
00086 };
00087 
00088 
00089 class random_effects_bounded_vector;
00090 
00094 class random_effects_bounded_vector : public param_init_bounded_vector
00095 {
00096   virtual void set_random_effects_active();
00097   virtual void set_random_effects_inactive();
00098   virtual void set_only_random_effects_active();
00099   virtual void set_only_random_effects_inactive();
00100 };
00101 
00105 class random_effects_matrix : public param_init_matrix
00106 {
00107   virtual void set_random_effects_active();
00108   virtual void set_random_effects_inactive();
00109   virtual void set_only_random_effects_active();
00110   virtual void set_only_random_effects_inactive();
00111 };
00112 
00116 class random_effects_bounded_matrix : public param_init_bounded_matrix
00117 {
00118   virtual void set_random_effects_active();
00119   virtual void set_random_effects_inactive();
00120   virtual void set_only_random_effects_active();
00121   virtual void set_only_random_effects_inactive();
00122 };
00123 
00124 class gauss_hermite_stuff;
00125 
00126 class nested_calls_shape;
00127 //class sparse_symbolic;
00128 
00132 class nested_calls_indices
00133 {
00134   imatrix * ptr1;
00135   i3_array * ptr2;
00136   i4_array * ptr3;
00137   i5_array * ptr4;
00138 public:
00139   nested_calls_indices(void) : ptr1(0),ptr2(0),ptr3(0),ptr4(0){}
00140   ivector & level1(int i) {return (*ptr1)(i);}
00141   ivector & level2(int i,int j) {return (*ptr2)(i,j);}
00142   ivector & level3(int i,int j,int k) {return (*ptr3)(i,j,k);}
00143   ivector & level4(int i,int j,int k,int l) {return (*ptr4)(i,j,k,l);}
00144   void allocate(const nested_calls_shape& nsc);
00145 };
00146 
00150 class nested_calls_shape
00151 {
00152   ivector * ptr1;
00153   imatrix * ptr2;
00154   i3_array * ptr3;
00155   i4_array * ptr4;
00156 public:
00157   nested_calls_shape(void) : ptr1(0),ptr2(0),ptr3(0),ptr4(0){}
00158   ivector * get_ptr1(void);
00159   imatrix * get_ptr2(void);
00160   i3_array * get_ptr3(void);
00161   i4_array * get_ptr4(void);
00162   void trim(void);
00163   void allocate(int);
00164   void allocate(int,int);
00165   void allocate(int,int,int);
00166   void allocate(int,int,int,int);
00167   void initialize(void);
00168   ~nested_calls_shape();
00169   int & level1(int i) {return (*ptr1)(i);}
00170   int & operator () (int i) {return (*ptr1)(i);}
00171   int & level2(int i,int j) {return (*ptr2)(i,j);}
00172   int & operator () (int i,int j) {return (*ptr2)(i,j);}
00173   int & level3(int i,int j,int k) {return (*ptr3)(i,j,k);}
00174   int & operator () (int i,int j,int k) {return (*ptr3)(i,j,k);}
00175   int & level4(int i,int j,int k,int l) {return (*ptr4)(i,j,k,l);}
00176   int & operator () (int i,int j,int k,int l) {return (*ptr4)(i,j,k,l);}
00177 };
00178 
00182 class laplace_approximation_calculator
00183 {
00184 public:
00185   int init_switch;
00186   nested_calls_indices nested_indices;
00187   nested_calls_shape nested_shape;
00188   int separable_call_level;
00189   int dd_nr_flag;
00190   int no_re_ders_flag;
00191   dmatrix * antiepsilon;
00192   i3_array * triplet_information;
00193   imatrix * compressed_triplet_information;
00194   imatrix * calling_set;
00195   dvector *  importance_sampling_values;
00196   dvector *  importance_sampling_weights;
00197   int is_diagnostics_flag;
00198   static int saddlepointflag;
00199   static int alternative_user_function_flag;
00200   static int sparse_hessian_flag;
00201   static int antiflag;
00202   int rseed;
00203   static int print_importance_sampling_weights_flag;
00204   static dvar_vector * variance_components_vector;
00205   // 1 in inner opt 2 in newton-raphson 3 in laplace approximation
00206   static int where_are_we_flag;
00207   // 1 in inner opt 2 in newton-raphson 3 in
00208   int num_separable_calls;
00209   // 1 in inner opt 2 in newton-raphson 3 in
00210   int separable_calls_counter;
00211   ivector nested_separable_calls_counter;
00212   ivector nested_tree_position;
00213   // 1 in inner opt 2 in newton-raphson 3 in
00214   ivector* num_local_re_array;
00215   // 1 in inner opt 2 in newton-raphson 3 in
00216   ivector* num_local_fixed_array;
00217   int use_outliers;
00218   int use_gauss_hermite;
00219   int in_gauss_hermite_phase;
00220   int multi_random_effects;
00221   int no_function_component_flag;
00222   dvector * separable_function_difference;
00223   gauss_hermite_stuff * gh;
00224   dvar_matrix * importance_sampling_components;
00225   int importance_sampling_counter;
00226   d3_array * block_diagonal_hessian;
00227   d3_array * block_diagonal_ch;
00228   dvar3_array * block_diagonal_vch;
00229   d3_array * block_diagonal_Dux;
00230   imatrix * block_diagonal_re_list;
00231   imatrix * block_diagonal_fe_list;
00232   d3_array * block_diagonal_vhessianadjoint;
00233   dvar3_array * block_diagonal_vhessian;
00234   dvector local_dtemp;
00235   double inner_crit;
00236   double max_separable_g;
00237   double nr_crit;
00238   ivector used_flags;
00239   int have_users_hesstype;
00240   int inner_maxfn;
00241   int nr_debug;
00242   int inner_lmnflag;
00243   int inner_lmnsteps;
00244   int inner_iprint;
00245   int inner_noprintx;
00246   function_minimizer * pmin;
00247   int block_diagonal_flag;
00248   int bw;
00249   int xsize;
00250   int usize;
00251   int nvariables;
00252   unsigned int nvar;
00253   ivector minder;
00254   ivector maxder;
00255   int num_der_blocks;
00256   int num_nr_iters;
00257   int num_importance_samples;
00258   dmatrix epsilon;
00259   int hesstype;
00260   int var_flag;
00261   int have_bounded_random_effects;
00262   int isfunnel_flag;
00263   int nfunnelblocks;
00264   int bHess_pd_flag;
00265   banded_symmetric_dmatrix * bHess;
00266   banded_symmetric_dmatrix * bHessadjoint;
00267   dcompressed_triplet * sparse_triplet;
00268   ivector * sparse_iterator;
00269   int sparse_count;
00270   int sparse_count_adjoint;
00271   dcompressed_triplet * sparse_triplet2;
00272   dvar_compressed_triplet * vsparse_triplet;
00273   dcompressed_triplet * vsparse_triplet_adjoint;
00274   hs_symbolic * sparse_symbolic;
00275   hs_symbolic * sparse_symbolic2;
00276 
00277   void make_sparse_triplet(void);
00278   void check_for_need_to_reallocate(int ip);
00279   void get_newton_raphson_info(function_minimizer * pmin);
00280   void get_newton_raphson_info_slave(function_minimizer * pmin);
00281   void get_newton_raphson_info_master(function_minimizer * pmin);
00282   dvector get_newton_raphson_info_block_diagonal(function_minimizer * pmin);
00283   dvector get_newton_raphson_info_banded(function_minimizer * pmin);
00284   double do_one_feval(const dvector& x,function_minimizer * pfmin);
00285   void test_trust_region_method(function_minimizer * pmin);
00286   void get_complete_hessian(dmatrix& H,function_minimizer * pfmin);
00287   void get_complete_hessian(dmatrix& H,dvector& g,function_minimizer * pfmin);
00288   dvector lincg(dvector& x,dvector& c, dmatrix& H,double tol,double Delta,
00289     function_minimizer * pfmin,double& truef,double& e,double& f,
00290     double& fbest,int& iflag,int& iter,int maxfn);
00291   dvector test_trust_region_method(const dvector& _x,const double& _f,
00292     function_minimizer * pfmin);
00293   dmatrix get_gradient_for_hessian_calcs(const dmatrix& local_Hess,
00294     double & f);
00295   fmm fmc1;
00296   //fmmt1 fmc1;
00297   fmm fmc;
00298   dvector scale;
00299   dvector curv;
00300   dvector xadjoint;
00301   dvector check_local_uadjoint;
00302   dvector check_local_uadjoint2;
00303   dvector check_local_xadjoint;
00304   dvector check_local_xadjoint2;
00305   dvector uadjoint;
00306   dvector uhat;
00307   dvector ubest;
00308   dvector grad;
00309   dvector * grad_x_u;
00310   dvector * grad_x;
00311   dvector step;
00312   dmatrix Hess;
00313   d3_array * Hess_components;
00314   dmatrix * pHess_non_quadprior_part;
00315   imatrix * derindex;
00316   dmatrix Hessadjoint;
00317   dmatrix Dux;
00318   init_df1b2vector y;
00319   dvector get_uhat_quasi_newton_block_diagonal(const dvector& x,
00320     function_minimizer * pfmin);
00321   dvector get_uhat_quasi_newton(const dvector& x,function_minimizer * pfmin);
00322   dvector get_uhat_quasi_newton_qd(const dvector& x,function_minimizer * pfmin);
00323   void set_u_dot(int i);
00324 
00325   void do_separable_stuff_laplace_approximation_banded_adjoint
00326     (const df1b2variable& ff);
00327 
00328   dvector get_uhat_lm_newton2(const dvector& x,function_minimizer * pfmin);
00329 
00330   dvector get_uhat_lm_newton(const dvector& x,function_minimizer * pfmin);
00331   laplace_approximation_calculator(int _xsize,int _usize,int _minder,
00332     int _maxder,function_minimizer * pfmin);
00333 
00334   laplace_approximation_calculator(int _xsize,int _usize,ivector _minder,
00335     ivector _maxder, function_minimizer * pfmin);
00336 
00337   dvector operator () (const dvector& _x,const double& _f,
00338     function_minimizer * pfmin);
00339 
00340   dvector get_uhat(const dvector& x,function_minimizer * pfmin);
00341 
00342   ~laplace_approximation_calculator();
00343 
00344   void  do_separable_stuff(void);
00345   void  do_separable_stuff_newton_raphson_block_diagonal(df1b2variable&);
00346   void  do_separable_stuff_newton_raphson_banded(df1b2variable&);
00347   void  do_separable_stuff_laplace_approximation_block_diagonal(df1b2variable&);
00348   void  do_separable_stuff_laplace_approximation_banded(df1b2variable&);
00349   dvector default_calculations_check_derivatives(const dvector& _x,
00350     function_minimizer * pfmin,const double& f);
00351   dvector default_calculations(const dvector& _x,const double& _f,
00352     function_minimizer * pfmin);
00353   dvector banded_calculations(const dvector& _x,const double& _f,
00354     function_minimizer * pfmin);
00355   dvector block_diagonal_calculations(const dvector& _x,const double& _f,
00356     function_minimizer * pfmin);
00357   dvector default_calculations_parallel_master(const dvector& _x,
00358     const double& _f,function_minimizer * pfmin);
00359   void default_calculations_parallel_slave(const dvector& _x,
00360     const double& _f,function_minimizer * pfmin);
00361   void pvm_slave_function_evaluation_random_effects(void);
00362   dvector banded_calculations_trust_region_approach(const dvector& _uhat,
00363     function_minimizer * pmin);
00364   void do_newton_raphson_banded(function_minimizer * pmin,double,int&);
00365   double inner_optimization_banded(/*dvector& uhat,*/ dvector& x,
00366     function_minimizer * pfmin,int& no_converge_flag);
00367   void set_default_hessian_type(void);
00368   void check_hessian_type(const dvector& _x,function_minimizer * );
00369   void check_hessian_type(function_minimizer * pfmin);
00370   void check_hessian_type2(function_minimizer * pfmin);
00371   void do_separable_stuff_hessian_type_information(void);
00372   void get_block_diagonal_hessian(df1b2variable&);
00373   void allocate_block_diagonal_stuff(void);
00374 
00375   void do_separable_stuff_laplace_approximation_importance_sampling_adjoint
00376     (df1b2variable&);
00377   void get_hessian_components_banded_lme(function_minimizer * pfmin);
00378   dvar_matrix get_hessian_from_components_lme(function_minimizer * pfmin);
00379   dvector banded_calculations_lme(const dvector& _x,const double& _f,
00380     function_minimizer * pfmin);
00381   dvector get_gradient_lme(const dvector& x,function_minimizer * pfmin);
00382   dvector get_gradient_lme(function_minimizer * pfmin);
00383   dvector get_gradient_lme_hp(const double& x,function_minimizer * pfmin);
00384   void check_pool_size(void);
00385   void check_derivatives(const dvector&,function_minimizer * pfmin,
00386     double fval1);
00387   dvector local_minimization_routine(dvector& s,dmatrix& Hess,
00388     dvector& grad,double lambda);
00389   dvector local_minimization(dvector& s,dmatrix& Hess,
00390     dvector& grad,double lambda);
00391   imatrix check_sparse_matrix_structure(void);
00392   double get_fx_fu(function_minimizer * pfmin);
00393   void do_separable_stuff_x_u_block_diagonal(df1b2variable& ff);
00394   void generate_antithetical_rvs();
00395   double standard_type3_loop(int no_converge_flag);
00396   void do_newton_raphson_state_space
00397     (function_minimizer * pfmin,double f_from_1,int& no_converge_flag);
00398   void begin_separable_call_stuff(void);
00399   void end_separable_call_stuff(void);
00400   void build_up_nested_shape(void);
00401 };
00402 
00406 class gauss_hermite_stuff
00407 {
00408 public:
00409   dvar_matrix gauss_hermite_values;
00410   dvector x;
00411   dvector w;
00412   int is;
00413   multi_index * mi;
00414 
00415   gauss_hermite_stuff(laplace_approximation_calculator * lapprox,
00416     int use_gauss_hermite,int num_separable_calls ,const ivector& itmp);
00417 
00418   friend class laplace_approximation_calculator;
00419 };
00420 
00421  ostream& operator <<  (const ostream &,const nested_calls_shape &);
00422 
00423 #endif  //  #if !defined(__RANDEFFECTS__)