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