ADMB Documentation  11.5.3197
 All Classes Files Functions Variables Typedefs Friends Defines
scbound.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 #ifdef __TURBOC__
00013   #pragma hdrstop
00014 #endif
00015 
00016 #include <stdlib.h>
00017 #include <stdio.h>
00018 #include <math.h>
00019 
00020 double dmin(double,double);
00021 double dmax(double, double);
00022 
00023 #define USE_BARD_PEN
00024 
00029 dvariable boundp(dvariable xx, double fmin, double fmax,
00030   const prevariable& _fpen, const double& s)
00031 {
00032   prevariable& fpen=(prevariable&) _fpen;
00033 
00034   dvariable t,y,x;
00035   x=xx/s;
00036   const double l4=log(4.0);
00037   dvariable ss=(sin(x*1.57079632679489661)+1.)/2.;
00038   double diff=fmax-fmin;
00039   t=fmin + diff*ss;
00040 #ifdef USE_BARD_PEN
00041   //cout << "xxxxx" << endl;
00042   double pen=.00001/diff;
00043   fpen=fpen-pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00044 #else
00045 
00046   if (x < -.9999)
00047   {
00048     fpen+=(x+0.9999)*(x+0.9999);
00049   }
00050 
00051   if (x > 0.9999)
00052   {
00053     fpen+=(x-0.9999)*(x-0.9999);
00054   }
00055 
00056   if (x < -1.)
00057   {
00058     fpen+=1000.*(x+1.)*(x+1.);
00059   }
00060 
00061   if (x > 1.)
00062   {
00063     fpen+=1000.*(x-1.)*(x-1.);
00064   }
00065 #endif
00066   return(t);
00067 }
00068 
00073 double boundp( double xx, double fmin, double fmax, const double& _fpen,
00074   const double& s)
00075 {
00076   double& fpen = (double&)_fpen;
00077   double t,x;
00078   x=xx/s;
00079   const double l4=log(4.0);
00080   double ss=(sin(x*1.57079632679489661)+1.)/2.;
00081   double diff=fmax-fmin;
00082   t=fmin + diff*ss;
00083 #ifdef USE_BARD_PEN
00084   //cout << "xxxxx" << endl;
00085   double pen=.00001/diff;
00086   fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00087 #else
00088   if (x < -.9999)
00089   {
00090     fpen+=(x+0.9999)*(x+0.9999);
00091   }
00092 
00093   if (x > 0.9999)
00094   {
00095     fpen+=(x-0.9999)*(x-0.9999);
00096   }
00097 
00098   if (x < -1)
00099   {
00100     fpen+=1000*(x+1)*(x+1);
00101   }
00102 
00103   if (x > 1)
00104   {
00105     fpen+=1000*(x-1)*(x-1);
00106   }
00107 #endif
00108   return(t);
00109 }
00110 
00115 double boundpin(double x, double fmin, double fmax,const double& s)
00116 {
00117   double tinv;
00118 
00119   if (x < fmin)
00120   {
00121     if (ad_printf)
00122       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00123     if (ad_printf) (*ad_printf)("; min = %lg", fmin);
00124     if (ad_printf) (*ad_printf)("; max = %lg\n", fmax);
00125 
00126     x=dmin(fmin+.001,fmin+.01*(fmax-fmin));
00127   }
00128 
00129   if (x > fmax)
00130   {
00131     if (ad_printf)
00132       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00133     if (ad_printf) (*ad_printf)("; min = %lg", fmin);
00134     if (ad_printf) (*ad_printf)("; max = %lg\n", fmax);
00135 
00136     x=dmax(fmax-.001,fmax-.01*(fmax-fmin));
00137   }
00138 
00139   tinv=::asin(2.*(x-fmin)/(fmax-fmin)-1.)/1.57079632679489661;
00140   return(s*tinv);
00141 }
00142 
00143 /*
00144 double boundpin(const prevariable& xx, double fmin, double fmax,
00145   const double& s)
00146 {
00147   double tinv;
00148   double x=value(xx);
00149 
00150   if (x < fmin)
00151   {
00152     if (ad_printf)
00153       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00154     if (ad_printf) (*ad_printf)("; min = %lg", fmin);
00155     if (ad_printf) (*ad_printf)("; max = %lg\n", fmax);
00156 
00157     x=dmin(fmin+.001,fmin+.01*(fmax-fmin));
00158   }
00159 
00160   if (x > fmax)
00161   {
00162     if (ad_printf)
00163       (*ad_printf)("variable out of bounds in boundpin: variable = %lg", x);
00164     if (ad_printf) (*ad_printf)("; min = %lg", fmin);
00165     if (ad_printf) (*ad_printf)("; max = %lg\n", fmax);
00166 
00167     x=dmax(fmax-.001,fmax-.01*(fmax-fmin));
00168   }
00169 
00170   tinv=::asin(2.*(x-fmin)/(fmax-fmin)-1.)/1.57079632679489661;
00171   return(s*tinv);
00172 }
00173 */