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