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