$search
00001 #ifndef _GP_HET_REGRESSION_HPP_ 00002 #define _GP_HET_REGRESSION_HPP_ 00003 00004 #include "gpReg/types.hpp" 00005 #include "gpReg/global.h" 00006 #include "gpRegression.hpp" 00007 00008 #include <assert.h> 00009 #include "covarianceFunction.hpp" 00010 #ifdef DO_PROFILING 00011 #include "profiler.hpp" 00012 #endif 00013 00014 00015 00016 // ----------------------------------------------------------------------- 00017 class CovNonStatSE : public CovFunc<TDoubleVector> 00018 { 00019 public: 00020 00021 CovNonStatSE( GPReg<TDoubleVector> &ellGP, double ellMean, double sf2 ) : 00022 m_ellGP(ellGP), m_ellMean(ellMean), m_sf2(sf2) 00023 { 00024 m_ell = 0.13; 00025 } 00026 //Cov1dSEPolar( ConfigFile &cfg, const char *cfgSection ); 00027 00028 virtual double getCov( const TDoubleVector &x, const TDoubleVector &y ) 00029 { 00030 double dist = fabs(x[0]-y[0]); 00032 double covBlock = 0.0; 00033 double a, b; 00034 if (x[0]<y[0]) { 00035 a = x[0]; 00036 b = y[0]; 00037 } else { 00038 a = y[0]; 00039 b = x[0]; 00040 } 00041 double temp; 00042 TDoubleVector z(1); 00043 for (z[0]=a; z[0]<=b; z[0]+=0.01) { 00044 m_ellGP.evalGP( z, temp ); 00045 temp += m_ellMean; 00046 covBlock += 1.0 / (1.0+exp(-temp)); 00047 } 00048 double f = (1.0-covBlock/0.01); 00049 if (f<0) f=0.0; 00050 //*/ 00051 00052 /* 00053 double covBlock = 1.0; 00054 double f = 1.0; 00055 double blockPoint = 0.5; 00056 if ( (x[0]<blockPoint && y[0]>blockPoint) || (x[0]>blockPoint && y[0]<blockPoint) ) { 00057 f=0.0; 00058 } 00059 */ 00060 00061 double cov = (m_sf2*m_sf2) * f * exp( -0.5 * pow(dist,2.0)/m_ell ); 00062 00063 //fprintf( stderr, "[%f %f] (%f %f) %f\n", x[0], y[0], ell_x, ell_y, cov ); 00064 FILE *covFile = fopen( "cov.dat", "a" ); 00065 fprintf( covFile, "%f %f %f %f %f\n", x[0], y[0], cov, covBlock, f ); 00066 fclose( covFile ); 00067 return cov; 00068 //return (m_sf2*m_sf2) * f1 * f2 * f3 * exp( -0.5 * dist*dist/(ell_x-ell_y) ); 00069 } 00070 00071 // paciorek 1d 00072 /* 00073 virtual double getCov( const TDoubleVector &x, const TDoubleVector &y ) 00074 { 00075 double ell_x, ell_y; 00076 m_ellGP.evalGP( x, ell_x ); 00077 m_ellGP.evalGP( y, ell_y ); 00078 double dist = x[0]-y[0]; 00079 double f1 = pow( ell_x, 0.25 ); 00080 double f2 = pow( ell_y, 0.25 ); 00081 double f3 = pow( (ell_x+ell_y)/2.0, -0.5 ); 00082 00083 double cov = (m_sf2*m_sf2) * f1 * f2 * f3 * exp( -0.5 * pow(dist,2.0)/((ell_x+ell_y)/2.0) ); 00084 //double cov = (m_sf2*m_sf2) * exp( -0.5 * pow(dist,2.0)/((ell_x+ell_y)/2.0) ); 00085 00086 //fprintf( stderr, "[%f %f] (%f %f) %f\n", x[0], y[0], ell_x, ell_y, cov ); 00087 FILE *covFile = fopen( "cov.dat", "a" ); 00088 fprintf( covFile, "%f %f %f %f %f %f\n", x[0], y[0], ell_x, ell_y, f1*f2*f3, cov ); 00089 fclose( covFile ); 00090 return cov; 00091 //return (m_sf2*m_sf2) * f1 * f2 * f3 * exp( -0.5 * dist*dist/(ell_x-ell_y) ); 00092 } 00093 */ 00094 00095 00096 00097 private: 00098 ConfigFile m_cfg; 00099 00100 public: 00101 GPReg<TDoubleVector> &m_ellGP; 00102 double m_ellMean; 00103 double m_ell; 00104 double m_sf2; 00105 }; 00106 00107 00108 00109 // ----------------------------------------------------------------------- 00110 template <class TInput> 00111 class GPNonStatReg 00112 { 00113 00114 public: 00115 00116 00117 00118 // ----------------------------------------------------------------------- 00119 GPNonStatReg( 00120 double ell, double ellMean, double sf2, double noise ) : 00121 m_numDataPoints(0), m_dataPoints(0), m_targets(0), 00122 m_ellPriorCovFunc(ell,0.0,sf2), m_ellPriorGP(m_ellPriorCovFunc,noise), 00123 m_covFunc(m_ellPriorGP,ellMean,sf2), m_dataGP(m_covFunc,noise) 00124 {} 00125 00126 00127 00128 // ----------------------------------------------------------------------- 00129 ~GPNonStatReg() 00130 {} 00131 00132 00133 00134 // ----------------------------------------------------------------------- 00135 void setDataPoints( TVector<TInput> &dataPoints, TVector<double> &targets, TVector<double> &ellTargets ) 00136 { 00137 m_dataPoints = dataPoints; 00138 m_numDataPoints = dataPoints.size(); 00139 m_targets = targets; 00140 m_ellTargets = ellTargets; 00141 } 00142 00143 00144 00145 // ----------------------------------------------------------------------- 00146 void buildGP() 00147 { 00148 // dummy ellPrior (all targets = ell) 00149 m_ellPriorGP.setDataPoints( m_dataPoints, m_ellTargets ); 00150 m_ellPriorGP.buildGP(); 00151 00152 m_dataGP.setDataPoints( m_dataPoints, m_targets ); 00153 m_dataGP.buildGP(); 00154 } 00155 00156 00157 00158 // ----------------------------------------------------------------------- 00159 void evalGP( const TInput &x, double &mean, double &var ) 00160 { 00161 m_dataGP.evalGP( x, mean, var ); 00162 } 00163 00164 00165 00166 // ----------------------------------------------------------------------- 00167 void evalGP( const TInput &x, double &mean ) 00168 { 00169 m_dataGP.evalGP( x, mean ); 00170 } 00171 00172 00173 00174 // ----------------------------------------------------------------------- 00175 public: 00176 int m_numDataPoints; 00177 TVector<TInput> m_dataPoints; 00178 TVector<double> m_targets; 00179 00180 Cov1dSE m_ellPriorCovFunc; 00181 GPReg<TInput> m_ellPriorGP; 00182 CovNonStatSE m_covFunc; 00183 GPReg<TInput> m_dataGP; 00184 00185 TVector<double> m_ellTargets; 00186 double m_noise; 00187 00188 }; 00189 00190 00191 00192 #endif //_GP_HET_REGRESSION_HPP_