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