00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #pragma once
00014
00015 #ifndef __FRACTBROWNIANPATH__
00016 #define __FRACTBROWNIANPATH__
00017
00018
00019 #include "DynamicallyConstructedPath.h"
00020 #include "NormalRandomNumberGenerator.h"
00021 #include <algorithm>
00022 #include <vector>
00023 #include <cmath>
00024 #include <math.h>
00025
00026
00027 #ifdef MKL
00028
00029 #define daxpy_ DAXPY
00030 #define ddot_ DDOT
00031 #define dgelsd_ DGELSD
00032 #define dgelss_ DGELSS
00033 #define dnrm2_ DNRM2
00034 #define dgesv_ DGESV
00035 #define dcopy_ DCOPY
00036 #define dgemm_ DGEMM
00037
00038 #endif
00039
00040
00041
00050
00051
00052 extern "C" double dtptrs_(char *UPLO,
00053 char *tran,
00054 char *diag,
00055 int *order,
00056 int *nrhs,
00057 double *matrovectro,
00058 double *vectro,
00059 int *order3,
00060 int *info);
00061
00063 template <typename my_alg_type>
00064 class FractBrownianPath :
00065 public DynamicallyConstructedPath<my_alg_type>
00066 {
00067 typedef std::vector<SCA> MYVEC;
00068 typedef std::vector<LIE> MYVECLIE;
00069 typedef double doublereal;
00070 typedef int integer;
00071
00073 const double Hurst;
00074
00075 static NormalRandomNumberGenerator vgCommonNormal;
00076
00077 virtual LIE GaussianIncrement( double t, int n ) const
00078 {
00079 LIE coord(n,S(vgNormal(t)));
00080 return coord;
00081 }
00082
00083 protected:
00084 NormalRandomNumberGenerator & vgNormal;
00085
00086 public:
00087 FractBrownianPath(const double H=0.8) : vgNormal( FractBrownianPath::vgCommonNormal ), Hurst(H) {};
00088
00089 FractBrownianPath(NormalRandomNumberGenerator & rand, const double H=0.8) : vgNormal( rand ), Hurst(H) {};
00090
00091 ~FractBrownianPath(void){};
00093 mutable MYVECLIE vectro;
00095 mutable MYVECLIE lievalues;
00097 mutable MYVECLIE vectro2;
00099 mutable MYVECLIE lievalues2;
00101 mutable MYVEC matrovectro;
00103 mutable MYVEC themidpoints;
00104
00105 mutable double newleftmidpoint;
00106
00107 mutable double newrightmidpoint;
00108
00109 void IncreaseVectors( const double &midpoint) const
00110 {
00111 char UPLO('U');
00112 char tran('T');
00113 char diag('N');
00114 int nrhs = 1;
00115 int info = 0;
00116 int sizeofvectro = vectro.size();
00117 int order = sizeofvectro;
00118 int order2 = sizeofvectro;
00120 themidpoints.push_back(midpoint);
00121 MYVEC ltemp;
00123 for (int i = 0; i < order; i++)
00124 {
00125 double tempcoeff = .5 *
00126 ( pow(themidpoints[i],2 * Hurst) +
00127 pow(midpoint,2 * Hurst) -
00128 pow( fabs( (themidpoints[i]-midpoint) ) , 2 * Hurst)) ;
00129 ltemp.push_back(tempcoeff);
00130 }
00131 double *pntrvec;
00132 pntrvec = <emp[0];
00133 dtptrs_( &UPLO,
00134 &tran,
00135 &diag,
00136 &order,
00137 &nrhs,
00138 &matrovectro[0],
00139 pntrvec,
00140 &order2,
00141 &info );
00142 double squaredsum(0);
00143 for (int i=0; i< order; i++)
00144 {
00145 squaredsum += pow( ltemp[i], 2 );
00146 matrovectro.push_back(ltemp[i]);
00147 }
00148 double thebeta = pow ( fabs( pow(midpoint, 2*Hurst) - squaredsum ) , .5 );
00149 ltemp.push_back(thebeta);
00150 matrovectro.push_back(thebeta);
00151 LIE temp = GaussianIncrement(1,1);
00152 LIE temp2 = GaussianIncrement(1,2);
00153 vectro.push_back(temp);
00154 LIE tobeaddedtolievalues(vectro[0] * ltemp[0]);
00155 for (int i=1; i< order + 1; i++)
00156 {
00157 tobeaddedtolievalues += ( vectro[i] * ltemp[i] );
00158 }
00159 lievalues.push_back(tobeaddedtolievalues);
00160 vectro2.push_back(temp2);
00161 LIE tobeaddedtolievalues2(vectro2[0] * ltemp[0]);
00162 for (int i=1; i< order + 1; i++)
00163 {
00164 tobeaddedtolievalues2 += ( vectro2[i] * ltemp[i] );
00165 }
00166 lievalues2.push_back(tobeaddedtolievalues2);
00167 };
00168
00169 LIE FindExistingLieValue(const double &point, const int n) const
00170 {
00171
00172 int findplace(-2);
00173 int size(themidpoints.size());
00174 for (int i=0; i< size; i++) if (themidpoints[i] == point) findplace = i;
00175 LIE tobereturned;
00176 if (n==1) tobereturned=lievalues[findplace];
00177 else
00178 {
00179 tobereturned=lievalues2[findplace];
00180 }
00181 if (findplace==-2) std::cout<< "Oh dear, FindExistingLieValue failed";
00182 return tobereturned;
00183 };
00184
00185 void ComputeChildLieIncrements( Increment<my_alg_type> & nvLeft, Increment<my_alg_type> & nvRight, ConstIterator itLeafAbove ) const
00186 {
00187 const dyadic_interval & diAbove = itLeafAbove->first;
00188 const Increment<my_alg_type> & nvAbove = itLeafAbove->second;
00189 newleftmidpoint = ((diAbove.sup() + diAbove.inf()) *.5);
00190 FractBrownianPath::IncreaseVectors(newleftmidpoint);
00191 newleftmidpoint = diAbove.inf();
00192 LIE latest = lievalues[lievalues.size() -1];
00193 LIE latest2 = lievalues2[lievalues2.size() -1];
00194 if (newleftmidpoint == 0) nvLeft.LieValue(latest +latest2);
00195 else
00196 {
00197 LIE leftlievalue = FractBrownianPath::FindExistingLieValue(newleftmidpoint,1);
00198 LIE leftlievalue2 = FractBrownianPath::FindExistingLieValue(newleftmidpoint,2);
00199 nvLeft.LieValue( latest - leftlievalue + latest2 - leftlievalue2);
00200 }
00201 newrightmidpoint = diAbove.sup();
00202 LIE rightlievalue = FractBrownianPath::FindExistingLieValue(newrightmidpoint,1);
00203 LIE rightlievalue2 = FractBrownianPath::FindExistingLieValue(newrightmidpoint,2);
00204 nvRight.LieValue( rightlievalue - latest + rightlievalue2 - latest2);
00205 };
00206
00207 LIE MakeRootLieIncrement( const dyadic_interval &increment ) const
00208 {
00209
00210 assert(
00211 vectro.empty()
00212 );
00213 assert(
00214 increment.inf() == 0
00215 );
00216 newrightmidpoint = increment.sup();
00217 newleftmidpoint = (increment.sup()* 0.5) + (increment.inf()*.5);
00218 double coefficient = pow(newrightmidpoint, Hurst);
00219 LIE tobeaddedtovectro = GaussianIncrement(1,1);
00220 LIE tobeaddedtovectro2 = GaussianIncrement(1,2);
00221 LIE tobeaddedtolievalues = tobeaddedtovectro * coefficient;
00222 LIE tobeaddedtolievalues2 = tobeaddedtovectro2 * coefficient;
00223 vectro.push_back(tobeaddedtovectro);
00224 vectro2.push_back(tobeaddedtovectro2);
00225 themidpoints.push_back(newrightmidpoint);
00226 matrovectro.push_back(coefficient);
00227 lievalues.push_back(tobeaddedtolievalues);
00228 lievalues2.push_back(tobeaddedtolievalues2);
00229 return (tobeaddedtolievalues + tobeaddedtolievalues2);
00230 };
00231
00232 void MakeNeighborRootLieIncrement( LIE & ans, const Iterator & OldRoot ) const
00233 {
00234
00235
00236 dyadic_interval diIncrement = OldRoot->first;
00237 diIncrement.flip_interval();
00238 newrightmidpoint = diIncrement.sup();
00239 IncreaseVectors(newrightmidpoint);
00240 newleftmidpoint = diIncrement.inf();
00241 LIE latest = lievalues[lievalues.size() -1];
00242 LIE latest2 = lievalues2[lievalues2.size() -1];
00243 LIE leftlievalue = FractBrownianPath::FindExistingLieValue(newleftmidpoint,1);
00244 LIE leftlievalue2 = FractBrownianPath::FindExistingLieValue(newleftmidpoint,2);
00245 ans = ( latest - leftlievalue + latest2 - leftlievalue2);
00246 };
00247 };
00248
00249 #endif // __FRACTBROWNIANPATH__