ADMB Documentation  11.2.2881
 All Classes Files Functions Variables Typedefs Friends Defines
xgradclc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xgradclc.cpp 2641 2014-11-13 20:28:20Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 #include <sys/stat.h>
00014 #include <fcntl.h>
00015 #include <string.h>
00016 
00017 #ifdef __TURBOC__
00018   #pragma hdrstop
00019   #include <iostream.h>
00020 #endif
00021 
00022 #ifdef __ZTC__
00023   #include <iostream.hpp>
00024 #endif
00025 
00026 #if defined (__WAT32__)
00027 #  include <io.h>
00028 #endif
00029 
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 
00033 #ifdef _MSC_VER
00034   #define lseek _lseek
00035   #define  read _read
00036   #define write _write
00037 #else
00038   #include <iostream>
00039   using namespace std;
00040   #include <sys/stat.h>
00041   #include <sys/types.h>
00042   #include <unistd.h>
00043 #endif
00044 
00045 #ifdef __SUN__
00046   #include <iostream.h>
00047   #include <sys/stat.h>
00048   #include <sys/types.h>
00049   #include <unistd.h>
00050 #endif
00051 
00052 #if defined(__NDPX__ )
00053   extern "C" {
00054     int lseek(int, int, int);
00055     int read(int, char*, int);
00056   };
00057 #endif
00058 
00059 #include <math.h>
00060 
00061 #ifndef OPT_LIB
00062   #include <cassert>
00063   #include <climits>
00064 #endif
00065 
00066 #ifdef ISZERO
00067   #undef ISZERO
00068 #endif
00069 #define ISZERO(d) ((d)==0.0)
00070 
00071 void funnel_derivatives(void);
00072 
00073 #if defined(__ZTC__)
00074   void _far * _cdecl _farptr_norm(void _far *);
00075   void _far * _cdecl _farptr_fromlong(unsigned long);
00076   long _cdecl _farptr_tolong(void _far *);
00077 #endif
00078 
00083 void funnel_gradcalc(void)
00084 {
00085 #  ifdef NO_DERIVS
00086      if (gradient_structure::no_derivatives)
00087      {
00088        return;
00089      }
00090 #  endif
00091   gradient_structure::TOTAL_BYTES = 0;
00092   gradient_structure::PREVIOUS_TOTAL_BYTES=0;
00093   if(!gradient_structure::instances)
00094   {
00095     return;
00096   }
00097 
00098    gradient_structure::GRAD_STACK1->_GRADFILE_PTR =
00099               gradient_structure::GRAD_STACK1->gradfile_handle();
00100 
00101   int& _GRADFILE_PTR=gradient_structure::GRAD_STACK1->_GRADFILE_PTR;
00102 
00103   off_t lpos = lseek(_GRADFILE_PTR,0L,SEEK_CUR);
00104 
00105   if(gradient_structure::GRAD_STACK1->ptr
00106        <= gradient_structure::GRAD_STACK1->ptr_first)
00107   {
00108 #ifdef SAFE_ALL
00109       cerr <<
00110         "warning -- calling funnel_gradcalc when no calculations generating"
00111            << endl << "derivative information have occurred" << endl;
00112 #endif
00113     return;
00114   }    // current is one past the end so -- it
00115 
00116   //if (gradient_structure::save_var_flag)
00117   {
00118     gradient_structure::save_arrays();
00119     gradient_structure::save_variables();
00120   }
00121 
00122   gradient_structure::GRAD_STACK1->ptr--;
00123 
00124   for (unsigned int i=0; i<gradient_structure::GRAD_LIST->nlinks; i++)
00125   {
00126     * (double*) (gradient_structure::GRAD_LIST->dlink_addresses[i]) = 0;
00127   }
00128 
00129   double_and_int* tmp =
00130     (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00131 
00132      unsigned long int max_last_offset =
00133        gradient_structure::ARR_LIST1->get_max_last_offset();
00134 
00135      unsigned int size = sizeof(double_and_int );
00136 
00137   double * zptr;
00138 
00139    for (unsigned int i = 0; i < (max_last_offset/size); i++)
00140    {
00141      tmp->x = 0;
00142 #if defined (__ZTC__)
00143      tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00144 #else
00145      tmp++;
00146 #endif
00147    }
00148 
00149     *gradient_structure::GRAD_STACK1->ptr->dep_addr = 1;
00150     zptr = gradient_structure::GRAD_STACK1->ptr->dep_addr;
00151 
00152 int break_flag=1;
00153 int funnel_flag=0;
00154 
00155 do
00156 {
00157   gradient_structure::GRAD_STACK1->ptr++;
00158   #ifdef FAST_ASSEMBLER
00159     gradloop();
00160   #else
00161     grad_stack_entry * grad_ptr_first=
00162       gradient_structure::GRAD_STACK1->ptr_first;
00163     while (gradient_structure::GRAD_STACK1->ptr-- >
00164            grad_ptr_first)
00165     {
00166       if (!gradient_structure::GRAD_STACK1->ptr->func)
00167       {
00168         funnel_flag=1;
00169         break;
00170       }
00171       else
00172         (*(gradient_structure::GRAD_STACK1->ptr->func))();
00173     }
00174 
00175   #endif
00176    if (funnel_flag) break;
00177 
00178   // back up the file one buffer size and read forward
00179   off_t offset = (off_t)(sizeof(grad_stack_entry)
00180     * gradient_structure::GRAD_STACK1->length);
00181   lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00182     -offset, SEEK_CUR);
00183 
00184   break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00185 }  while (break_flag); // do
00186 
00187  {
00188    if (lpos<0)
00189    {
00190      #ifdef GRAD_DIAG
00191       off_t ttmp =
00192      #endif
00193       lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00194 
00195      #ifdef GRAD_DIAG
00196       cout << "Offset in file at end of gradcalc is " << ttmp
00197            << " bytes from the beginning\n";
00198      #endif
00199    }
00200  }
00201 
00202   //if (gradient_structure::save_var_flag)
00203   {
00204     unsigned long bytes_needed = min(
00205       gradient_structure::ARR_LIST1->get_last_offset() + 1,
00206       gradient_structure::ARRAY_MEMBLOCK_SIZE);
00207     size_t _dsize = bytes_needed/sizeof(double);
00208 #ifndef OPT_LIB
00209     assert(_dsize <= INT_MAX);
00210 #endif
00211     int dsize = (int)_dsize;
00212 
00213     //dvector dtmp(0,dsize-1);
00214     //memcpy((char*)&(dtmp(0)),(char*)gradient_structure::ARRAY_MEMBLOCK_BASE,
00215       //dsize*sizeof(double));
00216 
00217     double* dptr=(double*) gradient_structure::ARRAY_MEMBLOCK_BASE;
00218     dptr-=1;
00219     int ii=0;
00220     int nzero=0;
00221     int nnzero=0;
00222     int dcount=0;
00223     int zero_flag=0;
00224     ivector offset(0,dsize-1);
00225     save_identifier_string("ue");
00226     if (!ISZERO(*(++dptr)))
00227     {
00228       save_double_value(*dptr);
00229       dcount++;
00230       zero_flag=0;
00231       offset(ii++)=0;
00232       nnzero++;
00233     }
00234     else
00235     {
00236       zero_flag=1;
00237       nzero++;
00238     }
00239 
00240     for (int i1=1;i1<dsize;i1++)
00241     {
00242       if (!ISZERO(*(++dptr)))
00243       {
00244         save_double_value(*dptr);
00245         dcount++;
00246         nnzero++;
00247         if (zero_flag)
00248         {
00249           offset(ii++)=nzero;
00250           nzero=0;
00251           zero_flag=0;
00252         }
00253       }
00254       else
00255       {
00256         nzero++;
00257         if (!zero_flag)
00258         {
00259           offset(ii++)=nnzero;
00260           nnzero=0;
00261           zero_flag=1;
00262         }
00263       }
00264     }
00265     save_int_value(dcount);
00266 
00267     for (int i=0;i<ii;i++)
00268     {
00269       save_int_value(offset(i));
00270     }
00271     save_int_value(ii);
00272 
00273     unsigned int ssize=gradient_structure::GRAD_LIST->nlinks;
00274 #ifndef OPT_LIB
00275     assert(ssize > 0);
00276     assert(ssize <= INT_MAX);
00277 #endif
00278     dvector stmp(0,(int)(ssize-1));
00279 
00280 #ifndef OPT_LIB
00281     assert(gradient_structure::GRAD_LIST->nlinks <= INT_MAX);
00282 #endif
00283     for (int i=0; i < (int)gradient_structure::GRAD_LIST->nlinks; i++)
00284     {
00285       memcpy((char*)&(stmp(i)),
00286         gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00287     }
00288     //dtmp.save_dvector_value();
00289     //dtmp.save_dvector_position();
00290     stmp.save_dvector_value();
00291     stmp.save_dvector_position();
00292 
00293     // save the address of the dependent variable for the funnel
00294     int wsize=sizeof(double_and_int*);
00295     gradient_structure::get_fp()->fwrite(&zptr,size_t(wsize));
00296     save_identifier_string("ae");
00297 
00298     gradient_structure::GRAD_STACK1->set_gradient_stack(funnel_derivatives);
00299     gradient_structure::restore_arrays();
00300     gradient_structure::restore_variables();
00301   }
00302 }
00303 
00308 void funnel_derivatives(void)
00309 {
00310   verify_identifier_string("ae");
00311   prevariable_position deppos=restore_prevariable_position();
00312   dvector_position stmp_pos=restore_dvector_position();
00313   dvector stmp=restore_dvector_value(stmp_pos);
00314   //dvector_position dtmp_pos=restore_dvector_position();
00315   //dvector dtmp=restore_dvector_value(dtmp_pos);
00316   int ii=restore_int_value();
00317   int i;
00318   int ip=ii;
00319   if (!ip) ip=1;
00320   ivector offset(0,ip);
00321   offset(ip)=0;
00322   //ivector offset(0,ip-1);
00323   for (i=ii-1;i>=0;i--)
00324   {
00325     offset(i)=restore_int_value();
00326   }
00327   int dcount=restore_int_value();
00328   int dc=dcount;
00329   if (!dc) dc=1;
00330   dvector dx(0,dc-1);
00331   for (i=dcount-1;i>=0;i--)
00332   {
00333     dx(i)=restore_double_value();
00334   }
00335 
00336   verify_identifier_string("ue");
00337 
00338   double df = restore_prevariable_derivative(deppos);
00339   double * dptr= (double *) gradient_structure::ARRAY_MEMBLOCK_BASE;
00340 
00341   //double * dd = &(dx(1));
00342   ii=0;
00343   int ic=0;
00344   dptr+=offset(ii++);
00345   for (i=0;i<dcount;i++)
00346   {
00347     *(dptr++)+=dx(i)*df;
00348     if (++ic==offset(ii))
00349     {
00350       if (i==dcount-1)
00351       {
00352         break;
00353       }
00354       dptr+=offset(ii+1);
00355       ii+=2;
00356       ic=0;
00357     }
00358   }
00359 
00360   int smax=stmp.indexmax();
00361   for (i=0;i<smax;i++)
00362   {
00363     if (!ISZERO(stmp(i)))
00364     {
00365       *(double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00366         +=stmp(i)*df;
00367     }
00368   }
00369 }
00370 
00375 dvariable& funnel_dvariable::operator=(const prevariable& t)
00376 {
00377   dvariable::operator = (t);
00378   funnel_gradcalc();
00379   return *this;
00380 }
00381 
00386 void ad_begin_funnel(void)
00387 {
00388   gradient_structure::GRAD_STACK1->set_gradient_stack(NULL);
00389 }
00390