ADMB Documentation  11.4.2891
 All Classes Files Functions Variables Typedefs Friends Defines
xgradclc.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 #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   gradient_structure::GRAD_LIST->initialize();
00125 
00126   double_and_int* tmp =
00127     (double_and_int*)gradient_structure::ARRAY_MEMBLOCK_BASE;
00128 
00129      unsigned long int max_last_offset =
00130        gradient_structure::ARR_LIST1->get_max_last_offset();
00131 
00132   size_t size = sizeof(double_and_int);
00133 
00134   double * zptr;
00135 
00136    for (unsigned int i = 0; i < (max_last_offset/size); i++)
00137    {
00138      tmp->x = 0;
00139 #if defined (__ZTC__)
00140      tmp = (double_and_int*)_farptr_norm((void*)(++tmp));
00141 #else
00142      tmp++;
00143 #endif
00144    }
00145 
00146     *gradient_structure::GRAD_STACK1->ptr->dep_addr = 1;
00147     zptr = gradient_structure::GRAD_STACK1->ptr->dep_addr;
00148 
00149 int break_flag=1;
00150 int funnel_flag=0;
00151 
00152 do
00153 {
00154   gradient_structure::GRAD_STACK1->ptr++;
00155   #ifdef FAST_ASSEMBLER
00156     gradloop();
00157   #else
00158     grad_stack_entry * grad_ptr_first=
00159       gradient_structure::GRAD_STACK1->ptr_first;
00160     while (gradient_structure::GRAD_STACK1->ptr-- >
00161            grad_ptr_first)
00162     {
00163       if (!gradient_structure::GRAD_STACK1->ptr->func)
00164       {
00165         funnel_flag=1;
00166         break;
00167       }
00168       else
00169         (*(gradient_structure::GRAD_STACK1->ptr->func))();
00170     }
00171 
00172   #endif
00173    if (funnel_flag) break;
00174 
00175   // back up the file one buffer size and read forward
00176   off_t offset = (off_t)(sizeof(grad_stack_entry)
00177     * gradient_structure::GRAD_STACK1->length);
00178   lpos = lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR,
00179     -offset, SEEK_CUR);
00180 
00181   break_flag=gradient_structure::GRAD_STACK1->read_grad_stack_buffer(lpos);
00182 }  while (break_flag); // do
00183 
00184  {
00185    if (lpos<0)
00186    {
00187      #ifdef GRAD_DIAG
00188       off_t ttmp =
00189      #endif
00190       lseek(gradient_structure::GRAD_STACK1->_GRADFILE_PTR, 0,SEEK_CUR);
00191 
00192      #ifdef GRAD_DIAG
00193       cout << "Offset in file at end of gradcalc is " << ttmp
00194            << " bytes from the beginning\n";
00195      #endif
00196    }
00197  }
00198 
00199   //if (gradient_structure::save_var_flag)
00200   {
00201     unsigned long bytes_needed = min(
00202       gradient_structure::ARR_LIST1->get_last_offset() + 1,
00203       gradient_structure::ARRAY_MEMBLOCK_SIZE);
00204     size_t _dsize = bytes_needed/sizeof(double);
00205 #ifndef OPT_LIB
00206     assert(_dsize <= INT_MAX);
00207 #endif
00208     int dsize = (int)_dsize;
00209 
00210     //dvector dtmp(0,dsize-1);
00211     //memcpy((char*)&(dtmp(0)),(char*)gradient_structure::ARRAY_MEMBLOCK_BASE,
00212       //dsize*sizeof(double));
00213 
00214     double* dptr=(double*) gradient_structure::ARRAY_MEMBLOCK_BASE;
00215     dptr-=1;
00216     int ii=0;
00217     int nzero=0;
00218     int nnzero=0;
00219     int dcount=0;
00220     int zero_flag=0;
00221     ivector offset(0,dsize-1);
00222     save_identifier_string("ue");
00223     if (!ISZERO(*(++dptr)))
00224     {
00225       save_double_value(*dptr);
00226       dcount++;
00227       zero_flag=0;
00228       offset(ii++)=0;
00229       nnzero++;
00230     }
00231     else
00232     {
00233       zero_flag=1;
00234       nzero++;
00235     }
00236 
00237     for (int i1=1;i1<dsize;i1++)
00238     {
00239       if (!ISZERO(*(++dptr)))
00240       {
00241         save_double_value(*dptr);
00242         dcount++;
00243         nnzero++;
00244         if (zero_flag)
00245         {
00246           offset(ii++)=nzero;
00247           nzero=0;
00248           zero_flag=0;
00249         }
00250       }
00251       else
00252       {
00253         nzero++;
00254         if (!zero_flag)
00255         {
00256           offset(ii++)=nnzero;
00257           nnzero=0;
00258           zero_flag=1;
00259         }
00260       }
00261     }
00262     save_int_value(dcount);
00263 
00264     for (int i=0;i<ii;i++)
00265     {
00266       save_int_value(offset(i));
00267     }
00268     save_int_value(ii);
00269 
00270     unsigned int ssize=gradient_structure::GRAD_LIST->nlinks;
00271 #ifndef OPT_LIB
00272     assert(ssize > 0);
00273     assert(ssize <= INT_MAX);
00274 #endif
00275     dvector stmp(0,(int)(ssize-1));
00276 
00277 #ifndef OPT_LIB
00278     assert(gradient_structure::GRAD_LIST->nlinks <= INT_MAX);
00279 #endif
00280     for (int i=0; i < (int)gradient_structure::GRAD_LIST->nlinks; i++)
00281     {
00282       memcpy((char*)&(stmp(i)),
00283         gradient_structure::GRAD_LIST->dlink_addresses[i],sizeof(double));
00284     }
00285     //dtmp.save_dvector_value();
00286     //dtmp.save_dvector_position();
00287     stmp.save_dvector_value();
00288     stmp.save_dvector_position();
00289 
00290     // save the address of the dependent variable for the funnel
00291     size_t wsize=sizeof(double_and_int*);
00292     gradient_structure::get_fp()->fwrite(&zptr, wsize);
00293     save_identifier_string("ae");
00294 
00295     gradient_structure::GRAD_STACK1->set_gradient_stack(funnel_derivatives);
00296     gradient_structure::restore_arrays();
00297     gradient_structure::restore_variables();
00298   }
00299 }
00300 
00305 void funnel_derivatives(void)
00306 {
00307   verify_identifier_string("ae");
00308   prevariable_position deppos=restore_prevariable_position();
00309   dvector_position stmp_pos=restore_dvector_position();
00310   dvector stmp=restore_dvector_value(stmp_pos);
00311   //dvector_position dtmp_pos=restore_dvector_position();
00312   //dvector dtmp=restore_dvector_value(dtmp_pos);
00313   int ii=restore_int_value();
00314   int i;
00315   int ip=ii;
00316   if (!ip) ip=1;
00317   ivector offset(0,ip);
00318   offset(ip)=0;
00319   //ivector offset(0,ip-1);
00320   for (i=ii-1;i>=0;i--)
00321   {
00322     offset(i)=restore_int_value();
00323   }
00324   int dcount=restore_int_value();
00325   int dc=dcount;
00326   if (!dc) dc=1;
00327   dvector dx(0,dc-1);
00328   for (i=dcount-1;i>=0;i--)
00329   {
00330     dx(i)=restore_double_value();
00331   }
00332 
00333   verify_identifier_string("ue");
00334 
00335   double df = restore_prevariable_derivative(deppos);
00336   double * dptr= (double *) gradient_structure::ARRAY_MEMBLOCK_BASE;
00337 
00338   //double * dd = &(dx(1));
00339   ii=0;
00340   int ic=0;
00341   dptr+=offset(ii++);
00342   for (i=0;i<dcount;i++)
00343   {
00344     *(dptr++)+=dx(i)*df;
00345     if (++ic==offset(ii))
00346     {
00347       if (i==dcount-1)
00348       {
00349         break;
00350       }
00351       dptr+=offset(ii+1);
00352       ii+=2;
00353       ic=0;
00354     }
00355   }
00356 
00357   int smax=stmp.indexmax();
00358   for (i=0;i<smax;i++)
00359   {
00360     if (!ISZERO(stmp(i)))
00361     {
00362       *(double*)(gradient_structure::GRAD_LIST->dlink_addresses[i])
00363         +=stmp(i)*df;
00364     }
00365   }
00366 }
00367 
00372 dvariable& funnel_dvariable::operator=(const prevariable& t)
00373 {
00374   dvariable::operator = (t);
00375   funnel_gradcalc();
00376   return *this;
00377 }
00378 
00383 void ad_begin_funnel(void)
00384 {
00385   gradient_structure::GRAD_STACK1->set_gradient_stack(NULL);
00386 }
00387