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