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