ADMB Documentation  11.1.2370
 All Classes Files Functions Variables Typedefs Friends Defines
derch.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: derch.cpp 1906 2014-04-17 20:52: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 int derchflag=0;
00025 double derch_stepsize=0.0;
00026 
00027 static ofstream* pofs=0;
00028 
00033 void derch(const double& _f, const independent_variables & _x,
00034  const dvector& _gg, int n, const int & _ireturn)
00035 {
00036   dvector& gg=(dvector&) _gg;
00037   double& f=(double&) _f;
00038   int& ireturn = (int&) _ireturn;
00039   independent_variables& x= (independent_variables&) _x;
00040   static long int i, n1 ,n2,ii;
00041   static double fsave;
00042   static int order_flag;
00043   static double s, f1, f2, g2, xsave;
00044   static long int j = 1;
00045   static int si;
00046   si=gg.indexmax();
00047   static dvector g(1,si);
00048   static dvector index(1,si);
00049 
00050   if (ireturn == 4 ) goto label4;
00051   else if (ireturn == 5) goto label5;
00052   g=gg;
00053   {
00054     int maxind=0;
00055     double  maxg;
00056     maxind=g.indexmin();
00057     maxg=fabs(g(maxind));
00058     dmatrix tmp(1,3,g.indexmin(),g.indexmax());
00059     tmp(1).fill_seqadd(1,1);
00060     tmp(2)=g;
00061     tmp(3)=-fabs(g);
00062     ofstream ofs("derivatives");
00063     dmatrix dtmp=sort(trans(tmp),3);
00064     ofs << dtmp << endl;
00065     for (int i=g.indexmin();i<=g.indexmax();i++)
00066     {
00067       if (fabs(g(i))>maxg)
00068       {
00069         maxg=fabs(g(i));
00070         maxind=i;
00071       }
00072     }
00073     cout << "maxind = " << maxind << " maxg = " << maxg << endl;
00074     index=ivector(column(dtmp,1));
00075   }
00076   while (j > 0)
00077   {
00078     cout << "\nEntering derivative checker.\n";
00079     cout << "   Enter index (1 ... "<< n <<") of derivative to check.\n";
00080     cout << "     To check all derivatives, enter 0;\n";
00081     cout << "     To quit  enter -1: ";
00082     flush(cout);
00083     cin >> j;
00084     if (j == -1)
00085     {
00086       ireturn = -1;
00087       return;
00088     }
00089     if (j == 0)
00090     {
00091       cout << "\n   Checking all derivatives. Press X to terminate checking.\n";
00092       flush(cout);
00093       pofs=new ofstream("ders.dat");
00094       (*pofs) << "      index  analytical        finite     % error " << endl;
00095       (*pofs) << "                              difference " << endl;
00096       n1 = 1;
00097       n2 = n;
00098     }
00099     else
00100     {
00101       n1 = j;
00102       n2 = j;
00103     }
00104     cout << "   Enter step size.\n";
00105     cout << "      To quit derivative checker enter -1;\n";
00106     cout << "      to check for unitialized variables enter 0): ";
00107     flush(cout);
00108 #   if defined(__BORLANDC__)
00109       char mystring[1000];
00110       cin >> mystring;
00111       s=atof(mystring);
00112 #   else
00113       cin >> s;
00114 #   endif
00115 
00116     if (s < 0) ad_exit(0);
00117     if (j==0)
00118     {
00119       cout << endl << "   If you want the derivatives approximated in order"
00120         << endl << "   of decreasing magnitude enter 1"
00121         << endl << "   Else enter 0" << endl;
00122       int tmpint=0;
00123       cin >> tmpint;
00124       if (tmpint==1)
00125         order_flag=1;
00126       else
00127         order_flag=0;
00128     }
00129     else
00130     {
00131 cout << "\n       X           Function     Analytical     Finite Diff;  Index"
00132          << endl;
00133     }
00134 
00135     for (ii=n1; ii<=n2; ii++)
00136     {
00137       if (order_flag==1)
00138         i=index(ii);
00139       else
00140         i=ii;
00141 
00142       derch_stepsize=s;
00143       derchflag=1;
00144       xsave=x(i);
00145       x(i)=xsave+s;
00146       fsave = f;
00147       ireturn = 4; // fgcomp(&f1,x,g1,n, params, vars);
00148       return;
00149 
00150     label4:
00151       derch_stepsize=s;
00152       derchflag=-1;
00153       f1 = f;
00154       x(i)=xsave-s;
00155       ireturn= 5; //fgcomp(&f2,x,g1,n, params, vars);
00156       return;
00157 
00158     label5:
00159       f2 = f;
00160       f = fsave;
00161       x(i)=xsave;
00162       if (s==0.0)
00163       {
00164         if (fabs(f1-f2)>0.0)
00165         {
00166           cout << "There appear to be uninitialized variables in "
00167                << "the objective function "  << endl
00168                << "    f1 = " << setprecision(15) << f1
00169                << " f2 = " << setprecision(15) << f2 << endl;
00170         }
00171         else
00172         {
00173           cout << "There do not appear to be uninitialized variables in" << endl
00174                << "the objective function " << endl;
00175         }
00176       }
00177       else
00178       {
00179         g2=(f1-f2)/(2.*s);
00180         derchflag=0;
00181         double perr= fabs(g2-g(i))/(fabs(g(i))+1.e-6);
00182 
00183         if (pofs)
00184         {
00185           if (perr > 1.e-3)
00186             (*pofs) << " ** ";
00187           else
00188             (*pofs) << "    ";
00189           (*pofs) << "  " << setw(4) << i
00190                   << "  " <<  setw(12) << g(i)
00191                   << "  " <<  setw(12) << g2
00192                   << "  " <<  setw(12) <<  perr
00193                   << endl;
00194         }
00195         if (ad_printf)
00196         {
00197           (*ad_printf)("  %12.5e  %12.5e  %12.5e  %12.5e ; %5d \n",
00198                 x(i), f, g(i), g2, i);
00199           fflush(stdout);
00200         }
00201       }
00202 #if defined(_MSC_VER)
00203       if ( kbhit() )
00204       {
00205         int c = toupper(getch());
00206         if ( c == 'X')
00207         {
00208           cout << "   Exiting derivative checker by user interrupt.\n";
00209           ad_exit(0);
00210         }
00211       }
00212 #endif
00213     } // for loop
00214   } // while (j > 0)
00215 //  ireturn = 2;
00216   ad_exit(0);
00217 }