ADMB Documentation  11.1x.2738
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fun.h
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2fun.h 2637 2014-11-13 05:58:08Z 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 #if defined(OPT_LIB)
01051   df1b2variable& operator () (int i,int j) const
01052     { return (df1b2variable&)(v[i][j]);}
01053   df1b2vector& operator [] (int i) const {return (df1b2vector&)(v[i]);}
01054   df1b2vector& operator () (int i) const {return (df1b2vector&)(v[i]);}
01055 #else
01056   df1b2variable& operator () (int i,int j) const;
01057   df1b2vector& operator [] (int i) const;
01058   df1b2vector& operator () (int i) const;
01059 #endif
01060   int colmin(void) const {return (*(df1b2matrix*)(this))(index_min).indexmin();}
01061   int colmax(void) const {return (*(df1b2matrix*)(this))(index_min).indexmax();}
01062   int colsize(void)const {return colmax()-colmin()+1;}
01063   df1b2matrix& operator *= (const df1b2variable& x);
01064   df1b2matrix& operator /= (const df1b2variable& x);
01065   df1b2matrix& operator *= (double x);
01066   df1b2matrix& operator /= (double x);
01067 };
01068 
01073 class df1b23array
01074 {
01075   int index_min;
01076   int index_max;
01077   df1b2matrix * v;
01078   vector_shapex * shape;
01079 public:
01080   int allocated(void){return v!=0;}
01081   void initialize(void);
01082   ~df1b23array();
01083   int indexmin(void){return index_min;}
01084   int indexmax(void){return index_max;}
01085   int size(void){return index_max-index_min+1;}
01086   df1b23array(int nrl,int nrh,int ncl,int nch,int,int);
01087   df1b23array(int nrl,int nrh);
01088   df1b23array(int nrl,int nrh,int,int);
01089   df1b23array(const df1b23array&);
01090   //df1b23array& operator =(const df3_one_matrix&);
01091 
01092   df1b23array& operator =(const df1b23array&);
01093   df1b23array& operator =(const df1b2variable&);
01094   df1b23array& operator =(double);
01095   df1b23array(void);
01096   df1b23array& operator +=(const df1b23array& M);
01097   df1b23array& operator -=(const df1b23array& M);
01098   void allocate(int nrl,int nrh,int ncl,int nch,int,int);
01099   void allocate(int nrl,int nrh,int,int);
01100   void allocate(int nrl,int nrh);
01101   void allocate(int nrl,int nrh,int ncl,int nch,int,int,const char *);
01102 
01103   void allocate(int nrl,int nrh,
01104     const index_type& ncl, const index_type& nch);
01105 
01106 
01107   void allocate(int nrl,int nrh,
01108     const index_type& ncl, const index_type& nch,
01109     const index_type& nxl, const index_type& nxh);
01110 
01111   void allocate(int nrl,int nrh,
01112     const index_type& ncl,const index_type& nch,
01113     const index_type& nxl,const index_type& nxh,
01114     const char *);
01115   void allocate(void);
01116   void deallocate(void);
01117 #  if defined(OPT_LIB)
01118   df1b2variable& operator () (int i,int j,int k){return v[i](j,k);}
01119   df1b2vector& operator () (int i,int j){return v[i][j];}
01120   df1b2matrix& operator () (int i){return v[i];}
01121   df1b2matrix& operator [] (int i){return v[i];}
01122 #else
01123   df1b2variable& operator () (int i,int j,int k);
01124   df1b2vector& operator () (int i,int j);
01125   df1b2matrix& operator () (int i);
01126   df1b2matrix& operator [] (int i);
01127 #endif
01128 };
01129 
01130 
01131 // **************************************
01132 // **************************************
01133 df1b2vector fabs(const df1b2vector& t1);
01134 df1b2vector mfexp(const df1b2vector& x);
01135 df1b2vector exp(const df1b2vector& x);
01136 df1b2vector sqrt(const df1b2vector& x);
01137 df1b2vector sin(const df1b2vector& x);
01138 df1b2vector tan(const df1b2vector& x);
01139 df1b2vector cos(const df1b2vector& x);
01140 df1b2vector log(const df1b2vector& x);
01141 
01142 df1b2matrix mfexp(const df1b2matrix& M);
01143 df1b2matrix log(const df1b2matrix& M);
01144 df1b2matrix trans(const df1b2matrix& M);
01145 df1b2matrix choleski_decomp(const df1b2matrix& M);
01146 df1b2matrix choleski_decomp_extra(const df1b2matrix& M);
01147 df1b2matrix exp(const df1b2matrix& M);
01148 df1b2vector column(const df1b2matrix&,int i);
01149 df1b2vector colsum(const df1b2matrix&);
01150 df1b2vector rowsum(const df1b2matrix&);
01151 
01152 df1b2vector operator * (double,const df1b2vector& y);
01153 df1b2vector operator * (const df1b2vector& y,const double);
01154 df1b2variable operator * (const df1b2vector& y,const dvector&);
01155 df1b2variable operator * (const dvector&,const df1b2vector& y);
01156 df1b2vector operator * (const df1b2vector& y,const df1b2variable&);
01157 df1b2vector operator * (const df1b2variable&,const df1b2vector& y);
01158 
01159 df1b2matrix outer_prod(const df1b2vector& x,const df1b2vector& y);
01160 df1b2matrix operator + (const df1b2matrix& M,const df1b2variable& x);
01161 
01162 df1b2matrix operator + (const df1b2matrix& M,const df1b2matrix& N);
01163 df1b2matrix operator + (const df1b2matrix& M,const dmatrix& N);
01164 
01165 df1b2matrix operator * (const df1b2matrix& M,const df1b2variable& x);
01166 df1b2matrix operator * (const dmatrix& M,const df1b2variable& x);
01167 df1b2vector operator * (const dmatrix& M,const df1b2vector& x);
01168 df1b2matrix operator * (const df1b2variable& x,const dmatrix& M);
01169 
01170 df1b2matrix operator * (const df1b2matrix& M,const df1b2matrix& x);
01171 df1b2matrix operator * (const dmatrix& M,const df1b2matrix& x);
01172 df1b2matrix operator + (const dmatrix& M,const df1b2matrix& x);
01173 df1b2matrix operator + (const df1b2variable&, const df1b2matrix& M);
01174 df1b2vector operator + (const df1b2variable&, const dvector& M);
01175 df1b2matrix operator * (const df1b2variable&, const df1b2matrix& M);
01176 df1b2matrix operator + (const df1b2matrix& M,const double x);
01177 df1b2matrix operator * (const df1b2matrix& M,const double x);
01178 
01179 df1b2matrix operator - (const df1b2matrix& M,const dmatrix& x);
01180 df1b2matrix operator - (const df1b2matrix& M,const double x);
01181 df1b2matrix operator - (const df1b2matrix& M,const df1b2variable& x);
01182 df1b2matrix operator - (const df1b2matrix& M,const df1b2matrix& x);
01183 df1b2matrix operator - (const dmatrix& M,const df1b2matrix& x);
01184 df1b2matrix operator - (const df1b2variable&, const df1b2matrix& M);
01185 df1b2matrix operator - (const double, const df1b2matrix& M);
01186 
01187 df1b2matrix operator + (const df1b2variable&, const df1b2matrix& M);
01188 
01189 df1b2matrix operator + (const double, const df1b2matrix& M);
01190 df1b2matrix operator * (const double, const df1b2matrix& M);
01191 df1b2matrix elem_prod(const df1b2matrix& x,const df1b2matrix& y);
01192 df1b2matrix elem_prod(const dmatrix& x,const df1b2matrix& y);
01193 df1b2matrix elem_prod(const df1b2matrix& x,const dmatrix& y);
01194 df1b2matrix elem_div(const df1b2matrix& x,const df1b2matrix& y);
01195 df1b2matrix elem_div(const dmatrix& x,const df1b2matrix& y);
01196 df1b2matrix elem_div(const df1b2matrix& x,const dmatrix& y);
01197 
01198 df1b2vector elem_prod(const df1b2vector& x,const df1b2vector& y);
01199 df1b2vector elem_prod(const df1b2vector& x,const dvector& y);
01200 df1b2vector elem_prod(const dvector& x,const df1b2vector& y);
01201 df1b2vector elem_div(const df1b2vector& x,const df1b2vector& y);
01202 df1b2vector elem_div(const dvector& x,const df1b2vector& y);
01203 df1b2vector elem_div(const df1b2vector& x,const dvector& y);
01204 
01205 df1b2vector operator + (const df1b2variable& x,const df1b2vector& y);
01206 df1b2vector operator + (double x,const df1b2vector& y);
01207 df1b2vector operator + (const df1b2vector& x,const df1b2vector& y);
01208 df1b2vector operator + (const df1b2vector& x,const df1b2variable& y);
01209 df1b2vector operator + (const dvector& x,const df1b2vector& y);
01210 df1b2vector operator + (const df1b2vector& y,const dvector& x);
01211 
01212 df1b2vector operator - (const df1b2vector& x,const df1b2vector& y);
01213 df1b2vector operator - (const dvector& x,const df1b2vector& y);
01214 df1b2vector operator - (const df1b2variable& x,const df1b2vector& y);
01215 df1b2vector operator - (const df1b2vector& x,const df1b2vector& y);
01216 df1b2vector operator - (const df1b2variable& x,const df1b2vector& y);
01217 df1b2vector operator - (const df1b2variable& x,const dvector& y);
01218 df1b2vector operator - (const df1b2vector& x,const df1b2variable& y);
01219 df1b2vector operator - (const df1b2vector& x,const dvector& y);
01220 
01221 df1b2vector operator * (const df1b2variable& y,const dvector&);
01222 
01223 df1b2vector operator * (const df1b2vector& x,const df1b2variable& y);
01224 df1b2vector operator * (const df1b2vector& x,double y);
01225 
01226 df1b2vector operator / (const df1b2variable& y,const df1b2vector& x);
01227 df1b2vector operator / (const double& y,const df1b2vector& x);
01228 
01229 df1b2vector operator / (const df1b2vector& x,const df1b2variable& y);
01230 df1b2vector operator / (const df1b2vector& x,double y);
01231 df1b2vector pow(const df1b2vector& x,double y);
01232 df1b2vector pow(const df1b2vector& v,const df1b2variable & x);
01233 df1b2vector pow(const df1b2vector& v,const df1b2vector & x);
01234 df1b2vector pow(const df1b2variable& v,const df1b2vector & x);
01235 df1b2vector pow(df1b2vector const& _x,dvector const& v);
01236 df1b2vector pow(double v,const df1b2vector & x);
01237 df1b2vector pow(const dvector& x,  const df1b2vector& a);
01238 df1b2vector pow(const dvector& x,  const df1b2variable& a);
01239 
01240 df1b2vector operator / (const dvector& x,const df1b2variable& y);
01241 df1b2vector operator + (const dvector& x,const df1b2variable& y);
01242 df1b2vector operator - (const dvector& x,const df1b2variable& y);
01243 df1b2vector operator * (const dvector& x,const df1b2variable& y);
01244 // **************************************
01245 // **************************************
01246 df1b2variable sum(const df1b2vector& x);
01247 df1b2variable mean(const df1b2vector& x);
01248 df1b2variable norm2(const df1b2vector& x);
01249 df1b2variable sumsq(const df1b2vector& x);
01250 df1b2variable norm(const df1b2vector& x);
01251 
01252 df1b2variable norm2(const df1b2matrix& x);
01253 df1b2variable sumsq(const df1b2matrix& x);
01254 df1b2variable norm(const df1b2matrix& x);
01255 df1b2variable mean(const df1b2matrix& x);
01256 df1b2variable sum(const df1b2matrix& x);
01257 df1b2matrix square(const df1b2matrix& x);
01258 df1b2vector square(const df1b2vector& x);
01259 df1b2vector cube(const df1b2vector& x);
01260 df1b2vector solve(const df1b2matrix& x,const dvector&);
01261 df1b2vector solve(const df1b2matrix& x,const df1b2vector&);
01262 
01263 int size_count(const df1b2matrix& x);
01264 int size_count(const df1b2vector& x);
01265 
01266 df1b2vector first_difference(const df1b2vector& x);
01267 
01268 df1b2vector operator * (const dvector& ,const df1b2matrix& );
01269 df1b2vector operator * (const df1b2vector& ,const df1b2matrix& );
01270 df1b2vector operator * (const df1b2vector& ,const dmatrix& );
01271 df1b2variable operator * (const df1b2vector& x,const df1b2vector& y);
01272 df1b2vector operator * (const df1b2matrix& M,const df1b2vector& x);
01273 df1b2vector operator * (const df1b2matrix& M,const dvector& x);
01274 void check_shape(const df1b2vector & x,const df1b2vector & y,const char * s);
01275 void check_shape(const df1b2vector & x,const dvector & y,const char * s);
01276 void check_shape(const dvector & x,const df1b2vector & y,const char * s);
01277 void check_shape(const df1b2vector & x,const df1b2matrix & y,const char * s);
01278 
01279 void checkidentiferstring(const char * ids,test_smartlist& list);
01280 
01281 double& value(const df1b2variable& _x);
01282 dvector value(const df1b2vector& _x);
01283 dmatrix value(const df1b2matrix& _x);
01284 
01289 class initial_df1b2params
01290 {
01291 public:
01292   static int stddev_scale(const dvector& d,const dvector& x);
01293   static int current_phase;
01294   static int separable_flag;
01295   static int separable_calculation_type;
01296   static int have_bounded_random_effects;
01297   int phase_start;
01298   int phase_save;
01299   int ind_index;
01300   double scalefactor;
01301 
01302   static void save_varsptr(void);
01303   static double cobjfun;
01304   static void restore_varsptr(void);
01305 #if defined(__x86_64) || (defined(_MSC_VER) && defined(_M_X64))
01306   static lmatrix* pointer_table;
01307 #else
01308   static imatrix* pointer_table;
01309 #endif
01310   static initial_df1b2params **  varsptr;  // this should be a resizeable
01311   static initial_df1b2params **  varsptr_sav;  // this should be a resizeable
01312   static int num_initial_df1b2params;         // array
01313   static int num_initial_df1b2params_sav;         // array
01314   virtual void set_value(const dvector&,const int& ii)=0;
01315   virtual void set_index(int aflag,const int& ii)=0;
01316   static int set_index(void);
01317   virtual void set_value(const init_df1b2vector&,const int& ii,
01318     const df1b2variable&)=0;
01319   static void reset(const init_df1b2vector&,const df1b2variable&);
01320   static void reset_all(const dvector&);
01321   static void reset(const df1b2vector&,const df1b2variable&);
01322   virtual void add_to_list(void);
01323   initial_df1b2params(void);
01324   void set_phase_start(int n){phase_start=n; phase_save=n;}
01325   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii)=0;
01326   double get_scalefactor();
01327   void set_scalefactor(const double);
01328 };
01329 
01330 typedef initial_df1b2params * P_INITIAL_DF1B2PARAMS;
01331 
01336 class do_naught_kludge
01337 {
01338 };
01339 
01344 class df1b2_init_vector : public df1b2vector, public initial_df1b2params
01345 {
01346 public:
01347   virtual ~df1b2_init_vector() {;}
01348 
01349   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01350   void set_phase_start(int n){phase_start=n; phase_save=n;}
01351   virtual void set_value(const dvector& x,const int& _ii);
01352   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01353     const df1b2variable& pen);
01354   void allocate(int min, int max,int phase_start, const adstring& s="unnamed");
01355   void allocate(int min, int max,const adstring& s="unnamed");
01356   virtual void set_index(int aflag,const int& ii);
01357 };
01358 
01363 class df1b2_init_matrix : public df1b2matrix, public initial_df1b2params
01364 {
01365 public:
01366   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01367 
01368   inline void allocate(int rmin,int rmax,const index_type& cmin,
01369     const index_type& cmax,int ps,const char * s)
01370   {
01371     set_phase_start(ps);
01372     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01373   }
01374 
01375   inline void allocate(int rmin,int rmax,const index_type& cmin,
01376     const index_type& cmax,const char * s)
01377   {
01378     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01379   }
01380 
01381   inline void allocate(int rmin,int rmax,int cmin,int cmax)
01382   {
01383     df1b2matrix::allocate(rmin,rmax,cmin,cmax);
01384   }
01385   inline void allocate(int rmin,int rmax,int cmin,int cmax,int ps,
01386     const char * s)
01387   {
01388     set_phase_start(ps);
01389     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01390   }
01391   inline void allocate(int rmin,int rmax,int cmin,int cmax,const char * s)
01392   {
01393     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01394   }
01395   inline void allocate(int rmin, int rmax, int cmin, int cmax, double, double,
01396     const char* s)
01397   {
01398     df1b2matrix::allocate(rmin,rmax,cmin,cmax,s);
01399   }
01400   void set_phase_start(int n){phase_start=n; phase_save=n;}
01401   virtual void set_value(const dvector& x,const int& _ii);
01402   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01403     const df1b2variable& pen);
01404   virtual void set_index(int aflag,const int& ii);
01405 };
01406 
01407 
01412 class df1b2_init_number : public df1b2variable, public initial_df1b2params
01413 {
01414   static do_naught_kludge do_naught_kludge_a;
01415 public:
01416   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01417   virtual int & get_ind_index(void) { return ind_index; }
01418   void set_phase_start(int n){phase_start=n; phase_save=n;}
01419   virtual void set_value(const dvector& x,const int& _ii);
01420   virtual void set_value(const init_df1b2vector&,const int& ii,
01421     const df1b2variable&);
01422   virtual void set_index(int aflag,const int& ii);
01423 
01424   void allocate(const adstring&);
01425   void allocate(int,const adstring&);
01426   df1b2_init_number();
01427   void operator = (const df1b2variable& _x);
01428 };
01429 
01434 class df1b2_init_number_vector
01435 {
01436   df1b2_init_number * v;
01437   int index_min;
01438   int index_max;
01439   double_index_type * it;
01440 public:
01441   df1b2_init_number_vector();
01442 
01443 #if defined(OPT_LIB)
01444    df1b2_init_number& operator [] (int i) { return v[i];}
01445    df1b2_init_number& operator () (int i) { return v[i];}
01446 #else
01447    df1b2_init_number& operator [] (int i);
01448    df1b2_init_number& operator () (int i);
01449 #endif
01450   void allocate(int min1,int max1,const index_type& phase_start,
01451     const char * s);
01452 
01453   void allocate(int min1,int max1,const char * s);
01454 
01455   int allocated(void) { return (v!=NULL); }
01456   int indexmin(void) {return (index_min);}
01457   int indexmax(void) {return (index_max);}
01458   ~df1b2_init_number_vector();
01459   void deallocate(void);
01460 };
01461 
01462 class df1b2_init_bounded_number;
01463 
01468 class df1b2_init_bounded_number :public df1b2_init_number
01469 {
01470   double minb;
01471   double maxb;
01472 public:
01473   void allocate(double _minb,double _maxb,int _phase_start,
01474     const char * s="UNNAMED");
01475   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01476 
01477   void allocate(double _minb,double _maxb,const char * s="UNNAMED");
01478   //void allocate(const ad_double& _minb,const ad_double& _maxb,
01479    // const char * s="UNNAMED");
01480 
01481   virtual void set_value(const dvector& x,const int& _ii);
01482   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01483     const df1b2variable& pen);
01484 };
01485 
01490 class df1b2_init_bounded_number_vector
01491 {
01492   df1b2_init_bounded_number * v;
01493   int index_min;
01494   int index_max;
01495   double_index_type * it;
01496 public:
01497   //virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01498   df1b2_init_bounded_number_vector();
01499 
01500 #if defined(OPT_LIB)
01501    df1b2_init_bounded_number& operator [] (int i) { return v[i];}
01502    df1b2_init_bounded_number& operator () (int i) { return v[i];}
01503 #else
01504    df1b2_init_bounded_number& operator [] (int i);
01505    df1b2_init_bounded_number& operator () (int i);
01506 #endif
01507   void allocate(int min1,int max1, const double_index_type& bmin,
01508    const double_index_type& bmax, const char * s);
01509 
01510   void allocate(int min1,int max1, const double_index_type& bmin,
01511    const double_index_type& bmax, const index_type& phase_start,
01512     const char * s);
01513 
01514 
01515   int allocated(void) { return (v!=NULL); }
01516   int indexmin(void) {return (index_min);}
01517   int indexmax(void) {return (index_max);}
01518   ~df1b2_init_bounded_number_vector();
01519   void deallocate(void);
01520 };
01521 
01526 class df1b2_init_bounded_vector : public df1b2_init_vector
01527 {
01528 protected:
01529   double minb;
01530   double maxb;
01531 public:
01532   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01533   void allocate(int,int,double,double,int,const char *);
01534   void allocate(int,int,double,double,const char *);
01535   void allocate(int,int,double,double);
01536   //void allocate(int,int);
01537   virtual void set_value(const dvector& x,const int& _ii);
01538   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01539     const df1b2variable& pen);
01540   inline double getminb(void){return minb;}
01541   inline double getmaxb(void){return maxb;}
01542 };
01543 
01544 class re_df1b2_init_bounded_vector;
01545 
01550 class random_effects_bounded_vector_info
01551 {
01552   re_df1b2_init_bounded_vector * pv;
01553   int i;
01554 public:
01555   random_effects_bounded_vector_info
01556     ( re_df1b2_init_bounded_vector * _pv,int _i) : pv(_pv), i(_i) {}
01557   friend class funnel_init_df1b2variable;
01558   friend class df1b2variable;
01559   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01560     const df1b2variable& pen);
01561 };
01562 
01567 class re_df1b2_init_bounded_vector : public df1b2_init_bounded_vector
01568 {
01569 public:
01570   random_effects_bounded_vector_info operator () (int i)
01571   {
01572     return random_effects_bounded_vector_info(this,i);
01573   }
01574   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01575     const df1b2variable& pen);
01576   virtual void set_value(const dvector& x,const int& _ii);
01577   re_df1b2_init_bounded_vector(void);
01578 };
01579 
01584 class df1b2_init_bounded_matrix : public df1b2_init_matrix
01585 {
01586 protected:
01587   double minb;
01588   double maxb;
01589 public:
01590   virtual void sd_scale(const dvector& d,const dvector& x,const int& ii);
01591   void allocate(int,int,int,int,double,double,int,const char *);
01592   void allocate(int,int,const index_type& ,const index_type& ,double,double,int,
01593     const char *);
01594   void allocate(int,int,int,int,double,double,const char *);
01595   void allocate(int,int,const index_type& ,const index_type& ,double,double,
01596     const char *);
01597   void allocate(int,int,int,int,double,double);
01598   //void allocate(int,int);
01599   virtual void set_value(const dvector& x,const int& _ii);
01600   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01601     const df1b2variable& pen);
01602 };
01603 
01608 class df1b2_init_bounded_dev_vector : public df1b2_init_bounded_vector
01609 {
01610 public:
01611   virtual void set_value(const dvector& x,const int& _ii);
01612   virtual void set_value(const init_df1b2vector& x,const int& _ii,
01613     const df1b2variable& pen);
01614 };
01615 
01616 void set_value(const df1b2variable& u,const init_df1b2vector& x,const int& _ii,
01617   double fmin,double fmax,const df1b2variable& fpen);
01618 
01619 void set_value(const df1b2variable& u,const init_df1b2vector& x,const int& _ii,
01620   double fmin,double fmax,const df1b2variable& fpen,double s);
01621 
01622 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax);
01623 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax,
01624   const df1b2variable& _fpen);
01625 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax,
01626   const df1b2variable& _fpen,double s);
01627 
01628   ostream& operator << (const ostream& _os, const df1b2variable& _x);
01629   ostream& operator << (const ostream& _os, const df1b2vector& _x);
01630   ostream& operator << (const ostream& _os, const df1b2matrix& _x);
01631 
01632   class re_objective_function_value : public df1b2variable
01633   {
01634   public:
01635     static re_objective_function_value * pobjfun;
01636     static double fun_without_pen;
01637     re_objective_function_value(void);
01638     ~re_objective_function_value();
01639     re_objective_function_value& operator = (const df1b2variable& v);
01640     re_objective_function_value& operator = (double v);
01641     void allocate(void);
01642     void allocate(const char *);
01643   };
01644 
01645 df1b2variable posfun(const df1b2variable&x,const double eps);
01646 df1b2variable posfun2(const df1b2variable&x,const double eps,
01647   const df1b2variable& _pen);
01648 df1b2variable posfun(const df1b2variable&x,const double eps,
01649   const df1b2variable& _pen);
01650 
01655 inline df1b2variable operator - (const df1b2variable& x)
01656 {
01657   return -double(1.0)*x;
01658 }
01659 
01664 inline df1b2vector operator - (const df1b2vector& x)
01665 {
01666   return -double(1.0)*x;
01667 }
01668 
01669 df1b2variable gamma_density(const df1b2variable& _x,double r, double mu);
01670 df1b2variable gamma_density(const df1b2variable& _x,const df1b2variable& _r,
01671   const  df1b2variable& _mu);
01672 
01673 df1b2variable log_gamma_density(const df1b2variable& _x,double r, double mu);
01674 df1b2variable log_gamma_density(const df1b2variable& _x,const df1b2variable& _r,
01675   const  df1b2variable& _mu);
01676 
01677 
01678 df1b2variable log_negbinomial_density(double x,const df1b2variable& mu,
01679   const df1b2variable& tau);
01680 
01681 df1b2variable log_density_poisson(double x,const df1b2variable& mu);
01682 
01683 ostream& operator << (const ostream& ostr,const df1b2variable& z);
01684 ostream& operator << (const ostream& ostr,const df1b2vector& z);
01685 ostream& operator << (const ostream& ostr,const df1b2matrix& z);
01686 ostream& operator << (const ostream& ostr,const df1b2_init_number_vector& z);
01687 ostream& operator << (const ostream& ostr,const init_df1b2vector& z);
01688 ostream& operator << (const ostream& ostr,
01689   const df1b2_init_bounded_number_vector& z);
01690 
01691 class adkludge1;
01692 
01693 
01698   class df1b2_header_ptr_vector
01699   {
01700     int index_min;
01701     int index_max;
01702     df1b2_header ** v;
01703   public:
01704     df1b2_header_ptr_vector(int mmin,int mmax);
01705     ~df1b2_header_ptr_vector();
01706     inline df1b2_header * & operator () (int i) { return *(v+i); }
01707     inline df1b2_header * & operator [] (int i) { return *(v+i); }
01708     int indexmin(void) {return index_min;}
01709     int indexmax(void) {return index_max;}
01710   };
01711 
01716   class double_ptr_vector
01717   {
01718     int index_min;
01719     int index_max;
01720     double ** v;
01721   public:
01722     double_ptr_vector(int mmin,int mmax);
01723     ~double_ptr_vector();
01724     inline double * & operator () (int i) { return *(v+i); }
01725     inline double * & operator [] (int i) { return *(v+i); }
01726     int indexmin(void) {return index_min;}
01727     int indexmax(void) {return index_max;}
01728   };
01729 
01730 int active(const df1b2_init_vector& x);
01731 int active(const df1b2_init_number& x);
01732 int active(const df1b2_init_matrix& x);
01733 //int active(const df1b2_init_bounded_vector& x);
01734 
01735 //int active (const df1b2_init_vector &);
01736 //int active(const initial_df1b2params& x);
01737 //int active(const df1b2vector& x);
01738 //int active(const df1b2matrix& x);
01739 
01740 double evaluate_function_quiet(const dvector& x,function_minimizer * pfmin);
01741 double evaluate_function(const dvector& x,function_minimizer * pfmin);
01742 double evaluate_function(double& fval,const dvector& x,
01743   function_minimizer * pfmin);
01744 double evaluate_function(double& fval,const dvector& x,const dvector&,
01745   function_minimizer * pfmin);
01746 
01747 
01748 double evaluate_function_no_derivatives(const dvector& x,
01749   function_minimizer * pfmin);
01750 
01751 #include <df1b2fnl.h>
01752 
01753 int allocated(const df1b2_init_vector&);
01754 int allocated(const df1b2vector&);
01755 int allocated(const df1b2_init_matrix&);
01756 #include <df3fun.h>
01757 
01758 /*
01759 df1b2variable betai(const df1b2variable & _a, const df1b2variable & _b,
01760          const df1b2variable & _x);
01761 df1b2variable betai(const df1b2variable & _a, const df1b2variable & _b,
01762          double _x);
01763 */
01764 
01765 df1b2variable betai(const df1b2variable& a, const df1b2variable& b,
01766   double x, int maxit=100);
01767 
01768 double do_gauss_hermite_block_diagonal(const dvector& x,
01769   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01770   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01771   function_minimizer * pmin);
01772 
01773 double do_gauss_hermite_block_diagonal_multi(const dvector& x,
01774   const dvector& u0,const dmatrix& Hess,const dvector& _xadjoint,
01775   const dvector& _uadjoint,const dmatrix& _Hessadjoint,
01776   function_minimizer * pmin);
01777 
01778 double calculate_importance_sample_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 calculate_importance_sample_block_diagonal_option2(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_option_antithetical
01789   (const dvector& x,const dvector& u0,const dmatrix& Hess,
01790   const dvector& _xadjoint,const dvector& _uadjoint,
01791   const dmatrix& _Hessadjoint,function_minimizer * pmin);
01792 
01793 double calculate_importance_sample_block_diagonal_funnel(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(const dvector& x,const dvector& u0,
01799   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01800   const dmatrix& _Hessadjoint,function_minimizer * pmin);
01801 
01802 double calculate_importance_sample_shess(const dvector& x,const dvector& u0,
01803   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
01804   const dmatrix& _Hessadjoint,function_minimizer * pmin);
01805 double calculate_importance_sample(const dvector& x,const dvector& u0,
01806   const banded_symmetric_dmatrix& bHess,const dvector& _xadjoint,
01807   const dvector& _uadjoint,
01808   const banded_symmetric_dmatrix& _bHessadjoint,function_minimizer * pmin);
01809 
01810 dvector evaluate_function_with_quadprior(const dvector& x,int usize,
01811   function_minimizer * pfmin);
01812 
01817 class style_flag_class
01818 {
01819 public:
01820   int old_style_flag;
01821   virtual void set_old_style_flag(void)=0;
01822   //style_flag_class(void) {set_old_style_flag();}
01823 };
01824 
01829 class quadratic_prior : public style_flag_class
01830 {
01831   dmatrix * pMinv;
01832   dvar_matrix * dfpMinv;
01833   dvar_vector * pu;
01834   int xmyindex;
01835 public:
01836   static int qflag;
01837   static quadratic_prior * ptr[]; // this should be a resizeable array
01838   static void get_M_calculations(void);
01839   static void cleanup_pMinv();
01840   static void cleanup_dfpMinv();
01841   static int num_quadratic_prior;
01842   static const int max_num_quadratic_prior;
01843   void add_to_list(void);
01844 public:
01845   static int in_qp_calculations;
01846   int get_offset(int xs);
01847   int get_myindex(void) { return xmyindex;}
01848   static quadratic_prior * get_ptr(int i){ return ptr[i];}
01849   void operator = (const dvar_matrix &);
01850   void operator = (const dmatrix &);
01851   static int get_in_qp_calculations() { return in_qp_calculations; }
01852   static int get_num_quadratic_prior(void) { return num_quadratic_prior;}
01853   static dvariable get_quadratic_priors(void);
01854   dvariable get_function(void);
01855   quadratic_prior(void);
01856   ~quadratic_prior(void);
01857   void allocate(const dvar_vector & _u,const char * s);
01858   void allocate(const dvar_vector & _u);
01859   void allocate(const dvar_matrix & _M, const dvar_vector & _u,const char * s);
01860   void allocate(const dvar_matrix & _M, const dvar_vector & _u);
01861   //dmatrix get_cHessian(void);
01862   void get_cHessian(dmatrix,int);
01863   void get_cHessian(dvar_matrix,int);
01864   void get_cHessian_from_vHessian(dmatrix ,int);
01865   void get_vHessian(dvar_matrix ,int);
01866   void get_cgradient(dvector,int);
01867   dvar_matrix get_Hessian(void);
01868   dvar_vector get_gradient(void);
01869   static dvar_vector get_gradient_contribution(void);
01870   static dvar_matrix get_Hessian_contribution(void);
01871   static void get_cgradient_contribution(dvector,int);
01872   static void get_cHessian_contribution(dmatrix,int);
01873   static void get_vHessian_contribution(dvar_matrix  ,int );
01874   static void get_cHessian_contribution_from_vHessian(dmatrix,int);
01875   friend class df1b2quadratic_prior;
01876   virtual void get_cM(void)=0;
01877 };
01878 
01883 class df1b2quadratic_prior : public style_flag_class
01884 {
01885   ivector * index;
01886   df1b2matrix * M;
01887   dmatrix * CM;
01888 public:
01889   dmatrix * Lxu;
01890   df1b2_init_vector * pu;
01891   int xmyindex;
01892   static df1b2quadratic_prior * ptr[]; // this should be a resizeable array
01893   static int num_quadratic_prior;
01894   static const int max_num_quadratic_prior;
01895   void add_to_list(void);
01896 public:
01897   static df1b2quadratic_prior * get_ptr(int i){ return ptr[i];}
01898   unsigned int num_active_parameters;
01899   unsigned int get_num_active_parameters(void) { return num_active_parameters; }
01900   int get_myindex(void) { return xmyindex;}
01901   void operator = (const df1b2matrix &);
01902   void operator = (const dmatrix &);
01903   static int get_num_quadratic_prior(void) { return num_quadratic_prior;}
01904   static dvariable get_quadratic_priors(void);
01905   df1b2variable get_function(void);
01906   df1b2quadratic_prior(void);
01907   ~df1b2quadratic_prior(void);
01908   void allocate(const df1b2_init_vector & _u,const char * s);
01909   void allocate(const df1b2_init_vector & _u);
01910   void allocate(const df1b2matrix & _M, const df1b2_init_vector & _u,
01911     const char * s);
01912   void allocate(const df1b2matrix & _M, const df1b2_init_vector & _u);
01913   void allocate(const dvar_matrix & _M, const dvar_vector & _u,const char * s);
01914   void allocate(const dvar_matrix & _M, const dvar_vector & _u);
01915   dmatrix get_cHessian(void);
01916   dvector get_cgradient(void);
01917   dvar_matrix get_Hessian(void);
01918   dvar_vector get_gradient(void);
01919   static dvar_vector get_gradient_contribution(void);
01920   static dvar_matrix get_Hessian_contribution(void);
01921   static dvector get_cgradient_contribution(void);
01922   static dmatrix get_cHessian_contribution(void);
01923   static void get_Lxu_contribution(dmatrix&);
01924   virtual void get_Lxu(dmatrix&)=0;
01925   friend class quadratic_prior;
01926   friend class df1b2_parameters;
01927 };
01928 
01933 class normal_quadratic_prior : public quadratic_prior
01934 {
01935   virtual void set_old_style_flag(void);
01936 public:
01937   normal_quadratic_prior(void);
01938   void operator = (const dvar_matrix & M);
01939 };
01940 
01945 class normal_df1b2quadratic_prior : public df1b2quadratic_prior
01946 {
01947   virtual void set_old_style_flag(void);
01948 public:
01949   void operator = (const df1b2matrix & M);
01950   normal_df1b2quadratic_prior(void);
01951 };
01952 
01957 class quadratic_re_penalty : public quadratic_prior
01958 {
01959   virtual void set_old_style_flag(void);
01960 public:
01961   quadratic_re_penalty(void);
01962   void operator = (const dvar_matrix & M);
01963   void operator = (const dmatrix & M);
01964 };
01965 
01970 class df1b2quadratic_re_penalty : public df1b2quadratic_prior
01971 {
01972   virtual void set_old_style_flag(void);
01973 public:
01974   void operator = (const df1b2matrix & M);
01975   void operator = (const dmatrix & M);
01976   df1b2quadratic_re_penalty(void);
01977 };
01978 
01979 dvar_vector solve(const named_dvar_matrix &,const random_effects_vector&);
01980 
01985 class constant_quadratic_re_penalty : public quadratic_prior
01986 {
01987   virtual void set_old_style_flag(void);
01988 public:
01989   constant_quadratic_re_penalty(void);
01990   void operator = (const dmatrix & M);
01991 };
01992 
01997 class constant_df1b2quadratic_re_penalty : public df1b2quadratic_prior
01998 {
01999   virtual void set_old_style_flag(void);
02000 public:
02001   void operator = (const dmatrix & M);
02002   constant_df1b2quadratic_re_penalty(void);
02003 };
02004 
02005 df1b2vector solve(df1b2matrix& M,df1b2vector& v,const df1b2variable& ln_det);
02006 
02007 df1b2vector solve(df1b2matrix& M,df1b2vector& v,const df1b2variable& ln_det,
02008   const int& sgn);
02009 df1b2vector lower_triangular_solve(const df1b2matrix& m,const df1b2vector& v);
02010 df1b2vector lower_triangular_solve_trans(const df1b2matrix& m,
02011   const df1b2vector& v);
02012 
02013 //df1b2variable ln_det(df1b2matrix& M);
02014 // line above replaced with line below based on issue #37
02015 df1b2variable ln_det(const df1b2matrix & m1);
02016 
02017 df1b2variable ln_det(df1b2matrix& M,int & sgn);
02018 
02019 //df1b2vector solve(df1b2matrix& M,df1b2vector& v);
02020 
02021 df1b2matrix expm(const df1b2matrix & A);
02022 df1b2matrix solve(const df1b2matrix& aa,const df1b2matrix& tz,
02023   df1b2variable ln_unsigned_det,df1b2variable& sign);
02024 df1b2matrix solve(const df1b2matrix& aa,const df1b2matrix& tz);
02025 df1b2vector solve(const df1b2matrix& aa,const df1b2vector& z,
02026   const df1b2variable& ld, df1b2variable& sign);
02027 
02028 void check_pool_depths();
02029 
02030 df1b2variable lower_triangular_ln_det(const df1b2matrix& m);
02031 df1b2variable lower_triangular_ln_det(const df1b2matrix& m,int& sgn);
02032 df1b2variable bounder(const df1b2variable&  x,double min,double max,
02033   double scale);
02034 
02035 df1b2variable inv_cumd_beta_stable(const df1b2variable& a,
02036   const df1b2variable& b,const df1b2variable& x,double eps=1.e-7);
02037 
02038 df1b2variable bounded_cumd_norm(const df1b2variable& _x,double beta);
02039 df1b2variable cumd_norm(const df1b2variable& _x);
02040 df1b2variable inv_cumd_exponential(const df1b2variable& y);
02041 df1b2variable cumd_exponential(const df1b2variable& x);
02042 extern int make_sure_df1b2fun_gets_linked;
02043 
02044 df1b2variable inv_cumd_normal_mixture(const df1b2variable& _x,double _a);
02045 df1b2variable inv_cumd_normal_logistic_mixture(const df1b2variable& _x,
02046   double _a);
02047 
02048 df1b2variable beta_deviate(const df1b2variable& _x,const df1b2variable& _a,
02049   const df1b2variable& _b,double eps=1.e-7);
02050 df1b2variable robust_normal_logistic_mixture_deviate(const df1b2variable& x,
02051   double spread=3.0);
02052 df1b2variable robust_normal_mixture_deviate(const df1b2variable& x,
02053   double spread=3.0);
02054 
02055 df1b2variable gamma_deviate(const df1b2variable& _x,const df1b2variable& _a);
02056 df1b2variable inv_cumd_gamma(const df1b2variable& _y,const df1b2variable& _a);
02057 double inv_cumd_gamma(double y,double _a);
02058 
02059 df1b2variable inv_cumd_cauchy(const df1b2variable& n);
02060 df1b2variable inv_cumd_t(const df1b2variable& n,const df1b2variable&  u,
02061   double eps=1.e-7);
02062 
02067   class df1b2function_tweaker
02068   {
02069     double mult;
02070     double eps;
02071     dvector coffs;
02072   public:
02073     df1b2function_tweaker(double eps,double mult);
02074     df1b2variable operator () (const df1b2variable&);
02075   };
02076 
02077 df1b2vector posfun(const df1b2vector& x,double eps,
02078   const df1b2variable& _pen);
02079 
02080 df1b2variable norm_to_gamma(const df1b2variable & v,const df1b2variable& alpha,
02081   double bound=0.999999);
02082 
02083 void print_is_diagnostics(laplace_approximation_calculator *lapprox);
02084 
02085   ofstream &  operator << (const ofstream& _s,const nested_calls_shape& m);
02086 
02087   df1b2variable log_der_logistic(double a,double b,const df1b2variable& x);
02088   df1b2variable logistic(double a,double b,const df1b2variable& x);
02089   df1b2variable dflogistic(double a,double b,const df1b2variable& x);
02090 
02091   df1b2variable dfposfun(const df1b2variable& x,const double eps);
02092   void ADMB_getcallindex(const df1b2variable& x);
02093   void ADMB_getcallindex(const df1b2vector& x);
02094   void ADMB_getcallindex(const df1b2matrix& x);
02095 
02096 df1b2variable asin(const df1b2variable& xx);
02097 #endif //!defined(__DF1B2FUN__)