ADMB Documentation  11.5.3254
 All Classes Files Functions Variables Typedefs Friends Defines
derch.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #ifdef __ZTC__
00012   #include <conio.h>
00013 #endif
00014 
00015 #include "fvar.hpp"
00016 
00017 #if defined(_WIN32)
00018   #include <conio.h>
00019 #endif
00020 
00021 #include <ctype.h>
00022 #include <stdlib.h>
00023 
00024 #ifndef OPT_LIB
00025   #include <cassert>
00026   #include <climits>
00027 #endif
00028 
00029 int derchflag=0;
00030 double derch_stepsize=0.0;
00031 
00032 static ofstream* pofs=0;
00033 
00038 void derch(const double& _f, const independent_variables & _x,
00039  const dvector& _gg, int n, const int & _ireturn)
00040 {
00041   dvector& gg=(dvector&) _gg;
00042   double& f=(double&) _f;
00043   int& ireturn = (int&) _ireturn;
00044   independent_variables& x= (independent_variables&) _x;
00045   static int i, n1 ,n2,ii;
00046   static double fsave;
00047   static int order_flag;
00048   static double s, f1, f2, g2, xsave;
00049   static int j = 1;
00050   static int si;
00051   si=gg.indexmax();
00052   static dvector g(1,si);
00053   static dvector index(1,si);
00054 
00055   if (ireturn == 4 ) goto label4;
00056   else if (ireturn == 5) goto label5;
00057   g=gg;
00058   {
00059     int maxind = g.indexmin();
00060     double maxg=fabs(g(maxind));
00061     dmatrix tmp(1,3,g.indexmin(),g.indexmax());
00062     tmp(1).fill_seqadd(1,1);
00063     tmp(2)=g;
00064     tmp(3)=-fabs(g);
00065     ofstream ofs("derivatives");
00066     dmatrix dtmp=sort(trans(tmp),3);
00067     ofs << dtmp << endl;
00068     for (int i=g.indexmin();i<=g.indexmax();i++)
00069     {
00070       if (fabs(g(i))>maxg)
00071       {
00072         maxg=fabs(g(i));
00073         maxind=i;
00074       }
00075     }
00076     cout << "maxind = " << maxind << " maxg = " << maxg << endl;
00077     index=column(dtmp,1);
00078   }
00079   while (j > 0)
00080   {
00081     cout << "\nEntering derivative checker.\n";
00082     cout << "   Enter index (1 ... "<< n <<") of derivative to check.\n";
00083     cout << "     To check all derivatives, enter 0;\n";
00084     cout << "     To quit  enter -1: ";
00085     flush(cout);
00086     cin >> j;
00087     if (j == -1)
00088     {
00089       ireturn = -1;
00090       return;
00091     }
00092     else if (j == 0)
00093     {
00094       cout << "\n   Checking all derivatives. Press X to terminate checking.\n";
00095       flush(cout);
00096       pofs=new ofstream("ders.dat");
00097       (*pofs) << "      index  analytical        finite     % error " << endl;
00098       (*pofs) << "                              difference " << endl;
00099       n1 = 1;
00100       n2 = n;
00101     }
00102     else if (j >= 1 && j <= n)
00103     {
00104       n1 = j;
00105       n2 = j;
00106     }
00107     else
00108     {
00109       cout << "Error: Invalid derivative index \"" << j
00110            << "\" entered which is not in range (1 ... "<< n <<").\n";
00111       ad_exit(1);
00112     }
00113     cout << "\n   Enter step size.\n";
00114     cout << "      To quit derivative checker enter -1;\n";
00115     cout << "      to check for uninitialized variables enter 0): ";
00116     flush(cout);
00117 #if defined(__BORLANDC__)
00118     char mystring[1000];
00119     cin >> mystring;
00120     s=atof(mystring);
00121 #else
00122     cin >> s;
00123 #endif
00124     if (s < 0) ad_exit(0);
00125 
00126     if (j==0)
00127     {
00128       cout << "\n   If you want the derivatives approximated in order\n"
00129            << "   of decreasing magnitude enter 1;\n"
00130            << "   else enter 0: ";
00131       flush(cout);
00132       int tmpint=0;
00133       cin >> tmpint;
00134       if (tmpint==1)
00135         order_flag=1;
00136       else
00137         order_flag=0;
00138     }
00139     if (s != 0)
00140     {
00141       cout <<
00142         "\n   X             Function      Analytical    Finite Diff ; Index"
00143         << endl;
00144     }
00145 
00146     for (ii=n1; ii<=n2; ii++)
00147     {
00148       i = ii;
00149       if (order_flag==1)
00150       {
00151         int count = 0;
00152         for (int _ii = index.indexmin(); _ii <= index.indexmax(); ++_ii)
00153         {
00154 #ifdef OPT_LIB
00155           int idx = (int)index(_ii);
00156 #else
00157           double _idx = index(_ii);
00158           assert(_idx <= (double)INT_MAX);
00159           int idx = (int)_idx;
00160 #endif
00161           if (x.indexmin() <= idx && idx <= x.indexmax())
00162           {
00163             ++count;
00164           }
00165           if (count == ii)
00166           {
00167             i = idx;
00168             break;
00169           }
00170         }
00171       }
00172       derch_stepsize=s;
00173       derchflag=1;
00174       xsave=x(i);
00175       x(i)=xsave+s;
00176       fsave = f;
00177       ireturn = 4; // fgcomp(&f1,x,g1,n, params, vars);
00178       return;
00179 
00180     label4:
00181       derch_stepsize=s;
00182       derchflag=-1;
00183       f1 = f;
00184       x(i)=xsave-s;
00185       ireturn= 5; //fgcomp(&f2,x,g1,n, params, vars);
00186       return;
00187 
00188     label5:
00189       f2 = f;
00190       f = fsave;
00191       x(i)=xsave;
00192       if (s==0.0)
00193       {
00194         if (fabs(f1-f2)>0.0)
00195         {
00196           cout << "There appear to be uninitialized variables in "
00197                << "the objective function "  << endl
00198                << "    f1 = " << setprecision(15) << f1
00199                << " f2 = " << setprecision(15) << f2 << endl;
00200         }
00201         else
00202         {
00203           cout << "There do not appear to be uninitialized variables in" << endl
00204                << "the objective function " << endl;
00205         }
00206       }
00207       else
00208       {
00209         g2=(f1-f2)/(2.*s);
00210         derchflag=0;
00211         double perr= fabs(g2-g(i))/(fabs(g(i))+1.e-6);
00212 
00213         if (pofs)
00214         {
00215           if (perr > 1.e-3)
00216             (*pofs) << " ** ";
00217           else
00218             (*pofs) << "    ";
00219           (*pofs) << "  " << setw(4) << i
00220                   << "  " <<  setw(12) << g(i)
00221                   << "  " <<  setw(12) << g2
00222                   << "  " <<  setw(12) <<  perr
00223                   << endl;
00224         }
00225         if (ad_printf)
00226         {
00227           (*ad_printf)("  %12.5e  %12.5e  %12.5e  %12.5e ; %5d \n",
00228                 x(i), f, g(i), g2, i);
00229           fflush(stdout);
00230         }
00231       }
00232 #if defined(_MSC_VER)
00233       if ( kbhit() )
00234       {
00235         int c = toupper(getch());
00236         if ( c == 'X')
00237         {
00238           cout << "   Exiting derivative checker by user interrupt.\n";
00239           ad_exit(0);
00240         }
00241       }
00242 #endif
00243     } // for loop
00244   } // while (j > 0)
00245 //  ireturn = 2;
00246   ad_exit(0);
00247 }