ADMB Documentation  11.4.2891
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fun.h
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(__DF1B2FUN__)
00047 #  define __DF1B2FUN__
00048 
00049 #if defined(__GNUC__) && (__GNUC__ < 3)
00050   #pragma interface
00051 #endif
00052 
00053 #include <adpool.h>
00054 #ifndef FVAR_HPP
00055 #  include <fvar.hpp>
00056 #endif
00057 #include <sys/stat.h>
00058 #if !defined(_MSC_VER)
00059   #include <stddef.h>
00060   #include <fcntl.h>
00061 #endif
00062 
00063 #define USE_BARD_PEN
00064 
00065 void ncount_checker(int ncount,int ncount_check);
00066 
00067 class named_dvar_matrix;
00068 
00069 int withinbound(int lb,int n,int ub);
00070 
00071 class do_naught_kludge;
00072 class newadkludge;
00073 
00074 
00079 class twointsandptr
00080 {
00081 public:
00082   short int ncopies;
00083   short int nvar;
00084 #if defined (INCLUDE_BLOCKSIZE)
00085   int blocksize;
00086 #endif
00087   adpool * ptr;
00088 };
00089 
00090 
00091 #if defined(__DERCHECK__)
00092 
00096 class dercheck_info
00097 {
00098   int ind_index;
00099 public:
00100   virtual void get_ind_index(void)
00101   {
00102     cerr << "cannot use this type here" << endl; ad_exit(1);
00103   }
00104 
00105   int node_number;
00106   int counter;
00107   int pass_number;
00108   int vartype;
00109   int dotflag;
00110   double der_value;
00111   double f1;
00112   double f2;
00113   double f3;
00114   double fd;
00115   double delta;
00116   dercheck_info(int _node_number,double _delta,int _index);
00117   void set_delta(double d){delta=d;}
00118   void set_node_number(int n){node_number=n;}
00119   void set_pass_number(int n){pass_number=n;}
00120   void set_index(int n){index=n;}
00121   void set_dotflag(int n){dotflag=n;}
00122   void set_vartype(int n){vartype=n;}
00123   void set_counter(int n){counter=n;}
00124   void set_f1(double d){f1=d;}
00125   void set_f2(double d){f2=d;}
00126   void set_f3(double d){f3=d;}
00127   void set_fd(void){fd=(f2-f3)/(2.0*delta);}
00128 };
00129 extern dercheck_info * derchecker;
00130 #endif
00131 
00132 typedef void * &  vreference;
00133 
00134 /*
00135 #if !defined(_MSC_VER)
00136 inline void increment_pointer(vreference p,int n)
00137 {
00138   char * cs =(char*)(p);
00139   cs+=n;
00140 }
00141 #endif
00142 */
00143 
00144 #include <df32fun.h>
00145 class df1b2_gradlist;
00146 extern df1b2_gradlist* f1b2gradlist;
00147 extern df1b2_gradlist* localf1b2gradlist;
00148 extern df1b2_gradlist* globalf1b2gradlist;
00149 extern df1b2_gradlist** f1b2gradlist_stack;
00150 extern int max_num_init_df1b2variable;
00151 
00152 void df1b2_gradcalc1(void);
00153 
00154 extern char AD_allocation_error_message[];
00155 
00156 #if defined(__BORLANDC__)
00157 int adptr_diff(void * x, void * y) { return int(x)-int(y); }
00158 #else
00159 int adptr_diff(void* x, void* y);
00160 #endif
00161 
00162 void read_pass1_1(void);
00163 void read_pass1_2(void);
00164 
00165 #undef AD_ALLOCATE
00166 #define AD_ALLOCATE(ptr,type,n,classname) \
00167   if ( (ptr = new type[n])==NULL) \
00168   { \
00169     cerr << AD_allocation_error_message << "classname" << endl; \
00170     ad_exit(1); \
00171   }
00172 
00173 
00174 #undef ADUNCONST
00175 #define ADUNCONST(type,obj) type & obj = (type&) _##obj;
00176 
00181 struct df1b2_header
00182 {
00183   //double* ptr;
00184   double* u;
00185   double* u_dot;
00186   double* u_bar;
00187   double* u_dot_bar;
00188   double* u_tilde;
00189   double* u_dot_tilde;
00190   double* u_bar_tilde;
00191   double* u_dot_bar_tilde;
00192 #if defined(__x86_64) || (defined(_MSC_VER) && defined(_M_X64))
00193   long indindex;
00194 #else
00195   int indindex;
00196 #endif
00197 
00198   //double * get_ptr(void){return ptr;}
00199 
00200   double* get_u(void) const {return (double*)u;}
00201   double* get_u_dot(void) const {return (double*)u_dot;}
00202   double* get_u_bar(void) const {return (double*)u_bar;}
00203   double* get_u_dot_bar(void) const {return (double*)u_dot_bar;}
00204   double* get_u_tilde(void) const {return (double*)u_tilde;}
00205   double* get_u_dot_tilde(void) const {return (double*)u_dot_tilde;}
00206   double* get_u_bar_tilde(void) const {return (double*)u_bar_tilde;}
00207   double* get_u_dot_bar_tilde(void) const {return (double*)u_dot_bar_tilde;}
00208 };
00209   class adkludge1;
00210 
00211   class df1b2vector;
00212 
00217   class predf1b2vector
00218   {
00219     df1b2vector * p;
00220     int lb;
00221     int ub;
00222     inline predf1b2vector(df1b2vector * _p,int _lb,int _ub)
00223       {p=_p;lb=_lb,ub=_ub;}
00224     friend class df1b2vector;
00225   };
00226 
00227   class df3_one_variable;
00228   class df3_two_variable;
00229   class df3_three_variable;
00230   class random_effects_bounded_vector_info;
00231 
00236   class df1b2variable : public df1b2_header
00237   {
00238   public:
00239     double* ptr;
00240     int get_local_nvar(void) const {return int(((twointsandptr*)ptr)->nvar);}
00241     int allocated(void){return ptr!=0;}
00242     int unallocated(void){return ptr==0;}
00243     static adpool* pool;
00244     static int current_allocation_index;
00245     static int adpool_use_index[];
00246     static int pool_allocation_number[];
00247     static int allocation_counter;
00248     static const int adpool_vectorsize;
00249     static int adpool_counter;
00250     static void increment_adpool_counter(void);
00251     static const int adpool_stack_size;
00252     static adpool* adpool_vector[];
00253     static adpool* adpool_stack[];
00254     static unsigned int adpool_nvar_stack[];
00255     static int adpool_stack_pointer;
00256     static void save_adpool_pointer(void);
00257     static void restore_adpool_pointer(void);
00258     static unsigned int nvar_vector[];
00259     static int passnumber;
00260     static unsigned int nvar;  // the number of independent variables
00261     static int minder;  // the number of independent variables
00262     static int maxder;  // the number of independent variables
00263     static unsigned int blocksize;
00264     static int noallocate;
00265 
00266     static int get_passnumber(void){return passnumber;}
00267     static void set_nvar(unsigned int n) { nvar = n;}
00268     static unsigned int get_nvar() { return nvar; }
00269     static void set_minder(int n){minder=n;}
00270     static void set_maxder(int n){maxder=n;}
00271     static void set_blocksize(void);
00272     static unsigned int get_blocksize(void);
00273     static unsigned int get_blocksize(const unsigned int n);
00274 #if defined(__x86_64) || (defined(_MSC_VER) && defined(_M_X64))
00275     long& get_ind_index() { return indindex; }
00276     const long& get_ind_index() const { return indindex; }
00277 #else
00278     int& get_ind_index() { return indindex; }
00279     const int& get_ind_index() const { return indindex; }
00280 #endif
00281     short int* ncopies;
00282     // for fixed size n whole thing is 6n+2
00283     void initialize(void);
00284     void initialize(const unsigned int n);
00285 
00286     df1b2variable(const do_naught_kludge&){ ptr = 0; }
00287 
00288 #if defined(USE_DDOUBLE)
00289 #undef double
00290     df1b2variable(double d);
00291 #define double dd_real
00292 #endif
00293     df1b2variable(double d);
00294     df1b2variable(void);
00295     void allocate(void);
00296     void allocate(const char *);
00297 
00298     df1b2variable(adkludge1 * );
00299     df1b2variable(const random_effects_bounded_vector_info& rebv);
00300 
00301     df1b2variable(const newadkludge * );
00302 
00303     df1b2variable(const df1b2variable& v);
00304 
00305     df1b2variable& operator = (const df3_one_variable& v);
00306     df1b2variable& operator = (const df3_two_variable& v);
00307     df1b2variable& operator = (const df3_three_variable& v);
00308     df1b2variable& operator = (const df1b2variable& v);
00309     df1b2variable& operator += (const df1b2variable& v);
00310     df1b2variable& operator -= (const df1b2variable& v);
00311     df1b2variable& operator *= (const df1b2variable& v);
00312     df1b2variable& operator /= (const df1b2variable& v);
00313 
00314     void operator = (double x);
00315 
00316     virtual ~df1b2variable();
00317     void deallocate(void);
00318   };
00319 
00320 void print_dot_derivatives(df1b2_header * px,const char * s);
00321 void print_derivatives(const adstring&, double f, double df,
00322   double d2f,double d3f,int bflag=0);
00323 void print_derivatives(const adstring&, double f, double df1,
00324   double df2,double df11,double df12, double df22,
00325   double df111, double df112, double df122, double df222,int bflag=0);
00326 
00327 void print_derivatives(df1b2_header * px,const char * s,
00328   int bflag=0);
00329 
00334   class init_df1b2variable : public df1b2variable
00335   {
00336   public:
00337     int ind_index;
00338     int get_index(void){return ind_index;}
00339     static int num_variables;
00340     static init_df1b2variable ** list;
00341     init_df1b2variable(double v=0.0);
00342     void operator = (double x);
00343     void set_u_dot(void);
00344   };
00345 
00350   class init_df1b2vector
00351   {
00352     int index_min;
00353     int index_max;
00354     int * ncopies;
00355     init_df1b2variable * trueptr;
00356     init_df1b2variable * ptr;
00357   public:
00358     void deallocate(void);
00359     void reallocate(void);
00360     void allocate(int lib,int ub);
00361     void allocate(void);
00362     init_df1b2vector(int lib,int ub);
00363     init_df1b2vector(void);
00364     init_df1b2vector(const init_df1b2vector& v);
00365     int indexmin(void) const {return index_min;}
00366     int indexmax(void) const {return index_max;}
00367 
00368 #if defined(OPT_LIB)
00369     init_df1b2variable& operator () (int i) { return ptr[i];}
00370     init_df1b2variable& operator [] (int i) { return ptr[i];}
00371 #else
00372     init_df1b2variable& operator () (int i);
00373     init_df1b2variable& operator [] (int i);
00374 #endif
00375     void set_value(const dvector&);
00376     ~init_df1b2vector();
00377   };
00378 
00384 class df1b2function1
00385 {
00386 public:
00387   double (*f)(double);
00388   double (*df)(double);
00389   double (*d2f)(double);
00390   double (*d3f)(double);
00391   adstring funname;
00392 
00393   df1b2function1(double (*_f)(double),double (*_df)(double),
00394       double (*d2f)(double),double (*_d3f)(double),
00395       const adstring& s="unnamed");
00396 
00397   df1b2variable operator () (const df1b2variable& x);
00398   //df1b2variable& operator () (const df1b2variable& z,const df1b2variable& x);
00399   df1b2variable& operator () (const df1b2variable& z,const df1b2variable& x,
00400     int s); // s is a switch for picking special function for simple
00401                // operations like +=
00402 };
00403 
00408   class df1b2function2
00409   {
00410   public:
00411     double (*f)(double,double);
00412     double (*df1)(double,double);
00413     double (*df2)(double,double);
00414     double (*d2f11)(double,double);
00415     double (*d2f12)(double,double);
00416     double (*d2f22)(double,double);
00417     double (*d3f111)(double,double);
00418     double (*d3f112)(double,double);
00419     double (*d3f122)(double,double);
00420     double (*d3f222)(double,double);
00421     adstring funname;
00422 
00423     df1b2function2
00424     (
00425       double (*_f)(double,double),
00426       double (*_df1)(double,double),double (*_df2)(double,double),
00427       double (*d2f11)(double,double),
00428       double (*d2f12)(double,double),
00429       double (*d2f22)(double,double),
00430       double (*_d3f111)(double,double),
00431       double (*_d3f112)(double,double),
00432       double (*_d3f122)(double,double),
00433       double (*_d3f222)(double,double),
00434       const adstring & funame="unnamed"
00435     );
00436 
00437     df1b2variable operator () (const df1b2variable& x,const df1b2variable& y);
00438   };
00439 
00444   class mydf1b2function2
00445   {
00446   public:
00447     double (*f)(double,double);
00448     double (*df1)(double,double);
00449     double (*df2)(double,double);
00450     double (*d2f11)(double,double);
00451     double (*d2f12)(double,double);
00452     double (*d2f22)(double,double);
00453     double (*d3f111)(double,double);
00454     double (*d3f112)(double,double);
00455     double (*d3f122)(double,double);
00456     double (*d3f222)(double,double);
00457 
00458     mydf1b2function2
00459     (
00460       double (*_f)(double,double),
00461       double (*_df1)(double,double),double (*_df2)(double,double),
00462       double (*d2f11)(double,double),
00463       double (*d2f12)(double,double),
00464       double (*d2f22)(double,double),
00465       double (*_d3f111)(double,double),
00466       double (*_d3f112)(double,double),
00467       double (*_d3f122)(double,double),
00468       double (*_d3f222)(double,double)
00469     );
00470 
00471     df1b2variable operator () (const df1b2variable& x,const df1b2variable& y);
00472   };
00473 
00478 class smartlist
00479 {
00480 public:
00481   char* buffer;
00482   char* buffend;
00483   char* bptr;
00484   char* sbptr;
00485   unsigned int bufsize;
00486   adstring filename;
00487   int fp;
00488   void saveposition(void){sbptr=bptr;}
00489   void reset(void){bptr=buffer;}
00490   void restoreposition(void){bptr=sbptr;}
00491   void restoreposition(int offset){bptr=sbptr+offset;}
00492   smartlist(unsigned int bufsize,const adstring& filename);
00493   double report_usage(void)
00494   {
00495     return double(size_t(bptr)-size_t(buffer))
00496            / double(size_t(buffend)-size_t(buffer));
00497   }
00498 };
00499 
00504   class test_smartlist
00505   {
00506   public:
00507     int direction;
00508     int written_flag;
00509     int noreadflag;
00510     void save_end(void);
00511     void restore_end(void);
00512     int eof_flag;
00513     int end_saved;
00514     char * true_buffer;
00515     char * true_buffend;
00516     char * recend;
00517     char * buffer;
00518     char * buffend;
00519     char * bptr;
00520     char * sbptr;
00521     size_t bufsize;
00522     adstring filename;
00523     int fp;
00524 public:
00525     test_smartlist(void);
00526     ~test_smartlist();
00527 
00528     void saveposition(void){sbptr=bptr;}
00529     void set_recend(){bptr=recend+1;} // one passed the end so that when
00530                                       // you back up n bytes it points to the
00531                                       // beginning of an n byte record
00532     void reset(void);
00533     int eof(void){ return eof_flag;}
00534     int get_noreadflag(void){ return noreadflag; }
00535     void set_noreadflag(int n){ noreadflag=n; }
00536     void restoreposition(void){bptr=sbptr;}
00537     void restoreposition(int offset){bptr=sbptr+offset;}
00538     test_smartlist(const size_t bufsize,const adstring& filename);
00539     void allocate(const size_t bufsize,const adstring& filename);
00540     void read_buffer(void);
00541     void set_forward(void){direction=0;}
00542     void set_backword(void){direction=-1;}
00543     void set_reverse(void){direction=-1;}
00544     void rewind(void);
00545     void initialize(void);
00546     void operator-=(const int);
00547     void operator+=(const int);
00548     double report_usage(void)
00549     {
00550       return double(size_t(bptr)-size_t(buffer))
00551              / double(size_t(buffend)-size_t(buffer));
00552     }
00553     //void write(void * p,const size_t n);
00554     void write(const size_t n);
00555     void write_buffer(void);
00556     void check_buffer_size(const size_t);
00557     void add_buffer_fringe(int n){buffend-=n;}
00558     int written(void){return written_flag;}
00559   };
00560 
00561   typedef void (*ADrfptr)(void);
00562 
00567   class fixed_list_entry
00568   {
00569   public:
00570     int numbytes;
00571     ADrfptr pf;
00572   };
00573 
00578   class fixed_smartlist
00579   {
00580   public:
00581     fixed_smartlist();
00582     fixed_smartlist(const size_t bufsize,const adstring& filename);
00583     ~fixed_smartlist();
00584 
00585     void operator ++ (void);
00586     void operator -- (void);
00587     size_t nentries;
00588     int direction;
00589     off_t endof_file_ptr;
00590     int written_flag;
00591     int noreadflag;
00592     void save_end(void);
00593     void restore_end(void);
00594     int eof_flag;
00595     int end_saved;
00596     fixed_list_entry * true_buffer;
00597     fixed_list_entry * true_buffend;
00598     fixed_list_entry * recend;
00599     fixed_list_entry * buffer;
00600     fixed_list_entry * buffend;
00601     fixed_list_entry * bptr;
00602     fixed_list_entry * sbptr;
00603     size_t bufsize;
00604     adstring filename;
00605     int fp;
00606     void saveposition(void){sbptr=bptr;}
00607     void set_recend(){bptr=recend+1;} // one passed the end so that when
00608                                       // you back up n bytes it points to the
00609                                       // beginning of an n byte record
00610     void reset(void);
00611     int eof(void){ return eof_flag;}
00612     void read_file(void);
00613     int get_noreadflag(void){ return noreadflag; }
00614     void set_noreadflag(int n){ noreadflag=n; }
00615     void restoreposition(void){bptr=sbptr;}
00616     void restoreposition(int offset){bptr=sbptr+offset;}
00617     void allocate(const size_t bufsize,const adstring& filename);
00618     void read_buffer(void);
00619     void set_forward(void){direction=0;}
00620     void set_backword(void){direction=-1;}
00621     void set_reverse(void){direction=-1;}
00622     void rewind(void);
00623     void initialize(void);
00624     void operator -= (int);
00625     void operator += (int);
00626     double report_usage(void)
00627     {
00628       return double(size_t(bptr)-size_t(buffer))
00629              / double(size_t(buffend)-size_t(buffer));
00630     }
00631     //void write(void * p,int n);
00632     void write(const size_t);
00633     void write_buffer(void);
00634     void write_buffer_one_less(void);
00635     void check_buffer_size(const size_t);
00636     void add_buffer_fringe(int n){buffend-=n;}
00637     int written(void){return written_flag;}
00638   };
00639 
00644   class fixed_smartlist2
00645   {
00646   public:
00647     fixed_smartlist2(void);
00648     fixed_smartlist2(const size_t bufsize,const adstring& filename);
00649     ~fixed_smartlist2();
00650 
00651     void operator ++ (void);
00652     void operator -- (void);
00653     size_t nentries;
00654     int direction;
00655     off_t endof_file_ptr;
00656     int written_flag;
00657     int noreadflag;
00658     void save_end(void);
00659     void restore_end(void);
00660     int eof_flag;
00661     int end_saved;
00662     int * true_buffer;
00663     int * true_buffend;
00664     int * recend;
00665     int * buffer;
00666     int * buffend;
00667     int * bptr;
00668     int * sbptr;
00669     size_t bufsize;
00670     adstring filename;
00671     int fp;
00672     void saveposition(void){sbptr=bptr;}
00673     void set_recend(){bptr=recend+1;} // one passed the end so that when
00674                                       // you back up n bytes it points to the
00675                                       // beginning of an n byte record
00676     void reset(void){bptr=buffer;}
00677     int eof(void){ return eof_flag; /*eof_flag=0;*/}
00678 
00679     void read_file(void);
00680     int get_noreadflag(void){ return noreadflag; }
00681     void set_noreadflag(int n){ noreadflag=n; }
00682     void restoreposition(void){bptr=sbptr;}
00683     void restoreposition(int offset){bptr=sbptr+offset;}
00684     void allocate(const size_t bufsize,const adstring& filename);
00685     void read_buffer(void);
00686     void set_forward(void){direction=0;}
00687     void set_backword(void){direction=-1;}
00688     void set_reverse(void){direction=-1;}
00689     void rewind(void);
00690     void initialize(void);
00691     void operator -= (int);
00692     void operator += (int);
00693     double report_usage(void)
00694     {
00695       return double(size_t(bptr)-size_t(buffer)) /
00696              double(size_t(buffend)-size_t(buffer));
00697     }
00698     //void write(void * p,int n);
00699     void write(const size_t n);
00700     void write_buffer(void);
00701     void write_buffer_one_less(void);
00702     void check_buffer_size(const size_t);
00703     void add_buffer_fringe(int n){buffend-=n;}
00704     int written(void){return written_flag;}
00705   };
00706 
00707 
00708   void write(const test_smartlist &,void *,int nsize);
00709   void read(const test_smartlist &,void *,int nsize);
00710   void memcpy(test_smartlist&, void*, const size_t nsize);
00711   void memcpy(void*, test_smartlist&, const size_t nsize);
00712 
00713   class df1b2function2c;
00714 
00719   class df1b2_gradlist
00720   {
00721   public:
00722 #if defined(CHECK_COUNT)
00723     static int ncount_check;
00724     static set_ncount_check(int n){ncount_check=n;}
00725 #endif
00726     int ncount;
00727     test_smartlist list;
00728     fixed_smartlist nlist;
00729     test_smartlist list2;
00730     fixed_smartlist2 nlist2;
00731     test_smartlist list3;
00732     fixed_smartlist2 nlist3;
00733     static int no_derivatives;
00734     static void set_no_derivatives(void) {no_derivatives=1;}
00735     static void set_yes_derivatives(void) {no_derivatives=0;}
00736     df1b2_gradlist(unsigned int bufsize,unsigned int nbufsize,
00737      unsigned int bufsize1, unsigned int nbufsize1,
00738      unsigned int bufsize2, unsigned int nbufsize2,const adstring& filename);
00739     int mywrite_pass1(const df1b2variable * px, const df1b2variable * py,
00740       df1b2variable * pz, mydf1b2function2 * pf);
00741     int write_pass1_sum(double x, const df1b2variable * py,
00742       df1b2variable * pz);
00743     int write_pass1_sum(const df1b2variable * px, const df1b2variable * py,
00744       df1b2variable * pz);
00745     int write_pass1_prod(double x, const df1b2variable * py,
00746       df1b2variable * pz);
00747     int write_pass1_prod(const df1b2variable * px, double py,
00748       df1b2variable * pz);
00749     int write_pass1x(const df1b2variable * _px, df1b2variable * pz,
00750       df1b2function1 * pf);
00751 
00752     int write_pass1_prod(const df1b2vector * px, const df1b2vector * py,
00753       df1b2variable * pz);
00754 
00755     int write_pass1_prod(const df1b2variable * px, const df1b2variable * py,
00756       df1b2variable * pz);
00757     int write_pass1_minuscv(const df1b2variable * py,df1b2variable * pz);
00758     int write_pass1_minusvc(const df1b2variable * py,df1b2variable * pz);
00759     int write_pass1_minus(const df1b2variable * px, const df1b2variable * py,
00760       df1b2variable * pz);
00761     int write_pass1c(const df1b2variable * px, double y, df1b2variable * pz,
00762       df1b2function2c * pf);
00763     int write_pass1c(double x,const df1b2variable * py, df1b2variable * pz,
00764       df1b2function2c * pf);
00765     int write_pass1(const df1b2variable * px, const df1b2variable * py,
00766       df1b2variable * pz, df1b2function2 * pf);
00767     int write_pass1_pluseq(const df1b2variable * px,df1b2variable * pz);
00768     int write_pass1_minuseq(const df1b2variable * px,df1b2variable * pz);
00769     int write_pass1_initialize(df1b2variable * pz);
00770     int write_pass1_eq(const df1b2variable * px,df1b2variable * pz);
00771     int write_pass1(const df1b2variable * px,
00772       df1b2variable * pz, df1b2function1 * pf);
00773     int write_pass1(const df1b2variable * px,
00774       df1b2variable * pz, double df, double d2f, double d3f );
00775 
00776     int write_pass1(const df1b2variable * _px,
00777       const df1b2variable * _py,df1b2variable * pz,double df_x,
00778       double df_y,
00779       double df_xx,
00780       double df_xy,
00781       double df_yy,
00782       double df_xxx,
00783       double df_xxy,
00784       double df_xyy,
00785       double df_yyy);
00786 
00787  int write_pass1(const df1b2variable * _px,
00788    const df1b2variable * _py,const df1b2variable * pw,
00789    const df1b2variable * pz,
00790    double df_x, double df_y, double df_z,
00791    double df_xx, double df_xy,double df_xz,double df_yy,
00792    double df_yz,double df_zz,
00793    double df_xxx,
00794    double df_xxy,
00795    double df_xxz,
00796    double df_xyy,
00797    double df_xyz,
00798    double df_xzz,
00799    double df_yyy,
00800    double df_yyz,
00801    double df_yzz,
00802    double df_zzz);
00803 
00804 
00805 
00806     int write_save_pass2_tilde_values(const df1b2variable * px);
00807     void reset(void);
00808   };
00809 
00814 class ad_dstar
00815 {
00816 public:
00817   static unsigned int n;
00818   static void allocate(const unsigned int _n);
00819   double* p;
00820   ad_dstar(void);
00821   ~ad_dstar(void);
00822   operator double* (){return p;}
00823 };
00824 
00825 typedef void (**ADprfptr)(void);
00826 typedef void (*ADrfptr)(void);
00827 
00828 void set_dependent_variable(const df1b2variable& _x);
00829 
00830 dmatrix get_hessian(const init_df1b2vector& x);
00831 //df1b2variable user_function(const init_df1b2vector& x);
00832 
00833 dmatrix check_second_derivatives(const init_df1b2vector& x);
00834 d3_array check_third_derivatives(const init_df1b2vector& x);
00835 
00836 // ************************************************************
00837 // ************************************************************
00838 
00839 df1b2variable mfexp(const df1b2variable& x);
00840 df1b2variable mfexp(const df1b2variable& x,double b);
00841 
00842 // ************************************************************
00843 // ************************************************************
00844 
00845 df1b2variable fabs(const df1b2variable& x);
00846 df1b2variable max(const df1b2vector& t1);
00847 df1b2variable sfabs(const df1b2variable& x);
00848 df1b2variable pow(const df1b2variable& x,const df1b2variable& y);
00849 df1b2variable pow(const df1b2variable& x,double y);
00850 df1b2variable mypow(const df1b2variable& x,double y);
00851 df1b2variable pow(double x,const df1b2variable& y);
00852 df1b2variable sin(const df1b2variable& x);
00853 df1b2variable atan(const df1b2variable& x);
00854 df1b2variable tan(const df1b2variable& x);
00855 df1b2variable inv(const df1b2variable& x);
00856 df1b2variable sqrt(const df1b2variable& x);
00857 df1b2variable cos(const df1b2variable& x);
00858 df1b2variable exp(const df1b2variable& x);
00859 df1b2variable log(const df1b2variable& x);
00860 df1b2variable square(const df1b2variable& x);
00861 df1b2variable cube(const df1b2variable& x);
00862 df1b2variable fourth(const df1b2variable& x);
00863 df1b2variable operator * (const df1b2variable& x,double y);
00864 df1b2variable operator * (const df1b2variable& x,const df1b2variable& y);
00865 df1b2variable operator * (double x,const df1b2variable& y);
00866 df1b2variable mult_add(const df1b2variable& x,const df1b2variable& y);
00867 df1b2variable operator + (const df1b2variable& x,const df1b2variable& y);
00868 df1b2variable div(const df1b2variable& x,const df1b2variable& y);
00869 df1b2variable operator + (double x,const df1b2variable& y);
00870 df1b2variable operator + (const df1b2variable& x,double y);
00871 
00876 inline df1b2variable operator + (const df1b2variable& x,double y)
00877 {
00878   return y+x;
00879 }
00880 
00881 df1b2variable operator - (const df1b2variable& x,const df1b2variable& y);
00882 df1b2variable operator - (double x,const df1b2variable& y);
00883 df1b2variable operator / (const df1b2variable& x,const df1b2variable& y);
00884 df1b2variable operator / (const df1b2variable& x,double y);
00885 df1b2variable operator - (const df1b2variable& x,double y);
00886 df1b2variable operator / (double x,const df1b2variable& y);
00887 
00888 df1b2variable lgamma2(const df1b2variable& _x);//new log gamma using forward AD
00889 df1b2variable gammln(const df1b2variable& _xx);
00890 df1b2vector gammln(const df1b2vector&  _xx);
00891 df1b2variable log_comb(const df1b2variable& n, double k);
00892 df1b2variable log_comb(const df1b2variable& n, const df1b2variable& k);
00893 df1b2variable log_comb(double n, const df1b2variable& k);
00894 df1b2variable factln(const df1b2variable& n);
00895 
00896 // vector and matrix stuff
00897 
00902 class df1b2vector_header
00903 {
00904 protected:
00905   int offset;
00906   int index_min;
00907   int index_max;
00908 };
00909 
00910 class df3_one_vector;
00911 class funnel_init_df1b2vector;
00912 class df3_two_vector;
00913 
00914  // class predf1b2vector
00915  // {
00916  //   df1b2vector * p;
00917  //   int lb;
00918  //   int ub;
00919  //inline predf1b2vector(df1b2vector * _p,int _lb,int _ub) {p=_p;lb=_lb,ub=_ub;}
00920  //   friend class df1b2vector;
00921  // };
00922 
00927 class df1b2vector : public df1b2vector_header
00928 {
00929   df1b2variable * v;
00930   vector_shapex * shape;
00931 public:
00932   inline df1b2vector& operator -- (void)
00933   {
00934     index_min--;index_max--;v++; return *this;
00935   }
00936   inline df1b2vector& operator ++ (void)
00937   {
00938     index_min++;index_max++;v--; return *this;
00939   }
00940   int pointersize() const { return (int)sizeof(df1b2variable); }
00941   inline df1b2variable * getv(void) {return v;}
00942   int allocated(void){return v!=0;}
00943   int indexmin(void)const {return index_min;}
00944   int indexmax(void)const {return index_max;}
00945   int size(){return index_max-index_min+1;}
00946   df1b2vector(int lb,int ub);
00947   df1b2vector(const dvector& v);
00948   ~df1b2vector();
00949   df1b2vector(const df1b2vector&);
00950   void copy(const df1b2vector&);
00951   df1b2vector(void);
00952   void allocate(int lb,int ub,const char *);
00953   void allocate(int lb,int ub);
00954   void allocate(const ad_integer&,const ad_integer& ub);
00955   void noallocate(int lib,int ub);
00956   df1b2vector& shift(int);
00957   void allocate(void);
00958   void initialize(void);
00959   void deallocate(void);
00960   int get_offset(void){return offset;}
00961   df1b2vector& operator = (const df3_one_vector&);
00962   df1b2vector& operator = (const df1b2vector&);
00963   df1b2vector& operator = (const dvector&);
00964   df1b2vector& operator = (double);
00965   df1b2vector& operator = (const df1b2variable&);
00966   df1b2vector& operator = (const df3_two_vector& M);
00967 #  if defined(OPT_LIB)
00968   df1b2variable& operator () (int i) const
00969   {
00970    return *((df1b2variable*)((char*)(v)+i*pointersize()));
00971   }
00972   df1b2variable& operator [] (int i) const
00973   {
00974    return *((df1b2variable*)((char*)(v)+i*pointersize()));
00975   }
00976  //  df1b2vector operator () (int i,int j)
00977  //  {
00978  //    return predf1b2vector(this,i,j);
00979  //  }
00980  // df1b2vector(const predf1b2vector&);
00981 #  else
00982   //df1b2vector operator () (int i,int j);
00983   //df1b2variable& operator () const (int i);
00984   df1b2variable& operator () (int i) const;
00985   df1b2variable& operator [] (int i) const;
00986 #  endif
00987   df1b2vector operator() (const ivector & iv);
00988   df1b2vector& operator += (const df1b2vector& x);
00989   df1b2vector& operator += (const dvector& x);
00990   df1b2vector& operator += (double x);
00991   df1b2vector& operator -= (const df1b2vector& x);
00992   df1b2vector& operator -= (const dvector& x);
00993   df1b2vector& operator /= (const df1b2vector& x);
00994   df1b2vector& operator *= (const df1b2vector& x);
00995   df1b2vector& operator *= (const df1b2variable& x);
00996   df1b2vector& operator -= (const df1b2variable& x);
00997   df1b2vector& operator += (const df1b2variable& x);
00998   df1b2vector& operator /= (const df1b2variable& x);
00999   df1b2vector(const predf1b2vector&);
01000   inline df1b2vector operator () (int lb,int ub)
01001   {
01002     return predf1b2vector(this,lb,ub);
01003   }
01004   friend class random_effects_bounded_vector_info;
01005   friend class df1b2variable;
01006   friend class xf_df1b2vector;
01007   //df1b2vector(const funnel_init_df1b2vector& _s);
01008 };
01009 class df3_one_matrix;
01010 class df3_two_matrix;
01011 
01016 class df1b2matrix
01017 {
01018   int index_min;
01019   int index_max;
01020   df1b2vector * v;
01021   mat_shapex * shape;
01022 public:
01023   int allocated(void){return v!=0;}
01024   void initialize(void);
01025   ~df1b2matrix();
01026   void colfill(const int j, const df1b2vector& v);
01027   int rowmin(void) const {return index_min;}
01028   int indexmin(void) const {return index_min;}
01029   int indexmax(void) const {return index_max;}
01030   int rowmax(void) const {return index_max;}
01031   int size(void) const {return index_max-index_min+1;}
01032   df1b2matrix(int nrl,int nrh,int ncl,int nch);
01033   df1b2matrix(int nrl,int nrh);
01034   df1b2matrix(const df1b2matrix&);
01035   df1b2matrix(int nrl,int nrh,const index_type& ncl,
01036     const index_type& nch);
01037   df1b2matrix& operator =(const df3_one_matrix&);
01038 
01039   df1b2matrix& operator =(const df1b2matrix&);
01040   df1b2matrix& operator =(const dmatrix&);
01041   df1b2matrix& operator =(const df1b2variable&);
01042   df1b2matrix& operator =(double);
01043   df1b2matrix& operator = (const df3_two_matrix& M);
01044   df1b2matrix(void);
01045   df1b2matrix& operator +=(const df1b2matrix& M);
01046   df1b2matrix& operator -=(const df1b2matrix& M);
01047   void allocate(int nrl,int nrh,int ncl,int nch);
01048   void allocate(int nrl,int nrh);
01049   void allocate(int nrl,int nrh,int ncl,int nch,const char *);
01050   void allocate(int nrl,int nrh,const index_type& ncl,
01051     const index_type& nch);
01052   void allocate(int nrl,int nrh,const index_type& ncl,
01053     const index_type& nch,const char *);
01054   void allocate(void);
01055   void deallocate(void);
01056   df1b2variable& operator()(int i,int j) const;
01057   df1b2vector& operator[](int i) const;
01058   df1b2vector& operator()(int i) const;
01059   int colmin(void) const {return (*(df1b2matrix*)(this))(index_min).indexmin();}
01060   int colmax(void) const {return (*(df1b2matrix*)(this))(index_min).indexmax();}
01061   int colsize(void)const {return colmax()-colmin()+1;}
01062   df1b2matrix& operator *= (const df1b2variable& x);
01063   df1b2matrix& operator /= (const df1b2variable& x);
01064   df1b2matrix& operator *= (double x);
01065   df1b2matrix& operator /= (double x);
01066 };
01067 #if defined(OPT_LIB)
01068 inline df1b2variable& df1b2matrix::operator()(int i,int j) const
01069 {
01070   return (df1b2variable&)(v[i][j]);
01071 }
01072 inline df1b2vector& df1b2matrix::operator[](int i) const
01073 {
01074   return (df1b2vector&)(v[i]);
01075 }
01076 inline df1b2vector& df1b2matrix::operator()(int i) const
01077 {
01078   if (!v)
01079     throw std::bad_alloc();
01080   else
01081     return v[i];
01082 }
01083 #endif
01084 
01089 class df1b23array
01090 {
01091   int index_min;
01092   int index_max;
01093   df1b2matrix * v;
01094   vector_shapex * shape;
01095 public:
01096   int allocated(void){return v!=0;}
01097   void initialize(void);
01098   ~df1b23array();
01099   int indexmin(void){return index_min;}
01100   int indexmax(void){return index_max;}
01101   int size(void){return index_max-index_min+1;}
01102   df1b23array(int nrl,int nrh,int ncl,int nch,int,int);
01103   df1b23array(int nrl,int nrh);
01104   df1b23array(int nrl,int nrh,int,int);
01105   df1b23array(const df1b23array&);
01106   //df1b23array& operator =(const df3_one_matrix&);
01107 
01108   df1b23array& operator =(const df1b23array&);
01109   df1b23array& operator =(const df1b2variable&);
01110   df1b23array& operator =(double);
01111   df1b23array(void);
01112   df1b23array& operator +=(const df1b23array& M);
01113   df1b23array& operator -=(const df1b23array& M);
01114   void allocate(int nrl,int nrh,int ncl,int nch,int,int);
01115   void allocate(int nrl,int nrh,int,int);
01116   void allocate(int nrl,int nrh);
01117   void allocate(int nrl,int nrh,int ncl,int nch,int,int,const char *);
01118 
01119   void allocate(int nrl,int nrh,
01120     const index_type& ncl, const index_type& nch);
01121 
01122 
01123   void allocate(int nrl,int nrh,
01124     const index_type& ncl, const index_type& nch,
01125     const index_type& nxl, const index_type& nxh);
01126 
01127   void allocate(int nrl,int nrh,
01128     const index_type& ncl,const index_type& nch,
01129     const index_type& nxl,const index_type& nxh,
01130     const char *);
01131   void allocate(void);
01132   void deallocate(void);
01133 #  if defined(OPT_LIB)
01134   df1b2variable& operator () (int i,int j,int k){return v[i](j,k);}
01135   df1b2vector& operator () (int i,int j){return v[i][j];}
01136   df1b2matrix& operator () (int i){return v[i];}
01137   df1b2matrix& operator [] (int i){return v[i];}
01138 #else
01139   df1b2variable& operator () (int i,int j,int k);
01140   df1b2vector& operator () (int i,int j);
01141   df1b2matrix& operator () (int i);
01142   df1b2matrix& operator [] (int i);
01143 #endif
01144 };
01145 
01146 
01147 // **************************************
01148 // **************************************
01149 df1b2vector fabs(const df1b2vector& t1);
01150 df1b2vector mfexp(const df1b2vector& x);
01151 df1b2vector exp(const df1b2vector& x);
01152 df1b2vector sqrt(const df1b2vector& x);
01153 df1b2vector sin(const df1b2vector& x);
01154 df1b2vector tan(const df1b2vector& x);
01155 df1b2vector cos(const df1b2vector& x);
01156 df1b2vector log(const df1b2vector& x);
01157 
01158 df1b2matrix mfexp(const df1b2matrix& M);
01159 df1b2matrix log(const df1b2matrix& M);
01160 df1b2matrix trans(const df1b2matrix& M);
01161 df1b2matrix choleski_decomp(const df1b2matrix& M);
01162 df1b2matrix choleski_decomp_extra(const df1b2matrix& M);
01163 df1b2matrix exp(const df1b2matrix& M);
01164 df1b2vector column(const df1b2matrix&,int i);
01165 df1b2vector colsum(const df1b2matrix&);
01166 df1b2vector rowsum(const df1b2matrix&);
01167 
01168 df1b2vector operator * (double,const df1b2vector& y);
01169 df1b2vector operator * (const df1b2vector& y,const double);
01170 df1b2variable operator * (const df1b2vector& y,const dvector&);
01171 df1b2variable operator * (const dvector&,const df1b2vector& y);
01172 df1b2vector operator * (const df1b2vector& y,const df1b2variable&);
01173 df1b2vector operator * (const df1b2variable&,const df1b2vector& y);
01174 
01175 df1b2matrix outer_prod(const df1b2vector& x,const df1b2vector& y);
01176 df1b2matrix operator + (const df1b2matrix& M,const df1b2variable& x);
01177 
01178 df1b2matrix operator + (const df1b2matrix& M,const df1b2matrix& N);
01179 df1b2matrix operator + (const df1b2matrix& M,const dmatrix& N);
01180 
01181 df1b2matrix operator * (const df1b2matrix& M,const df1b2variable& x);
01182 df1b2matrix operator * (const dmatrix& M,const df1b2variable& x);
01183 df1b2vector operator * (const dmatrix& M,const df1b2vector& x);
01184 df1b2matrix operator * (const df1b2variable& x,const dmatrix& M);
01185 
01186 df1b2matrix operator * (const df1b2matrix& M,const df1b2matrix& x);
01187 df1b2matrix operator * (const dmatrix& M,const df1b2matrix& x);
01188 df1b2matrix operator + (const dmatrix& M,const df1b2matrix& x);
01189 df1b2matrix operator + (const df1b2variable&, const df1b2matrix& M);
01190 df1b2vector operator + (const df1b2variable&, const dvector& M);
01191 df1b2matrix operator * (const df1b2variable&, const df1b2matrix& M);
01192 df1b2matrix operator + (const df1b2matrix& M,const double x);
01193 df1b2matrix operator * (const df1b2matrix& M,const double x);
01194 
01195 df1b2matrix operator - (const df1b2matrix& M,const dmatrix& x);
01196 df1b2matrix operator - (const df1b2matrix& M,const double x);
01197 df1b2matrix operator - (const df1b2matrix& M,const df1b2variable& x);
01198 df1b2matrix operator - (const df1b2matrix& M,const df1b2matrix& x);
01199 df1b2matrix operator - (const dmatrix& M,const df1b2matrix& x);
01200 df1b2matrix operator - (const df1b2variable&, const df1b2matrix& M);
01201 df1b2matrix operator - (const double, const df1b2matrix& M);
01202 
01203 df1b2matrix operator + (const df1b2variable&, const df1b2matrix& M);
01204 
01205 df1b2matrix operator + (const double, const df1b2matrix& M);
01206 df1b2matrix operator * (const double, const df1b2matrix& M);
01207 df1b2matrix elem_prod(const df1b2matrix& x,const df1b2matrix& y);
01208 df1b2matrix elem_prod(const dmatrix& x,const df1b2matrix& y);
01209 df1b2matrix elem_prod(const df1b2matrix& x,const dmatrix& y);
01210 df1b2matrix elem_div(const df1b2matrix& x,const df1b2matrix& y);
01211 df1b2matrix elem_div(const dmatrix& x,const df1b2matrix& y);
01212 df1b2matrix elem_div(const df1b2matrix& x,const dmatrix& y);
01213 
01214 df1b2vector elem_prod(const df1b2vector& x,const df1b2vector& y);
01215 df1b2vector elem_prod(const df1b2vector& x,const dvector& y);
01216 df1b2vector elem_prod(const dvector& x,const df1b2vector& y);
01217 df1b2vector elem_div(const df1b2vector& x,const df1b2vector& y);
01218 df1b2vector elem_div(const dvector& x,const df1b2vector& y);
01219 df1b2vector elem_div(const df1b2vector& x,const dvector& y);
01220 
01221 df1b2vector operator + (const df1b2variable& x,const df1b2vector& y);
01222 df1b2vector operator + (double x,const df1b2vector& y);
01223 df1b2vector operator + (const df1b2vector& x,const df1b2vector& y);
01224 df1b2vector operator + (const df1b2vector& x,const df1b2variable& y);
01225 df1b2vector operator + (const dvector& x,const df1b2vector& y);
01226 df1b2vector operator + (const df1b2vector& y,const dvector& x);
01227 
01228 df1b2vector operator - (const df1b2vector& x,const df1b2vector& y);
01229 df1b2vector operator - (const dvector& x,const df1b2vector& y);
01230 df1b2vector operator - (const df1b2variable& x,const df1b2vector& y);
01231 df1b2vector operator - (const df1b2vector& x,const df1b2vector& y);
01232 df1b2vector operator - (const df1b2variable& x,const df1b2vector& y);
01233 df1b2vector operator - (const df1b2variable& x,const dvector& y);
01234 df1b2vector operator - (const df1b2vector& x,const df1b2variable& y);
01235 df1b2vector operator - (const df1b2vector& x,const dvector& y);
01236 
01237 df1b2vector operator * (const df1b2variable& y,const dvector&);
01238 
01239 df1b2vector operator * (const df1b2vector& x,const df1b2variable& y);
01240 df1b2vector operator * (const df1b2vector& x,double y);
01241 
01242 df1b2vector operator / (const df1b2variable& y,const df1b2vector& x);
01243 df1b2vector operator / (const double& y,const df1b2vector& x);
01244 
01245 df1b2vector operator / (const df1b2vector& x,const df1b2variable& y);
01246 df1b2vector operator / (const df1b2vector& x,double y);
01247 df1b2vector pow(const df1b2vector& x,double y);
01248 df1b2vector pow(const df1b2vector& v,const df1b2variable & x);
01249 df1b2vector pow(const df1b2vector& v,const df1b2vector & x);
01250 df1b2vector pow(const df1b2variable& v,const df1b2vector & x);
01251 df1b2vector pow(df1b2vector const& _x,dvector const& v);
01252 df1b2vector pow(double v,const df1b2vector & x);
01253 df1b2vector pow(const dvector& x,  const df1b2vector& a);
01254 df1b2vector pow(const dvector& x,  const df1b2variable& a);
01255 
01256 df1b2vector operator / (const dvector& x,const df1b2variable& y);
01257 df1b2vector operator + (const dvector& x,const df1b2variable& y);
01258 df1b2vector operator - (const dvector& x,const df1b2variable& y);
01259 df1b2vector operator * (const dvector& x,const df1b2variable& y);
01260 // **************************************
01261 // **************************************
01262 df1b2variable sum(const df1b2vector& x);
01263 df1b2variable mean(const df1b2vector& x);
01264 df1b2variable norm2(const df1b2vector& x);
01265 df1b2variable sumsq(const df1b2vector& x);
01266 df1b2variable norm(const df1b2vector& x);
01267 
01268 df1b2variable norm2(const df1b2matrix& x);
01269 df1b2variable sumsq(const df1b2matrix& x);
01270 df1b2variable norm(const df1b2matrix& x);
01271 df1b2variable mean(const df1b2matrix& x);
01272 df1b2variable sum(const df1b2matrix& x);
01273 df1b2matrix square(const df1b2matrix& x);
01274 df1b2vector square(const df1b2vector& x);
01275 df1b2vector cube(const df1b2vector& x);
01276 df1b2vector solve(const df1b2matrix& x,const dvector&);
01277 df1b2vector solve(const df1b2matrix& x,const df1b2vector&);
01278 
01279 int size_count(const df1b2matrix& x);
01280 int size_count(const df1b2vector& x);
01281 
01282 df1b2vector first_difference(const df1b2vector& x);
01283 
01284 df1b2vector operator * (const dvector& ,const df1b2matrix& );
01285 df1b2vector operator * (const df1b2vector& ,const df1b2matrix& );
01286 df1b2vector operator * (const df1b2vector& ,const dmatrix& );
01287 df1b2variable operator * (const df1b2vector& x,const df1b2vector& y);
01288 df1b2vector operator * (const df1b2matrix& M,const df1b2vector& x);
01289 df1b2vector operator * (const df1b2matrix& M,const dvector& x);
01290 void check_shape(const df1b2vector & x,const df1b2vector & y,const char * s);
01291 void check_shape(const df1b2vector & x,const dvector & y,const char * s);
01292 void check_shape(const dvector & x,const df1b2vector & y,const char * s);
01293 void check_shape(const df1b2vector & x,const df1b2matrix & y,const char * s);
01294 
01295 void checkidentiferstring(const char * ids,test_smartlist& list);
01296 
01297 double& value(const df1b2variable& _x);
01298 dvector value(const df1b2vector& _x);
01299 dmatrix value(const df1b2matrix& _x);
01300 
01305 class initial_df1b2params
01306 {
01307 public:
01308   static int stddev_scale(const dvector& d,const dvector& x);
01309   static int current_phase;
01310   static int separable_flag;
01311   static int separable_calculation_type;
01312   static int have_bounded_random_effects;
01313   int phase_start;
01314   int phase_save;
01315   int ind_index;
01316   double scalefactor;
01317 
01318   static void save_varsptr(void);
01319   static double cobjfun;
01320   static void restore_varsptr(void);
01321 #if defined(__x86_64) || (defined(_MSC_VER) && defined(_M_X64))
01322   static lmatrix* pointer_table;
01323 #else
01324   static imatrix* pointer_table;
01325 #endif
01326   static initial_df1b2params **  varsptr;  // this should be a resizeable
01327   static initial_df1b2params **  varsptr_sav;  // this should be a resizeable
01328   static int num_initial_df1b2params;         // array
01329   static int num_initial_df1b2params_sav;         // array
01330   virtual void set_value(const dvector&,const int& ii)=0;
01331   virtual void set_index(int aflag,const int& ii)=0;
01332   static int set_index(void);
01333   virtual void set_value(const init_df1b2vector&,const int& ii,
01334     const df1b2variable&)=0;
01335   static void reset(const init_df1b2vector&,const df1b2variable&);
01336   static void reset_all(const dvector&);
01337   static void reset(const df1b2vector&,const df1b2variable&);
01338   virtual void add_to_list(void);
01339   initial_df1b2params(void);
01340   void set_phase_start(int n){phase_start=n; phase_save=n;}
01341   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii)=0;
01342   double get_scalefactor();
01343   void set_scalefactor(const double);
01344 };
01345 
01346 typedef initial_df1b2params * P_INITIAL_DF1B2PARAMS;
01347 
01352 class do_naught_kludge
01353 {
01354 };
01355 
01360 class df1b2_init_vector : public df1b2vector, public initial_df1b2params
01361 {
01362 public:
01363   virtual ~df1b2_init_vector() {;}
01364 
01365   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01366   void set_phase_start(int n){phase_start=n; phase_save=n;}
01367   virtual void set_value(const dvector& x,const int& _ii);
01368   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01369     const df1b2variable& pen);
01370   void allocate(int min, int max,int phase_start, const adstring& s="unnamed");
01371   void allocate(int min, int max,const adstring& s="unnamed");
01372   virtual void set_index(int aflag,const int& ii);
01373 };
01374 
01379 class df1b2_init_matrix : public df1b2matrix, public initial_df1b2params
01380 {
01381 public:
01382   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01383 
01384   inline void allocate(int rmin,int rmax,const index_type& cmin,
01385     const index_type& cmax,int ps,const char * s)
01386   {
01387     set_phase_start(ps);
01388     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01389   }
01390 
01391   inline void allocate(int rmin,int rmax,const index_type& cmin,
01392     const index_type& cmax,const char * s)
01393   {
01394     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01395   }
01396 
01397   inline void allocate(int rmin,int rmax,int cmin,int cmax)
01398   {
01399     df1b2matrix::allocate(rmin,rmax,cmin,cmax);
01400   }
01401   inline void allocate(int rmin,int rmax,int cmin,int cmax,int ps,
01402     const char * s)
01403   {
01404     set_phase_start(ps);
01405     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01406   }
01407   inline void allocate(int rmin,int rmax,int cmin,int cmax,const char * s)
01408   {
01409     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01410   }
01411   inline void allocate(int rmin, int rmax, int cmin, int cmax, double, double,
01412     const char* s)
01413   {
01414     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01415   }
01416   void set_phase_start(int n){phase_start=n; phase_save=n;}
01417   virtual void set_value(const dvector& x,const int& _ii);
01418   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01419     const df1b2variable& pen);
01420   virtual void set_index(int aflag,const int& ii);
01421 };
01422 
01423 
01428 class df1b2_init_number : public df1b2variable, public initial_df1b2params
01429 {
01430   static do_naught_kludge do_naught_kludge_a;
01431 public:
01432   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01433   virtual int & get_ind_index(void) { return ind_index; }
01434   void set_phase_start(int n){phase_start=n; phase_save=n;}
01435   virtual void set_value(const dvector& x,const int& _ii);
01436   virtual void set_value(const init_df1b2vector&,const int& ii,
01437     const df1b2variable&);
01438   virtual void set_index(int aflag,const int& ii);
01439 
01440   void allocate(const adstring&);
01441   void allocate(int,const adstring&);
01442   df1b2_init_number();
01443   void operator = (const df1b2variable& _x);
01444 };
01445 
01450 class df1b2_init_number_vector
01451 {
01452   df1b2_init_number * v;
01453   int index_min;
01454   int index_max;
01455   double_index_type * it;
01456 public:
01457   df1b2_init_number_vector();
01458 
01459 #if defined(OPT_LIB)
01460    df1b2_init_number& operator [] (int i) { return v[i];}
01461    df1b2_init_number& operator () (int i) { return v[i];}
01462 #else
01463    df1b2_init_number& operator [] (int i);
01464    df1b2_init_number& operator () (int i);
01465 #endif
01466   void allocate(int min1,int max1,const index_type& phase_start,
01467     const char * s);
01468 
01469   void allocate(int min1,int max1,const char * s);
01470 
01471   int allocated(void) { return (v!=NULL); }
01472   int indexmin(void) {return (index_min);}
01473   int indexmax(void) {return (index_max);}
01474   ~df1b2_init_number_vector();
01475   void deallocate(void);
01476 };
01477 
01478 class df1b2_init_bounded_number;
01479 
01484 class df1b2_init_bounded_number :public df1b2_init_number
01485 {
01486   double minb;
01487   double maxb;
01488 public:
01489   void allocate(double _minb,double _maxb,int _phase_start,
01490     const char * s="UNNAMED");
01491   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01492 
01493   void allocate(double _minb,double _maxb,const char * s="UNNAMED");
01494   //void allocate(const ad_double& _minb,const ad_double& _maxb,
01495    // const char * s="UNNAMED");
01496 
01497   virtual void set_value(const dvector& x,const int& _ii);
01498   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01499     const df1b2variable& pen);
01500 };
01501 
01506 class df1b2_init_bounded_number_vector
01507 {
01508   df1b2_init_bounded_number * v;
01509   int index_min;
01510   int index_max;
01511   double_index_type * it;
01512 public:
01513   //virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01514   df1b2_init_bounded_number_vector();
01515 
01516 #if defined(OPT_LIB)
01517    df1b2_init_bounded_number& operator [] (int i) { return v[i];}
01518    df1b2_init_bounded_number& operator () (int i) { return v[i];}
01519 #else
01520    df1b2_init_bounded_number& operator [] (int i);
01521    df1b2_init_bounded_number& operator () (int i);
01522 #endif
01523   void allocate(int min1,int max1, const double_index_type& bmin,
01524    const double_index_type& bmax, const char * s);
01525 
01526   void allocate(int min1,int max1, const double_index_type& bmin,
01527    const double_index_type& bmax, const index_type& phase_start,
01528     const char * s);
01529 
01530 
01531   int allocated(void) { return (v!=NULL); }
01532   int indexmin(void) {return (index_min);}
01533   int indexmax(void) {return (index_max);}
01534   ~df1b2_init_bounded_number_vector();
01535   void deallocate(void);
01536 };
01537 
01542 class df1b2_init_bounded_vector : public df1b2_init_vector
01543 {
01544 protected:
01545   double minb;
01546   double maxb;
01547 public:
01548   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01549   void allocate(int,int,double,double,int,const char *);
01550   void allocate(int,int,double,double,const char *);
01551   void allocate(int,int,double,double);
01552   //void allocate(int,int);
01553   virtual void set_value(const dvector& x,const int& _ii);
01554   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01555     const df1b2variable& pen);
01556   inline double getminb(void){return minb;}
01557   inline double getmaxb(void){return maxb;}
01558 };
01559 
01560 class re_df1b2_init_bounded_vector;
01561 
01566 class random_effects_bounded_vector_info
01567 {
01568   re_df1b2_init_bounded_vector * pv;
01569   int i;
01570 public:
01571   random_effects_bounded_vector_info
01572     ( re_df1b2_init_bounded_vector * _pv,int _i) : pv(_pv), i(_i) {}
01573   friend class funnel_init_df1b2variable;
01574   friend class df1b2variable;
01575   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01576     const df1b2variable& pen);
01577 };
01578 
01583 class re_df1b2_init_bounded_vector : public df1b2_init_bounded_vector
01584 {
01585 public:
01586   random_effects_bounded_vector_info operator () (int i)
01587   {
01588     return random_effects_bounded_vector_info(this,i);
01589   }
01590   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01591     const df1b2variable& pen);
01592   virtual void set_value(const dvector& x,const int& _ii);
01593   re_df1b2_init_bounded_vector(void);
01594 };
01595 
01600 class df1b2_init_bounded_matrix : public df1b2_init_matrix
01601 {
01602 protected:
01603   double minb;
01604   double maxb;
01605 public:
01606   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01607   void allocate(int,int,int,int,double,double,int,const char *);
01608   void allocate(int,int,const index_type& ,const index_type& ,double,double,int,
01609     const char *);
01610   void allocate(int,int,int,int,double,double,const char *);
01611   void allocate(int,int,const index_type& ,const index_type& ,double,double,
01612     const char *);
01613   void allocate(int,int,int,int,double,double);
01614   //void allocate(int,int);
01615   virtual void set_value(const dvector& x,const int& _ii);
01616   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01617     const df1b2variable& pen);
01618 };
01619 
01624 class df1b2_init_bounded_dev_vector : public df1b2_init_bounded_vector
01625 {
01626 public:
01627   virtual void set_value(const dvector& x,const int& _ii);
01628   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01629     const df1b2variable& pen);
01630 };
01631 
01632 void set_value(const df1b2variable& u,const init_df1b2vector& x,const int& _ii,
01633   double fmin,double fmax,const df1b2variable& fpen);
01634 
01635 void set_value(const df1b2variable& u,const init_df1b2vector& x,const int& _ii,
01636   double fmin,double fmax,const df1b2variable& fpen,double s);
01637 
01638 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax);
01639 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax,
01640   const df1b2variable& _fpen);
01641 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax,
01642   const df1b2variable& _fpen,double s);
01643 
01644   ostream& operator << (const ostream& _os, const df1b2variable& _x);
01645   ostream& operator << (const ostream& _os, const df1b2vector& _x);
01646   ostream& operator << (const ostream& _os, const df1b2matrix& _x);
01647 
01648   class re_objective_function_value : public df1b2variable
01649   {
01650   public:
01651     static re_objective_function_value * pobjfun;
01652     static double fun_without_pen;
01653     re_objective_function_value(void);
01654     ~re_objective_function_value();
01655     re_objective_function_value& operator = (const df1b2variable& v);
01656     re_objective_function_value& operator = (double v);
01657     void allocate(void);
01658     void allocate(const char *);
01659   };
01660 
01661 df1b2variable posfun(const df1b2variable&x,const double eps);
01662 df1b2variable posfun2(const df1b2variable&x,const double eps,
01663   const df1b2variable& _pen);
01664 df1b2variable posfun(const df1b2variable&x,const double eps,
01665   const df1b2variable& _pen);
01666 
01671 inline df1b2variable operator - (const df1b2variable& x)
01672 {
01673   return -double(1.0)*x;
01674 }
01675 
01680 inline df1b2vector operator - (const df1b2vector& x)
01681 {
01682   return -double(1.0)*x;
01683 }
01684 
01685 df1b2variable gamma_density(const df1b2variable& _x,double r, double mu);
01686 df1b2variable gamma_density(const df1b2variable& _x,const df1b2variable& _r,
01687   const  df1b2variable& _mu);
01688 
01689 df1b2variable log_gamma_density(const df1b2variable& _x,double r, double mu);
01690 df1b2variable log_gamma_density(const df1b2variable& _x,const df1b2variable& _r,
01691   const  df1b2variable& _mu);
01692 
01693 
01694 df1b2variable log_negbinomial_density(double x,const df1b2variable& mu,
01695   const df1b2variable& tau);
01696 
01697 df1b2variable log_density_poisson(double x,const df1b2variable& mu);
01698 
01699 ostream& operator << (const ostream& ostr,const df1b2variable& z);
01700 ostream& operator << (const ostream& ostr,const df1b2vector& z);
01701 ostream& operator << (const ostream& ostr,const df1b2matrix& z);
01702 ostream& operator << (const ostream& ostr,const df1b2_init_number_vector& z);
01703 ostream& operator << (const ostream& ostr,const init_df1b2vector& z);
01704 ostream& operator << (const ostream& ostr,
01705   const df1b2_init_bounded_number_vector& z);
01706 
01707 class adkludge1;
01708 
01709 
01714   class df1b2_header_ptr_vector
01715   {
01716     int index_min;
01717     int index_max;
01718     df1b2_header ** v;
01719   public:
01720     df1b2_header_ptr_vector(int mmin,int mmax);
01721     ~df1b2_header_ptr_vector();
01722     inline df1b2_header * & operator () (int i) { return *(v+i); }
01723     inline df1b2_header * & operator [] (int i) { return *(v+i); }
01724     int indexmin(void) {return index_min;}
01725     int indexmax(void) {return index_max;}
01726   };
01727 
01732   class double_ptr_vector
01733   {
01734     int index_min;
01735     int index_max;
01736     double ** v;
01737   public:
01738     double_ptr_vector(int mmin,int mmax);
01739     ~double_ptr_vector();
01740     inline double * & operator () (int i) { return *(v+i); }
01741     inline double * & operator [] (int i) { return *(v+i); }
01742     int indexmin(void) {return index_min;}
01743     int indexmax(void) {return index_max;}
01744   };
01745 
01746 int active(const df1b2_init_vector& x);
01747 int active(const df1b2_init_number& x);
01748 int active(const df1b2_init_matrix& x);
01749 //int active(const df1b2_init_bounded_vector& x);
01750 
01751 //int active (const df1b2_init_vector &);
01752 //int active(const initial_df1b2params& x);
01753 //int active(const df1b2vector& x);
01754 //int active(const df1b2matrix& x);
01755 
01756 double evaluate_function_quiet(const dvector& x,function_minimizer * pfmin);
01757 double evaluate_function(const dvector& x,function_minimizer * pfmin);
01758 double evaluate_function(double& fval,const dvector& x,
01759   function_minimizer * pfmin);
01760 double evaluate_function(double& fval,const dvector& x,const dvector&,
01761   function_minimizer * pfmin);
01762 
01763 
01764 double evaluate_function_no_derivatives(const dvector& x,
01765   function_minimizer * pfmin);
01766 
01767 #include <df1b2fnl.h>
01768 
01769 int allocated(const df1b2_init_vector&);
01770 int allocated(const df1b2vector&);
01771 int allocated(const df1b2_init_matrix&);
01772 #include <df3fun.h>
01773 
01774 /*
01775 df1b2variable betai(const df1b2variable & _a, const df1b2variable & _b,
01776          const df1b2variable & _x);
01777 df1b2variable betai(const df1b2variable & _a, const df1b2variable & _b,
01778          double _x);
01779 */
01780 
01781 df1b2variable betai(const df1b2variable& a, const df1b2variable& b,
01782   double x, int maxit=100);
01783 
01784 double do_gauss_hermite_block_diagonal(const dvector& x,
01785   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01786   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01787   function_minimizer * pmin);
01788 
01789 double do_gauss_hermite_block_diagonal_multi(const dvector& x,
01790   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01791   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01792   function_minimizer * pmin);
01793 
01794 double calculate_importance_sample_block_diagonal(const dvector& x,
01795   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01796   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01797   function_minimizer * pmin);
01798 
01799 double calculate_importance_sample_block_diagonal_option2(const dvector& x,
01800   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01801   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01802   function_minimizer * pmin);
01803 
01804 double calculate_importance_sample_block_diagonal_option_antithetical
01805   (const dvector& x,const dvector& u0,const dmatrix& Hess,
01806   const dvector& _xadjoint,const dvector& _uadjoint,
01807   const dmatrix& _Hessadjoint,function_minimizer * pmin);
01808 
01809 double calculate_importance_sample_block_diagonal_funnel(const dvector& x,
01810   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01811   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01812   function_minimizer * pmin);
01813 
01814 double calculate_importance_sample(const dvector& x,const dvector& u0,
01815   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01816   const dmatrix& _Hessadjoint,function_minimizer * pmin);
01817 
01818 double calculate_importance_sample_shess(const dvector& x,const dvector& u0,
01819   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01820   const dmatrix& _Hessadjoint,function_minimizer * pmin);
01821 double calculate_importance_sample(const dvector& x,const dvector& u0,
01822   const banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
01823   const dvector& _uadjoint,
01824   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin);
01825 
01826 dvector evaluate_function_with_quadprior(const dvector& x,int usize,
01827   function_minimizer * pfmin);
01828 
01833 class style_flag_class
01834 {
01835 public:
01836   int old_style_flag;
01837   virtual void set_old_style_flag(void)=0;
01838   //style_flag_class(void) {set_old_style_flag();}
01839 };
01840 
01845 class quadratic_prior : public style_flag_class
01846 {
01847   dmatrix * pMinv;
01848   dvar_matrix * dfpMinv;
01849   dvar_vector * pu;
01850   int xmyindex;
01851 public:
01852   static int qflag;
01853   static quadratic_prior * ptr[]; // this should be a resizeable array
01854   static void get_M_calculations(void);
01855   static void cleanup_pMinv();
01856   static void cleanup_dfpMinv();
01857   static int num_quadratic_prior;
01858   static const int max_num_quadratic_prior;
01859   void add_to_list(void);
01860 public:
01861   static int in_qp_calculations;
01862   int get_offset(int xs);
01863   int get_myindex(void) { return xmyindex;}
01864   static quadratic_prior * get_ptr(int i){ return ptr[i];}
01865   void operator = (const dvar_matrix &);
01866   void operator = (const dmatrix &);
01867   static int get_in_qp_calculations() { return in_qp_calculations; }
01868   static int get_num_quadratic_prior(void) { return num_quadratic_prior;}
01869   static dvariable get_quadratic_priors(void);
01870   dvariable get_function(void);
01871   quadratic_prior(void);
01872   ~quadratic_prior(void);
01873   void allocate(const dvar_vector & _u,const char * s);
01874   void allocate(const dvar_vector & _u);
01875   void allocate(const dvar_matrix & _M, const dvar_vector & _u,const char * s);
01876   void allocate(const dvar_matrix & _M, const dvar_vector & _u);
01877   //dmatrix get_cHessian(void);
01878   void get_cHessian(dmatrix,int);
01879   void get_cHessian(dvar_matrix,int);
01880   void get_cHessian_from_vHessian(dmatrix ,int);
01881   void get_vHessian(dvar_matrix ,int);
01882   void get_cgradient(dvector,int);
01883   dvar_matrix get_Hessian(void);
01884   dvar_vector get_gradient(void);
01885   static dvar_vector get_gradient_contribution(void);
01886   static dvar_matrix get_Hessian_contribution(void);
01887   static void get_cgradient_contribution(dvector,int);
01888   static void get_cHessian_contribution(dmatrix,int);
01889   static void get_vHessian_contribution(dvar_matrix  ,int );
01890   static void get_cHessian_contribution_from_vHessian(dmatrix,int);
01891   friend class df1b2quadratic_prior;
01892   virtual void get_cM(void)=0;
01893 };
01894 
01899 class df1b2quadratic_prior : public style_flag_class
01900 {
01901   ivector * index;
01902   df1b2matrix * M;
01903   dmatrix * CM;
01904 public:
01905   dmatrix * Lxu;
01906   df1b2_init_vector * pu;
01907   int xmyindex;
01908   static df1b2quadratic_prior * ptr[]; // this should be a resizeable array
01909   static int num_quadratic_prior;
01910   static const int max_num_quadratic_prior;
01911   void add_to_list(void);
01912 public:
01913   static df1b2quadratic_prior * get_ptr(int i){ return ptr[i];}
01914   unsigned int num_active_parameters;
01915   unsigned int get_num_active_parameters(void) { return num_active_parameters; }
01916   int get_myindex(void) { return xmyindex;}
01917   void operator = (const df1b2matrix &);
01918   void operator = (const dmatrix &);
01919   static int get_num_quadratic_prior(void) { return num_quadratic_prior;}
01920   static dvariable get_quadratic_priors(void);
01921   df1b2variable get_function(void);
01922   df1b2quadratic_prior(void);
01923   ~df1b2quadratic_prior(void);
01924   void allocate(const df1b2_init_vector & _u,const char * s);
01925   void allocate(const df1b2_init_vector & _u);
01926   void allocate(const df1b2matrix & _M, const df1b2_init_vector & _u,
01927     const char * s);
01928   void allocate(const df1b2matrix & _M, const df1b2_init_vector & _u);
01929   void allocate(const dvar_matrix & _M, const dvar_vector & _u,const char * s);
01930   void allocate(const dvar_matrix & _M, const dvar_vector & _u);
01931   dmatrix get_cHessian(void);
01932   dvector get_cgradient(void);
01933   dvar_matrix get_Hessian(void);
01934   dvar_vector get_gradient(void);
01935   static dvar_vector get_gradient_contribution(void);
01936   static dvar_matrix get_Hessian_contribution(void);
01937   static dvector get_cgradient_contribution(void);
01938   static dmatrix get_cHessian_contribution(void);
01939   static void get_Lxu_contribution(dmatrix&);
01940   virtual void get_Lxu(dmatrix&)=0;
01941   friend class quadratic_prior;
01942   friend class df1b2_parameters;
01943 };
01944 
01949 class normal_quadratic_prior : public quadratic_prior
01950 {
01951   virtual void set_old_style_flag(void);
01952 public:
01953   normal_quadratic_prior(void);
01954   void operator = (const dvar_matrix & M);
01955 };
01956 
01961 class normal_df1b2quadratic_prior : public df1b2quadratic_prior
01962 {
01963   virtual void set_old_style_flag(void);
01964 public:
01965   void operator = (const df1b2matrix & M);
01966   normal_df1b2quadratic_prior(void);
01967 };
01968 
01973 class quadratic_re_penalty : public quadratic_prior
01974 {
01975   virtual void set_old_style_flag(void);
01976 public:
01977   quadratic_re_penalty(void);
01978   void operator = (const dvar_matrix & M);
01979   void operator = (const dmatrix & M);
01980 };
01981 
01986 class df1b2quadratic_re_penalty : public df1b2quadratic_prior
01987 {
01988   virtual void set_old_style_flag(void);
01989 public:
01990   void operator = (const df1b2matrix & M);
01991   void operator = (const dmatrix & M);
01992   df1b2quadratic_re_penalty(void);
01993 };
01994 
01995 dvar_vector solve(const named_dvar_matrix &,const random_effects_vector&);
01996 
02001 class constant_quadratic_re_penalty : public quadratic_prior
02002 {
02003   virtual void set_old_style_flag(void);
02004 public:
02005   constant_quadratic_re_penalty(void);
02006   void operator = (const dmatrix & M);
02007 };
02008 
02013 class constant_df1b2quadratic_re_penalty : public df1b2quadratic_prior
02014 {
02015   virtual void set_old_style_flag(void);
02016 public:
02017   void operator = (const dmatrix & M);
02018   constant_df1b2quadratic_re_penalty(void);
02019 };
02020 
02021 df1b2vector solve(df1b2matrix& M,df1b2vector& v,const df1b2variable& ln_det);
02022 
02023 df1b2vector solve(df1b2matrix& M,df1b2vector& v,const df1b2variable& ln_det,
02024   const int& sgn);
02025 df1b2vector lower_triangular_solve(const df1b2matrix& m,const df1b2vector& v);
02026 df1b2vector lower_triangular_solve_trans(const df1b2matrix& m,
02027   const df1b2vector& v);
02028 
02029 //df1b2variable ln_det(df1b2matrix& M);
02030 // line above replaced with line below based on issue #37
02031 df1b2variable ln_det(const df1b2matrix & m1);
02032 
02033 df1b2variable ln_det(df1b2matrix& M,int & sgn);
02034 
02035 //df1b2vector solve(df1b2matrix& M,df1b2vector& v);
02036 
02037 df1b2matrix expm(const df1b2matrix & A);
02038 df1b2matrix solve(const df1b2matrix& aa,const df1b2matrix& tz,
02039   df1b2variable ln_unsigned_det,df1b2variable& sign);
02040 df1b2matrix solve(const df1b2matrix& aa,const df1b2matrix& tz);
02041 df1b2vector solve(const df1b2matrix& aa,const df1b2vector& z,
02042   const df1b2variable& ld, df1b2variable& sign);
02043 
02044 void check_pool_depths();
02045 
02046 df1b2variable lower_triangular_ln_det(const df1b2matrix& m);
02047 df1b2variable lower_triangular_ln_det(const df1b2matrix& m,int& sgn);
02048 df1b2variable bounder(const df1b2variable&  x,double min,double max,
02049   double scale);
02050 
02051 df1b2variable inv_cumd_beta_stable(const df1b2variable& a,
02052   const df1b2variable& b,const df1b2variable& x,double eps=1.e-7);
02053 
02054 df1b2variable bounded_cumd_norm(const df1b2variable& _x,double beta);
02055 df1b2variable cumd_norm(const df1b2variable& _x);
02056 df1b2variable inv_cumd_exponential(const df1b2variable& y);
02057 df1b2variable cumd_exponential(const df1b2variable& x);
02058 extern int make_sure_df1b2fun_gets_linked;
02059 
02060 df1b2variable inv_cumd_normal_mixture(const df1b2variable& _x,double _a);
02061 df1b2variable inv_cumd_normal_logistic_mixture(const df1b2variable& _x,
02062   double _a);
02063 
02064 df1b2variable beta_deviate(const df1b2variable& _x,const df1b2variable& _a,
02065   const df1b2variable& _b,double eps=1.e-7);
02066 df1b2variable robust_normal_logistic_mixture_deviate(const df1b2variable& x,
02067   double spread=3.0);
02068 df1b2variable robust_normal_mixture_deviate(const df1b2variable& x,
02069   double spread=3.0);
02070 
02071 df1b2variable gamma_deviate(const df1b2variable& _x,const df1b2variable& _a);
02072 df1b2variable inv_cumd_gamma(const df1b2variable& _y,const df1b2variable& _a);
02073 double inv_cumd_gamma(double y,double _a);
02074 
02075 df1b2variable inv_cumd_cauchy(const df1b2variable& n);
02076 df1b2variable inv_cumd_t(const df1b2variable& n,const df1b2variable&  u,
02077   double eps=1.e-7);
02078 
02083   class df1b2function_tweaker
02084   {
02085     double mult;
02086     double eps;
02087     dvector coffs;
02088   public:
02089     df1b2function_tweaker(double eps,double mult);
02090     df1b2variable operator () (const df1b2variable&);
02091   };
02092 
02093 df1b2vector posfun(const df1b2vector& x,double eps,
02094   const df1b2variable& _pen);
02095 
02096 df1b2variable norm_to_gamma(const df1b2variable & v,const df1b2variable& alpha,
02097   double bound=0.999999);
02098 
02099 void print_is_diagnostics(laplace_approximation_calculator *lapprox);
02100 
02101   ofstream &  operator << (const ofstream& _s,const nested_calls_shape& m);
02102 
02103   df1b2variable log_der_logistic(double a,double b,const df1b2variable& x);
02104   df1b2variable logistic(double a,double b,const df1b2variable& x);
02105   df1b2variable dflogistic(double a,double b,const df1b2variable& x);
02106 
02107   df1b2variable dfposfun(const df1b2variable& x,const double eps);
02108   void ADMB_getcallindex(const df1b2variable& x);
02109   void ADMB_getcallindex(const df1b2vector& x);
02110   void ADMB_getcallindex(const df1b2matrix& x);
02111 
02112 df1b2variable asin(const df1b2variable& xx);
02113 #endif //!defined(__DF1B2FUN__)