ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
derch.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: derch.cpp 2473 2014-10-04 00:16:49Z johnoel $
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 long int i, n1 ,n2,ii;
00046   static double fsave;
00047   static int order_flag;
00048   static double s, f1, f2, g2, xsave;
00049   static long 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     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
00103     {
00104       n1 = j;
00105       n2 = j;
00106     }
00107     cout << "\n   Enter step size.\n";
00108     cout << "      To quit derivative checker enter -1;\n";
00109     cout << "      to check for uninitialized variables enter 0): ";
00110     flush(cout);
00111 #if defined(__BORLANDC__)
00112     char mystring[1000];
00113     cin >> mystring;
00114     s=atof(mystring);
00115 #else
00116     cin >> s;
00117 #endif
00118     if (s < 0) ad_exit(0);
00119 
00120     if (j==0)
00121     {
00122       cout << "\n   If you want the derivatives approximated in order\n"
00123            << "   of decreasing magnitude enter 1;\n"
00124            << "   else enter 0: ";
00125       flush(cout);
00126       int tmpint=0;
00127       cin >> tmpint;
00128       if (tmpint==1)
00129         order_flag=1;
00130       else
00131         order_flag=0;
00132     }
00133     if (s != 0)
00134     {
00135       cout <<
00136         "\n   X             Function      Analytical    Finite Diff ; Index"
00137         << endl;
00138     }
00139 
00140     for (ii=n1; ii<=n2; ii++)
00141     {
00142       i = ii;
00143       if (order_flag==1)
00144       {
00145         int count = 0;
00146         for (int _ii = index.indexmin(); _ii <= index.indexmax(); ++_ii)
00147         {
00148 #ifdef OPT_LIB
00149           int idx = (int)index(_ii);
00150 #else
00151           double _idx = index(_ii);
00152           assert(_idx <= (double)INT_MAX);
00153           int idx = (int)_idx;
00154 #endif
00155           if (x.indexmin() <= idx && idx <= x.indexmax())
00156           {
00157             ++count;
00158           }
00159           if (count == ii)
00160           {
00161             i = idx;
00162             break;
00163           }
00164         }
00165       }
00166       derch_stepsize=s;
00167       derchflag=1;
00168       xsave=x(i);
00169       x(i)=xsave+s;
00170       fsave = f;
00171       ireturn = 4; // fgcomp(&f1,x,g1,n, params, vars);
00172       return;
00173 
00174     label4:
00175       derch_stepsize=s;
00176       derchflag=-1;
00177       f1 = f;
00178       x(i)=xsave-s;
00179       ireturn= 5; //fgcomp(&f2,x,g1,n, params, vars);
00180       return;
00181 
00182     label5:
00183       f2 = f;
00184       f = fsave;
00185       x(i)=xsave;
00186       if (s==0.0)
00187       {
00188         if (fabs(f1-f2)>0.0)
00189         {
00190           cout << "There appear to be uninitialized variables in "
00191                << "the objective function "  << endl
00192                << "    f1 = " << setprecision(15) << f1
00193                << " f2 = " << setprecision(15) << f2 << endl;
00194         }
00195         else
00196         {
00197           cout << "There do not appear to be uninitialized variables in" << endl
00198                << "the objective function " << endl;
00199         }
00200       }
00201       else
00202       {
00203         g2=(f1-f2)/(2.*s);
00204         derchflag=0;
00205         double perr= fabs(g2-g(i))/(fabs(g(i))+1.e-6);
00206 
00207         if (pofs)
00208         {
00209           if (perr > 1.e-3)
00210             (*pofs) << " ** ";
00211           else
00212             (*pofs) << "    ";
00213           (*pofs) << "  " << setw(4) << i
00214                   << "  " <<  setw(12) << g(i)
00215                   << "  " <<  setw(12) << g2
00216                   << "  " <<  setw(12) <<  perr
00217                   << endl;
00218         }
00219         if (ad_printf)
00220         {
00221           (*ad_printf)("  %12.5e  %12.5e  %12.5e  %12.5e ; %5d \n",
00222                 x(i), f, g(i), g2, i);
00223           fflush(stdout);
00224         }
00225       }
00226 #if defined(_MSC_VER)
00227       if ( kbhit() )
00228       {
00229         int c = toupper(getch());
00230         if ( c == 'X')
00231         {
00232           cout << "   Exiting derivative checker by user interrupt.\n";
00233           ad_exit(0);
00234         }
00235       }
00236 #endif
00237     } // for loop
00238   } // while (j > 0)
00239 //  ireturn = 2;
00240   ad_exit(0);
00241 }