ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
boundfun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: boundfun.cpp 1712 2014-02-28 23:52:33Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00014 #include "fvar.hpp"
00015 
00021 #define USE_BARD_PEN
00022 
00023 #include <stdlib.h>
00024 #include <stdio.h>
00025 #include <math.h>
00026 
00027 // function prototypes duplicated in fvar.hpp
00028 // double dmin(double, double);
00029 // double dmax(double, double);
00030 
00040 dvariable dfatan1(dvariable x, double fmin, double fmax,
00041   const prevariable&  _fpen)
00042 {
00043   prevariable&  fpen=(prevariable&)  _fpen;
00044   dvariable t;
00045 
00046   t= (atan(x)/PI);
00047   t=( t +.5 );
00048   t= t *( fmax-fmin ) + fmin;
00049   t=( (atan(x)/PI) +.5 )*( fmax-fmin ) + fmin;
00050 
00051   if (x < -12.)
00052   {
00053     fpen+=.1*(x+12.)*(x+12.);
00054   }
00055 
00056   if (x > 12.)
00057   {
00058     fpen+=.1*(x-12.)*(x-12.);
00059   }
00060   return(t);
00061 }
00069 double dftinv(double x, double fmin, double fmax)
00070 {
00071   if (x <= fmin)
00072   {
00073     if (ad_printf)
00074     {
00075       (*ad_printf)("variable out of bounds in dftinv\nvariable = %lg", x);
00076       (*ad_printf)("lower bound = %lg", fmin);
00077       (*ad_printf)("upper bound = %lg\n", fmax);
00078     }
00079 
00080     x=dmin(fmin+.001,fmin+.01*(fmax-fmin));
00081   }
00082 
00083   double tinv = tan( ((x-fmin)/(fmax-fmin) -.5) * PI);
00084   return tinv;
00085 }
00095 dvariable boundp(const prevariable& x, double fmin, double fmax,
00096   const prevariable& _fpen,double s)
00097 {
00098   return boundp(x/s,fmin,fmax,_fpen);
00099 }
00108 dvariable boundp(const prevariable& x, double fmin, double fmax,
00109   const prevariable& _fpen)
00110 {
00111   if (gradient_structure::Hybrid_bounded_flag==0)
00112   {
00113     prevariable&  fpen=(prevariable&)  _fpen;
00114     dvariable t,y;
00115     double diff=fmax-fmin;
00116     const double l4=log(4.0);
00117     dvariable ss=0.4999999999999999*sin(x*1.57079632679489661)+0.50;
00118     t=fmin + diff*ss;
00119 
00120   #ifdef USE_BARD_PEN
00121     double pen=.000001/diff;
00122     fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00123   #else
00124     if (x < -.9999)
00125     {
00126       fpen+=cube(-0.9999-x);
00127       if (x < -1.)
00128       {
00129         fpen+=1.e+6*cube(-1.0-x);
00130         if (x < -1.02)
00131         {
00132           fpen+=1.e+10*cube(-1.02-x);
00133         }
00134       }
00135     }
00136     if (x > 0.9999)
00137     {
00138       fpen+=cube(x-0.9999);
00139       if (x > 1.)
00140       {
00141         fpen+=1.e+6*cube(x-1.);
00142         if (x > 1.02)
00143         {
00144           fpen+=1.e+10*cube(x-1.02);
00145         }
00146       }
00147     }
00148   #endif
00149     return(t);
00150   }
00151   else
00152   {
00153     double diff=fmax-fmin;
00154     dvariable t,y;
00155     if (x>-20)
00156     {
00157       y=1.0/(1+exp(-x));
00158     }
00159     else
00160     {
00161       dvariable u=exp(x);
00162       y=u/(1.0+u);
00163     }
00164     t=fmin + diff*y;
00165     return(t);
00166   }
00167 }
00176 dvariable dfboundp(const prevariable& x, double fmin,double fmax)
00177 {
00178   if (gradient_structure::Hybrid_bounded_flag==0)
00179   {
00180     return (fmax-fmin)*0.499999999999999*1.57079632679489661
00181       *cos(x*1.57079632679489661);
00182   }
00183   else
00184   {
00185     double diff=fmax-fmin;
00186     dvariable dfy;
00187     if (x>-20)
00188     {
00189       dvariable u=exp(-x);
00190       //y=1.0/(1+u);
00191       dfy=u/square(1.0+u);
00192     }
00193     else
00194     {
00195       dvariable u=exp(x);
00196       //y=u/(1.0+u);
00197       dfy=u/square(1.0+u);
00198     }
00199     if (dfy==0)
00200     {
00201       cout << "error in dfboundp" << endl;
00202     }
00203     return diff*dfy;
00204   }
00205 }
00206 
00216 double ndfboundp( double x, double fmin, double fmax,const double& fpen)
00217 {
00218   if (gradient_structure::Hybrid_bounded_flag==0)
00219   {
00220     return (fmax-fmin)*0.499999999999999*1.57079632679489661
00221       *cos(x*1.57079632679489661);
00222   }
00223   else
00224   {
00225     double diff=fmax-fmin;
00226     double dfy;
00227     if (x>-20)
00228     {
00229       double u=exp(-x);
00230       //y=1.0/(1+u);
00231       dfy=u/square(1.0+u);
00232     }
00233     else
00234     {
00235       double u=exp(x);
00236       //y=u/(1.0+u);
00237       dfy=u/square(1.0+u);
00238     }
00239     return diff*dfy;
00240   }
00241 }
00242 
00250 double boundp(double x, double fmin, double fmax)
00251 {
00252   if (gradient_structure::Hybrid_bounded_flag==0)
00253   {
00254     double t;
00255     double diff=fmax-fmin;
00256     double ss=0.49999999999999999*sin(x*1.57079632679489661)+0.50;
00257     t=fmin + diff*ss;
00258     return(t);
00259   }
00260   else
00261   {
00262     double diff=fmax-fmin;
00263     double t,y;
00264     if (x>-20)
00265     {
00266       y=1.0/(1+exp(-x));
00267     }
00268     else
00269     {
00270       double u=exp(x);
00271       y=u/(1.0+u);
00272     }
00273     t=fmin + diff*y;
00274     return(t);
00275   }
00276 }
00277 
00288 double nd2fboundp( double x, double fmin, double fmax,const double& fpen)
00289 {
00290   if (x<-0.99999)
00291   {
00292     return (boundp(x,fmin,fmax,fpen)-2.*boundp(x+1.e-6,fmin,fmax,fpen)
00293       +boundp(x+2.e-6,fmin,fmax,fpen))/1.e-12;
00294   }
00295   else if (x>0.99999)
00296   {
00297     return (boundp(x-2.e-6,fmin,fmax,fpen)-2.*boundp(x-1.e-6,fmin,fmax,fpen)
00298       +boundp(x,fmin,fmax,fpen))/1.e-12;
00299   }
00300   else
00301   {
00302     return (boundp(x+1.e-6,fmin,fmax,fpen)-2.*boundp(x,fmin,fmax,fpen)
00303       +boundp(x-1.e-6,fmin,fmax,fpen))/1.e-12;
00304   }
00305 }
00316 double boundp( double x, double fmin, double fmax,const double& _fpen)
00317 {
00318   if (gradient_structure::Hybrid_bounded_flag==0)
00319   {
00320     double t;
00321     double& fpen=(double&) _fpen;
00322     double diff=fmax-fmin;
00323     const double l4=log(4.0);
00324     double ss=0.499999999999999*sin(x*1.57079632679489661)+0.50;
00325     t=fmin + diff*ss;
00326   #ifdef USE_BARD_PEN
00327     double pen=.001/diff;
00328     fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00329   #else
00330     if (x < -.9999)
00331     {
00332       fpen+=(x+0.9999)*(x+0.9999);
00333       if (x < -1.)
00334       {
00335         fpen+=1.e+6*(x+1.)*(x+1.);
00336         if (x < -1.02)
00337         {
00338           fpen+=1.e+10*(x+1.02)*(x+1.02);
00339         }
00340       }
00341     }
00342     if (x > 0.9999)
00343     {
00344       fpen+=(x-0.9999)*(x-0.9999);
00345       if (x > 1.)
00346       {
00347         fpen+=1.e+6*(x-1.)*(x-1.);
00348         if (x > 1.02)
00349         {
00350           fpen+=1.e+10*(x-1.02)*(x-1.02);
00351         }
00352       }
00353     }
00354   #endif
00355     return(t);
00356   }
00357   else
00358   {
00359     double diff=fmax-fmin;
00360     double t,y;
00361     if (x>-20)
00362     {
00363       y=1.0/(1+exp(-x));
00364     }
00365     else
00366     {
00367       double u=exp(x);
00368       y=u/(1.0+u);
00369     }
00370     t=fmin + diff*y;
00371     return(t);
00372   }
00373 }
00374 
00384 double boundpin(double x, double fmin, double fmax,double s)
00385 {
00386   return s*boundpin(x,fmin,fmax);
00387 }
00388 
00397 double boundpin(double x, double fmin, double fmax)
00398 {
00399   if (x < fmin)
00400   {
00401     if (ad_printf)
00402     {
00403       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00404       (*ad_printf)("; min = %lg", fmin);
00405       (*ad_printf)("; max = %lg\n", fmax);
00406     }
00407     x=dmin(fmin+.001,fmin+.01*(fmax-fmin));
00408   }
00409 
00410   if (x > fmax)
00411   {
00412     if (ad_printf)
00413       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00414     if (ad_printf) (*ad_printf)("; min = %lg", fmin);
00415     if (ad_printf) (*ad_printf)("; max = %lg\n", fmax);
00416 
00417     x=dmax(fmax-.001,fmax-.01*(fmax-fmin));
00418   }
00419 
00420   double tinv;
00421   if (gradient_structure::Hybrid_bounded_flag==0)
00422   {
00423     tinv=::asin(2.*(x-fmin)/(fmax-fmin)-1.)/1.57079632679489661;
00424   }
00425   else
00426   {
00427     //double y=(x-fmin)/(fmax-fmin);
00428     //double u=1.e-20+y/(1.e-20+(1.0-y));
00429     double y=1.e-20+(fmax-x)/(1.e-20+(x-fmin));
00430     tinv=-log(y);
00431   }
00432   return(tinv);
00433 }
00444 double boundpin(const prevariable& x, double fmin, double fmax,double s)
00445 {
00446   return s*boundpin(x,fmin,fmax);
00447 }
00448 
00457 double boundpin(const prevariable& xx, double fmin, double fmax)
00458 {
00459   double tinv;
00460   double x=value(xx);
00461 
00462   if (x < fmin)
00463   {
00464     if (ad_printf)
00465       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00466     if (ad_printf) (*ad_printf)("; min = %lg", fmin);
00467     if (ad_printf) (*ad_printf)("; max = %lg\n", fmax);
00468 
00469     x=dmin(fmin+.001,fmin+.01*(fmax-fmin));
00470   }
00471 
00472   if (x > fmax)
00473   {
00474     if (ad_printf)
00475     {
00476       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00477       (*ad_printf)("; min = %lg", fmin);
00478       (*ad_printf)("; max = %lg\n", fmax);
00479     }
00480 
00481     x=dmax(fmax-.001,fmax-.01*(fmax-fmin));
00482   }
00483   if (gradient_structure::Hybrid_bounded_flag==0)
00484   {
00485     tinv=::asin(2.*(x-fmin)/(fmax-fmin)-1.)/1.57079632679489661;
00486   }
00487   else
00488   {
00489     //double y=(x-fmin)/(fmax-fmin);
00490     //double u=1.e-20+y/(1.e-20+(1.0-y));
00491     //tinv=-log(u);
00492     double y=1.e-20+(fmax-x)/(1.e-20+(x-fmin));
00493     tinv=-log(y);
00494   }
00495 
00496   return(tinv);
00497 }
00503 double dmin(double x, double y)
00504 {
00505   if (x<y)
00506   {
00507     return (x);
00508   }
00509   else
00510   {
00511     return(y);
00512   }
00513 }
00519 double dmax(double x, double y)
00520 {
00521   if (x>y)
00522   {
00523     return (x);
00524   }
00525   else
00526   {
00527     return(y);
00528   }
00529 }