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