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