@@ -125,7 +125,8 @@ void GP::ComputeAlpha()
125125 timer.tic ();
126126#endif
127127 this ->alpha .set_size (y.n_cols );
128- cholbacksub (alpha, this ->L , trans (this ->y - this ->meanvals ));
128+ // cholbacksub(alpha, this->L, trans(this->y - this->meanvals));
129+ solve (alpha, trimatu (this ->L ), trans (this ->y - this ->meanvals ));
129130 this ->need_to_compute_alpha = false ;
130131#ifdef BENCHMARK
131132 printf (" %f seconds to compute alpha\n " , timer.toc ());
@@ -138,7 +139,8 @@ void GP::ComputeW()
138139 if (this ->need_to_compute_w ) {
139140 ComputeChol ();
140141 this ->W .set_size (this ->L .n_rows , this ->L .n_cols );
141- cholbacksub (W, this ->L , eye<Mat<REAL> >(K.n_rows , K.n_cols ));
142+ // cholbacksub(W, this->L, eye<Mat<REAL> >(K.n_rows, K.n_cols));
143+ solve (W, trimatu (this ->L ), eye<Mat<REAL> >(K.n_rows , K.n_cols ));
142144 this ->need_to_compute_w = false ;
143145 }
144146}
@@ -321,7 +323,8 @@ void GP::Predict(const Col<REAL> &Xs, REAL &mu, REAL &var)
321323
322324 // trans(L) is LOWER TRIANGULAR
323325 Col<REAL> v (kstar.n_elem );
324- solve_tri (v, trans (L), kstar, false );
326+ // solve_tri(v, trans(L), kstar, false);
327+ solve (v, trimatu (trans (L)), kstar);
325328 REAL kxsxs = this ->kernel ->Eval (Xs, Xs);
326329 var = kxsxs - dot (v, v) + this ->s2_n ;
327330}
@@ -375,12 +378,15 @@ void GP::PredictGradient(const Col<REAL> &Xs, Col<REAL> &grad, Col<REAL> &vargra
375378 Col<REAL> v (kstar.n_elem );
376379 Col<REAL> tmp (kstar.n_elem );
377380 Mat<REAL> dv (kstar.n_elem , grad.n_elem );
378- solve_tri (v, trans (L), kstar, false );
381+ // solve_tri(v, trans(L), kstar, false);
382+ solve (v, trimatu (trans (L)), kstar);
379383
380384 // printf("About to solve for dv\n");
381- solve_tri (tmp, trans (L), trans (dkstar.row (0 )), false );
385+ // solve_tri(tmp, trans(L), trans(dkstar.row(0)), false);
386+ solve (tmp, trimatu (trans (L)), trans (dkstar.row (0 )));
382387 dv.col (0 ) = tmp;
383- solve_tri (tmp, trans (L), trans (dkstar.row (1 )), false );
388+ // solve_tri(tmp, trans(L), trans(dkstar.row(1)), false);
389+ solve (tmp, trimatu (trans (L)), trans (dkstar.row (1 )));
384390 dv.col (1 ) = tmp;
385391 // dv.print("dv: ");
386392
@@ -408,7 +414,7 @@ void GP::OptimizeNoiseParam(REAL &noise_param, int max_iterations)
408414
409415 Col<REAL> noise_params (1 );
410416 noise_params (0 ) = noise_param;
411-
417+
412418 opt.Initialize (noise_params, &f_eval_noise, &df_eval_noise, &fdf_eval_noise, 0.5 , 0.1 );
413419 opt.Optimize (noise_params, max_iterations);
414420}
@@ -432,7 +438,7 @@ void GP::OptimizeKernelParam(Col<REAL> &kernel_param, int max_iterations)
432438double f_eval_mean (const gsl_vector *x, void *param)
433439{
434440 GP *gp_obj = reinterpret_cast <GP*>(param);
435-
441+
436442 Col<REAL> mean_param (gp_obj->GetMeanFunction ()->GetParamDim ());
437443 for (unsigned int i=0 ; i < mean_param.n_elem ; ++i) {
438444 mean_param (i) = gsl_vector_get (x, i);
@@ -441,14 +447,14 @@ double f_eval_mean(const gsl_vector *x, void *param)
441447 gp_obj->SetMeanFuncParams (mean_param);
442448
443449 double ret = -gp_obj->ComputeLikelihood ();
444-
450+
445451 return ret;
446452}
447453
448454void df_eval_mean (const gsl_vector *x, void *param, gsl_vector *g)
449455{
450456 GP *gp_obj = reinterpret_cast <GP*>(param);
451-
457+
452458 Col<REAL> mean_param (gp_obj->GetMeanFunction ()->GetParamDim ());
453459 for (unsigned int i=0 ; i < mean_param.n_elem ; ++i) {
454460 mean_param (i) = gsl_vector_get (x, i);
@@ -472,26 +478,26 @@ void fdf_eval_mean(const gsl_vector *x, void *param, double *f, gsl_vector *g)
472478double f_eval_kernel (const gsl_vector *x, void *param)
473479{
474480 GP *gp_obj = reinterpret_cast <GP*>(param);
475-
481+
476482 Col<REAL> kernel_param (gp_obj->GetKernelFunction ()->GetParamDim ());
477483 for (unsigned int i=0 ; i < kernel_param.n_elem ; ++i) {
478484 kernel_param (i) = gsl_vector_get (x, i);
479485 if (kernel_param (i) < 1e-6 ) {
480486 return 1e6 ;
481487 }
482488 }
483-
489+
484490 gp_obj->SetKernelFuncParams (kernel_param);
485491
486492 double ret = -gp_obj->ComputeLikelihood ();
487-
493+
488494 return ret;
489495}
490496
491497void df_eval_kernel (const gsl_vector *x, void *param, gsl_vector *g)
492498{
493499 GP *gp_obj = reinterpret_cast <GP*>(param);
494-
500+
495501 Col<REAL> kernel_param (gp_obj->GetKernelFunction ()->GetParamDim ());
496502 for (unsigned int i=0 ; i < kernel_param.n_elem ; ++i) {
497503 kernel_param (i) = gsl_vector_get (x, i);
@@ -515,17 +521,17 @@ void fdf_eval_kernel(const gsl_vector *x, void *param, double *f, gsl_vector *g)
515521double f_eval_noise (const gsl_vector *x, void *param)
516522{
517523 GP *gp_obj = reinterpret_cast <GP*>(param);
518-
524+
519525 REAL noise_param;
520526 noise_param = gsl_vector_get (x, 0 );
521527 if (noise_param < 1e-6 ) {
522528 return 1e6 ;
523529 }
524-
530+
525531 gp_obj->SetNoise (noise_param);
526532
527533 double ret = -gp_obj->ComputeLikelihood ();
528-
534+
529535 return ret;
530536}
531537
0 commit comments