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