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