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