ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
fvar.hpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar.hpp 1919 2014-04-22 22:02:01Z 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 "New 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  */
00041 #ifndef FVAR_HPP
00042 #define FVAR_HPP
00043 
00050 #include <math.h>
00051 // Borrow definition of M_PI from GCC
00052 #ifndef M_PI
00053 #   define M_PI 3.14159265358979323846
00054 #endif
00055 #ifndef PI
00056 #   define PI 3.14159265358979323846
00057 #endif
00058 
00059 #if defined(__GNUC__) && (__GNUC__ < 3)
00060   #pragma interface
00061 #endif
00062 
00063 #if defined(THREAD_EXPERIMENT)
00064 #   define THREAD_SAFE
00065 #endif
00066 
00067 #if defined(THREAD_SAFE)
00068 #   include <pthread.h>
00069 #endif
00070 
00071 #define USE_VECTOR_SHAPE_POOL
00072 
00073 #if defined(USE_DDOUBLE)
00074 #   include <qd/qd.h>
00075 #   define  double dd_real
00076     dd_real operator ++(dd_real & d)
00077     {
00078        return d += 1.0;
00079     }
00080 
00081     dd_real pow(const dd_real & d, int i);
00082     dd_real pow(const dd_real & d, const dd_real &);
00083 #endif
00084 
00085 #if defined(__TURBOC__)
00086 #   define __STL_PTHREADS
00087 #   include <k3stl.h>
00088 #   include <pthread_alloc>
00089 #endif
00090 
00091 #define  __NUMBERVECTOR__
00092 
00096 #define ADUNCONST(type,obj) type & obj = (type&) _##obj;
00097 
00098 #define my_off_t long int
00099 
00100 #ifdef __BORLANDC__
00101 #   define my_u_off_t unsigned long int
00102 #else
00103 #   define my_u_off_t unsigned long long int
00104 #endif
00105 
00106 #define  MFCL2_CONSTRUCTORS
00107 
00108 #if defined(_ADEXEP)
00109 #   define USE_EXECPTIONS
00110 #endif
00111 
00112 #define __USE_IOSTREAM__
00113 
00114 #if defined(__BORLANDC__)
00115   #if (__BORLANDC__  >= 0x0550)
00116     #include <fcntl.h>
00117   #endif
00118 #endif
00119 
00120 #define USE_HIGHER_ARRAYS
00121 
00122 #if defined(__BORLANDC__) || defined (_MSC_VER) || defined(__WAT32__)
00123 #   include <io.h>
00124 #endif
00125 
00126 #if !defined(AD_LONG_INT)
00127   #define AD_LONG_INT long int
00128 #endif
00129 
00130 #if !defined(_MSC_VER)
00131   #include <unistd.h>
00132   #include <fcntl.h> // to get fstreambase
00133   #if !defined(O_BINARY)
00134     #define O_BINARY 0
00135   #endif
00136   #define __GNU__
00137 #endif
00138 
00139 #if defined(__DJGPP__)
00140   #include <dos.h>
00141   #include <pc.h>
00142 #endif
00143 
00144 #if !defined(CGNU_DOUBLE)
00145   #define CGNU_DOUBLE double
00146 #endif
00147 
00148 #ifdef _SNGL_
00149    #define double float
00150 #endif
00151 
00152 #ifndef NO_DERIVS
00153   #define  NO_DERIVS
00154 #endif
00155 
00156 #ifndef OPT_LIB
00157   #ifndef SAFE_ARRAYS
00158   #define SAFE_ARRAYS
00159   #endif
00160 #endif
00161 
00162 // C language function prototypes
00163 extern "C"
00164 {
00165    typedef int (*fptr) (const char *format, ...);
00166    extern fptr ad_printf;
00167    typedef void (*exitptr) (int);
00168    extern exitptr ad_exit;
00169 
00170    void spdll_exit(int);
00171 }
00172 
00177 class smart_counter
00178 {
00179   int *ncopies;
00180 public:
00181   int *get_ncopies();
00182   smart_counter();
00183   smart_counter(const smart_counter& sc);
00184   ~smart_counter();
00185 };
00186 
00191 class double_and_int
00192 {
00193  public:
00195    double x;
00199    inline double &xvalue(void)
00200    {
00201       return x;
00202    }
00203 };
00204 
00205 
00206 // "forward" class definitions
00207 class banded_symmetric_dvar_matrix;
00208 class banded_symmetric_dmatrix;
00209 class banded_lower_triangular_dmatrix;
00210 class banded_lower_triangular_dvar_matrix;
00211 class random_number_generator;
00212 class ad_integer;
00213 class index_type;
00214 class dlink;
00215 class dvector;
00216 class dmatrix;
00217 class tdmatrix;
00218 class imatrix;
00219 class prevariable;
00220 class dvariable;
00221 class dvar_vector;
00222 class ivector;
00223 class lvector;
00224 class dvar_matrix;
00225 class uostream;
00226 class uistream;
00227 class arr_link;
00228 class arr_list;
00229 class d3_array;
00230 class d4_array;
00231 class d5_array;
00232 class d6_array;
00233 class d7_array;
00234 class dvar3_array;
00235 class dvar4_array;
00236 class dvar5_array;
00237 class dvar6_array;
00238 class dvar7_array;
00239 class dlist;
00240 class indvar_offset_list;
00241 class dvar_matrix_position;
00242 class DF_FILE;
00243 class dmatrix_position;
00244 class d3_array_position;
00245 class dvector_position;
00246 class vector_shape;
00247 class ivector_position;
00248 class kkludge_object;
00249 class dvar_vector_position;
00250 dvector restore_dvar_vector_derivatives(void);
00251 class gradient_structure;
00252 class dependent_variables_information;
00253 class vector_shapex;
00254 class predvar_vector;
00255 class independent_variables;
00256 
00257 #if defined(__GNUC__)
00258 #   if (__GNUC__  >= 3)
00259 #      include <fstream>
00260 #   else
00261 #      include <fstream.h>
00262 #   endif
00263 #elif defined(_MSC_VER)
00264 #   if (_MSC_VER >= 1300)
00265 #      include <fstream>
00266 #   else
00267 #      include <fstream.h>
00268 #   endif
00269 #else
00270 #   include <fstream.h>
00271 #endif
00272 
00273 #include <stdio.h>
00274 
00275 #if defined(__BORLANDC__)
00276 #   if (__BORLANDC__  < 0x0550)
00277 #      include <fcntl.h> // to get fstreambase
00278 #      include <strstrea.h>
00279 #   else
00280 #      include <strstream>
00281 #   endif
00282 #endif
00283 
00284 #if defined(_MSC_VER)
00285 #   if (_MSC_VER < 1300)
00286 #     include <iostream.h>
00287 #     include <strstrea.h>
00288 #   else
00289 #     include <iostream>
00290 #   endif
00291 #   include <stddef.h>
00292 #   include <fcntl.h> // to get fstreambase
00293 #   include <dos.h>
00294 #   undef __ZTC__
00295 #   undef __SUN__
00296 #endif
00297 
00298 #ifdef __ZTC__
00299 #   include <fstream.hpp> // to get fstreambase
00300 #   include <iostream.hpp>
00301 #   include <dos.h>
00302 #   undef __SUN__
00303 #endif
00304 
00305 #ifdef __SUN__
00306 #   undef __NDPX__
00307 #   include <fstream.h>
00308 #   include <iostream.h>
00309 #   ifndef _FPOS_T_DEFINED
00310 #      ifdef __GNUDOS__
00311 #         if defined(linux) || defined(__ADSGI__)|| defined(__CYGWIN32__)
00312              typedef long int fpos_t;
00313 #         else
00314              typedef unsigned long int fpos_t;
00315 #         endif
00316 #         undef __SUN__
00317           typedef int fpos_t;
00318 #      endif
00319 #      define _FPOS_T_DEFINED
00320 #   endif
00321 #endif
00322 
00323 #if defined(__BORLANDC__)
00324 #   if (__BORLANDC__  > 0x0520)
00325        using std::ofstream;
00326        using std::ifstream;
00327 #   endif
00328 #endif
00329 
00330 #if defined(__cplusplus)
00331   #include <iostream>
00332   #include <iomanip>
00333   #include <sstream>
00334   #include <istream>
00335   #include <sstream>
00336   using std::ofstream;
00337   using std::ostream;
00338   using std::ifstream;
00339   using std::istream;
00340   using std::istringstream;
00341   using std::streampos;
00342   using std::streambuf;
00343   using std::setw;
00344   using std::setprecision;
00345   using std::istringstream;
00346   using std::ios;
00347   using std::cerr;
00348   using std::cin;
00349   using std::cout;
00350   using std::endl;
00351 #else
00352   #include <iomanip.h>
00353   #include <strstream.h>
00354 #endif
00355 
00356 #define BEGIN_MINIMIZATION(nvar, objective_function, ind_vars, gradient, \
00357 cntrl) \
00358 gradient_structure gs; \
00359 while (cntrl.ireturn >= 0) \
00360 { \
00361   cntrl.fmin(objective_function,ind_vars,gradient ); \
00362   if (cntrl.ireturn > 0) \
00363   {
00364 #define END_MINIMIZATION(nvar,gradient) \
00365     gradcalc(nvar, gradient); \
00366   } \
00367 }
00368 
00369 void default_evaluation4ind(void);
00370 void default_evaluation(void);
00371 void default_evaluation0(void);
00372 void default_evaluation1(void);
00373 void default_evaluation1m(void);
00374 void default_evaluation2(void);
00375 void default_evaluation3(void);
00376 void default_evaluation4(void);
00377 void default_evaluation4m(void);
00378 
00379 void myheapcheck(char *);
00380 
00381 void RETURN_ARRAYS_INCREMENT(void);
00382 void RETURN_ARRAYS_DECREMENT(void);
00383 
00384 void *farptr_norm(void *);
00385 long int farptr_tolong(void *);
00386 long int _farptr_tolong(void *);
00387 
00388 class i3_array;
00389 
00390 ostream & operator<<(const ostream & ostr, const dmatrix & z);
00391 istream & operator>>(const istream & istr, const dmatrix & z);
00392 uostream & operator<<(const uostream & ostr, const dmatrix & z);
00393 uistream & operator>>(const uistream & istr, const dmatrix & z);
00394 
00395 ostream & operator<<(const ostream & ostr, const d3_array & z);
00396 istream & operator>>(const istream & istr, const d3_array & z);
00397 uostream & operator<<(const uostream & ostr, const d3_array & z);
00398 uistream & operator>>(const uistream & istr, const d3_array & z);
00399 
00400 ostream & operator<<(const ostream & ostr, const dvar3_array & z);
00401 istream & operator>>(const istream & istr, const dvar3_array & z);
00402 uostream & operator<<(const uostream & ostr, const dvar3_array & z);
00403 uistream & operator>>(const uistream & istr, const dvar3_array & z);
00404 
00405 ostream & operator<<(const ostream & ostr, const ivector & z);
00406 istream & operator>>(const istream & istr, const ivector & z);
00407 uostream & operator<<(const uostream & ostr, const ivector & z);
00408 uistream & operator>>(const uistream & istr, const ivector & z);
00409 
00410 ostream & operator<<(const ostream & ostr, const lvector & z);
00411 istream & operator>>(const istream & istr, const lvector & z);
00412 uostream & operator<<(const uostream & ostr, const lvector & z);
00413 uistream & operator>>(const uistream & istr, const lvector & z);
00414 
00415 ostream & operator<<(const ostream & ostr, const dvector & z);
00416 istream & operator>>(const istream & istr, const dvector & z);
00417 uostream & operator<<(const uostream & ostr, const dvector & z);
00418 uistream & operator>>(const uistream & istr, const dvector & z);
00419 
00420 ostream & operator<<(const ostream & ostr, const dvar_vector & z);
00421 istream & operator>>(const istream & istr, const dvar_vector & z);
00422 uostream & operator<<(const uostream & ostr, const dvar_vector & z);
00423 uistream & operator>>(const uistream & istr, const dvar_vector & z);
00424 
00425 ostream & operator<<(const ostream & ostr, const dvar_matrix & z);
00426 istream & operator>>(const istream & istr, const dvar_matrix & z);
00427 uostream & operator<<(const uostream & ostr, const dvar_matrix & z);
00428 uistream & operator>>(const uistream & istr, const dvar_matrix & z);
00429 
00430 
00431 ostream & operator<<(const ostream & ostr, const prevariable & z);
00432 istream & operator>>(const istream & istr, const prevariable & z);
00433 uostream & operator<<(const uostream & ostr, const prevariable & z);
00434 uistream & operator>>(const uistream & istr, const prevariable & z);
00435 
00436 ostream & setscientific(const ostream & s);
00437 
00438 //ostream& setshowpoint(const ostream& s);
00439 
00440 class preshowpoint {};
00441 preshowpoint setshowpoint(void);
00442 ostream & operator <<(const ostream &, preshowpoint);
00443 #if (__BORLANDC__  >= 0x0560)
00444   #define setfixed() std::fixed
00445 #else
00446 ostream & setfixed(const ostream & s);
00447 class prefixed {};
00448 prefixed setfixed(void);
00449 ostream & operator<<(const ostream &, prefixed);
00450 #endif
00451 
00452 #if (__BORLANDC__  >= 0x0560)
00453   #define setscientific() std::scientific
00454 #else
00455 class prescientific {};
00456 prescientific setscientific(void);
00457 ostream & operator<<(const ostream &, prescientific);
00458 #endif
00459 
00460 istream & operator>>(const istream & istr, const imatrix & z);
00461 ostream & operator<<(const ostream & istr, const imatrix & z);
00462 istream & operator>>(const istream & istr, const i3_array & z);
00463 ostream & operator<<(const ostream & istr, const i3_array & z);
00464 class grad_stack;
00465 
00471 class kkludge_object{};
00472 
00477 class humungous_pointer
00478 {
00479  public:
00480    char *ptr;
00481    int adjustment;
00482    void free(void);
00483    ~humungous_pointer()
00484    {;
00485    }
00486    void adjust(int);
00487    humungous_pointer();
00488    humungous_pointer operator +(unsigned long int &offset);
00489    humungous_pointer & operator +=(unsigned long int &offset);
00490    humungous_pointer & operator =(void *p);
00491    int operator ==(void *p);
00492    int operator !=(void *p);
00493    operator  char *();
00494    operator  void *();
00495    operator  double_and_int *();
00496    operator  double *();
00497 };
00498 
00499 #ifndef _VECTOR_SHAPE
00500 #define _VECTOR_SHAPE
00501 #include <dfpool.h>
00502 
00507 class vector_shape_pool:public dfpool
00508 {
00509  public:
00510    vector_shape_pool(void);
00511    vector_shape_pool(int);
00512 };
00513 
00518 class ts_vector_shape_pool:public tsdfpool
00519 {
00520  public:
00521    ts_vector_shape_pool(void);
00522    ts_vector_shape_pool(int);
00523 };
00524 #endif
00525 
00530 class vector_shape
00531 {
00532  public:
00533 #if defined(USE_VECTOR_SHAPE_POOL)
00534    static vector_shape_pool *xpool;
00535    void *operator new(size_t);
00536    void operator delete(void *ptr, size_t)
00537    {
00538       xpool->free(ptr);
00539    }
00540 #endif
00541    unsigned int ncopies;
00542    void shift(int min);
00543    int index_min;
00544    int index_max;
00545  private:
00546    friend class dvector;
00547    //friend class tdvector;
00548    friend class subdvector;
00549    friend class dvar_vector;
00550    friend class ivector;
00551    friend class lvector;
00552    friend class ptr_vector;
00553  public:
00554    int decr_ncopies(void)
00555    {
00556       return --ncopies;
00557    }
00558    int get_ncopies(void)
00559    {
00560       return ncopies;
00561    }
00562    int incr_ncopies(void)
00563    {
00564       return ++ncopies;
00565    }
00566    vector_shape(int lb, int lu)
00567    {
00568       index_min = lb;
00569       index_max = lu;
00570       ncopies = 0;
00571    }
00572    int indexmin() const
00573    {
00574       return index_min;
00575    }
00576    int indexmax() const
00577    {
00578       return index_max;
00579    }
00580 };
00581 
00582 
00587 class ptr_vector
00588 {
00589   void **v;
00590   vector_shape *shape;
00591 
00592 public:
00593   ptr_vector();
00594   ptr_vector(const ptr_vector & t);
00595   ptr_vector(int ncl, int ncu);
00596   ~ptr_vector();
00597 
00598   // returns the minimum allowable index
00599   int indexmin() const
00600   {
00601     return shape->index_min;
00602   }
00603   // returns the maximum allowable index
00604   int indexmax() const
00605   {
00606     return shape->index_max;
00607   }
00608   // returns the maximum allowable index
00609   int size() const
00610   {
00611     return shape->index_max - shape->index_min + 1;
00612   }
00613 
00614   void shift(int min);
00615   void allocate(int, int);
00616   //operator void ** ();
00617   void *&operator[] (int i);
00618   void *&operator() (int i);
00619   //void*& elem(int i);
00620   void *&elem(int i)
00621   {
00622     return (*(v + i));
00623   }
00624   int operator!(void) const
00625   {
00626     return (shape == NULL);
00627   }
00628   int operator() (void) const
00629   {
00630     return (shape != NULL);
00631   }
00632 
00633   ptr_vector& operator=(const ptr_vector& t);
00634   void initialize();
00635 };
00636 
00641 class preivector
00642 {
00643    ivector *p;
00644    int lb;
00645    int ub;
00646    inline preivector(ivector * _p, int _lb, int _ub)
00647    {
00648       p = _p;
00649       lb = _lb, ub = _ub;
00650    }
00651    friend class ivector;
00652 };
00653 
00654 
00659 class ivector
00660 {
00661  protected:
00662    int index_min;
00663    int index_max;
00664    int *v;
00665    vector_shapex *shape;
00666 
00667  public:
00668    int allocated(void) const
00669    {
00670       return (shape != NULL);
00671    }
00672 
00673    ivector & operator --(void)
00674    {
00675       index_min--;
00676       index_max--;
00677       v++;
00678       return *this;
00679    }
00680    void reallocate(double size);
00681 
00682    ivector & operator ++(void)
00683    {
00684       index_min++;
00685       index_max++;
00686       v--;
00687       return *this;
00688    }
00689    friend class vector_index;
00690    int operator!(void) const
00691    {
00692       return (shape == NULL);
00693    }
00694 
00695    int &elem(int i)
00696    {
00697       return (v[i]);
00698    }
00699    const int &elem(int i) const
00700    {
00701       return v[i];
00702    }
00703    // returns the minimum allowable index
00704    int indexmin() const
00705    {
00706       return index_min;
00707    }
00708    // returns the maximum allowable index
00709    int indexmax() const
00710    {
00711       return index_max;
00712    }
00713    // returns the maximum allowable index
00714    int size() const
00715    {
00716       return index_max - index_min + 1;
00717    }
00718    int *get_v(void) const
00719    {
00720       return v;
00721    }
00722    ivector & shift(int min);
00723    //ivector(unsigned int sz); //makes an array [0..sz]
00724    //inline allocated(void) {return ::allocated(this);}
00725    ivector(const lvector &);
00726    ivector(void);
00727 
00728    ivector(const dvector &);
00729 
00730    void fill(const char *s);
00731    void fill_seqadd(int, int);
00732    void fill_multinomial(const int &seed, const dvector & p);
00733    void fill_multinomial(const random_number_generator & rng,
00734      const dvector & p);
00735    ivector(int ncl, int ncu);
00736    ivector(const ad_integer & ncl, const index_type & ncu);
00737    void allocate(const ad_integer & ncl, const index_type & ncu);
00738    // makes an array [ncl..ncu]
00739 
00740    ivector(unsigned int sz, long int *x);
00741 
00742    void safe_allocate(int ncl, int ncu);
00743    void allocate(int ncl, int ncu);
00744    void allocate(const ivector & iv);
00745    void allocate(const dvector & dv);
00746    void allocate(const char *);
00747    void allocate(void);
00748    void deallocate(void);
00749    void safe_deallocate(void);
00750 
00751    //operator int* () {return v;}
00752 
00753    ~ivector();
00754 
00755    void save_ivector_position(void) const;
00756 
00757    void save_ivector_value(void) const;
00758 
00759    void write_on(const ostream & s) const;
00760 
00761    void read_from(const istream & s);
00762 
00763    void write_on(const uostream & s) const;
00764 
00765    void read_from(const uistream & s);
00766 
00767 #ifdef OPT_LIB
00768    inline int &operator[] (int i)
00769    {
00770       return (*(v + i));
00771    }
00772    inline int &operator() (int i)
00773    {
00774       return (*(v + i));
00775    }
00776 #else
00777    int &operator[] (int i);
00778    int &operator() (int i);
00779 #endif
00780 
00781 #ifdef OPT_LIB
00782    inline const int &operator[] (int i) const
00783    {
00784       return (*(v + i));
00785    }
00786    inline const int &operator() (int i) const
00787    {
00788       return (*(v + i));
00789    }
00790 #else
00791    const int &operator[] (int i) const;
00792    const int &operator() (int i) const;
00793 #endif
00794 
00795    ivector sub(int lb, int ub)
00796    {
00797       return preivector(this, lb, ub);
00798    }
00799    ivector operator () (int lb, int ub)
00800    {
00801       return preivector(this, lb, ub);
00802    }
00803    ivector(const preivector & pd);
00804 
00805    ivector operator() (const ivector & u);
00806 
00807    ivector(const ivector & t);
00808    void shallow_copy(const ivector & t);
00809 
00810    ivector & operator=(const ivector & t);
00811    ivector & operator =(int t);
00812 
00813    ivector & operator=(const lvector & t);
00814 
00815    friend class imatrix;
00816    friend class i3_array;
00817    void initialize(void);
00818    ivector & operator+=(const ivector & t);
00819    ivector & operator +=(int t);
00820 };
00821 
00822 ivector operator+(const ivector & v, const ivector & w);
00823 ivector operator-(const ivector & v, const ivector & w);
00824 ivector operator+(const ivector & v, int w);
00825 ivector operator-(const ivector & v, int w);
00826 ivector operator+(int v, const ivector & w);
00827 ivector operator-(int v, const ivector & w);
00828 
00833 class lvector_position
00834 {
00835    int min;
00836    int max;
00837    int *v;
00838  public:
00839    lvector_position(void);
00840    lvector_position(const lvector & v);
00841    lvector_position(const lvector_position & dvp);
00842 };
00843 
00848 class lvector
00849 {
00850    AD_LONG_INT *v;
00851    vector_shape *shape;
00852 
00853  public:
00854    int operator!(void) const
00855    {
00856       return (shape == NULL);
00857    }
00858 
00859    inline AD_LONG_INT & elem(int i)
00860    {
00861       return (v[i]);
00862    }
00863    inline const AD_LONG_INT & elem(int i) const
00864    {
00865       return v[i];
00866    }
00867    // returns the minimum allowable index
00868    int indexmin() const
00869    {
00870       return shape->index_min;
00871    }
00872    // returns the maximum allowable index
00873    int indexmax() const
00874    {
00875       return shape->index_max;
00876    }
00877    // returns the maximum allowable index
00878    int size() const
00879    {
00880       return shape ? shape->index_max - shape->index_min + 1 : 0;
00881    }
00882    void shift(int min);
00883 
00884    void fill(const char *s);
00885    void fill_seqadd(const AD_LONG_INT &, const AD_LONG_INT &);
00886    void fill_multinomial(const int &seed, const dvector & p);
00887    void fill_multinomial(const random_number_generator & rng,
00888      const dvector & p);
00889 
00890    //lvector(unsigned int sz); //makes an array [0..sz]
00891 
00892    lvector(const lvector &);
00893 
00894    lvector(const dvector &);
00895 
00896    lvector(const ivector &);
00897 
00898    lvector(void);
00899 
00900    lvector(int ncl, int ncu);
00901    void allocate(int ncl, int ncu);
00902    void allocate(const lvector &);
00903    void allocate(void);
00904    // makes an array [ncl..ncu]
00905 
00906    lvector(unsigned int sz, AD_LONG_INT * x);
00907 
00908    operator   AD_LONG_INT *();
00909 
00910    ~lvector();
00911 
00912    void write_on(const ostream & s) const;
00913 
00914    void read_from(const istream & s);
00915 
00916    void write_on(const uostream & s) const;
00917 
00918    void read_from(const uistream & s);
00919 
00920 #ifdef OPT_LIB
00921    inline AD_LONG_INT & operator[] (int i)
00922    {
00923       return v[i];
00924    }
00925    inline AD_LONG_INT & operator() (int i)
00926    {
00927       return v[i];
00928    }
00929 #else
00930    AD_LONG_INT & operator[](int i);
00931    AD_LONG_INT & operator()(int i);
00932 #endif
00933 
00934 #ifdef OPT_LIB
00935    inline const AD_LONG_INT & operator[] (int i) const
00936    {
00937       return v[i];
00938    }
00939    inline const AD_LONG_INT & operator() (int i) const
00940    {
00941       return v[i];
00942    }
00943 #else
00944    const AD_LONG_INT & operator[] (int i) const;
00945    const AD_LONG_INT & operator() (int i) const;
00946 #endif
00947 
00948    lvector operator() (const lvector & u);
00949 
00950    lvector & operator=(const lvector & t);
00951 
00952    void initialize(void);
00953    friend class lmatrix;
00954 };
00955 
00956 AD_LONG_INT sum(const lvector &);
00957 
00962 class dependent_variables_information
00963 {
00964    int max_num_dependent_variables;
00965    int depvar_count;
00966    ptr_vector grad_buffer_position;
00967    lvector cmpdif_buffer_position;
00968    lvector grad_file_position;
00969    lvector cmpdif_file_position;
00970    lvector grad_file_count;
00971    lvector cmpdif_file_count;
00972    dependent_variables_information(int ndv);
00973    friend class gradient_structure;
00974 };
00975 dvar_vector_position restore_dvar_vector_position(void);
00976 dvector restore_dvar_vector_value(const dvar_vector_position & tmp);
00977 void arr_free(double_and_int *);
00978 double_and_int *arr_new(unsigned int sz);
00979 
00984 class gradient_structure
00985 {
00986    static char cmpdif_file_name[61];
00987    static DF_FILE *fp;
00988  public:
00989 #if defined(NO_DERIVS)
00990    static int no_derivatives;
00991 #endif
00992  private:
00993    static long int USE_FOR_HESSIAN;
00994    static long int NVAR;
00995    static int NUM_RETURN_ARRAYS;
00996    static dvariable **RETURN_ARRAYS;
00997    static int RETURN_ARRAYS_PTR;
00998    static dvariable **RETURN_PTR_CONTAINER;
00999    static long int TOTAL_BYTES;
01000    static long int PREVIOUS_TOTAL_BYTES;
01001    static unsigned long ARRAY_MEMBLOCK_SIZE;//js
01002    static humungous_pointer ARRAY_MEMBLOCK_BASE;
01003    static humungous_pointer ARRAY_MEMBLOCK_BASEA;
01004    static humungous_pointer ARRAY_MEMBLOCK_SAVE;
01005  public:
01006    static double *get_ARRAY_MEMBLOCK_BASE()
01007    {
01008       return (double *) ARRAY_MEMBLOCK_BASE;
01009    }
01010  private:
01011    static double *variables_save;
01012 #ifdef __BORLANDC__
01013    static long int CMPDIF_BUFFER_SIZE;
01014    static long int GRADSTACK_BUFFER_SIZE;
01015 #else
01016    static long long int CMPDIF_BUFFER_SIZE;
01017    static long long int GRADSTACK_BUFFER_SIZE;
01018 #endif
01019    static unsigned int MAX_NVAR_OFFSET;
01020    static int save_var_file_flag;
01021    static int save_var_flag;
01022 
01023    static unsigned int MAX_DLINKS;
01024    static indvar_offset_list *INDVAR_LIST;
01025    static int NUM_DEPENDENT_VARIABLES;
01026    static dependent_variables_information *DEPVARS_INFO;
01027 
01028    // this needs to be a static member function so other static
01029    // member functions can call it
01030    static void check_set_error(const char *variable_name);
01031 
01032    static int instances;
01033    int x;
01034 
01035  public:
01036    // exception class
01037    class arrmemblerr
01038    {
01039    };
01040 
01041    static int Hybrid_bounded_flag;
01042    static double *hessian_ptr;
01043    static int get_USE_FOR_HESSIAN(void)
01044    {
01045       return USE_FOR_HESSIAN;
01046    }
01047    static void set_USE_FOR_HESSIAN(int i)
01048    {
01049       USE_FOR_HESSIAN = i;
01050    }
01051    friend class dfsdmat;
01052    gradient_structure(long int size = 100000L);// constructor
01053    ~gradient_structure(void);// destructor
01054    static void save_variables(void);
01055    static void restore_variables(void);
01056    static void save_arrays(void);
01057    static void restore_arrays(void);
01058    static long int totalbytes(void);
01059    friend dvector restore_dvar_vector_value(
01060      const dvar_vector_position& tmp);
01061    friend void cleanup_temporary_files();
01062    //friend dvector restore_dvar_vector_value(const dvar_vector_position&);
01063    friend dvector restore_dvar_vector_derivatives(void);
01064    friend dmatrix restore_dvar_matrix_derivatives(void);
01065    friend dmatrix restore_dvar_matrix_value(void);
01066    //friend dmatrix restore_derivatives(void);
01067    friend void gradfree(dlink * v);
01068    friend double_and_int *arr_new(unsigned int sz);//js
01069    friend void arr_free(double_and_int *);
01070    friend void RETURN_ARRAYS_DECREMENT(void);
01071    friend void RETURN_ARRAYS_INCREMENT(void);
01072    friend void make_indvar_list(const dvar_vector & t);
01073    //friend void gradcalc( int , double *);
01074    friend void gradcalc(int nvar, const dvector & g);
01075    friend void slave_gradcalc(void);
01076    friend void funnel_gradcalc(void);
01077    friend void allocate_dvariable_space(void);
01078    friend void wide_funnel_gradcalc(void);
01079    friend dvar_vector_position restore_dvar_vector_position(void);
01080    static grad_stack *GRAD_STACK1;
01081    friend double_and_int *gradnew();
01082    static dlist *GRAD_LIST;
01083    static int RETURN_ARRAYS_SIZE;
01084    //static int RETURN_INDEX;
01085    static dvariable *RETURN_PTR;
01086    static dvariable *MIN_RETURN;
01087    static dvariable *MAX_RETURN;
01088    static arr_list *ARR_LIST1;
01089    static arr_list *ARR_FREE_LIST1;
01090    //static void funnel_jacobcalc(void);
01091    static void jacobcalc(int nvar, const dmatrix & jac);
01092    static void jacobcalc(int nvar, const ofstream & jac);
01093    static void jacobcalc(int nvar, const uostream & jac);
01094 
01095    friend void default_evaluation(void);
01096    //access functions
01097 
01098    friend class DF_FILE;
01099    static DF_FILE *get_fp(void)
01100    {
01101       return fp;
01102    }
01103    static void set_NUM_RETURN_ARRAYS(int i);
01104 #if defined(NO_DERIVS)
01105    static void set_NO_DERIVATIVES(void);
01106    static void set_YES_DERIVATIVES(void);
01107 #endif
01108    static void set_YES_SAVE_VARIABLES_VALUES();
01109    static void set_NO_SAVE_VARIABLES_VALUES();
01110    //static int _GRADFILE_PTR; // should be int gradfile_handle;
01111    //static int _GRADFILE_PTR1; // should be int gradfile_handle;
01112    //static int _GRADFILE_PTR2; // should be int gradfile_handle;
01113    //static int _VARSSAV_PTR; // should be int gradfile_handle;
01114    static void set_NUM_DEPENDENT_VARIABLES(int i);
01115    static void set_RETURN_ARRAYS_SIZE(int i);
01116    static void set_ARRAY_MEMBLOCK_SIZE(unsigned long i);
01117 #ifdef __BORLANDC__
01118    static void set_CMPDIF_BUFFER_SIZE(long int i);
01119    static void set_GRADSTACK_BUFFER_SIZE(long int i);
01120    static void set_GRADSTACK_BUFFER_BYTES(long int i);
01121 #else
01122    static void set_CMPDIF_BUFFER_SIZE(long long int i);
01123    static void set_GRADSTACK_BUFFER_SIZE(long long int i);
01124    static void set_GRADSTACK_BUFFER_BYTES(long long int i);
01125 #endif
01126    static void set_MAX_NVAR_OFFSET(unsigned int i);
01127    static void set_MAX_DLINKS(int i);
01128    static long int NUM_GRADSTACK_BYTES_WRITTEN(void);
01129    friend class dlist;
01130    friend class grad_stack;
01131    static void save_dependent_variable_position(const prevariable & v1);
01132    static long int max_last_offset;
01133    friend class function_minimizer;
01134    friend void funnel_derivatives(void);
01135 };
01136 
01137 void jacobcalc(int nvar, const dmatrix & g);
01138 void jacobcalc(int nvar, const ofstream & ofs);
01139 void jacobcalc(int nvar, const uostream & ofs);
01140 
01141 #if defined(__BORLANDC__ )
01142 #if defined(__GNUC__)
01143 #if (__GNUC__ < 3)
01144 #pragma interface
01145 #endif
01146 #else
01147 #pragma interface
01148 #endif
01149 #endif
01150 
01151 //class dvect_ptr_ptr { dvector **m; };
01152 
01157 class dlink
01158 {
01159    double_and_int di;
01160    dlink *prev;
01161  public:// comments
01162    dlink * previous();
01163    //access function
01164    inline double_and_int *get_address()
01165    {
01166       return &di;
01167    }
01168 
01169    //friend tempvar();
01170    //friend class prevariable;
01171    //friend class tempvar;
01172    friend class dlist;
01173    friend void gradcalc(int nvar, const dvector & g);
01174    friend void slave_gradcalc(void);
01175    friend void gradloop();
01176    friend double_and_int *gradnew();
01177    friend void allocate_dvariable_space(void);
01178 };
01179 
01183 class dlist
01184 {
01185   dlink* last;
01186   unsigned int last_offset;
01187   unsigned int nlinks;
01188   dlink** dlink_addresses;
01189   friend double_and_int *gradnew();
01190   friend void df_check_derivative_values(void);
01191   friend void df_check_derivative_values_indexed(void);
01192   friend void df_check_derivative_values_indexed_break(void);
01193 
01194 public:
01195   // constructor
01196   dlist();
01197   // destructor
01198   ~dlist();
01199   //create a new link
01200   dlink* create();
01201   // add a link
01202   dlink* append(dlink*);
01203   dlink* last_remove();
01204 
01205   // check list integrity
01206   void check_list(void);
01207   size_t total_addresses() const;
01208 
01209   friend void funnel_gradcalc(void);
01210   friend void slave_gradcalc(void);
01211   friend void gradcalc(int nvar, const dvector& g);
01212   friend void gradloop();
01213   friend void gradient_structure::restore_variables();
01214   friend void gradient_structure::save_variables();
01215   friend void gradient_structure::jacobcalc(int nvar,
01216                                             const dmatrix& jac);
01217    friend void allocate_dvariable_space(void);
01218    //friend void gradient_structure::funnel_jacobcalc(void);
01219    friend void gradient_structure::jacobcalc(int nvar,
01220      const ofstream& jac);
01221    friend void gradient_structure::jacobcalc(int nvar,
01222      const uostream& jac);
01223    friend void funnel_derivatives(void);
01224 };
01225 
01226 class indvar_offset_list;
01227 
01234 class grad_stack_entry
01235 {
01236  public:
01238    void (*func) (void);
01239    double *dep_addr;    
01240    double *ind_addr1;   
01241    double mult1;        
01242    double *ind_addr2;   
01243    double mult2;        
01244  public:
01245    friend void gradcalc(int nvar, const dvector & g);
01246    friend void slave_gradcalc(void);
01247    friend void gradloop();
01248    friend void default_evaluation(void);
01249    friend class grad_stack;
01250    friend void gradient_structure::jacobcalc(int nvar,
01251      const dmatrix & jac);
01252    //friend void gradient_structure::funnel_jacobcalc(void);
01253 };
01254 void default_evaluation3ind(void);
01255 
01260 class grad_stack
01261 {
01262    grad_stack_entry *true_ptr_first;
01263    grad_stack_entry *ptr_first;
01264    grad_stack_entry *ptr_last;
01265 #ifdef __BORLANDC__
01266    long int length;
01267    long int true_length;
01268 #else
01269    long long int length;
01270    long long int true_length;
01271 #endif
01272  public:
01273    grad_stack_entry * ptr;
01274  private:
01275    //lvector * table;
01276    // js
01277    int _GRADFILE_PTR;  // should be int gradfile_handle;
01278    int _GRADFILE_PTR1;  // should be int gradfile_handle;
01279    int _GRADFILE_PTR2;  // should be int gradfile_handle;
01280    int _VARSSAV_PTR;  // should be int gradfile_handle;
01281    char gradfile_name[61];
01282    char gradfile_name1[61];
01283    char gradfile_name2[61];
01284    char var_store_file_name[61];
01285    void create_gradfile();
01286 #ifdef __BORLANDC__
01287    long end_pos;
01288    long end_pos1;
01289    long end_pos2;
01290 #else
01291    long long end_pos;
01292    long long end_pos1;
01293    long long end_pos2;
01294 #endif
01295    dmatrix *table;
01296  public:
01297    friend void gradcalc(int nvar, const dvector & g);
01298    friend void slave_gradcalc(void);
01299    friend void funnel_gradcalc(void);
01300    friend void default_evaluation(void);
01301    friend void default_evaluation3ind(void);
01302    friend void reset_gradient_stack(void);
01303    friend void default_evaluation4ind(void);
01304    friend void grad_chk(void);
01305    friend void gradloop();
01306    friend void cleanup_temporary_files();
01307    ostream & operator  <<(grad_stack);
01308    void print();
01309    grad_stack();
01310    ~grad_stack();
01311    void write_grad_stack_buffer(void);
01312 
01313    void set_gradient_stack(void (*func) (void),
01314      double *dep_addr, double *ind_addr1 = NULL, double mult1 = 0,
01315      double *ind_addr2 = NULL, double mult2 = 0);
01316 
01317    void set_gradient_stack(void (*func) (void),
01318      double *dep_addr, double *ind_addr1,
01319      double *ind_addr2);
01320 
01321 
01322    void set_gradient_stack0(void (*func) (void), double *dep_addr);
01323 
01324    void set_gradient_stack1(void (*func) (void),
01325      double *dep_addr, double *ind_addr1);
01326 
01327    void set_gradient_stack2(void (*func) (void),
01328      double *dep_addr, double *ind_addr1, double mult1);
01329 
01330    void set_gradient_stack4(void (*func) (void),
01331      double *dep_addr, double *ind_addr1,
01332      double *ind_addr2);
01333 
01334    void set_gradient_stack(void (*func) (void),
01335      double *dep_addr, double *ind_addr1,
01336      double mult1, double *ind_addr2, double mult2,
01337      double *ind_addr3, double mult3,
01338      double *ind_addr4, double mult4);
01339 
01340    void set_gradient_stack(void (*func) (void),
01341      double *dep_addr, double *ind_addr1,
01342      double mult1, double *ind_addr2, double mult2,
01343      double *ind_addr3, double mult3);
01344 
01345    int read_grad_stack_buffer(my_off_t & lpos);
01346    void set_gradient_stack(void (*ptr) (void));
01347    void set_gbuffer_pointers(void);
01348    //js
01349    void increment_current_gradfile_ptr(void);
01350    int decrement_current_gradfile_ptr(void);
01351    //void open_gradfile();
01352    //void close_gradfile();
01353 #ifdef _MSC_VER
01354    int gradfile_handle();
01355 #else
01356    int &gradfile_handle();
01357 #endif
01358    char *get_gradfile_name();
01359    friend class gradient_structure;
01360    //int get_ngradfiles();
01361 };
01362 
01363 
01364 #ifdef OPT_LIB
01365 
01369 inline void grad_stack::set_gradient_stack(void (*func) (void),
01370   double *dep_addr, double *ind_addr1, double mult1, double *ind_addr2,
01371   double mult2)
01372 {
01373 #ifdef NO_DERIVS
01374    if (!gradient_structure::no_derivatives)
01375    {
01376 #endif
01377 #     if defined(MYDEBUG)
01378       int wrote_buffer = 0;
01379 #     endif
01380       if (ptr > ptr_last)
01381       {
01382         // current buffer is full -- write it to disk and reset pointer
01383         // and counter
01384         this->write_grad_stack_buffer();
01385 #     if defined(MYDEBUG)
01386         wrote_buffer = 1;
01387 #     endif
01388       }
01389 #     if defined(MYDEBUG)
01390       if (wrote_buffer == 1)
01391       {
01392         cout << "WROTE BUFFER" << endl;
01393       }
01394 #     endif
01395       ptr->func = func;
01396       ptr->dep_addr = dep_addr;
01397       ptr->ind_addr1 = ind_addr1;
01398       ptr->mult1 = mult1;
01399       ptr->ind_addr2 = ind_addr2;
01400       ptr->mult2 = mult2;
01401       ptr++;
01402 #ifdef NO_DERIVS
01403    }
01404 #endif
01405 }
01406 #else
01407   //  void grad_stack::set_gradient_stack(void (* func)(void),
01408    //   double * dep_addr,double * ind_addr1, double mult1, double * ind_addr2,
01409     //  double mult2);
01410 #endif
01411 
01416 inline void grad_stack::set_gradient_stack(void (*func) (void),
01417   double *dep_addr, double *ind_addr1, double mult1, double *ind_addr2,
01418   double mult2, double *ind_addr3, double mult3)
01419 {
01420 #ifdef NO_DERIVS
01421    if (!gradient_structure::no_derivatives)
01422    {
01423 #endif
01424       if (ptr > ptr_last)
01425       {
01426          // current buffer is full -- write it to disk and reset pointer
01427          // and counter
01428          this->write_grad_stack_buffer();
01429       }
01430       ptr->func = NULL;
01431       ptr->dep_addr = dep_addr;
01432       ptr->ind_addr1 = ind_addr1;
01433       ptr->mult1 = mult1;
01434       ptr->ind_addr2 = ind_addr2;
01435       ptr->mult2 = mult2;
01436       ptr++;
01437       if (ptr > ptr_last)
01438       {
01439          // current buffer is full -- write it to disk and reset pointer
01440          // and counter
01441          this->write_grad_stack_buffer();
01442       }
01443       ptr->func = func;
01444       ptr->ind_addr1 = ind_addr3;
01445       ptr->mult1 = mult3;
01446       ptr++;
01447 #ifdef NO_DERIVS
01448    }
01449 #endif
01450 }
01451 
01456 inline void grad_stack::set_gradient_stack(void (*func) (void),
01457   double *dep_addr, double *ind_addr1, double mult1, double *ind_addr2,
01458   double mult2, double *ind_addr3, double mult3, double *ind_addr4,
01459   double mult4)
01460 {
01461 #ifdef NO_DERIVS
01462    if (!gradient_structure::no_derivatives)
01463    {
01464 #endif
01465       if (ptr > ptr_last)
01466       {
01467          // current buffer is full -- write it to disk and reset pointer
01468          // and counter
01469          this->write_grad_stack_buffer();
01470       }
01471       ptr->func = NULL;
01472       ptr->dep_addr = dep_addr;
01473       ptr->ind_addr1 = ind_addr1;
01474       ptr->mult1 = mult1;
01475       ptr->ind_addr2 = ind_addr2;
01476       ptr->mult2 = mult2;
01477       ptr++;
01478       if (ptr > ptr_last)
01479       {
01480          // current buffer is full -- write it to disk and reset pointer
01481          // and counter
01482          this->write_grad_stack_buffer();
01483       }
01484       ptr->func = func;
01485       ptr->ind_addr1 = ind_addr3;
01486       ptr->mult1 = mult3;
01487       ptr->ind_addr2 = ind_addr4;
01488       ptr->mult2 = mult4;
01489       ptr++;
01490 #ifdef NO_DERIVS
01491    }
01492 #endif
01493 }
01494 
01499 inline void grad_stack::set_gradient_stack(void (*func) (void),
01500   double *dep_addr, double *ind_addr1, double *ind_addr2)
01501 {
01502 #ifdef NO_DERIVS
01503    if (!gradient_structure::no_derivatives)
01504    {
01505 #endif
01506       if (ptr > ptr_last)
01507       {
01508          // current buffer is full -- write it to disk and reset pointer
01509          // and counter
01510          this->write_grad_stack_buffer();
01511       }
01512       ptr->func = func;
01513       ptr->dep_addr = dep_addr;
01514       ptr->ind_addr1 = ind_addr1;
01515       ptr->ind_addr2 = ind_addr2;
01516       ptr++;
01517 #ifdef NO_DERIVS
01518    }
01519 #endif
01520 }
01521 
01526 inline void grad_stack::set_gradient_stack2(void (*func) (void),
01527   double *dep_addr, double *ind_addr1, double mult1)
01528 {
01529 #ifdef NO_DERIVS
01530    if (!gradient_structure::no_derivatives)
01531    {
01532 #endif
01533       if (ptr > ptr_last)
01534       {
01535          // current buffer is full -- write it to disk and reset pointer
01536          // and counter
01537          this->write_grad_stack_buffer();
01538       }
01539       ptr->func = func;
01540       ptr->dep_addr = dep_addr;
01541       ptr->ind_addr1 = ind_addr1;
01542       ptr->mult1 = mult1;
01543       ptr++;
01544 #ifdef NO_DERIVS
01545    }
01546 #endif
01547 }
01548 
01553 inline void grad_stack::set_gradient_stack4(void (*func) (void),
01554   double *dep_addr, double *ind_addr1, double *ind_addr2)
01555 {
01556 #ifdef NO_DERIVS
01557    if (!gradient_structure::no_derivatives)
01558    {
01559 #endif
01560       if (ptr > ptr_last)
01561       {
01562          // current buffer is full -- write it to disk and reset pointer
01563          // and counter
01564          this->write_grad_stack_buffer();
01565       }
01566       ptr->func = func;
01567       ptr->dep_addr = dep_addr;
01568       ptr->ind_addr1 = ind_addr1;
01569       ptr->ind_addr2 = ind_addr2;
01570       ptr++;
01571 #ifdef NO_DERIVS
01572    }
01573 #endif
01574 }
01575 
01580 inline void grad_stack::set_gradient_stack(void (*func) (void))
01581 {
01582 #ifdef NO_DERIVS
01583    if (!gradient_structure::no_derivatives)
01584    {
01585 #endif
01586       if (ptr > ptr_last)
01587       {
01588          // current buffer is full -- write it to disk and reset pointer
01589          // and counter
01590          this->write_grad_stack_buffer();
01591       }
01592 
01593       ptr->dep_addr = NULL;
01594       ptr->func = func;
01595       // want to put a long int into the memory space of a double
01596       ptr->ind_addr2 = NULL;
01597       ptr->mult2 = 0;
01598       ptr++;
01599 #ifdef NO_DERIVS
01600    }
01601 #endif
01602 }
01603 
01608 class indvar_offset_list
01609 {
01610    // The number of independent variables
01611    int nvar;
01612    double **address;
01613 
01614  public:
01615    friend class gradient_structure;
01616    inline double *get_address(const int &i)
01617    {
01618       return address[i];
01619    }
01620    void put_address(unsigned int &i, double *iaddress)
01621    {
01622       address[i] = iaddress;
01623       //  cerr << "In put_address i = " << i << "\n";
01624    }
01625 };
01626 
01627 void gradfree(dlink *);
01628 
01629 class prevariable_position;
01630 
01638 class prevariable
01639 {
01640  protected:
01641 #ifndef __SUN__
01642    prevariable(void)
01643    {
01644    }
01645 #endif
01646 #ifndef __NDPX__
01647    prevariable(double_and_int * u)
01648    {
01649       v = u;
01650    }
01651 #endif
01652 
01653  public:
01654    double_and_int * v; 
01655    friend class dvar_vector_iterator;
01656    friend class dvar_vector;
01657    friend class dvar_matrix;
01658    friend class dvar3_array;
01659    //shclass sc;
01660    friend class indvar_offset_list;
01661    friend class gradient_structure;
01662    friend double_and_int *gradnew();
01663    friend void make_indvar_list(int, dvariable *);
01664 
01665    friend class banded_symmetric_dvar_matrix;
01666    friend class banded_lower_triangular_dvar_matrix;
01667    friend double &value(const prevariable & v1);
01668 
01669    friend double *address(const prevariable & v1);
01670 
01671    //void gradfree(dlink * v)
01672 
01673    friend prevariable & operator*(const prevariable& v1, const prevariable& v2);
01674 
01675    friend prevariable & operator*(double v1, const prevariable & v2);
01676 
01677    friend prevariable & operator*(const prevariable & v1, double v2);
01678 
01679    friend prevariable & operator/(const prevariable& t1, const prevariable& t2);
01680 
01681    friend prevariable & operator/(double t1, const prevariable & t2);
01682 
01683    friend prevariable & operator/(const prevariable & t1, double t2);
01684 
01685    friend prevariable & sin(const prevariable & t1);
01686 
01687    friend prevariable & fabs(const prevariable & t1);
01688    friend prevariable & sigmoid(const prevariable & t1);
01689 
01690    //"smoothed absolute value function
01691    friend prevariable & sfabs(const prevariable & t1);
01692 
01693    friend prevariable & sqrt(const prevariable & t1);
01694    friend prevariable & sqr(const prevariable & t1);
01695 
01696    friend prevariable & exp(const prevariable & t1);
01697 
01698    friend prevariable & atan(const prevariable & t1);
01699 
01700    friend prevariable & tan(const prevariable & t1);
01701    friend prevariable & tanh(const prevariable & t1);
01702 
01703    friend prevariable & atan2(const prevariable & t1,
01704      const prevariable & t2);
01705    friend prevariable & atan2(const prevariable & t1, double t2);
01706    friend prevariable & atan2(double t1, const prevariable & t2);
01707 
01708    friend prevariable & acos(const prevariable & t1);
01709 
01710    friend prevariable & asin(const prevariable & t1);
01711 
01712    friend prevariable & cos(const prevariable & t1);
01713    friend prevariable & sinh(const prevariable & t1);
01714    friend prevariable & cosh(const prevariable & t1);
01715 
01716    friend prevariable & log(const prevariable & t1);
01717    friend prevariable & log10(const prevariable & t1);
01718 
01719    friend prevariable & ldexp(const prevariable &, const int &);
01720 
01721 
01722  public:
01723    void save_prevariable_position(void) const;
01724    prevariable_position restore_prevariable_position(void);
01725    void save_prevariable_value(void) const;
01726    double restore_prevariable_value(void);
01727    double restore_prevariable_derivative(void);
01728 
01729 
01730    inline double *xadr()
01731    {
01732       return (&(v->x));
01733    }
01734    inline double &xval()
01735    {
01736       return ((v->x));
01737    }
01738 
01739    inline double_and_int *&get_v()
01740    {
01741       return v;
01742    }
01743    inline double_and_int *get_v() const
01744    {
01745       return v;
01746    }
01747 
01748    prevariable & operator=(const prevariable &);
01749    prevariable & operator=(double);
01750 #if (__BORLANDC__  >= 0x0540)
01751    prevariable & operator=(const prevariable &) const;
01752    prevariable & operator=(double) const;
01753 #endif
01754 
01755    int operator==(const prevariable & v1) const;
01756    int operator<=(const prevariable & v1) const;
01757    int operator>=(const prevariable & v1) const;
01758    int operator>(const prevariable & v1) const;
01759    int operator<(const prevariable & v1) const;
01760    int operator!=(const prevariable & v1) const;
01761 
01762    int operator==(double v1) const;
01763    int operator<=(double v1) const;
01764    int operator>=(double v1) const;
01765    int operator>(double v1) const;
01766    int operator<(double v1) const;
01767    int operator!=(double v1) const;
01768 #if defined(USE_DDOUBLE)
01769    int operator==(int v1) const;
01770    int operator<=(int v1) const;
01771    int operator>=(int v1) const;
01772    int operator>(int v1) const;
01773    int operator<(int v1) const;
01774    int operator!=(int v1) const;
01775 #endif
01776 
01777  public:
01778 #ifdef __SUN__
01779    prevariable(void)
01780    {
01781    }
01782 #endif
01783 #ifdef __NDPX__
01784    prevariable(double_and_int * u)
01785    {
01786       v = u;
01787    }
01788 #endif
01789 
01790    void initialize(void);
01791 
01792    friend char *fform(const char *, const prevariable &);
01793 
01794    void operator+=(const prevariable & t1);
01795    void operator +=(double t1);
01796 
01797    void operator-=(const prevariable & t1);
01798    void operator -=(double t1);
01799 
01800    void operator/=(const prevariable & v1);
01801    void operator /=(double v1);
01802 
01803    void operator*=(const prevariable & v1);
01804    void operator *=(double v1);
01805 
01806    friend prevariable& operator+(const prevariable& v1, const prevariable& v2);
01807 
01808    friend prevariable & operator+(double v1, const prevariable & v2);
01809 
01810    friend prevariable & operator+(const prevariable & v1, double v2);
01811 
01812    friend prevariable & operator-(const prevariable & v1);
01813 
01814    friend prevariable & operator-(const prevariable& v1, const prevariable& v2);
01815 
01816    friend prevariable & operator-(double v1, const prevariable & v2);
01817 
01818    friend prevariable & operator-(const prevariable & v1, double v2);
01819 
01820    friend prevariable & pow(const prevariable & t1, const prevariable & t2);
01821 
01822    friend prevariable & pow(const prevariable & t1, double t2);
01823 
01824    friend prevariable & pow(double t1, const prevariable & t2);
01825 };
01826 
01827 inline double &value(const prevariable & v1)
01828 {
01829    return v1.v->x;
01830 }
01831 
01832 inline double *address(const prevariable & v1)
01833 {
01834    return (&(v1.v->x));
01835 }
01836 
01837 
01838 prevariable & operator<<(const prevariable & v1, const prevariable & v2);
01839 dvar_vector & operator<<(const dvar_vector & v1, const dvar_vector & v2);
01840 dvar_matrix & operator<<(const dvar_matrix & v1, const dvar_matrix & v2);
01841 
01842 class df1_one_variable;
01843 class df1_two_variable;
01844 class df1_three_variable;
01845 
01850 class dvariable:public prevariable
01851 {
01852  public:
01853    dvariable();
01854    ~dvariable();
01855    dvariable(double t);
01856    dvariable(const int &t);
01857    dvariable(kkludge_object);
01858    dvariable(const prevariable &);
01859    dvariable & operator=(const prevariable &);
01860    dvariable & operator =(const df1_one_variable & v);
01861    dvariable & operator =(const df1_two_variable & v);
01862    dvariable & operator =(const df1_three_variable & v);
01863    dvariable & operator=(double);
01864 #if defined(USE_DDOUBLE)
01865 #  undef double
01866    dvariable & operator=(double);
01867 #  define double dd_real
01868 #endif
01869    dvariable(const dvariable &);
01870 //#  if (__BORLANDC__  > 0x0520)
01871 //     dvariable& operator+=(const prevariable&);
01872 //#  endif
01873 };
01874 
01875 #if defined(max)
01876 #undef max
01877 #endif
01878 #if defined(min)
01879 #undef min
01880 #endif
01881 
01886 class funnel_dvariable:public dvariable
01887 {
01888  public:
01889    dvariable & operator=(const prevariable &);
01890 };
01891 
01892 dvar_vector operator*(const dvar_vector & t1, double x);
01893 dvar_vector operator/(double x, const dvar_vector & t1);
01894 dvar_vector operator/(const dvar_vector & t1, double x);
01895 dvar_vector operator+(double x, const dvar_vector & t1);
01896 dvar_vector operator+(const dvar_vector & t1, double x);
01897 dvar_vector operator-(double x, const dvar_vector & t1);
01898 dvar_vector operator-(const dvar_vector & t1, double x);
01899 dvar_vector operator-(const dvar_vector & t1);
01900 dvar_vector operator*(const dvar_vector & t1, const prevariable & x);
01901 dvar_vector operator/(const prevariable & x, const dvar_vector & t1);
01902 dvar_vector operator/(const dvar_vector & t1, const prevariable & x);
01903 dvar_vector operator+(const prevariable & x, const dvar_vector & t1);
01904 dvar_vector operator+(const dvar_vector & t1, const prevariable & x);
01905 dvar_vector operator-(const prevariable & x, const dvar_vector & t1);
01906 dvar_vector operator-(const dvar_vector & t1, const prevariable & x);
01907 dvar_vector operator-(const dvector & t1, const prevariable & x);
01908 dvar_vector operator*(const dvector & t1, const prevariable & x);
01909 dvar_vector operator*(const prevariable & x, const dvector & t1);
01910 
01911 dvector operator*(const dvector & t1, double x);
01912 dvector operator/(double x, const dvector & t1);
01913 dvector operator/(const dvector & t1, double x);
01914 dvector operator+(double x, const dvector & t1);
01915 dvector operator+(const dvector & t1, double x);
01916 dvector operator-(double x, const dvector & t1);
01917 dvector operator-(const dvector & t1, double x);
01918 dvector operator-(const dvector & t1);
01919 
01920 double min(const dmatrix &);
01921 double max(const dmatrix &);
01922 int max(const imatrix &);
01923 double max(const dvector &);
01924 dvariable max(const dvar_vector &);
01925 dvariable min(const dvar_vector &);
01926 
01927 dmatrix symmetrize(const dmatrix & m1);
01928 dvector eigenvalues(const dmatrix & m1);
01929 dmatrix eigenvectors(const dmatrix & m1);
01930 dmatrix eigenvectors(const dmatrix & m1, const dvector & eigenvalues);
01931 
01932 dvar_matrix symmetrize(const dvar_matrix & m1);
01933 dvar_vector eigenvalues(const dvar_matrix & m1);
01934 dvar_matrix eigenvectors(const dvar_matrix & m1);
01935 
01936 dmatrix outer_prod(const dvector & t1, const dvector & t2);
01937 dvar_matrix outer_prod(const dvar_vector & t1, const dvar_vector & t2);
01938 dvar_matrix outer_prod(const dvector & t1, const dvar_vector & t2);
01939 dvar_matrix outer_prod(const dvar_vector & t1, const dvector & t2);
01940 dmatrix operator*(double x, const dmatrix & m);
01941 dmatrix operator*(const dmatrix & m, double d);
01942 dmatrix operator/(const dmatrix & m, double d);
01943 dmatrix operator/(double d, const dmatrix & m);
01944 dmatrix operator+(double x, const dmatrix & m);
01945 dvar_matrix operator +(const dvariable & x, const dmatrix & m);
01946 dvar_matrix operator -(const dvariable & x, const dmatrix & m);
01947 dmatrix operator+(const dmatrix & m, double d);
01948 dmatrix operator-(double x, const dmatrix & m);
01949 dmatrix operator-(const dmatrix & m, double d);
01950 dvar_matrix operator/(const dvar_matrix & m, const prevariable & x);
01951 dvar_matrix operator/(const dmatrix & m, const prevariable & x);
01952 dvar_matrix operator/(const dvar_matrix & m, double x);
01953 dvar_matrix operator/(double x, const dvar_matrix & m);
01954 dvar_matrix operator/(const prevariable & x, const dvar_matrix & m);
01955 
01956 dvar_matrix operator*(const prevariable & x, const dmatrix & m);
01957 dvar_matrix operator*(const dvar_matrix & m, const prevariable & x);
01958 dvar_matrix operator*(const prevariable & x, const dvar_matrix & m);
01959 dvar_matrix operator*(double x, const dvar_matrix & m);
01960 
01961 dvector operator&(const dvector & t1, const dvector & t2);
01962 dvar_vector operator&(const dvar_vector & t1, const dvar_vector & t2);
01963 
01964 ivector column(const imatrix & m, int i);
01965 dvector extract_column(const dmatrix & m, int i);
01966 dvector column(const dmatrix & m, int i);
01967 dvector extract_row(const dmatrix & m, int j);
01968 dvector row(const dmatrix & m, int j);
01969 dvar_vector extract_column(const dvar_matrix & m, int i);
01970 dvar_vector column(const dvar_matrix & m, int i);
01971 dvector column_value(const dvar_matrix & m, int i);
01972 dvar_vector extract_row(const dvar_matrix & m, int j);
01973 dvar_vector row(const dvar_matrix & m, int j);
01974 
01975 // dvector mathematical functions
01976 dvector sin(const dvector & t1);
01977 dvector sqrt(const dvector & t1);
01978 dvector sqr(const dvector & t1);
01979 dvector exp(const dvector & t1);
01980 dvector mfexp(const dvector & t1);
01981 dvector mfexp(const dvector & t1, double d);
01982 dvector atan(const dvector & t1);
01983 dvector tan(const dvector & t1);
01984 dvector tanh(const dvector & t1);
01985 dvector atan2(const dvector & t1, const dvector & t2);
01986 dvector atan2(const dvector & t1, double t2);
01987 dvector atan2(double t1, const dvector & t2);
01988 dvector acos(const dvector & t1);
01989 dvector asin(const dvector & t1);
01990 dvector cos(const dvector & t1);
01991 dvector sinh(const dvector & t1);
01992 dvector cosh(const dvector & t1);
01993 dvector log(const dvector & t1);
01994 dvector log10(const dvector & t1);
01995 dvector pow(const dvector & t1, double);
01996 dvector pow(const dvector & t1, int);
01997 dvector pow(double, const dvector & t1);
01998 ivector pow(const ivector & v1, int x);
01999 ivector pow(int x, const ivector & v1);
02000 // end of dvector mathematical functions
02001 
02002 // dvar_vector mathematical functions
02003 dvar_vector sin(const dvar_vector & t1);
02004 dvar_vector sqrt(const dvar_vector & t1);
02005 dvar_vector sqr(const dvar_vector & t1);
02006 dvar_vector exp(const dvar_vector & t1);
02007 dvar_vector mfexp(const dvar_vector & t1);
02008 dvar_vector mfexp(const dvar_vector & t1, double d);
02009 dvar_vector atan(const dvar_vector & t1);
02010 dvar_vector tan(const dvar_vector & t1);
02011 dvar_vector tanh(const dvar_vector & t1);
02012 dvar_vector atan2(const dvar_vector & t1, const dvar_vector & t2);
02013 dvar_vector atan2(const dvar_vector & t1, double t2);
02014 dvar_vector atan2(double t1, const dvar_vector & t2);
02015 dvar_vector acos(const dvar_vector & t1);
02016 dvar_vector asin(const dvar_vector & t1);
02017 dvar_vector cos(const dvar_vector & t1);
02018 dvar_vector sinh(const dvar_vector & t1);
02019 dvar_vector cosh(const dvar_vector & t1);
02020 dvar_vector log(const dvar_vector & t1);
02021 dvar_vector log10(const dvar_vector & t1);
02022 dvar_vector pow(const dvar_vector &, const dvar_vector & t1);
02023 dvar_vector pow(const dvar_vector &, const dvector & t1);
02024 dvar_vector pow(const dvector &, const dvar_vector & t1);
02025 dvector pow(const dvector &, const dvector & t1);
02026 dvar_vector pow(const dvar_vector & t1, double);
02027 dvar_vector pow(const dvar_vector & t1, int);
02028 dvar_vector pow(const dvar_vector & t1, const prevariable &);
02029 dvar_vector pow(const dvector & t1, const prevariable &);
02030 dvar_vector pow(const prevariable &, const dvar_vector & t1);
02031 dvar_vector pow(const dvector & x, const dvar_vector & a);
02032 // end of dvar_vector mathematical functions
02033 
02034 // dmatrix mathematical functions
02035 dmatrix exp(const dmatrix & m);
02036 dmatrix mfexp(const dmatrix & m);
02037 dmatrix mfexp(const dmatrix & m, double d);
02038 dmatrix sqrt(const dmatrix & m);
02039 dmatrix sqr(const dmatrix & m);
02040 dmatrix pow(const dmatrix & m, double e);
02041 dmatrix pow(const dmatrix & m, int e);
02042 dmatrix log(const dmatrix & m);
02043 dmatrix sin(const dmatrix & m);
02044 dmatrix cos(const dmatrix & m);
02045 dmatrix tan(const dmatrix & m);
02046 dmatrix elem_div(const dmatrix & m, const dmatrix & m2);
02047 dmatrix elem_prod(const dmatrix & m, const dmatrix & m2);
02048 // end of dmatrix mathematical functions
02049 
02050 //  dvar_matrix mathematical functions
02051 dvar_matrix exp(const dvar_matrix & m);
02052 dvar_matrix mfexp(const dvar_matrix & m);
02053 dvar_matrix mfexp(const dvar_matrix & m, double d);
02054 dvar_matrix sqrt(const dvar_matrix & m);
02055 dvar_matrix sqr(const dvar_matrix & m);
02056 dvar_matrix log(const dvar_matrix & m);
02057 dvar_matrix sin(const dvar_matrix & m);
02058 dvar_matrix cos(const dvar_matrix & m);
02059 dvar_matrix tan(const dvar_matrix & m);
02060 dvar_matrix pow(const dvar_matrix & m, double e);
02061 dvar_matrix pow(const dvar_matrix & m, const prevariable & e);
02062 dvar_matrix pow(const dmatrix & m, const prevariable & e);
02063 dvar_matrix pow(const dvar_matrix & m, int e);
02064 dvar_matrix elem_prod(const dvar_matrix & m, const dvar_matrix & m2);
02065 dvar_matrix elem_prod(const dvar_matrix & m, const dmatrix & m2);
02066 dvar_matrix elem_prod(const dmatrix & m, const dvar_matrix & m2);
02067 dvar_matrix elem_div(const dvar_matrix & m, const dvar_matrix & m2);
02068 dvar_matrix elem_div(const dvar_matrix & m, const dmatrix & m2);
02069 dvar_matrix elem_div(const dmatrix & m, const dvar_matrix & m2);
02070 //  end of dvar_matrix mathematical functions
02071 
02072 int min(const ivector & t1);
02073 int max(const ivector & t1);
02074 int Max(const ivector & t1);
02075 
02076 double mfexp(double);
02077 double mfexp(double, double bound);
02078 dvariable mfexp(const prevariable &);
02079 dvariable mfexp(const prevariable &, double bound);
02080 
02081 #if defined(THREAD_SAFE)
02082 
02085 class ts_vector_shapex
02086 {
02087  public:
02088    void *trueptr;
02089    ts_vector_shapex(int lb, int ub, void *p):index_min(lb),
02090       index_max(ub), ncopies(0), trueptr(p)
02091    {
02092    }
02093    void *get_truepointer(void)
02094    {
02095       return trueptr;
02096    }
02097    //friend class dvector;
02098    friend class ivector;
02099    //friend class tdvector;
02100    friend class dvar_vector;
02101 
02102 #  if defined(USE_VECTOR_SHAPE_POOL)
02103    static ts_vector_shape_pool **xpool;
02104    void *operator  new(size_t);
02105    void operator  delete(void *ptr, size_t n);
02106 #  endif
02107 
02108    unsigned int ncopies;
02109    void shift(int min);
02110    int index_min;
02111    int index_max;
02112  private:
02113    friend class subdvector;
02114    friend class lvector;
02115    friend class ptr_vector;
02116  public:
02117    int decr_ncopies(void)
02118    {
02119       return --ncopies;
02120    }
02121    int get_ncopies(void)
02122    {
02123       return ncopies;
02124    }
02125    int incr_ncopies(void)
02126    {
02127       return ++ncopies;
02128    }
02129    int indexmin()
02130    {
02131       return index_min;
02132    }
02133    int indexmax()
02134    {
02135       return index_max;
02136    }
02137 };
02138 #endif
02139 
02143 class vector_shapex
02144 {
02145  public:
02146    vector_shapex(int lb, int ub, void *p):index_min(lb),
02147       index_max(ub), ncopies(0), trueptr(p)
02148    {
02149    }
02151    void *get_truepointer(void)
02152    {
02153       return trueptr;
02154    }
02155    friend class dvector;
02156    friend class ivector;
02157    friend class dvar_vector;
02158 #if defined(USE_VECTOR_SHAPE_POOL)
02159    static vector_shape_pool *xpool;
02160    void *operator new(size_t);
02161    void operator delete(void *ptr, size_t)
02162    {
02163       xpool->free(ptr);
02164    }
02165 #endif
02166    void shift(int min);
02167    int index_min;  
02168    int index_max;  
02169 
02176    unsigned int ncopies;
02179    void *trueptr;
02180 
02181  private:
02182    friend class subdvector;
02183    friend class lvector;
02184    friend class ptr_vector;
02185  public:
02186    int decr_ncopies(void)  
02187    {
02188       return --ncopies;
02189    }
02190    int get_ncopies(void) 
02191    {
02192       return ncopies;
02193    }
02194    int incr_ncopies(void)  
02195    {
02196       return ++ncopies;
02197    }
02198    int indexmin() 
02199    {
02200       return index_min;
02201    }
02202    int indexmax() 
02203    {
02204       return index_max;
02205    }
02206 };
02207 
02212 class predvector
02213 {
02214    dvector *p;
02215    int lb;
02216    int ub;
02217    inline predvector(dvector * _p, int _lb, int _ub)
02218    {
02219       p = _p;
02220       lb = _lb, ub = _ub;
02221    }
02222    friend class dvector;
02223 };
02224 
02229 class predvar_vector
02230 {
02231    dvar_vector *p;
02232    int lb;
02233    int ub;
02234    inline predvar_vector(dvar_vector * _p, int _lb, int _ub)
02235    {
02236       p = _p;
02237       lb = _lb, ub = _ub;
02238    }
02239    friend class dvar_vector;
02240 };
02241 
02242 
02249 class dvector
02250 {
02251  protected:
02252    double *v; 
02253    int index_min; 
02254    int index_max; 
02255 
02256 #if defined(THREAD_SAFE)
02257    ts_vector_shapex *shape;
02258 #else
02259    vector_shapex *shape;
02260 #endif
02261 
02262  public:
02263    dvector operator -();
02265    int operator!(void) const
02266    {
02267       return (shape == NULL);
02268    }
02270    int allocated(void) const
02271    {
02272       return (shape != NULL);
02273    }
02275    dvector & operator --(void)
02276    {
02277       index_min--;
02278       index_max--;
02279       v++;
02280       return *this;
02281    }
02283    dvector & operator ++(void)
02284    {
02285       index_min++;
02286       index_max++;
02287       v--;
02288       return *this;
02289    }
02290    void reallocate(double size);
02297    dvector sub(int lb, int ub)
02298    {
02299       return predvector(this, lb, ub);
02300    }
02307    dvector operator () (int lb, int ub)
02308    {
02309       return predvector(this, lb, ub);
02310    }
02311 
02312    dvector(const predvector & pd);
02313 
02314    void fill_randpoisson(double lambda, const random_number_generator & rng);
02315    void fill_randnegbinomial(double lambda, double tau,
02316      const random_number_generator & rng);
02317    void shallow_copy(const dvector &);
02318 
02319    friend class banded_symmetric_dvar_matrix;
02320    friend class banded_lower_triangular_dvar_matrix;
02321    friend class banded_symmetric_dmatrix;
02322    friend class banded_lower_triangular_dmatrix;
02323 
02324    void allocate(int ncl, int ncu);
02325    void allocate(const dvector & dv);
02326    void allocatec(const dvector & dv);
02327    void allocate(const dvar_vector &);
02328    void allocate(const char *);
02329    void allocate(void);
02330    void deallocate(void);
02331    void safe_deallocate(void);
02332    void safe_allocate(int, int);
02333 
02334    const double *address() const
02335    {
02336       return v;
02337    }
02338    double *&get_v(void)
02339    {
02340       return (v);
02341    }
02342    double &elem(int i)
02343    {
02344       return (v[i]);
02345    }
02346    double *initpointer(void)
02347    {
02348       return (v + indexmin());
02349    }
02350    const double *initpointer(void) const
02351    {
02352       return (v + indexmin());
02353    }
02354    double *get_v(void) const
02355    {
02356       return v;
02357    }
02358    const double &elem(int i) const
02359    {
02360       return v[i];
02361    }
02362 
02363    void fill(const char *s);
02364    void fill_randu_ni(long int &n);
02365    void fill_randn_ni(long int &n);
02366    void fill_randbi_ni(long int &n, double);
02367 
02368    void fill_randu(long int &n);
02369    void fill_randn(long int &n);
02370    void fill_randbi(long int &n, double);
02371 
02372 
02373    void fill_randbi(double p, const random_number_generator & rng);
02374    void fill_randu(const random_number_generator & rng);
02375    void fill_randn(const random_number_generator & rng);
02376    void fill_randcau(const random_number_generator & rng);
02377 
02378    void fill_seqadd(double, double);
02379    void fill_multinomial(const int &seed, const dvector & p);
02380    void fill_multinomial(const random_number_generator& rng, const dvector& p);
02381    void initialize(void);
02382 
02383    // returns the minimum allowable index
02384    int &testmin()
02385    {
02386       return shape->index_min;
02387    }
02388 
02389    // returns the minimum allowable index
02390    int &testmax()
02391    {
02392       return shape->index_max;
02393    }
02394 
02396    int indexmin() const
02397    {
02398       return index_min;
02399    }
02400 
02402    int indexmax() const
02403    {
02404       return index_max;
02405    }
02406 
02408    int size() const
02409    {
02410       return (index_max - index_min + 1);
02411    }
02412 
02413    dvector & shift(int min);
02414 
02415    dvector(const dvar_vector_position & dvp, const kkludge_object &);
02416 
02417    dvector(const ad_integer &, const index_type &);
02418    void allocate(const ad_integer &, const index_type &);
02419    dvector(void);
02420    dvector(const dvector &);
02421    //dvector(const dvector&,int lb,int ub);
02422 
02423    dvector(const ivector &);
02424 
02425    dvector(const lvector &);
02426 
02427    dvector(const char *);
02428 
02429    dvector(int ncl, int ncu);
02430    // makes an array [ncl..ncu]
02431 
02432    dvector(unsigned int sz, double *x);
02433 
02434    dvector(char *filename, const int &column);
02435 
02436    //operator double* () { return v;}
02437 
02438    ~dvector();
02439 
02440    void save_dvector_position(void) const;
02441    void save_dvector_derivatives(const dvar_vector_position & pos) const;
02442    void save_dvector_derivatives_na(const dvar_vector_position & pos)
02443       const;
02444    void save_dvector_derivatives(void) const;
02445    void save_dvector_value(void) const;
02446 
02447    //dvector operator()(int,int);
02448    dvector operator() (const lvector &);
02449    dvector operator() (const ivector & u);
02450 
02451    dvector & operator+=(const dvector & v1);
02452    dvector & operator-=(const dvector & v1);
02453    dvector & operator +=(double v1);
02454    dvector & operator /=(double v1);
02455    dvector & operator *=(double v1);
02456    dvector & operator -=(double v1);
02457 
02458    void read_from(const uistream &);
02459 
02460    friend class sdmatrix;
02461    friend double norm(const dvector &);
02462    friend double norm2(const dvector &);
02463    friend double sumsq(const dvector &);
02464    friend class dvar_vector;
02465    friend class dmatrix;
02466    friend class d3_array;
02467    friend char *fform(const char *, const dvector &);
02468 
02469    void write_on(const ostream &) const;
02470    void write_on(const uostream &) const;
02471    void read_from(const istream &);
02472    friend double operator*(const dvector &, const dvector &);
02473 
02474    friend dvariable operator*(const dvector &, const dvar_vector &);
02475 
02476    friend dvariable operator*(const dvar_vector &, const dvector &);
02477 
02478    friend dvar_vector operator*(const prevariable &, const dvector &);
02479 
02480    friend dvector operator*(double, const dvector &);
02481 
02482    friend dvector operator+(const dvector &, const dvector &);
02483 
02484    friend dvector elem_prod(const dvector &, const dvector &);
02485 
02486    friend dvector first_difference(const dvector &);
02487    friend dvector second_difference(const dvector &);
02488 
02489    friend dvector elem_div(const dvector &, const dvector &);
02490 
02491    friend dvar_vector elem_prod(const dvector &, const dvar_vector &);
02492 
02493    friend dvar_vector elem_div(const dvector &, const dvar_vector &);
02494 
02495    friend dvar_vector elem_prod(const dvar_vector &, const dvector &);
02496 
02497    friend dvar_vector elem_div(const dvar_vector &, const dvector &);
02498 
02499    friend dvar_vector operator+(const dvar_vector &, const dvector &);
02500 
02501    friend dvar_vector operator+(const dvector &, const dvar_vector &);
02502 
02503    friend dvector operator-(const dvector &, const dvector &);
02504 
02505    friend dvar_vector operator-(const dvar_vector &, const dvector &);
02506 
02507    friend dvar_vector operator-(const dvector &, const dvar_vector &);
02508 
02509    friend dvector operator*(const dvector & x, const dmatrix & m);
02510 
02511    friend dvector operator*(const dmatrix & x, const dvector & m);
02512 
02513    friend dvar_vector operator*(const dvector & x, const dvar_matrix & m);
02514 
02515    friend dvar_vector operator*(const dvar_matrix & x, const dvector & m);
02516 
02517 
02518 #ifdef OPT_LIB
02519    inline double &operator[] (register int i)
02520    {
02521       return (*(v + i));
02522    }
02523    inline double &operator() (register int i)
02524    {
02525       return (*(v + i));
02526    }
02527 #else
02528    double &operator[] (int i);
02529    double &operator() (int i);
02530 #endif
02531 
02532 #ifdef OPT_LIB
02533    inline const double &operator[] (register int i) const
02534    {
02535       return (*(v + i));
02536    }
02537    inline const double &operator() (register int i) const
02538    {
02539       return (*(v + i));
02540    }
02541 #else
02542    const double &operator[] (int i) const;
02543    const double &operator() (int i) const;
02544 #endif
02545 
02546    dvector & operator=(const dvector & t);
02547 
02548    dvector & operator =(double x);
02549 
02550    friend dvector exp(const dvector &);
02551 
02552    friend dvector log(const dvector &);
02553 
02554    friend dvector fabs(const dvector &);
02555 
02556    friend double max(const dvector &);
02557 
02558    friend double min(const dvector &);
02559 };
02560 
02565 class independent_variables:public dvector
02566 {
02567  public:
02568    independent_variables(const independent_variables & v):dvector(v)
02569    {
02570    }
02571 
02572    independent_variables(int ncl, int ncu):dvector(ncl, ncu)
02573    {
02574    }
02575    // makes an array [ncl..ncu]
02576 
02577    independent_variables(unsigned int sz, double *x):dvector(sz, x)
02578    {
02579    }
02580 
02581    independent_variables & operator=(const dvector & t);
02582 };
02583 
02584 
02585 
02586 dvariable dfatan1(dvariable, double, double, const prevariable & fpen);
02587 
02588 double boundp(double x, double fmin, double fmax, const double &fpen);
02589 double boundp(double x, double fmin, double fmax);
02590 
02591 dvariable boundp(const prevariable & x, double fmin, double fmax,
02592   const prevariable & fpen);
02593 dvariable boundp(const prevariable & x, double fmin, double fmax,
02594   const prevariable & fpen, double s);
02595 
02596 double nd2fboundp(double x, double minb, double maxb, const double &pen);
02597 double boundpin(double x, double fmin, double fmax);
02598 double boundpin(const prevariable & x, double fmin, double fmax);
02599 double boundpin(const prevariable & x, double fmin, double fmax, double s);
02600 
02601 double dmin(double, double);
02602 
02603 double dmax(double i, double j);
02604 
02605 #include <stdlib.h>
02606 #ifdef __TURBOC__
02607   #include <alloc.h>
02608 #endif
02609 
02610 double sigmoid(double t1);
02611 
02616 class mat_shape
02617 {
02618    unsigned int ncopies;
02619    unsigned int nrows;
02620    unsigned int ncols;
02621    int row_min;
02622    int row_max;
02623    int col_min;
02624    int col_max;
02625    mat_shape(int rl, int ru, int cl = 0, int cu = -1);
02626    mat_shape()
02627    {
02628    };
02629    void colshift(int min);
02630    void rowshift(int min);
02631 
02632    //friend class const_dmatrix;
02633    friend class dmatrix;
02634    friend class sdmatrix;
02635    friend class dvar_matrix;
02636    friend class imatrix;
02637    friend class lmatrix;
02638    friend class i3_array;
02639 };
02640 
02645 class mat_shapex
02646 {
02647  public:
02648    void *trueptr;
02649    unsigned int ncopies;
02650    mat_shapex(const void *m)
02651    {
02652       trueptr = (void *) m;
02653       ncopies = 0;
02654    }
02655    mat_shapex()
02656    {
02657       trueptr = NULL;
02658       ncopies = 0;
02659    };
02660 
02661    void *get_pointer(void)
02662    {
02663       return trueptr;
02664    }
02665    friend class dmatrix;
02666    friend class sdmatrix;
02667    friend class dvar_matrix;
02668    friend class imatrix;
02669    friend class i3_array;
02670    friend class lmatrix;
02671 };
02672 
02673 class arr_link;
02674 
02679 class arr_list
02680 {
02681    arr_link *last;
02682    arr_link *free_last;
02683    unsigned long int last_offset;
02684    unsigned long int max_last_offset;
02685  public:
02686    long int number_arr_links;
02687    friend class arr_link;
02688 
02689  public:
02690    arr_list(void)
02691    {
02692       last = 0;
02693       free_last = 0;
02694       last_offset = 0;
02695       max_last_offset = 0;
02696       number_arr_links = 0;
02697    }
02698    unsigned long int get_last_offset()
02699    {
02700       return last_offset;
02701    }
02702    unsigned long int get_number_arr_links()
02703    {
02704       return (number_arr_links);
02705    }
02706    unsigned long int get_max_last_offset()
02707    {
02708       return (max_last_offset);
02709    }
02710    void reset_max_last_offset()
02711    {
02712       max_last_offset = 0;
02713    }
02714    friend double_and_int *arr_new(unsigned int);
02715    friend void arr_free(double_and_int *);
02716    friend void arr_remove(arr_link **);
02717    friend void arr_free_list_remove(arr_link **);
02718    friend void arr_free_add(arr_link *);
02719    friend void arr_free_remove(arr_link *);
02720 };
02721 
02726 class arr_link
02727 {
02728 #if defined(USE_VECTOR_SHAPE_POOL)
02729    static vector_shape_pool *xpool;
02730    void *operator new(size_t);
02731    void operator delete(void* ptr, size_t)
02732    {
02733       xpool->free(ptr);
02734    }
02735 #endif
02736    arr_link *prev;
02737    arr_link *next;
02738    arr_link *free_prev;
02739    arr_link *free_next;
02740    unsigned int status;
02741    // unsigned int     free_list_status;
02742    unsigned int size;
02743    unsigned long int offset;
02744  public:
02745    arr_link();
02746 
02747    friend double_and_int *arr_new(unsigned int);
02748    friend void arr_free(double_and_int *);
02749    friend void arr_remove(arr_link **);
02750    friend void arr_free_remove(arr_link *);
02751    friend void arr_free_add(arr_link *);
02752 };
02753 
02754 #if defined(__NUMBERVECTOR__)
02755 class param_init_number_vector;
02756 class param_init_bounded_number_vector;
02757 class param_init_bounded_number_matrix;
02758 class param_init_vector_vector;
02759 class param_init_bounded_vector_vector;
02760 #endif
02761 
02766 class dvar_vector
02767 {
02768  public:
02769    double_and_int * va;
02770    int index_min;
02771    int index_max;
02772    arr_link *link_ptr;
02773    vector_shapex *shape;
02774  public:
02775    dvar_vector operator -();
02776 
02777    dvar_vector & operator --(void)
02778    {
02779       index_min--;
02780       index_max--;
02781       va++;
02782       return *this;
02783    }
02784    dvar_vector & operator ++(void)
02785    {
02786       index_min++;
02787       index_max++;
02788       va--;
02789       return *this;
02790    }
02791    dvar_vector sub(int lb, int ub)
02792    {
02793       return predvar_vector(this, lb, ub);
02794    }
02795    dvar_vector operator () (int lb, int ub)
02796    {
02797       return predvar_vector(this, lb, ub);
02798    }
02799    void shallow_copy(const dvar_vector &);
02800    int operator!(void) const
02801    {
02802       return (shape == NULL);
02803    }
02804    friend class dvar_matrix;
02805    friend class dvar3_array;
02806    friend class banded_symmetric_dvar_matrix;
02807    friend class banded_lower_triangular_dvar_matrix;
02808    friend class banded_symmetric_dmatrix;
02809    friend class banded_lower_triangular_dmatrix;
02810 
02811    void fill_randpoisson(double lambda, const random_number_generator & rng);
02812    void fill_randnegbinomial(double lambda, double tau,
02813      const random_number_generator & rng);
02814    prevariable elem(int i)
02815    {
02816       return (va + i);
02817    }
02818 
02819    double &elem_value(int i)
02820    {
02821       return (va[i].x);
02822    }
02823 
02824    double_and_int *get_va()
02825    {
02826       return va;
02827    }
02828 
02829    prevariable elem(int i) const
02830    {
02831       return (va + i);
02832    }
02833 
02834    double &elem_value(int i) const
02835    {
02836       return (va[i].x);
02837    }
02838 
02839    double_and_int *get_va() const
02840    {
02841       return va;
02842    }
02843 
02844    friend dvar_matrix operator*(const dvar_matrix & m1, const dmatrix & m2);
02845 
02846    void deallocate();
02847    dvar_vector(const dvar_vector &);
02848    dvar_vector(const predvar_vector &);
02849    dvar_vector();
02850    dvar_vector(int ncl, int ncu); // makes an array [ncl..ncu]
02851    dvar_vector(int ncl, int ncu, kkludge_object);
02852 
02853    //dvar_vector(const ad_integer&,const ad_integer&);
02854    dvar_vector(unsigned int sz, double *x);
02855    dvar_vector(const independent_variables &);
02856    friend char *fform(const char *, const dvar_vector &);
02857 #   if defined(__NUMBERVECTOR__)
02858    dvar_vector(const param_init_number_vector &);
02859    dvar_vector(const param_init_bounded_number_vector &);
02860 #   endif
02861    dvar_vector(const dvector &);
02862    dvar_vector(const char *);
02863    ~dvar_vector();
02864    void allocate(int, int);
02865    void allocate(void);
02866    void allocate(const dvector &);
02867    void allocatec(const dvar_vector &);
02868    void allocate(const dvar_vector &);
02869 
02870    void allocate(const ad_integer &, const ad_integer &);
02871    void initialize(const dvector & ww);
02872    void initialize(void);
02873    void save_dvar_vector_position(void) const;
02874    void save_dvar_vector_value(void) const;
02875    void write_on(const ostream &) const;
02876    void write_on(const uostream &) const;
02877    void read_from(const istream &);
02878    void read_from(const uistream &);
02879    // returns the minimum allowable index
02880    int indexmin() const
02881    {
02882       return index_min;
02883    }
02884    // returns the maximum allowable index
02885    int indexmax() const
02886    {
02887       return index_max;
02888    }
02889    // returns the number of elements
02890    int size() const
02891    {
02892       return (index_max - index_min + 1);
02893    }
02894    dvar_vector & shift(int min);
02895 
02896 #ifdef OPT_LIB
02897 #if defined(__NDPX__) || defined(__SUN__)
02898    inline prevariable operator() (register int i)
02899    {
02900       return (va + i);
02901    }
02902    inline prevariable operator[] (register int i)
02903    {
02904       return (va + i);
02905    }
02906    inline const prevariable operator() (int i) const
02907    {
02908       return (va + i);
02909    }
02910    inline const prevariable operator[] (int i) const
02911    {
02912       return (va + i);
02913    }
02914 #else
02915    inline prevariable operator() (int i)
02916    {
02917       return (va + i);
02918    }
02919    inline prevariable operator[] (int i)
02920    {
02921       return (va + i);
02922    }
02923    inline const prevariable operator() (int i) const
02924    {
02925       return (va + i);
02926    }
02927    inline const prevariable operator[] (int i) const
02928    {
02929       return (va + i);
02930    }
02931 #endif
02932 #else
02933    prevariable operator[] (int i);
02934    prevariable operator() (int i);
02935    const prevariable operator[] (int i) const;
02936    const prevariable operator() (int i) const;
02937 #endif
02938 
02939    double *initpointer(void)
02940    {
02941       return ((double *) (va + indexmin()));
02942    }
02943    const double *initpointer(void) const
02944    {
02945       return ((double *) (va + indexmin()));
02946    }
02947    dvar_vector operator() (const lvector &);
02948    //dvar_vector operator()(int,int);
02949    dvar_vector operator () (const ivector & u);
02950    dvar_vector & operator+=(const prevariable & d);
02951    dvar_vector & operator+=(double d);
02952    dvar_vector & operator/=(const prevariable & d);
02953    //dvar_vector& operator*=(const dvariable& d);
02954    dvar_vector & operator*=(const prevariable & d);
02955    dvar_vector & operator*=(double d);
02956    dvar_vector & operator/=(double d);
02957    dvar_vector & operator-=(const prevariable & d);
02958    dvar_vector & operator-=(double d);
02959    dvar_vector & operator+=(const dvector & v1);
02960    dvar_vector & operator-=(const dvector & v1);
02961    dvar_vector & operator+=(const dvar_vector & v1);
02962    dvar_vector & operator-=(const dvar_vector & v1);
02963    dvar_vector & operator=(const dvar_vector & t);
02964    dvar_vector & operator=(const dvector & t);
02965    dvar_vector & operator =(double t);
02966    dvar_vector & operator=(const prevariable & t);
02967    void fill(const char *);
02968    void fill_randu(long int &n);
02969    void fill_randn(long int &n);
02970    void fill_randbi(long int &n, double);
02971 
02972    void fill_randu_ni(long int &n);
02973    void fill_randn_ni(long int &n);
02974    void fill_randbi_ni(long int &n, double);
02975 
02976    void fill_seqadd(double, double);
02977    void fill_multinomial(const int &seed, const dvector & p);
02978    void fill_multinomial(const random_number_generator& rng, const dvector& p);
02979    friend dvar_vector operator+(const dvar_vector &, const dvar_vector &);
02980    friend dvar_vector operator+(const dvar_vector &, const dvector &);
02981    friend dvar_vector operator+(const dvector &, const dvar_vector &);
02982    friend dvar_vector operator-(const dvar_vector &, const dvar_vector &);
02983 
02984    friend dvar_vector operator-(const dvector &, const dvar_vector &);
02985 
02986    friend dvar_vector operator-(const dvar_vector &, const dvector &);
02987 
02988    friend dvar_vector sigmoid(const dvar_vector & t1);
02989 
02990    friend dvariable operator*(const dvar_vector &, const dvar_vector &);
02991 
02992    friend dvar_vector elem_div(const dvar_vector &, const dvar_vector &);
02993 
02994    friend dvariable operator*(const dvector &, const dvar_vector &);
02995 
02996    friend dvariable operator*(const dvar_vector &, const dvector &);
02997 
02998    friend dvar_vector operator*(const prevariable &, const dvar_vector &);
02999 
03000    friend dvar_vector operator*(const prevariable &, const dvector &);
03001 
03002    friend dvar_vector operator*(double, const dvar_vector &);
03003 
03004    friend dvar_vector operator*(const dvar_vector &, const dmatrix &);
03005 
03006    friend dvar_vector operator*(const dmatrix &, const dvar_vector &);
03007 
03008    friend dvar_vector operator*(const dvar_vector &, const dvar_matrix &);
03009 
03010    friend dvar_vector operator*(const dvar_matrix &, const dvar_vector &);
03011 
03012    friend dvar_matrix operator*(const dvar_matrix &, const dvar_matrix &);
03013 
03014    friend dvar_matrix operator*(const dmatrix &, const dvar_matrix &);
03015 
03016    friend dvar_vector elem_prod(const dvar_vector &, const dvar_vector &);
03017 
03018    friend dvar_vector first_difference(const dvar_vector &);
03019    friend dvar_vector second_difference(const dvar_vector &);
03020 
03021    // js, see above
03022    //friend dvar_vector elem_div(const dvar_vector&, const dvar_vector&);
03023 
03024    friend dvar_vector elem_prod(const dvector &, const dvar_vector &);
03025 
03026    friend dvar_vector elem_div(const dvector &, const dvar_vector &);
03027 
03028    friend dvar_vector elem_prod(const dvar_vector &, const dvector &);
03029 
03030    friend dvar_vector elem_div(const dvar_vector &, const dvector &);
03031 
03032    friend dvariable norm(const dvar_vector &);
03033    friend dvariable norm2(const dvar_vector &);
03034    friend dvariable sumsq(const dvar_vector &);
03035 
03036    friend void copy_status(const ostream & s, const dvar_vector & v);
03037 
03038    friend dvar_vector exp(const dvar_vector &);
03039 
03040    friend dvar_vector log(const dvar_vector &);
03041 
03042    friend dvar_vector sin(const dvar_vector &);
03043 
03044    friend dvar_vector fabs(const dvar_vector &);
03045 
03046    friend dvector value(const dvar_vector & v1);
03047 
03048    friend dvar_vector sfabs(const dvar_vector &);
03049 
03050    friend void make_indvar_list(const dvar_vector &);
03051    friend class array_size;
03052 };
03053 
03054  /*
03055     class funnel_dvar_vector : public dvar_vector
03056     {
03057     public:
03058     funnel_dvar_vector(int l,int u);
03059     dvar_vector& operator=(const dvar_vector&);
03060     };
03061   */
03062 
03063 //class fvar_ptr { dvar_vector *p; };
03064 
03070 class dvar_matrix
03071 {
03072    int index_min;
03073    int index_max;
03074    dvar_vector *m;
03075    mat_shapex *shape;
03076 
03077  public:
03078    dvar_matrix & operator --(void)
03079    {
03080       index_min--;
03081       index_max--;
03082       m++;
03083       return *this;
03084    }
03085    dvar_matrix & operator ++(void)
03086    {
03087       index_min++;
03088       index_max++;
03089       m--;
03090       return *this;
03091    }
03092 
03093    int operator!(void) const
03094    {
03095       return (shape == NULL);
03096    }
03097    inline dvar_vector & elem(int i)
03098    {
03099       return (m[i]);
03100    }
03101    inline prevariable elem(int i, int j)
03102    {
03103       return (elem(i).elem(j));
03104    }
03105    inline dvar_vector & elem(int i) const
03106    {
03107       return (m[i]);
03108    }
03109    inline prevariable elem(int i, int j) const
03110    {
03111       return (elem(i).elem(j));
03112    }
03113 
03114    friend class banded_symmetric_dvar_matrix;
03115    friend class banded_lower_triangular_dvar_matrix;
03116    friend class banded_symmetric_dmatrix;
03117    friend class banded_lower_triangular_dmatrix;
03118    friend class dvar3_array;
03119    void shallow_copy(const dvar_matrix &);
03120    dvar_matrix();
03121    void allocate(int nrl, int nrh, int ncl, int nch);
03122    void allocate(int nrl, int nrh);
03123    void allocate(ad_integer nrl, ad_integer nrh);
03124    void allocate(const dmatrix & m1);
03125    void allocate(const dvar_matrix & m1);
03126    void allocate(int nrl, int nrh, const ivector & ncl, const ivector & nch);
03127    void allocate(int nrl, int nrh, int ncl, const ivector & nch);
03128    void allocate(int nrl, int nrh, const ivector & ncl, int nch);
03129    void allocate(void);
03130    void deallocate();
03131    dvar_matrix(const banded_symmetric_dvar_matrix & v);
03132    dvar_matrix(const banded_lower_triangular_dvar_matrix & v);
03133 # if defined(__NUMBERVECTOR__)
03134    dvar_matrix(const param_init_vector_vector &);
03135    dvar_matrix(const param_init_bounded_vector_vector &);
03136    dvar_matrix(const param_init_bounded_number_matrix &);
03137 # endif
03138    dvar_matrix sub(int, int);
03139 
03140    double fill_seqadd(double, double);
03141 
03142    int colmin(void) const
03143    {
03144       return ((*this) (indexmin()).indexmin());
03145    }
03146    int colmax(void) const
03147    {
03148       return ((*this) (indexmin()).indexmax());
03149    }
03150    int rowmin(void) const
03151    {
03152       return (index_min);
03153    }
03154    int rowmax(void) const
03155    {
03156       return (index_max);
03157    }
03158    int indexmin(void) const
03159    {
03160       return (index_min);
03161    }
03162    int indexmax(void) const
03163    {
03164       return (index_max);
03165    }
03166    // returns the number of rows
03167    int rowsize() const
03168    {
03169       return (rowmax() - rowmin() + 1);
03170    }
03171    // returns the number of columns
03172    int colsize() const
03173    {
03174       return (colmax() - colmin() + 1);
03175    }
03176    void colshift(int min);
03177    void rowshift(int min);
03178 
03179    friend char *fform(const char *, const dvar_matrix &);
03180 
03181    friend class dvar_vector;
03182 
03183    dvar_matrix(const ad_integer & nrl, const ad_integer & nrh,
03184      const index_type & ncl, const index_type & nch);
03185 
03186    void allocate(const ad_integer & nrl, const ad_integer & nrh,
03187      const index_type & ncl, const index_type & nch);
03188 
03189    dvar_matrix(int, int, int, int);
03190    dvar_matrix(int, int);
03191    dvar_matrix(int, int, kkludge_object kk);
03192    // makes a matrix [nrl..nrh][ncl..nch]
03193 
03194    dvar_matrix(int, int, const ivector &, const ivector &);
03195    // makes a ragged dvar_matrix [nrl..nrh][ncl..nch]
03196 
03197    dvar_matrix(int, int, int, const ivector &);
03198    // makes a ragged dvar_matrix [nrl..nrh][ncl..nch]
03199 
03200    dvar_matrix(const dvar_matrix &);
03201    // copy initializer
03202    void initialize(void);
03203 
03204    dvar_matrix(const dmatrix &);
03205 
03206    //dvar_matrix(char *);
03207 
03208    ~dvar_matrix();
03209 
03210    void save_dvar_matrix_position(void) const;
03211    void save_dvar_matrix_value(void) const;
03212 
03213    void fill(const char *);
03214    //void colfill(const int&n,...);
03215    //void rowfill(const int&n,...);
03216 
03217    void colfill_randu(const int &j, long int &n);
03218    void rowfill_randu(const int &i, long int &n);
03219    void colfill_randn(const int &j, long int &n);
03220    void rowfill_randn(const int &i, long int &n);
03221    void fill_randn(long int &n);
03222    void fill_randu(long int &n);
03223 
03224    void colfill_seqadd_ni(const int &, double, double);
03225    void colfill_randu_ni(const int &j, long int &n);
03226    void rowfill_randu_ni(const int &i, long int &n);
03227    void colfill_randn_ni(const int &j, long int &n);
03228    void rowfill_randn_ni(const int &i, long int &n);
03229    void fill_randn_ni(long int &n);
03230    void fill_randu_ni(long int &n);
03231 
03232    void colfill_seqadd(const int &, double, double);
03233    void rowfill_seqadd(const int &, double, double);
03234    void colfill(int j, const dvar_vector & v);
03235    void rowfill(int j, const dvar_vector & v);
03236 
03237    void write_on(const ostream &) const;
03238    void write_on(const uostream &) const;
03239    void read_from(const istream &);
03240    void read_from(const uistream &);
03241 
03242 #ifdef OPT_LIB
03243    inline dvar_vector & operator[] (register int i)
03244    {
03245       return (m[i]);
03246    }
03247    inline dvar_vector & operator() (register int i)
03248    {
03249       return (m[i]);
03250    }
03251 #else
03252    dvar_vector & operator ()(int i);
03253    dvar_vector & operator[](int);
03254 #endif
03255 
03256 #ifdef OPT_LIB
03257    inline const dvar_vector & operator[] (register int i) const
03258    {
03259       return m[i];
03260    }
03261    inline const dvar_vector & operator() (register int i) const
03262    {
03263       return m[i];
03264    }
03265 #else
03266    const dvar_vector & operator() (int i) const;
03267    const dvar_vector & operator[] (int) const;
03268 #endif
03269 
03270 #ifdef OPT_LIB
03271 #ifdef __NDPX__
03272    prevariable operator () (register int i, register int j)
03273    {
03274       return (prevariable((m[i]).va + j));
03275    }
03276 #else
03277    inline prevariable operator () (register int i, register int j)
03278    {
03279       return ((m[i]).va + j);
03280    }
03281 #endif
03282 #else
03283    prevariable operator () (int i, int j);
03284 #endif
03285 
03286    inline double &elem_value(register int i, register int j)
03287    {
03288       return *(double *) ((m[i]).va + j);
03289    }
03290 
03291    inline const double &elem_value(register int i, register int j) const
03292    {
03293       return *(double *) ((m[i]).va + j);
03294    }
03295 #ifdef OPT_LIB
03296 #ifdef __NDPX__
03297    prevariable operator() (register int i, register int j) const
03298    {
03299       return (prevariable((m[i]).va + j));
03300    }
03301 #else
03302    inline prevariable operator() (register int i, register int j) const
03303    {
03304       return ((m[i]).va + j);
03305    }
03306 #endif
03307 #else
03308    const prevariable operator() (int i, int j) const;
03309 #endif
03310 
03311    dvar_matrix & operator+=(const dvar_matrix & x);
03312    dvar_matrix & operator-=(const dvar_matrix & x);
03313    dvar_matrix & operator+=(const dmatrix & x);
03314    dvar_matrix & operator-=(const dmatrix & x);
03315 
03316 
03317    dvar_matrix & operator=(const dvar_matrix &);
03318 
03319    dvar_matrix & operator=(const dmatrix &);
03320    dvar_matrix & operator =(double t);
03321    dvar_matrix & operator=(const prevariable & t);
03322 
03323    dvar_matrix & operator*=(const prevariable & t);
03324    dvar_matrix & operator *=(double t);
03325    dvar_matrix & operator/=(const prevariable & t);
03326    dvar_matrix & operator /=(double t);
03327 
03328    friend dvar_vector operator*(const dvar_vector &, const dvar_matrix &);
03329 
03330    friend dvar_vector operator*(const dvar_matrix &, const dvar_vector &);
03331 
03332    friend dvar_vector operator*(const dvector &, const dvar_matrix &);
03333 
03334    friend dvar_vector operator*(const dvar_matrix &, const dvector &);
03335 
03336    friend dvar_matrix operator*(const dvar_matrix &, const dvar_matrix &);
03337 
03338    friend dvar_matrix operator*(const dvar_matrix &, const dmatrix &);
03339 
03340    friend dvar_matrix operator*(const dmatrix &, const dvar_matrix &);
03341 
03342    friend dvar_matrix operator+(const dvar_matrix &, const dvar_matrix &);
03343    friend dvar_matrix operator+(const dvar_matrix &, const dmatrix &);
03344    friend dvar_matrix operator+(const dmatrix &, const dvar_matrix &);
03345 
03346    friend dvar_matrix operator+(double, const dvar_matrix &);
03347    friend dvar_matrix operator+(const dvar_matrix &, double);
03348    friend dvar_matrix operator-(double, const dvar_matrix &);
03349    friend dvar_matrix operator-(const dvar_matrix &, double);
03350 
03351    friend dvar_matrix operator+(const dvariable &, const dvar_matrix &);
03352    friend dvar_matrix operator+(const dvar_matrix &, const dvariable &);
03353    friend dvar_matrix operator-(const dvariable &, const dvar_matrix &);
03354    friend dvar_matrix operator-(const dvar_matrix &, const dvariable &);
03355 
03356    friend dvar_matrix operator-(const dvar_matrix &, const dvar_matrix &);
03357    friend dvar_matrix operator-(const dvar_matrix &, const dmatrix &);
03358    friend dvar_matrix operator-(const dmatrix &, const dvar_matrix &);
03359 
03360    friend dvar_matrix inv(const dvar_matrix &);
03361 
03362    friend dvariable det(const dvar_matrix &);
03363    friend dvariable ln_det(const dvar_matrix &, const int &sgn);
03364 
03365    //friend dvar_matrix testsub(dvar_matrix);
03366 
03367    friend dvar_matrix trans(const dvar_matrix &);
03368 
03369    friend dvariable norm(const dvar_matrix &);
03370    friend dvariable norm2(const dvar_matrix &);
03371    friend dvariable sumsq(const dvar_matrix &);
03372 
03373    friend void copy_status(const ostream & s, const dvar_matrix & m1);
03374 };
03375 
03376 dvariable ln_det(const dvar_matrix &);
03377 dvar_matrix operator *(const dvar_matrix & t1, double x);
03378 dmatrix value(const dvar_matrix & m);
03379 d3_array value(const dvar3_array & a);
03380 dvar_vector sort(const dvar_vector &, int NSTACK = 60);
03381 dvector sort(const dvector &, int NSTACK = 60);
03382 ivector sort(const ivector &, int NSTACK = 60);
03383 dvector sort(const dvector &, const ivector & index, int NSTACK = 60);
03384 ivector sort(const ivector &, const ivector & index, int NSTACK = 60);
03385 dmatrix sort(const dmatrix &, int column, int NSTACK = 60);
03386 imatrix sort(const imatrix &, int column, int NSTACK = 60);
03387 
03388 
03389 #include "factors.h"
03390 int count_factor(const dvector & v, const double &eps);
03391 ivector as_factor(const dvector & v, const double eps = 1.0e-6);
03392 int count_factor(const ivector & v);
03393 
03394  //void gradcalc( int , double *);
03395 void gradcalc(int nvar, const dvector & g);
03396 double gradcalc(int nvar, const dvector& g, dvariable& f);
03397 void slave_gradcalc(void);
03398 
03403 class dmatrix
03404 {
03405  protected:
03406    int index_min;
03407    int index_max;
03408    dvector *m;
03409    mat_shapex *shape;
03410    friend char *fform(const char *, const dmatrix &);
03411    friend class dvar_matrix;
03412  public:
03413    dmatrix & operator --(void)
03414    {
03415       index_min--;
03416       index_max--;
03417       m++;
03418       return *this;
03419    }
03420    dmatrix & operator ++(void)
03421    {
03422       index_min++;
03423       index_max++;
03424       m--;
03425       return *this;
03426    }
03427    void shallow_copy(const dmatrix &);
03428    int operator!(void) const
03429    {
03430       return (shape == NULL);
03431    }
03432 
03433    dmatrix sub(int, int);
03434    dmatrix(void);
03435    dmatrix(int, int, kkludge_object);
03436    dmatrix(int, int);
03437    void allocate(void);
03438    void allocate(const dmatrix & dm);
03439    void allocate(const dvar_matrix &);
03440    void allocate(int nrl, int nrh, int ncl, int nch);
03441    void allocate(int nrl, int nrh);
03442    void allocate(ad_integer nrl, ad_integer nrh);
03443    void allocate(int nrl, int nrh, int ncl, const ivector & nch);
03444    //void allocate(int nrl,int nrh,
03445    // const index_type& ncl,const index_type& nch);
03446    void allocate(int nrl, int nrh, const ivector & ncl, int nch);
03447    void deallocate();
03448    void allocate(const ad_integer & nrl, const ad_integer & nrh,
03449      const index_type & ncl, const index_type & nch);
03450    void allocate(int nrl, int nrh, const ivector & ncl,
03451      const ivector & nch);
03452    friend class banded_symmetric_dmatrix;
03453    friend class banded_lower_triangular_dmatrix;
03454 
03455    dmatrix(int, int, int, int);
03456    // makes a matrix [nrl..nrh][ncl..nch]
03457 
03458    dmatrix(const ad_integer & nrl, const ad_integer & nrh,
03459      const index_type & ncl, const index_type & nch);
03460 
03461    dmatrix(int, int, const ivector & coll, const ivector & colh);
03462    // makes a ragged dmatrix[nrl..nrh][ncl..nch]
03463 
03464    dmatrix(int, int, int coll, const ivector & colh);
03465    // makes a ragged dmatrix[nrl..nrh][ncl..nch]
03466 
03467    dmatrix(const dvar_matrix_position &);
03468 
03469    dmatrix(const dmatrix_position &);
03470 
03471    dmatrix(const dmatrix &);
03472    dmatrix(const banded_symmetric_dmatrix &);
03473    dmatrix(const banded_lower_triangular_dmatrix &);
03474    dmatrix(char *);
03475    void fill(const char *);
03476    double fill_seqadd(double, double);
03477    void initialize(void);
03478    // copy initializer
03479 
03480    ~dmatrix();
03481    void save_dmatrix_derivatives(const dvar_matrix_position & pos) const;
03482    void save_dmatrix_derivatives_na(const dvar_matrix_position & pos)
03483       const;
03484    void save_dmatrix_value(void) const;
03485    void save_dmatrix_position(void) const;
03486    //void save_dmatrix_derivatives(void);
03487 
03488    int indexmin(void) const
03489    {
03490       return index_min;
03491    }
03492    int indexmax(void) const
03493    {
03494       return index_max;
03495    }
03496    int rowmin(void) const
03497    {
03498       return index_min;
03499    }
03500    int rowmax(void) const
03501    {
03502       return index_max;
03503    }
03504    int colmin(void) const
03505    {
03506       return ((*this) (indexmin()).indexmin());
03507    }
03508    int colmax(void) const
03509    {
03510       return ((*this) (indexmin()).indexmax());
03511    }
03512    // returns the number of rows
03513    int rowsize() const
03514    {
03515       return (rowmax() - rowmin() + 1);
03516    }
03517    // returns the number of columns
03518    int colsize() const
03519    {
03520       return (colmax() - colmin() + 1);
03521    }
03522    void rowshift(int min);
03523    void colshift(int min);
03524 
03525    void write_on(const ostream &) const;
03526    void write_on(const uostream &) const;
03527    void read_from(const istream &);
03528    void read_from(const uistream &);
03529 
03530    //void colfill(const int&n,...);
03531    //void rowfill(const int&n,...);
03532 
03533    void colfill_randu(const int &j, long int &n);
03534    void rowfill_randu(const int &i, long int &n);
03535    void colfill_randn(const int &j, long int &n);
03536    void fill_randn(long int &n);
03537    void fill_randu(long int &n);
03538    void rowfill_randn(const int &i, long int &n);
03539 
03540 
03541 
03542    void colfill_randu(const int &j, const random_number_generator & rng);
03543    void rowfill_randu(const int &i, const random_number_generator & rng);
03544    void fill_randn(const random_number_generator & rng);
03545    void fill_randcau(const random_number_generator & rng);
03546    void fill_randu(const random_number_generator & rng);
03547    void colfill_randn(const int &j, const random_number_generator & rng);
03548    void rowfill_randn(const int &i, const random_number_generator & rng);
03549 
03550    void colfill_randu_ni(const int &j, long int &n);
03551    void rowfill_randu_ni(const int &i, long int &n);
03552    void colfill_randn_ni(const int &j, long int &n);
03553    void fill_randn_ni(long int &n);
03554    void fill_randu_ni(long int &n);
03555    void rowfill_randn_ni(const int &i, long int &n);
03556 
03557 
03558 
03559    void colfill_seqadd(const int &, const int &, const int &);
03560    void colfill_seqadd(const int &, double, double);
03561    void rowfill_seqadd(const int &, double, double);
03562    void colfill(int j, const dvector & v);
03563    void rowfill(int j, const dvector & v);
03564 
03565 #if defined(OPT_LIB)
03566    inline dvector & operator() (register int i)
03567    {
03568       return m[i];
03569    }
03570    inline dvector & operator[] (register int i)
03571    {
03572       return m[i];
03573    }
03574 #else
03575    dvector & operator()(int i);
03576    dvector & operator[](int);
03577 #endif
03578 
03579 #if defined(OPT_LIB)
03580    inline const dvector & operator() (register int i) const
03581    {
03582       return m[i];
03583    }
03584    inline const dvector & operator[] (register int i) const
03585    {
03586       return m[i];
03587    }
03588 #else
03589    const dvector & operator() (int i) const;
03590    const dvector & operator[] (int) const;
03591 #endif
03592 
03593 #if defined(OPT_LIB) && !defined(__INTEL_COMPILER)
03594    inline double &operator() (register int i, register int j)
03595    {
03596       return (*(m[i].v + j));
03597    }
03598    inline const double &operator() (register int i, register int j) const
03599    {
03600       return (*(m[i].v + j));
03601    }
03602 #else
03603    double &operator() (int i, int j);
03604    const double &operator() (int i, int j) const;
03605 #endif
03606 
03607    inline dvector & elem(int i)
03608    {
03609       return (*(m + i));
03610    }
03611    inline double &elem(int i, int j)
03612    {
03613       return (*((*(m + i)).v + j));
03614    }
03615    inline const dvector & elem(int i) const
03616    {
03617       return (*(m + i));
03618    }
03619    inline const double &elem(int i, int j) const
03620    {
03621       return (*((*(m + i)).v + j));
03622    }
03623    friend class d3_array;
03624    friend dvector operator*(const dvector &, const dmatrix &);
03625 
03626    friend dvector operator*(const dmatrix &, const dvector &);
03627 
03628    friend dvar_vector operator*(const dvar_vector &, const dmatrix &);
03629 
03630    friend dvar_vector operator*(const dmatrix &, const dvar_vector &);
03631 
03632    friend dmatrix operator*(const dmatrix &, const dmatrix &);
03633 
03634    friend dvar_matrix operator*(const dvar_matrix &, const dmatrix &);
03635 
03636    friend dvar_matrix operator*(const dmatrix &, const dvar_matrix &);
03637 
03638    friend dvar_matrix::dvar_matrix(const dmatrix &);
03639 
03640    friend dmatrix operator-(const dmatrix &, const dmatrix &);
03641    friend dmatrix operator+(const dmatrix &, const dmatrix &);
03642 
03643    friend dvar_matrix operator+(const dvar_matrix &, const dmatrix &);
03644 
03645    friend dvar_matrix operator+(const dmatrix &, const dvar_matrix &);
03646 
03647    friend dmatrix trans(const dmatrix & m1);
03648 
03649    friend dmatrix inv(const dmatrix &);
03650    friend dmatrix inv(const dmatrix & m1, const double &_ln_det,
03651      const int &_sgn);
03652 
03653    friend double det(const dmatrix &);
03654    friend double ln_det(const dmatrix & m1, const int &sgn);
03655 
03656    friend double norm(const dmatrix &);
03657    friend double norm2(const dmatrix &);
03658    friend double sumsq(const dmatrix &);
03659 
03660    dmatrix & operator+=(const dmatrix & t);
03661    dmatrix & operator-=(const dmatrix & t);
03662 
03663    dmatrix & operator =(const dmatrix & t);
03664    dmatrix & operator =(double t);
03665 
03666    dmatrix operator() (const ivector & t);
03667 
03668    friend dvar_matrix & dvar_matrix::operator=(const dmatrix &);
03669 
03670    dmatrix(const tdmatrix & t);
03671 
03672    dmatrix & operator /=(double d);
03673    dmatrix & operator *=(double d);
03674 };
03675 
03676 imatrix operator*(const imatrix &, const imatrix &);
03677 
03678 dmatrix trans(const dmatrix & m1);
03679 
03680 imatrix trans(const imatrix & m1);
03681 
03682 dvariable dfatan1(dvariable, double, double, double *);
03683 
03684 double dftinv(double, double, double);
03685 
03686 dvariable boundp(double, double, double, double *);
03687 
03688 dvariable dfboundp(double, double, double, double *);
03689 dvariable dfboundp(const prevariable &, double, double);
03690 
03691 double mean(const dvector &);
03692 double mean(const dmatrix &);
03693 double mean(const d3_array &);
03694 
03695 double std_dev(const dvector &);
03696 double var(const dvector &);
03697 
03698 dvariable mean(const dvar_vector &);
03699 dvariable mean(const dvar_matrix &);
03700 dvariable mean(const dvar3_array &);
03701 
03702 dvariable std_dev(const dvar_vector &);
03703 dvariable var(const dvar_vector &);
03704 
03705 dvariable sum(const dvar_vector &);
03706 double sum(const dvector &);
03707 int sum(const ivector &);
03708 
03709 dvar_vector rowsum(const dvar_matrix &);
03710 dvar_vector colsum(const dvar_matrix &);
03711 
03712 dvector colsum(const dmatrix &);
03713 dvector rowsum(const dmatrix &);
03714 
03715 ivector colsum(const imatrix &);
03716 ivector rowsum(const imatrix &);
03717 
03718 int colsum(const imatrix &, int column);
03719 double colsum(const dmatrix &, int column);
03720 dvariable colsum(const dvar_matrix &, int column);
03721 
03722 double sfabs(double t1); //"smoothed absolute value function
03723 
03724 dvector sfabs(const dvector & t1); //"smoothed absolute value function
03725 
03727 
03732 class imatrix_position
03733 {
03734  public:
03735    int row_min;
03736    int row_max;
03737    ivector lb;
03738    ivector ub;
03739    ptr_vector ptr;
03740    imatrix_position(const imatrix &);
03741    imatrix_position(int min, int max);
03742    imatrix_position(const imatrix_position &);
03743    ivector_position operator () (int i);
03744 };
03745 
03746 imatrix_position restore_imatrix_position(void);
03747 imatrix restore_imatrix_value(const imatrix_position &);
03748 
03753 class imatrix
03754 {
03755  protected:
03756    int index_min;
03757    int index_max;
03758    ivector *m;
03759    mat_shapex *shape;
03760    friend char *fform(const char *, const dmatrix &);
03761 
03762  public:
03763    int operator!(void) const
03764    {
03765       return (shape == NULL);
03766    }
03767 
03768    imatrix(int, int);
03769    // makes a matrix [0..nr][0..nc]
03770 
03771    imatrix(int nrl, int nrh, const ivector & iv);
03772    void allocate(int nrl, int nrh, const ivector & iv);
03773 
03774    imatrix(int, int, int, int);
03775    // makes a matrix [nrl..nrh][ncl..nch]
03776    imatrix(int, int, int, const ivector &);
03777    imatrix sub(int, int);
03778    imatrix(int, int, const ivector &, const ivector &);
03779    imatrix(const ad_integer & nrl, const ad_integer & nrh,
03780      const index_type & ncl, const index_type & nch);
03781 
03782    imatrix & operator=(const imatrix & t);
03783    imatrix & operator =(const int);
03784    imatrix(const imatrix &);
03785    // copy initializer
03786    imatrix(const imatrix_position &);
03787    imatrix(void);
03788 
03789    ~imatrix();
03790    void shallow_copy(const imatrix &);
03791 
03792    void save_imatrix_value(void);
03793    void save_imatrix_position(void);
03794    imatrix restore_imatrix_value(const imatrix_position & mpos);
03795    imatrix_position restore_imatrix_position(void);
03796 
03797    void allocate(void);
03798    void allocate(const imatrix & dm);
03799    void allocate(int nrl, int nrh, int ncl, int nch);
03800    void allocate(int nrl, int nrh);
03801    void allocate(int nrl, int nrh, int ncl, const ivector & nch);
03802    void allocate(int nrl, int nrh, const ivector & ncl, const ivector & nch);
03803    void allocate(const ad_integer & nrl, const ad_integer & nrh,
03804      const index_type & ncl, const index_type & nch);
03805    void deallocate();
03806 
03807 #ifdef OPT_LIB
03808    inline ivector & operator() (int i)
03809    {
03810       return m[i];
03811    }
03812    inline int &operator() (int i, int j)
03813    {
03814       return (*((*(m + i)).v + j));
03815    }
03816    inline ivector & operator[] (int i)
03817    {
03818       return m[i];
03819    }
03820 #else
03821    ivector & operator [] (int);
03822    ivector & operator () (int);
03823    int &operator () (int, int);
03824 #endif
03825 
03826 #ifdef OPT_LIB
03827    inline const ivector & operator() (int i) const
03828    {
03829       return m[i];
03830    }
03831    inline const int &operator() (int i, int j) const
03832    {
03833       return (*((*(m + i)).v + j));
03834    }
03835    inline const ivector & operator[] (int i) const
03836    {
03837       return m[i];
03838    }
03839 #else
03840    const ivector & operator[] (int) const;
03841    const ivector & operator() (int) const;
03842    const int &operator() (int, int) const;
03843 #endif
03844    int indexmin(void) const
03845    {
03846       return index_min;
03847    }
03848    int indexmax(void) const
03849    {
03850       return index_max;
03851    }
03852    int rowmin(void) const
03853    {
03854       return index_min;
03855    }
03856    int rowmax(void) const
03857    {
03858       return index_max;
03859    }
03860    int colmin(void) const
03861    {
03862       return ((*this) (indexmin()).indexmin());
03863    }
03864    int colmax(void) const
03865    {
03866       return ((*this) (indexmin()).indexmax());
03867    }
03868    // returns the number of rows
03869    int rowsize() const
03870    {
03871       return (rowmax() - rowmin() + 1);
03872    }
03873    // returns the number of columns
03874    int colsize() const
03875    {
03876       return (colmax() - colmin() + 1);
03877    }
03878    void rowshift(int min);
03879    inline ivector & elem(int i)
03880    {
03881       return (*(m + i));
03882    }
03883    inline int &elem(int i, int j)
03884    {
03885       return (*((*(m + i)).v + j));
03886    }
03887    inline const ivector & elem(int i) const
03888    {
03889       return (*(m + i));
03890    }
03891    inline const int &elem(int i, int j) const
03892    {
03893       return (*((*(m + i)).v + j));
03894    }
03895 
03896 
03897    void write_on(const ostream &) const;
03898    void write_on(const uostream &) const;
03899    void read_from(const istream &);
03900    void read_from(const uistream &);
03901    void initialize(void);
03902    friend class i3_array;
03903    void fill_seqadd(int, int);
03904    void colfill_seqadd(int, int, int);
03905 };
03906 
03907 dvariable regression(const dvector & obs, const dvar_vector & pred);
03908 double regression(const dvector & obs, const dvector & pred);
03909 
03910 dvariable robust_regression_fixed(const dvector& obs, const dvar_vector& pred,
03911   double a = 0.7);
03912 dvariable robust_regression(const dvector & obs, const dvar_vector & pred,
03913   double a = 0.7);
03914 
03915 dvariable robust_regression(const dvector & obs, const dvar_vector & pred,
03916   const dvariable & cutoff);
03917 
03918 dmatrix column_vector(const dvector &);
03919 dmatrix row_vector(const dvector &);
03920 
03921 dvar_matrix column_vector(const dvar_vector &);
03922 dvar_matrix row_vector(const dvar_vector &);
03923 
03924 dmatrix identity_matrix(int min, int max);
03925 
03926 istream & operator>>(const istream & s, const ptr_vector & v);
03927 ostream & operator<<(const ostream & s, const ptr_vector & v);
03928 
03932 class fmm_control
03933 {
03934  public:
03935    int noprintx;
03936    long maxfn;
03937    long iprint;
03938    double crit;
03939    double fringe;
03940    long imax;
03941    double dfn;
03942    long ifn;
03943    long iexit;
03944    long ialph;
03945    long ihflag;
03946    long ihang;
03947    long scroll_flag;
03948    int maxfn_flag;
03949    int quit_flag;
03950    double min_improve;
03951    int ireturn;
03952    int dcheck_flag;
03953    int use_control_c;
03954 
03955    void set_defaults();
03956    fmm_control();
03957    fmm_control(const fmm_control &);
03958    fmm_control(const lvector & ipar);
03959    void writeon(const ostream & s) const;
03960 };
03961 
03966 class sdmatrix:public dmatrix
03967 {
03968  public:
03969    void allocate(int);
03970    void allocate();
03971    sdmatrix(int);
03972    sdmatrix();
03973    ~sdmatrix();
03974    void deallocate();
03975 };
03976 
03977 class dfsdmat;
03978 
03979 uistream & operator>>(const uistream &, const dfsdmat &);
03980 uostream & operator<<(const uostream &, const dfsdmat &);
03981 
03986 class dfsdmat
03987 {
03988    int tmp_file;
03989    int disk_save_flag;
03990    double *ptr;
03991    double **m;
03992    double *minp;
03993    double *maxp;
03994    int shared_memory;
03995    int n;
03996  public:
03997    int disk_save(void)
03998    {
03999       return disk_save_flag;
04000    }
04001    void save(void);
04002    void restore(void);
04003    double *getminp(void)
04004    {
04005       return minp;
04006    }
04007    int operator!(void) const
04008    {
04009       return (ptr == NULL);
04010    }
04011    int size(void)
04012    {
04013       return n;
04014    }
04015    dfsdmat(int n);
04016    dfsdmat();
04017    dfsdmat(int n, const gradient_structure & gs);
04018    void allocate(int n);
04019    void allocate(int n, const gradient_structure & gs);
04020    void allocate(void);
04021    ~dfsdmat();
04022    void deallocate(void);
04023    friend uistream & operator>>(const uistream &, const dfsdmat &);
04024    friend uostream & operator<<(const uostream &, const dfsdmat &);
04025 
04026 #if defined(OPT_LIB) && !defined(__INTEL_COMPILER)
04027    double &elem(int i, int j)
04028    {
04029       return *(m[i] + j);
04030    }
04031    double &operator () (int i, int j)
04032    {
04033       return *(m[i] + j);
04034    }
04035 #else
04036    double &elem(int i, int j);
04037    double &operator () (int i, int j);
04038 #endif
04039 };
04040 
04044 class fmm:public fmm_control
04045 {
04046  private:
04047    dfsdmat h;
04048    dvector w;
04049    dvector funval;
04050  public:
04051    double dmin, fbest, df;
04052 
04053    long int llog, n1, ic, iconv, i1, xxlink;
04054    double z, zz, gys, gs, sig, gso, alpha, tot, fy, dgs;
04055    long int itn, icc, np, nn, is, iu, iv, ib;
04056    int i, j;
04057    double gmax; 
04058    double fsave;
04059    dvector xx;
04060    dvector gbest;
04061    dvector xsave;
04062    dvector gsave;
04063 
04064    int n;
04065    int disk_save;
04066 
04067  public:
04068    fmm(int nvar, int disk_save = 0);
04069    fmm(int nvar, const lvector & ipar, int disk_save = 0);
04070    double minimize(const independent_variables & x,
04071      double (*pf) (const dvar_vector &));
04072 
04074    double minimize(const independent_variables & x, const dvector & c,
04075      double (*pf) (const dvar_vector &, const dvector &));
04076 
04077  //void fmin(const double& f, const independent_variables &x, const dvector& g);
04078    void fmin(const double &f, const dvector & x, const dvector & g);
04079 
04080    dmatrix & hessian(); 
04081 };
04082 
04083 class function_minimizer;
04084 
04089 class fmmt1:public fmm_control
04090 {
04091  private:
04092    dvector w;
04093    dvector funval;
04094    int xm;
04095    dmatrix xstep;
04096    dvector xrho;
04097    dvector rrr;
04098    dmatrix xy;
04099    dvector xold;
04100    dvector gold;
04101  public:
04102    double dmin, fbest, df;
04103 
04104    long int llog, n1, ic, iconv, i1, link;
04105    double z, zz, gys, gs, sig, gso, alpha, tot, fy, dgs;
04106    long int itn, icc, np, nn, is, iu, iv, ib;
04107    int i, j;
04108    double gmax;
04109    double fsave;
04110    dvector xx;
04111    dvector gbest;
04112    dvector xsave;
04113    dvector gsave;
04114 
04115    int n;
04116 
04117  public:
04118    fmmt1(int nvar, int _xm = 7);
04119    fmmt1(int nvar, const lvector & ipar);
04120    double minimize(const independent_variables & x,
04121      double (*pf) (const dvar_vector &));
04122 
04123    double minimize(const independent_variables & x, const dvector & c,
04124      double (*pf) (const dvar_vector &, const dvector &));
04125 
04126    void fmin2(const double &f, const independent_variables & x,
04127      const dvector & g, function_minimizer *);
04128 
04129    void fmin(const double &f, const dvector & x, const dvector & g);
04130 
04131 //  dmatrix& hessian();
04132 };
04133 
04134 void derch(const double &f, const independent_variables & x,
04135   const dvector & g, int n, const int &ireturn);
04136 
04137 void fmin(double f, const independent_variables & x, const dvector & g,
04138   const int &n, const dvector & w, const dvector & h, const fmm_control & fmc);
04139 
04140 void fmmdisp(const dvector & x, const dvector & g, const int &nvar,
04141   int scroll_flag, int noprintx = 0);
04142 
04143 void fmmdisp(const double *x, const double *g, const int &nvar,
04144   int scroll_flag, int noprintx = 0);
04145 
04146 ostream & operator<<(const ostream & s, const fmm_control & fmc);
04147 
04152 class uostream:public ofstream
04153 {
04154  public:
04155 #if defined(__TURBOC__) && (__BORLANDC__  <= 0x0520)
04156   uostream(const char*, int = ios::out | ios::binary, int = filebuf::openprot);
04157   void open(const char*, int = ios::out | ios::binary, int = filebuf::openprot);
04158 #endif
04159 #if (__BORLANDC__  >= 0x0540)
04160    uostream(const char *, int = ios::out | ios::binary, int protection = 666);
04161    void open(const char *, int = ios::out | ios::binary, int protection = 666);
04162 #endif
04163 #if defined (_MSC_VER) || defined (__WAT32__)
04164   #if (_MSC_VER < 1300)
04165   uostream(const char*, int = ios::out | ios::binary, int = filebuf::openprot);
04166   void open(const char*, int = ios::out | ios::binary, int = filebuf::openprot);
04167   #else
04168   uostream(const char *, int = ios::out | ios::binary, int prot = 0664);
04169   void open(const char *, int = ios::out | ios::binary, int prot = 0664);
04170   #endif
04171 #endif
04172 
04173 #ifdef __ZTC__
04174    uostream(const char *, int = ios::out, int = filebuf::openprot);
04175    void open(const char *, int = ios::out, int = filebuf::openprot);
04176 #endif
04177 
04178 #ifdef __NDPX__
04179    uostream(const char *, int = ios::out, int = filebuf::openprot);
04180    void open(const char *, int = ios::out, int = filebuf::openprot);
04181 #endif
04182 
04183 #ifdef __SUN__
04184    //uostream(const char*, int = ios::out, int = openprot);
04185    //void open(const char*, int = ios::out, int = openprot);
04186 #endif
04187 
04188 #if !defined(_MSC_VER)
04189   #if defined(__ADSGI__)
04190   uostream(const char *name, int mode = ios::out, int prot = 0664);
04191   void open(const char *name, int mode = ios::out, int prot = 0664);
04192   #else
04193   uostream(const char *name, int mode = ios::out | ios::binary,
04194     int prot = 0664);
04195   void open(const char *name, int mode = ios::out | ios::binary,
04196     int prot = 0664);
04197   #endif
04198 #endif
04199 
04200    // insert character
04201 #ifndef __SUN__
04202    uostream & operator<<(signed char);
04203 #endif
04204    uostream & operator<<(unsigned char);
04205 
04206    // insert numeric value
04207    uostream & operator<<(short);
04208    uostream & operator<<(unsigned short);
04209    uostream & operator<<(int);
04210    uostream & operator<<(unsigned int);
04211    uostream & operator<<(long);
04212    uostream & operator<<(unsigned long);
04213    uostream & operator<<(float);
04214    uostream & operator<<(double);
04215    uostream & operator<<(const char *)
04216    {
04217       return *this;
04218    };
04219 #ifdef __TURBOC__
04220    uostream & operator<<(long double);
04221 #endif
04222 
04223 
04224    // insert pointer
04225    uostream & operator<<(void *);
04226 
04227    virtual void sss(void);
04228 };
04229 
04230 
04231  // inline void uostream::open(const char* name, int m, int prot)
04232  // {
04233  // #if defined (__TURBOC__) &&   (__BORLANDC__  <= 0x0520)
04234  //   fstreambase::open(name, m, prot);
04235  // #endif
04236  // #if (__BORLANDC__  >= 0x0540 && __BORLANDC__  <= 0x0550)
04237  //   ofstream::open(name, m, prot);
04238  // #else
04239  // #  if defined(linux)
04240  // #    if (__GNUC__  >= 3)
04241  //        ofstream::open(name, std::_Ios_Openmode(m));
04242  // #    else
04243  //        ofstream::open(name, m);
04244  // #    endif
04245  // #  else
04246  //      ofstream::open(name, m);
04247  // #  endif
04248  // #endif
04249  //
04250  // #ifdef _MSC_VER
04251  // #  if (_MSC_VER >= 1400)
04252  //   ofstream::open(name, m);
04253  // #  else
04254  //   //fstreambase::open(name, m, prot);
04255  //   ofstream::open(name, m, prot);
04256  // #  endif
04257  // #endif
04258  // #ifdef __ZTC__
04259  //   fstream_common::open(name, m, prot);
04260  // #endif
04261  // #ifdef __NDPX__
04262  //   ofstream::open(name, m, prot);
04263  // #endif
04264  // #ifdef __SUN__
04265  //   ofstream::open(name, m, prot);
04266  // #endif
04267  // }
04268 
04273 class uistream:public ifstream
04274 {
04275  public:
04276 #if defined (__TURBOC__) &&   (__BORLANDC__  <= 0x0520)
04277    uistream(const char *, int = ios::in | ios::binary, int = filebuf::openprot);
04278   void open(const char *, int = ios::in | ios::binary, int = filebuf::openprot);
04279 #endif
04280 #if (__BORLANDC__  >= 0x0540)
04281    uistream(const char *, int = ios::in | ios::binary, int protection = 666);
04282    void open(const char *, int = ios::in | ios::binary, int protection = 666);
04283 #endif
04284 #if defined (_MSC_VER) || defined (__WAT32__)
04285   #if (_MSC_VER < 1300)
04286   uistream(const char *, int = ios::in | ios::binary, int = filebuf::openprot);
04287   void open(const char *, int = ios::in | ios::binary, int = filebuf::openprot);
04288   #else
04289   uistream(const char *, int = ios::in | ios::binary, int prot = 0664);
04290   void open(const char *, int = ios::in | ios::binary, int prot = 0664);
04291   #endif
04292 #endif
04293 #ifdef __ZTC__
04294    uistream(const char *, int = ios::in, int = filebuf::openprot);
04295    void open(const char *, int = ios::in, int = filebuf::openprot);
04296 #endif
04297 
04298 #ifdef __NDPX__
04299    uistream(const char *, int = ios::in, int = filebuf::openprot);
04300    void open(const char *, int = ios::in, int = filebuf::openprot);
04301 #endif
04302 
04303 #ifdef __SUN__
04304    // uistream(const char* name, int mode = ios::in, int prot=0664);
04305    // void open(const char* name, int mode = ios::in, int prot=0664);
04306 #endif
04307 
04308 #if !defined(_MSC_VER)
04309   #if defined(__ADSGI__)
04310   uistream(const char *name, int mode = ios::in, int prot = 0664);
04311   void open(const char *name, int mode = ios::in, int prot = 0664);
04312   #else
04313   uistream(const char *name, int mode = ios::in | ios::binary,
04314     int prot = 0664);
04315   void open(const char *name, int mode = ios::in | ios::binary,
04316     int prot = 0664);
04317   #endif
04318 #endif
04319 
04320    // extract characters into an array
04321 #ifndef __SUN__
04322    uistream & get(signed char *, int, char = '\n');
04323 #endif
04324    uistream & get(unsigned char *, int, char = '\n');
04325 
04326    // extract a single character
04327    uistream & get(unsigned char &);
04328 #ifndef __SUN__
04329    uistream & get(signed char &);
04330 #endif
04331    int get();
04332 
04333 
04334    // extract and discard chars but stop at delim
04335    uistream & ignore(int = 1, int = EOF);
04336 
04337 #ifndef __SUN__
04338    uistream & operator>>(const signed char *);
04339 #endif
04340    uistream & operator>>(const unsigned char *);
04341    uistream & operator>>(const unsigned char &);
04342 #ifndef __SUN__
04343    uistream & operator>>(const signed char &);
04344 #endif
04345    uistream & operator>>(const short &);
04346    uistream & operator>>(const int &);
04347    uistream & operator>>(const long &);
04348    uistream & operator>>(const unsigned short &);
04349    uistream & operator>>(const unsigned int &);
04350    uistream & operator>>(const unsigned long &);
04351    uistream & operator>>(const float &);
04352    uistream & operator>>(const double &);
04353    uistream & operator>>(const char &);
04354 #if defined(__TURBOC__) || defined (_MSC_VER)
04355    uistream & operator>>(const long double &);
04356 #endif
04357    virtual void sss(void);
04358 };
04359 
04360   // inline void   uistream::open(const char* name, int m, int prot)
04361   // {
04362   // #if defined(__TURBOC__) && (__BORLANDC__  <= 0x0520)
04363   //   fstreambase::open(name, m, prot);
04364   // #endif
04365   // #ifdef __ZTC__
04366   //   fstream_common::open(name, m, prot);
04367   // #endif
04368   // #ifdef __NDPX__
04369   //   ifstream::open(name, m, prot);
04370   // #endif
04371   // #ifdef __SUN__
04372   //   ifstream::open(name, m, prot);
04373   // #endif
04374   // }
04375 
04376 class fmmc;
04377 
04378 void derch(const double &f, const dvector & x, const dvector & gg, int n,
04379   const int &ireturn);
04380 
04385 class fmmc
04386 {
04387  public:
04388    int maxfn;
04389    double crit;
04390    double min_improve;
04391    int iprint;
04392    int scroll_flag;
04393    int j;
04394    int J;
04395    long int ifn;
04396    long int iter;
04397    int imax;
04398    int ihang;
04399    int quit_flag;
04400    dvector *funval;
04401    dvector *left_bracket_gradient;
04402    dvector *right_bracket_gradient;
04403    dvector *g;
04404    dvector *h;
04405    dvector *xi;
04406    dvector *d;
04407    dvector *extx;
04408    dvector *g2;
04409    dvector *grad;
04410    dvector *extg;
04411    dvector *theta;
04412    dvector *xbest;
04413    dvector *gbest;
04414    int lin_flag;
04415    int ext_flag;
04416    int int_flag;
04417    int ifnex;
04418    int ireturn;
04419    int frp_flag;
04420    double gg;
04421    double gam;
04422    double fp;
04423    double dgg;
04424    double rho_min;
04425    double converge_flag;
04426    double gamma;
04427    double Psi_0;
04428    double extf;
04429    double crit1;
04430    double rho_1;
04431    double Psi_1;
04432    double dir_deriv;
04433    double rho_i;
04434    double left_bracket;
04435    double left_bracket_value;
04436    double right_bracket;
04437    double right_bracket_value;
04438    double rho_0;
04439    double fbest;
04440    fmmc(const int &n);
04441    ~fmmc();
04442    void fmin(const double &f, const dvector & p, const dvector & gg);
04443    double dfn;
04444    int maxfn_flag;
04445    int iexit;
04446    int ihflag;
04447 };
04448 
04449 class dd3_array;
04450 
04455 class three_array_shape
04456 {
04457    //unsigned int nslices;
04458    unsigned int ncopies;
04459    //unsigned int nrows;
04460    //unsigned int ncols;
04461    int slice_min;
04462    int slice_max;
04463    // int row_min;
04464    // int row_max;
04465    //int col_min;
04466    //int col_max;
04467    three_array_shape(int sl, int sh);
04468    //mat_shape(){};
04469 
04470    friend class i3_array;
04471    friend class d3_array;
04472    friend class dd3_array;
04473    friend class qd3_array;
04474    friend class dvar3_array;
04475 };
04476 
04477 //class dmatrix_ptr { dmatrix *p; };
04478 //class dvar_matrix_ptr { dvar_matrix *p; };
04479 
04484 class d3_array
04485 {
04486    dmatrix *t;
04487    three_array_shape *shape;
04488    friend class d4_array;
04489  public:
04490    int operator!(void) const
04491    {
04492       return (shape == NULL);
04493    }
04494    // conclass cgors
04495    d3_array(void);
04496    void save_d3_array_value(void) const;
04497    void shallow_copy(const d3_array &);
04498    d3_array sub(int, int);
04499    d3_array(int sl, int sh, int nrl, int nrh, int ncl, int nch);
04500    d3_array(int sl, int sh, int nrl, int nrh);
04501    d3_array(int sl, int sh, const index_type & nrl, const index_type & nrh);
04502    d3_array(int sl, int sh);
04503    d3_array(const d3_array_position &);
04504 
04505    void save_d3_array_position(void) const;
04506 
04507    d3_array(int sl, int sh, int nrl, int nrh, const ivector & ncl, int nch);
04508 
04509    d3_array(const ad_integer & sl, const ad_integer & sh,
04510      const index_type & nrl, const index_type & nrh,
04511      const index_type & ncl, const index_type & nch);
04512 
04513    void allocate(const ad_integer & sl, const ad_integer & sh,
04514      const index_type & nrl, const index_type & nrh,
04515      const index_type & ncl, const index_type & nch);
04516 
04517    d3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04518      const imatrix & ncl, const imatrix & nch);
04519    d3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04520      int ncl, const imatrix & nch);
04521    d3_array(int sl, int sh, int nrl, const ivector & nrh,
04522      int ncl, const imatrix & nch);
04523    d3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04524      const ivector & ncl, const ivector & nch);
04525    d3_array(int sl, int sh, int nrl, int nrh, const ivector & ncl,
04526      const ivector & nch);
04527    d3_array(int sl, int sh, int nrl, const ivector & nrh,
04528      int ncl, const ivector & nch);
04529    d3_array(int sl, int sh, int nrl, const ivector & nrh,
04530      int ncl, int nch);
04531    d3_array(const d3_array & m2);
04532    ~d3_array();
04533 
04534    void allocate(const dvar3_array &);
04535    void allocate(const d3_array & d3v);
04536    void allocate(int sl, int sh, int nrl, int nrh, int ncl, int nch);
04537    void allocate(int sl, int sh, int nrl, int nrh);
04538    void allocate(int sl, int sh, const index_type & nrl,
04539      const index_type & nrh);
04540    void allocate(int sl, int sh);
04541 
04542    void allocate(int sl, int sh, int nrl, int nrh, const ivector& ncl, int nch);
04543    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04544      const imatrix & ncl, const imatrix & nch);
04545    void allocate(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04546      const imatrix & nch);
04547    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04548      int ncl, const imatrix & nch);
04549    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04550      int ncl, int nch);
04551    void allocate(int sl, int sh, const ivector& nrl, int nrh, int ncl, int nch);
04552    void allocate(int sl, int sh, int nrl, const ivector& nrh, int ncl, int nch);
04553    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04554      const ivector & ncl, const ivector & nch);
04555    void allocate(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04556      const ivector & nch);
04557    void allocate(int sl, int sh, int nrl, int nrh, const ivector & ncl,
04558      const ivector & nch);
04559    void allocate(int sl, int sh, int nrl, int nrh, int ncl,
04560      const ivector & nch);
04561    void allocate(void);
04562    void deallocate(void);
04563    void initialize(int sl, int sh, int nrl, const ivector & nrh,
04564      int ncl, const ivector & nch);
04565 
04566    //access functions
04567    int indexmin(void) const
04568    {
04569       return shape->slice_min;
04570    }
04571    int indexmax(void) const
04572    {
04573       return shape->slice_max;
04574    }
04575    int slicemin(void) const
04576    {
04577       return shape->slice_min;
04578    }
04579    int slicemax(void) const
04580    {
04581       return shape->slice_max;
04582    }
04583    int colmin(void) const
04584    {
04585       return ((*this) (slicemin()).colmin());
04586    }
04587    int colmax(void) const
04588    {
04589       return ((*this) (slicemin()).colmax());
04590    }
04591    int rowmin(void) const
04592    {
04593       return ((*this) (slicemin()).rowmin());
04594    }
04595    int rowmax(void) const
04596    {
04597       return ((*this) (slicemin()).rowmax());
04598    }
04599    // returns the number of rows
04600    int slicesize(void) const
04601    {
04602       return (slicemax() - slicemin() + 1);
04603    }
04604    // returns the number of rows
04605    int rowsize(void) const
04606    {
04607       return (rowmax() - rowmin() + 1);
04608    }
04609    // returns the number of columns
04610    int colsize(void) const
04611    {
04612       return (colmax() - colmin() + 1);
04613    }
04614    void initialize(void);
04615 
04616    dmatrix & elem(int k)
04617    {
04618       return (t[k]);
04619    }
04620    const dmatrix & elem(int k) const
04621    {
04622       return t[k];
04623    }
04624 #ifdef OPT_LIB
04625    inline const double &operator() (int k, int i, int j) const
04626    {
04627       return (((t[k].m[i]).v)[j]);
04628    }
04629 
04630    inline const dvector & operator() (int k, int i) const
04631    {
04632       return (t[k].m[i]);
04633    }
04634 
04635    inline const dmatrix & operator() (int i) const
04636    {
04637       return (t[i]);
04638    }
04639 
04640    inline const dmatrix & operator[] (int i) const
04641    {
04642       return (t[i]);
04643    }
04644 
04645    inline double &operator () (int k, int i, int j)
04646    {
04647       return (((t[k].m[i]).v)[j]);
04648    }
04649 
04650    inline dvector & operator () (int k, int i)
04651    {
04652       return (t[k].m[i]);
04653    }
04654 
04655    inline dmatrix & operator() (int i)
04656    {
04657       return (t[i]);
04658    }
04659 
04660    inline dmatrix & operator[] (int i)
04661    {
04662       return (t[i]);
04663    }
04664 #else
04665    double &operator () (int k, int i, int j);
04666    dvector & operator ()(int k, int i);
04667    dmatrix & operator[](int i);
04668    dmatrix & operator()(int i);
04669    const double &operator() (int k, int i, int j) const;
04670    const dvector & operator() (int k, int i) const;
04671    const dmatrix & operator[] (int i) const;
04672    const dmatrix & operator() (int i) const;
04673 #endif
04674 
04675    d3_array & operator=(const d3_array & m1);
04676    d3_array & operator=(double x);
04677    friend d3_array value(const dvar3_array & ar);
04678 
04679    void fill_randn(const random_number_generator & rng);
04680    void fill_randcau(const random_number_generator & rng);
04681    void fill_randu(const random_number_generator & rng);
04682 
04683    void fill_randu(long int &n);
04684    void fill_randn(long int &n);
04685 
04686    void fill_randu_ni(long int &n);
04687    void fill_randn_ni(long int &n);
04688    double fill_seqadd(double, double);
04689    void operator /=(double d);
04690 };
04691 
04696 class i3_array
04697 {
04698    imatrix *t;
04699    three_array_shape *shape;
04700  public:
04701 #  if defined(MFCL2_CONSTRUCTORS)
04702    i3_array(int sl, int sh, int nrl, int nrh, const ivector & nc);
04703    void allocate(int sl, int sh, int nrl, int nrh, const ivector & nc);
04704 #  endif
04705    int operator!(void) const
04706    {
04707       return (shape == NULL);
04708    }
04709    // conclass cgors
04710    void shallow_copy(const i3_array &);
04711    i3_array(void);
04712    i3_array(int sl, int sh, const index_type & nrl, const index_type & nrh,
04713      const index_type & ncl, const index_type & nch);
04714 
04715    i3_array(int _sl, int _sh, const imatrix & m1);
04716 
04717    i3_array(int sl, int sh);
04718    i3_array(int sl, int sh, int nrl, int nrh, int ncl, int nch);
04719    i3_array(int sl, int sh, int nrl, int nrh, const ivector & ncl, int nch);
04720 
04721    i3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04722      const imatrix & ncl, const imatrix & nch);
04723    i3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04724      int ncl, const imatrix & nch);
04725    i3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04726      const ivector & ncl, const ivector & nch);
04727    i3_array(int sl, int sh, int nrl, int nrh, const ivector & ncl,
04728      const ivector & nch);
04729    i3_array(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04730      const ivector & nch);
04731    i3_array(int sl, int sh, int nrl, const ivector & nrh, int ncl, int nch);
04732    i3_array(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04733      const imatrix & nch);
04734    i3_array(const i3_array & m2);
04735    ~i3_array();
04736 
04737    void allocate(int sl, int sh, int nrl, const ivector& nrh, int ncl, int nch);
04738    void allocate(const dvar3_array &);
04739    void allocate(const i3_array & i3v);
04740    void allocate(int sl, int sh, int nrl, int nrh, int ncl, int nch);
04741    void allocate(int sl, int sh);
04742 
04743    void allocate(int sl, int sh, int nrl, int nrh, const ivector& ncl, int nch);
04744 
04745    void allocate(int sl, int sh, const index_type& nrl, const index_type& nrh,
04746      const index_type & ncl, const index_type & nch);
04747 
04748    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04749      const imatrix & ncl, const imatrix & nch);
04750    void allocate(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04751      const imatrix & nch);
04752    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04753      int ncl, const imatrix & nch);
04754    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04755      int ncl, int nch);
04756    void allocate(int sl, int sh, const ivector & nrl, int nrh,
04757      int ncl, int nch);
04758    //void allocate(int sl, int sh, int nrl, const ivector& nrh,
04759    //  int ncl,int nch);
04760    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04761      const ivector & ncl, const ivector & nch);
04762    void allocate(int sl, int sh, int nrl, const ivector & nrh,
04763      int ncl, const ivector & nch);
04764    void allocate(int sl, int sh, int nrl, int nrh, const ivector & ncl,
04765      const ivector & nch);
04766    void allocate(int sl, int sh, int nrl, int nrh, int ncl, const ivector& nch);
04767    void allocate(void);
04768    void deallocate(void);
04769    void initialize(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04770      const ivector & nch);
04771 
04772    //access functions
04773    int indexmin(void) const
04774    {
04775       return shape->slice_min;
04776    }
04777    int indexmax(void) const
04778    {
04779       return shape->slice_max;
04780    }
04781    int slicemin(void) const
04782    {
04783       return shape->slice_min;
04784    }
04785    int slicemax(void) const
04786    {
04787       return shape->slice_max;
04788    }
04789    int colmin(void) const
04790    {
04791       return ((*this) (slicemin()).colmin());
04792    }
04793    int colmax(void) const
04794    {
04795       return ((*this) (slicemin()).colmax());
04796    }
04797    int rowmin(void) const
04798    {
04799       return ((*this) (slicemin()).rowmin());
04800    }
04801    int rowmax(void) const
04802    {
04803       return ((*this) (slicemin()).rowmax());
04804    }
04805    // returns the number of rows
04806    int slicesize() const
04807    {
04808       return (slicemax() - slicemin() + 1);
04809    }
04810    // returns the number of rows
04811    int rowsize() const
04812    {
04813       return (rowmax() - rowmin() + 1);
04814    }
04815    // returns the number of columns
04816    int colsize() const
04817    {
04818       return (colmax() - colmin() + 1);
04819    }
04820    void initialize(void);
04821 
04822    imatrix & elem(int k)
04823    {
04824       return (t[k]);
04825    }
04826    const imatrix & elem(int k) const
04827    {
04828       return t[k];
04829    }
04830 #ifdef OPT_LIB
04831    inline const int &operator() (int k, int i, int j) const
04832    {
04833       return (((t[k].m[i]).v)[j]);
04834    }
04835 
04836    inline const ivector & operator() (int k, int i) const
04837    {
04838       return (t[k].m[i]);
04839    }
04840 
04841    inline const imatrix & operator() (int i) const
04842    {
04843       return (t[i]);
04844    }
04845 
04846    inline const imatrix & operator[] (int i) const
04847    {
04848       return (t[i]);
04849    }
04850 
04851    inline int &operator() (int k, int i, int j)
04852    {
04853       return (((t[k].m[i]).v)[j]);
04854    }
04855 
04856    inline ivector & operator () (int k, int i)
04857    {
04858       return (t[k].m[i]);
04859    }
04860 
04861    inline imatrix & operator() (int i)
04862    {
04863       return (t[i]);
04864    }
04865 
04866    inline imatrix & operator[] (int i)
04867    {
04868       return (t[i]);
04869    }
04870 #else
04871    int &operator () (int k, int i, int j);
04872    ivector & operator ()(int k, int i);
04873    imatrix & operator[](int i);
04874    imatrix & operator()(int i);
04875    const int &operator() (int k, int i, int j) const;
04876    const ivector & operator() (int k, int i) const;
04877    const imatrix & operator[] (int i) const;
04878    const imatrix & operator() (int i) const;
04879 #endif
04880 
04881    i3_array & operator=(const i3_array & m1);
04882    i3_array & operator=(int x);
04883 
04884    void fill_randu(long int &n);
04885    void fill_randn(long int &n);
04886    void fill_randu_ni(long int &n);
04887    void fill_randn_ni(long int &n);
04888 };
04889 
04890 #   if defined(__NUMBERVECTOR__)
04891 class param_init_matrix_vector;
04892 class param_init_bounded_matrix_vector;
04893 #   endif
04894 
04899 class dvar3_array
04900 {
04901    dvar_matrix *t;
04902    three_array_shape *shape;
04903 
04904  public:
04905    void shallow_copy(const dvar3_array &);
04906    dvar3_array sub(int, int);
04907    dvar3_array(int, int);
04908    int operator!(void) const
04909    {
04910       return (shape == NULL);
04911    }
04912    // conclass cgors
04913 
04914    void initialize(void);
04915    void allocate(int sl, int sh, int nrl, int nrh, int ncl, int nch);
04916    void allocate(int sl, int sh, int nrl, int nrh);
04917    void allocate(int sl, int sh, const index_type& nrl, const index_type& nrh);
04918    void allocate(int sl, int sh);
04919 
04920    void allocate(int sl, int sh, int nrl, int nrh, const ivector& ncl, int nch);
04921    void allocate(const d3_array & m1);
04922    void allocate(void);
04923    void allocate(const dvar3_array & m1);
04924    void allocate(int sl, int sh, int nrl, int nrh,
04925      const ivector & ncl, const ivector & nch);
04926    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04927      const ivector & ncl, const ivector & nch);
04928    void allocate(int sl, int sh, const ivector & nrl, const ivector & nrh,
04929      int ncl, int nch);
04930    void allocate(int sl, int sh, const ivector & nrl, int nrh,
04931      int ncl, int nch);
04932    void allocate(int sl, int sh, int nrl, const ivector & nrh,
04933      int ncl, int nch);
04934    void allocate(int sl, int sh, int nrl, const ivector & nrh,
04935      int ncl, const ivector & nch);
04936    void allocate(int sl, int sh, int nrl, int nrh,
04937      int ncl, const ivector & nch);
04938    void allocate(ad_integer sl, ad_integer sh, const index_type & nrl,
04939      const index_type & nrh, const index_type & ncl, const index_type & nch);
04940    void allocate(ad_integer sl, ad_integer sh, const index_type & nrl,
04941      const index_type & nrh);
04942    void allocate(ad_integer sl, ad_integer sh);
04943 
04944    void deallocate();
04945    dvar3_array(int sl, int sh, int nrl, int nrh, int ncl, int nch);
04946    dvar3_array(int sl, int sh, int nrl, int nrh, const ivector & ncl, int nch);
04947 
04948    dvar3_array(int sl, int sh, const ivector & nrl, const ivector & nrh,
04949      ivector & ncl, const ivector & nch);
04950 
04951    dvar3_array(int sl, int sh, int nrl, const ivector & nrh, int ncl,
04952      const ivector & nch);
04953 
04954    dvar3_array(int sl, int sh, int nrl, const ivector & nrh, int ncl, int nch);
04955 
04956    dvar3_array(ad_integer sl, ad_integer sh, const index_type & nrl,
04957      const index_type & nrh, const index_type & ncl, const index_type & nch);
04958 
04959    dvar3_array(const d3_array & m2);
04960 #   if defined(__NUMBERVECTOR__)
04961    dvar3_array(const param_init_matrix_vector & m2);
04962    dvar3_array(const param_init_bounded_matrix_vector & m2);
04963 #   endif
04964 
04965    dvar3_array(const dvar3_array & m2);
04966 
04967    dvar3_array(void);
04968 
04969    ~dvar3_array();
04970 
04971    d3_array value(const dvar3_array &);
04972 
04973    //access functions
04974    int indexmin(void) const
04975    {
04976       return shape->slice_min;
04977    }
04978    int indexmax(void) const
04979    {
04980       return shape->slice_max;
04981    }
04982    int slicemin(void) const
04983    {
04984       return shape->slice_min;
04985    }
04986    int slicemax(void) const
04987    {
04988       return (shape->slice_max);
04989    }
04990    int colmin(void) const
04991    {
04992       return ((*this) (slicemin()).colmin());
04993    }
04994    int colmax(void) const
04995    {
04996       return ((*this) (slicemin()).colmax());
04997    }
04998    int rowmin(void) const
04999    {
05000       return ((*this) (slicemin()).rowmin());
05001    }
05002    int rowmax(void) const
05003    {
05004       return ((*this) (slicemin()).rowmax());
05005    }
05006    // returns the number of rows
05007    int slicesize() const
05008    {
05009       return (slicemax() - slicemin() + 1);
05010    }
05011    // returns the number of rows
05012    int rowsize() const
05013    {
05014       return (rowmax() - rowmin() + 1);
05015    }
05016    // returns the number of columns
05017    int colsize() const
05018    {
05019       return (colmax() - colmin() + 1);
05020    }
05021    dvar_matrix & elem(int k)
05022    {
05023       return (t[k]);
05024    }
05025    prevariable elem(int i, int j, int k)
05026    {
05027       return (t[k].elem(i, j));
05028    }
05029    const dvar_matrix & elem(int k) const
05030    {
05031       return t[k];
05032    }
05033    const prevariable elem(int i, int j, int k) const
05034    {
05035       return t[k].elem(i, j);
05036    }
05037 
05038 #ifdef OPT_LIB
05039    inline const prevariable operator() (int k, int i, int j) const
05040    {
05041       return (((t[k].m[i]).va) + j);
05042    }
05043 
05044    inline const dvar_vector & operator() (int k, int i) const
05045    {
05046       return (t[k].m[i]);
05047    }
05048 
05049    inline const dvar_matrix & operator() (int i) const
05050    {
05051       return (t[i]);
05052    }
05053 
05054    inline const dvar_matrix & operator[] (int i) const
05055    {
05056       return (t[i]);
05057    }
05058 
05059    inline prevariable operator () (int k, int i, int j)
05060    {
05061       return (((t[k].m[i]).va) + j);
05062    }
05063 
05064    inline dvar_vector & operator () (int k, int i)
05065    {
05066       return (t[k].m[i]);
05067    }
05068 
05069    inline dvar_matrix & operator() (int i)
05070    {
05071       return (t[i]);
05072    }
05073 
05074    inline dvar_matrix & operator[] (int i)
05075    {
05076       return (t[i]);
05077    }
05078 #else
05079    prevariable operator () (int k, int i, int j);
05080    dvar_vector & operator ()(int k, int i);
05081    dvar_matrix & operator[](int i);
05082    dvar_matrix & operator()(int i);
05083    const prevariable operator() (int k, int i, int j) const;
05084    const dvar_vector & operator() (int k, int i) const;
05085    const dvar_matrix & operator[] (int i) const;
05086    const dvar_matrix & operator() (int i) const;
05087 #endif
05088 
05089    dvar3_array & operator=(const d3_array & m1);
05090    dvar3_array & operator=(double x);
05091    dvar3_array & operator=(const dvar3_array & m1);
05092 
05093    void fill_randu(long int &n);
05094    void fill_randn(long int &n);
05095 
05096    void fill_randu_ni(long int &n);
05097    void fill_randn_ni(long int &n);
05098    double fill_seqadd(double, double);
05099    void operator/=(const prevariable &);
05100    void operator /=(double);
05101 };
05102 
05103 dvariable inv_cumd_exponential(const prevariable & y);
05104 dvariable cumd_exponential(const prevariable & x);
05105 
05106 double inv_cumd_exponential(double y);
05107 double cumd_exponential(double x);
05108 
05109 double cumd_logistic(const double &x);
05110 double inv_cumd_logistic(const double &x);
05111 
05112 dvariable cumd_logistic(const prevariable & x);
05113 dvariable inv_cumd_logistic(const prevariable & x);
05114 double inv_cumd_norm(const double &x);
05115 double cumd_norm(const double &x);
05116 double cumd_norm(const double &x, double);
05117 dvariable inv_cumd_norm(const prevariable & x);
05118 dvar_vector inv_cumd_norm(const dvar_vector & x);
05119 prevariable & cumd_norm(const prevariable & x);
05120 prevariable & bounded_cumd_norm(const prevariable & x, double);
05121 double bounded_cumd_norm(double x, double);
05122 //dvariable& old_cumd_norm(const prevariable& x);
05123 double normal_tail_right(const double &x);
05124 
05125 dvariable inv_cumd_norm_logistic(const prevariable & x, double);
05126 prevariable & cumd_norm_logistic(const prevariable & x, double);
05127 double inv_cumd_norm_logistic(double x, double);
05128 double cumd_norm_logistic(double x, double);
05129 
05134 class prevariable_position
05135 {
05136    double_and_int *v;
05137  public:
05138    prevariable_position(const prevariable & x)
05139    {
05140       v = x.get_v();
05141    }
05142    prevariable_position(double_and_int * p)
05143    {
05144       v = p;
05145    }
05146    double &xval()
05147    {
05148       return ((v->x));
05149    }
05150 };
05151 
05152 double restore_prevariable_derivative(const prevariable_position & pre);
05153 double restore_prevariable_derivative(void);
05154 prevariable_position restore_prevariable_position(void);
05155 void save_double_derivative(double x, const prevariable_position & pos);
05156 double restore_prevariable_value(void);
05157 void save_double_value(double x);
05158 int sum(const imatrix &);
05159 double sum(const dmatrix &);
05160 double sum(const d3_array &);
05161 double sum(const d4_array &);
05162 double sum(const d5_array &);
05163 double sum(const d6_array &);
05164 double sum(const d7_array &);
05165 dvariable sum(const dvar_matrix &);
05166 dvariable sum(const dvar3_array &);
05167 dvariable sum(const dvar4_array &);
05168 dvariable sum(const dvar5_array &);
05169 dvariable sum(const dvar6_array &);
05170 dvariable sum(const dvar7_array &);
05171 
05172 dmatrix fabs(const dmatrix & m);
05173 //double& value(const  double& u);
05174    //const double& value(const double& u);
05175 double norm(const d3_array &);
05176 double norm2(const d3_array &);
05177 double sumsq(const d3_array &);
05178 d3_array exp(const d3_array & m);
05179 d3_array mfexp(const d3_array & m);
05180 d3_array mfexp(const d3_array & m, double d);
05181 d3_array log(const d3_array & m);
05182 d3_array fabs(const d3_array & m);
05183 d3_array sin(const d3_array & m);
05184 d3_array cos(const d3_array & m);
05185 d3_array tan(const d3_array & m);
05186 d3_array sqrt(const d3_array & m);
05187 d3_array sqr(const d3_array & m);
05188 d3_array elem_prod(const d3_array & m1, const d3_array & m2);
05189 d3_array elem_div(const d3_array & m1, const d3_array & m2);
05190 d3_array operator+(const d3_array & m1, const d3_array & m2);
05191 d3_array operator+(const d3_array & m1, double m2);
05192 d3_array operator/(const d3_array & m1, double m2);
05193 d3_array operator/(double m2, const d3_array & m1);
05194 d3_array operator+(double m1, const d3_array & m2);
05195 d3_array operator-(const d3_array & m1, const d3_array & m2);
05196 d3_array operator-(const d3_array & m1, double m2);
05197 d3_array operator-(double m1, const d3_array & m2);
05198 d3_array operator*(const d3_array & m1, const d3_array & m2);
05199 dmatrix operator *(const d3_array & m1, const dvector & m2);
05200 d3_array operator*(const d3_array & m1, double m2);
05201 d3_array operator*(double m1, const d3_array & m2);
05202 
05203 dvariable norm(const dvar3_array & m);
05204 dvariable norm2(const dvar3_array & m);
05205 dvariable sumsq(const dvar3_array & m);
05206 dvar3_array exp(const dvar3_array & m);
05207 dvar3_array mfexp(const dvar3_array & m);
05208 dvar3_array mfexp(const dvar3_array & m, double d);
05209 dvar3_array log(const dvar3_array & m);
05210 dvar3_array fabs(const dvar3_array & m);
05211 dvar3_array sin(const dvar3_array & m);
05212 dvar3_array cos(const dvar3_array & m);
05213 dvar3_array tan(const dvar3_array & m);
05214 dvar3_array sqrt(const dvar3_array & m);
05215 dvar3_array sqr(const dvar3_array & m);
05216 dvar3_array elem_prod(const dvar3_array & m1, const dvar3_array & m2);
05217 dvar3_array elem_div(const dvar3_array & m1, const dvar3_array & m2);
05218 dvar3_array operator+(const dvar3_array & m1, const dvar3_array & m2);
05219 dvar3_array operator-(const dvar3_array & m1, const dvar3_array & m2);
05220 dvar3_array elem_prod(const d3_array & m1, const dvar3_array & m2);
05221 dvar3_array elem_div(const d3_array & m1, const dvar3_array & m2);
05222 dvar3_array operator+(const d3_array & m1, const dvar3_array & m2);
05223 dvar3_array operator-(const d3_array & m1, const dvar3_array & m2);
05224 dvar3_array elem_prod(const dvar3_array & m1, const d3_array & m2);
05225 dvar3_array elem_div(const dvar3_array & m1, const d3_array & m2);
05226 dvar3_array operator+(const dvar3_array & m1, const d3_array & m2);
05227 dvar3_array operator+(const dvar3_array & m1, const dvariable & m2);
05228 dvar3_array operator+(const dvariable & d1, const dvar3_array & m1);
05229 
05230 dvar3_array operator/(const prevariable & m2, const dvar3_array & m1);
05231 dvar3_array operator/(const prevariable & m2, const d3_array & m1);
05232 dvar3_array operator/(double m2, const dvar3_array & m1);
05233 
05234 dvar3_array operator/(const dvar3_array & m1, const prevariable & m2);
05235 dvar3_array operator/(const d3_array & m1, const prevariable & m2);
05236 dvar3_array operator/(const dvar3_array & m1, double m2);
05237 
05238 dvar3_array operator+(const dvariable & m1, const d3_array & m2);
05239 dvar3_array operator+(double m1, const dvar3_array & m2);
05240 dvar3_array operator-(const dvar3_array & m1, const d3_array & m2);
05241 dvar3_array operator-(const dvar3_array & m1, const dvariable & m2);
05242 dvar3_array operator-(const dvariable & m1, const d3_array & m2);
05243 dvar3_array operator-(const dvariable & m1, const dvar3_array & m2);
05244 dvar3_array operator-(double m1, const dvar3_array & m2);
05245 dvar3_array operator*(const dvar3_array & m1, const d3_array & m2);
05246 dvar3_array operator*(const dvar3_array & m1, const dvariable & m2);
05247 dvar3_array operator*(const dvariable & m1, const d3_array & m2);
05248 dvar3_array operator*(const dvariable & m1, const dvar3_array & m2);
05249 dvar3_array operator*(double m1, const dvar3_array & m2);
05250 
05251 double square(double x);
05252 dvector square(const dvector & x);
05253 dmatrix square(const dmatrix & x);
05254 d3_array square(const d3_array & x);
05255 
05256 dvariable & square(const prevariable & x);
05257 dvar_vector square(const dvar_vector & x);
05258 dvar_matrix square(const dvar_matrix & x);
05259 dvar3_array square(const dvar3_array & x);
05260 
05261 double cube(double x);
05262 double