From dfd10bc6598d0fbc2d9753dccd4cacc256794731 Mon Sep 17 00:00:00 2001 From: Parikshit Ram Date: Sun, 20 Apr 2008 19:20:12 +0000 Subject: [PATCH] testing files added to svn --- fastlib2/contrib/pram/opt/build.py | 26 + fastlib2/contrib/pram/opt/main.cc | 201 +++ .../contrib/pram/opt/optimizers_reloaded.h | 1240 +++++++++++++++++ fastlib2/contrib/pram/opt/phi.h | 156 +++ 4 files changed, 1623 insertions(+) create mode 100644 fastlib2/contrib/pram/opt/build.py create mode 100755 fastlib2/contrib/pram/opt/main.cc create mode 100644 fastlib2/contrib/pram/opt/optimizers_reloaded.h create mode 100644 fastlib2/contrib/pram/opt/phi.h diff --git a/fastlib2/contrib/pram/opt/build.py b/fastlib2/contrib/pram/opt/build.py new file mode 100644 index 0000000000..8258311608 --- /dev/null +++ b/fastlib2/contrib/pram/opt/build.py @@ -0,0 +1,26 @@ + +#librule( + #name = "mog_em", # this line can be safely omitted + #sources = [], # files that must be compiled + #headers = [],# include files part of the 'lib' + #deplibs = ["fastlib:fastlib"], # depends on fastlib core + #tests = ["mog_em_tests.cc"] # this file contains a main with test functions + #) + +binrule( + name = "main", # the executable name + sources = ["main.cc"], # compile main.cc + headers = ["phi.h","optimizers_reloaded.h"], # no extra headers + deplibs = ["fastlib:fastlib"] # + ) + +# to build: +# 1. make sure have environment variables set up: +# $ source /full/path/to/fastlib/script/fl-env /full/path/to/fastlib +# (you might want to put this in bashrc) +# 2. fl-build main +# - this automatically will assume --mode=check, the default +# - type fl-build --help for help +# 3. ./main +# - to build same target again, type: make +# - to force recompilation, type: make clean diff --git a/fastlib2/contrib/pram/opt/main.cc b/fastlib2/contrib/pram/opt/main.cc new file mode 100755 index 0000000000..07a587a768 --- /dev/null +++ b/fastlib2/contrib/pram/opt/main.cc @@ -0,0 +1,201 @@ +#include +#include "phi.h" +#include "optimizers_reloaded.h" + +long double test_function(Vector& theta, const Matrix& data, Vector* grad) { + + // First we make a model out of the theta + // so mu <- theta[1,..,dim] + // U = [ theta[dim+1] 0 0 ...; + // theta[dim+2] theta[dim+3] 0 ...; + // theta[dim+4] theta[dim+5] theta[dim+6] 0 ...; + // ..... + // ]; + // Sigma = U' * U; + + index_t dim = data.n_rows(), n = data.n_cols(); + Matrix sigma, lower_triangle_matrix, upper_triangle_matrix; + ArrayList d_sigma; + Vector mu; + double s_min = 0.01; + double *temp_mu; + + // obtaining the mu values + temp_mu = (double*)malloc(dim * sizeof(double)); + + for(index_t i = 0; i < dim; i++) { + temp_mu[i] = theta.get(i); + } + mu.Copy(temp_mu, dim); + + //printf("Mu : ["); + //for(index_t i = 0; i < dim; i++) { + //printf("%lf ",mu.get(i)); + //} + //printf("\b];\n"); + + // obtaining the sigma and d_sigma values + d_sigma.Init(dim*(dim+1)/2); + + // the sigma values + lower_triangle_matrix.Init(dim, dim); + lower_triangle_matrix.SetAll(0.0); + for(index_t i = 0; i < dim; i++) { + for(index_t j = 0; j < i; j++) { + lower_triangle_matrix.set(i, j, theta[dim + i*(i+1)/2 + j]); + } + // adding small value to the diagonal of the + // covariance matrix to stop it from going to + // infinity by obtaining zero determinant of + // covariance + lower_triangle_matrix.set(i, i, theta[dim + i*(i+1)/2 + i] + s_min); + } + la::TransposeInit(lower_triangle_matrix, &upper_triangle_matrix); + la::MulInit(lower_triangle_matrix, upper_triangle_matrix, &sigma); + + // the d_sigma values + Matrix d_sigma_d_r, d_sigma_d_r_t, temp_matrix1, temp_matrix2; + d_sigma_d_r.Init(dim, dim); + d_sigma_d_r_t.Init(dim, dim); + temp_matrix1.Init(dim, dim); + temp_matrix2.Init(dim, dim); + + for(index_t i = 0; i < dim; i++) { + for(index_t j = 0; j < i+1; j++) { + d_sigma_d_r.SetAll(0.0); + d_sigma_d_r.set(i, j, 1.0); + la::TransposeOverwrite(d_sigma_d_r, &d_sigma_d_r_t); + + la::MulOverwrite(d_sigma_d_r, upper_triangle_matrix, &temp_matrix1); + la::MulOverwrite(lower_triangle_matrix, d_sigma_d_r_t, &temp_matrix2); + la::AddInit(temp_matrix1, temp_matrix2, &d_sigma[i*(i+1)/2 + j]); + } + } + + //printf("Sigma : ["); + //for(index_t i = 0; i < dim; i++) { + //for(index_t j = 0; j < dim; j++) { + // printf("%lf ",sigma.get(i,j)); + //} + //printf("\b;"); + //} + //printf("\b]\n"); + + // calculating the value of the function for each data point + // and adding it up + // f_theta(x_i) = -log phi(x_i, mu, sigma); + // g_mu(x_i) = (-1/phi(x_i, mu, sigma)) * d phi / d mu; + // g_sigma(x_i) = (-1/phi(x_i, mu, sigma)) * d phi / d sigma; + // l_theta = \sum_{i=1}^N f_theta(x_i); + // g_l_theta = grad = \sum_{i=1}^N [g_mu(x_i) g_sigma(x_i)] + + long double l_theta = 0.0, f_theta, tmp_val; + Vector g_mu, g_sigma, x; + x.Init(dim); + g_mu.Init(dim); + g_mu.SetZero(); + g_sigma.Init(dim*(dim+1)/2); + g_sigma.SetZero(); + + //printf("calc func val\n"); + //fflush(NULL); + for(index_t i = 0; i < n; i++) { + Vector d_phi_d_mu, d_phi_d_sigma; + x.CopyValues(data.GetColumnPtr(i)); + tmp_val = phi(x, mu, sigma, d_sigma, + &d_phi_d_mu, &d_phi_d_sigma); + f_theta = -log(tmp_val); + + double alpha = -1.0 / tmp_val; + la::AddExpert(alpha, d_phi_d_mu, &g_mu); + la::AddExpert(alpha, d_phi_d_sigma, &g_sigma); + + l_theta += f_theta; + } + + //printf("%Lf setting grad val\n", l_theta); + //fflush(NULL); + double *temp_grad; + temp_grad = (double*)malloc(theta.length() * sizeof(double)); + for(index_t i = 0; i < dim; i++) { + temp_grad[i] = g_mu.get(i); + } + + for(index_t i = 0; i < dim*(dim+1)/2; i++) { + temp_grad[dim+i] = g_sigma.get(i); + } + grad->CopyValues(temp_grad); + return l_theta; +} + +int main(int argc, char* argv[]) { + + fx_init(argc, argv); + + const char *datafile = fx_param_str_req(NULL, "data"); + + Matrix data_points; + data::Load(datafile, &data_points); + + double temp_array[] = {4, -2, 3, 1, 2}; + double *real_theta_array; + long double function_val = 4607.5; + Vector real_theta, grad; + index_t len = 5; + + real_theta_array = (double*)malloc(5 * sizeof(double)); + for(index_t i = 0; i < 5; i++) { + real_theta_array[i] = temp_array[i]; + } + real_theta.Copy(real_theta_array, len); + grad.Init(5); + + //printf("entering func\n"); + //fflush(NULL); + long double val = test_function(real_theta, data_points, &grad); + //printf("exiting func\n"); + //fflush(NULL); + printf("%Lf %Lf\n", val, function_val); + + datanode *opt_module = fx_submodule(NULL,"opt","opt"); + fx_param_int(opt_module,"param_space_dim", 5); + + //QuasiNewton opt; + SMD_SingleStep opt; + + // fx_param_str(opt_module, "method", "QuasiNewton"); + fx_param_str(opt_module, "method", "SMD_SingleStep"); + opt.Init(test_function, data_points, opt_module); + double *pt; + double p[] = {1, 2, 3, 1, 6}; + pt = (double*)malloc(5 * sizeof(double)); + for(index_t i = 0; i < 5; i++) { + pt[i] = p[i]; + } + fx_timer_start(opt_module,"opt_time"); + opt.Eval(pt); + fx_timer_stop(opt_module,"opt_time"); + + printf("theta : ["); + for(index_t i = 0; i < 5; i++) { + printf(" %lf,",pt[i]); + } + printf("\b ]\n"); + + Vector calc_theta; + calc_theta.Copy(pt, 5); + long double min_ob = test_function(calc_theta, data_points, &grad); + printf("%Lf\n",min_ob); + //fx_silence(); + fx_done(); + + return 1; +} + +/** + * The actual minimum obtained by the Quasi Newton + * method is 4599.772730 + * The minima is [3.870753, -1.952005, 2.819059, 0.831117, 2.048447] + * The time required was : 0.255378 sec + * Iterations through the data : 19 + */ diff --git a/fastlib2/contrib/pram/opt/optimizers_reloaded.h b/fastlib2/contrib/pram/opt/optimizers_reloaded.h new file mode 100644 index 0000000000..cefc9eb733 --- /dev/null +++ b/fastlib2/contrib/pram/opt/optimizers_reloaded.h @@ -0,0 +1,1240 @@ +/** + * @author Parikshit Ram (pram@cc.gatech.edu) + * @file optimizers.h + * + * Implements classes for two types of optimizer + * + */ + +#ifndef OPTIMIZERS_H +#define OPTIMIZERS_H + +#include + +/** + * An optimizer using the Nelder Mead method, + * also known as the polytope or the simplex + * method. + * + * It does multivariate minimization of an + * objective function. If it is optimizing in + * 'd' dimensions, it would require 'd+1' + * starting points. + * + * Example use: + * + * @code + * double init_pts[d+1][d]; + * index_t number_of_function_evaluations; + * struct datanode *opt_module = fx_submodule(NULL,"NelderMead","opt_module"); + * Matrix data; + * index_t dim_param_space; + * + * ... + * NelderMead opt; + * opt.Init(obj_function, data, dim_param_space, opt_module); + * ... + * opt.Eval(init_pts); + * // init_pts[0] contains the optimal point found + * @endcode + * + */ +class NelderMead { + + private: + index_t dimension_; + Matrix data_; + long double (*func_ptr_)(Vector&, const Matrix&); + datanode *opt_module_; + + public: + + NelderMead() { + } + + ~NelderMead() { + } + + void Init(long double (*fun)(Vector&, const Matrix&), + Matrix& data, datanode *opt_module) { + + data_.Copy(data); + func_ptr_ = fun; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix& data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double **pts) { + + index_t dim = dimension(), num_func_eval; + index_t i, j, ihi, ilo, inhi,mpts = dim + 1; + double sum, swap, *psum; + long double swap_y, rtol, ytry, ysave, TINY = 1.0e-10; + long double *y; + Vector param_passed; + long double tol = fx_param_double(opt_module_,"tolerance", 1.0e-7); + index_t NMAX = fx_param_int(opt_module_, "MAX_FUNC_EVAL", 50000); + + param_passed.Init(dim); + psum = (double*)malloc(dim * sizeof(double)); + num_func_eval = 0; + y = (long double*)malloc(mpts*sizeof(long double)); + for(i = 0; i < mpts; i++) { + param_passed.CopyValues(pts[i]); + y[i] = (*func_ptr_)(param_passed,data()); + + } + + + for(;;) { + ilo = 0; + ihi = y[0] > y[1] ? (inhi = 1,0) : (inhi = 0,1); + for( i = 0; i < mpts; i++ ) { + if(y[i] <= y[ilo]) ilo = i; + if(y[i] > y[ihi]) { + inhi = ihi; + ihi = i; + } + else if((y[i] > y[inhi])&&(i != ihi)) inhi = i; + } + + rtol = 2.0 * fabs(y[ihi] - y[ilo]) / ( fabs(y[ihi]) + fabs(y[ilo]) + TINY ) ; + if(rtol < tol) { + swap_y = y[0]; + y[0] = y[ilo]; + y[ilo] = swap_y; + for( i = 0; i < dim; i++ ) { + swap = pts[0][i]; + pts[0][i] = pts[ilo][i] ; + pts[ilo][i] = swap; + } + fx_format_result(opt_module_,"min_obtained","%Lf", y[0]); + break; + } + if(num_func_eval > NMAX){ + fx_format_result(opt_module_,"min_obtained","%Lf", y[ilo]); + NOTIFY("Maximum number of function evaluations exceeded"); + break; + } + num_func_eval += 2; + + // Beginning a new iteration. + // Extrapolating by a factor of -1.0 through the face of the simplex + // across from the high point, i.e, reflect the simplex from the high point + for( j = 0 ; j < dim ; j++ ){ + sum = 0.0; + for( i = 0 ; i < mpts ; i++ ) + if (i != ihi) + sum += pts[i][j]; + + psum[j] = sum / dim; + } + + ytry = ModSimplex_(pts, y, psum, ihi, -1.0); + if( ytry <= y[ilo] ) { + // result better than best point + // so additional extrapolation by a factor of 2 + ytry = ModSimplex_(pts, y, psum, ihi, 2.0); + } + else if( ytry >= y[ihi] ) { + // result worse than the worst point + // so there is a lower intermediate point, + // i.e., do a one dimensional contraction + ysave = y[ihi]; + + ytry = ModSimplex_(pts, y, psum, ihi, 0.5); + if( ytry > y[ihi] ) { + // Can't get rid of the high point, + // try to contract around the best point + for( i = 0; i < mpts; i++ ) { + if( i != ilo ) { + for( j = 0; j < dim; j++ ) { + pts[i][j] = psum[j] = 0.5 * ( pts[i][j] + pts[ilo][j] ); + } + param_passed.CopyValues(psum); + y[i] = (*func_ptr_)(param_passed, data()); + } + } + num_func_eval += dim; + for( j = 0 ; j < dim ; j++ ){ + sum = 0.0; + for( i = 0 ; i < mpts ; i++ ) + if (i != ihi) + sum += pts[i][j]; + psum[j] = sum / dim; + } + } + } + else --num_func_eval; + } + fx_format_result(opt_module_, "func_evals", "%d", num_func_eval); + return; + } + + long double ModSimplex_(double **pts, long double *y, double *psum, + index_t ihi, float fac) { + + + index_t j, dim = dimension(); + long double ytry; + double *ptry; + Vector param_passed; + + param_passed.Init(dim); + ptry = (double*) malloc (dim * sizeof(double)); + for (j = 0; j < dim; j++) { + ptry[j] = psum[j] * (1 - fac) + pts[ihi][j] * fac; + } + param_passed.CopyValues(ptry); + ytry = (*func_ptr_)(param_passed, data()); + + if (ytry < y[ihi]) { + y[ihi] = ytry; + for (j = 0; j < dim; j++) { + pts[ihi][j] = ptry[j]; + } + } + return ytry; + } + +}; + +/** + * An optimizer using the Quasi Newton method, + * also known as the variable metrics + * method. + * + * It does multivariate minimization of an + * objective function using only the function + * value and the gradients. + * + * Example use: + * + * @code + * double init_pt[d]; + * index_t number_of_iters; + * struct datanode *opt_module = fx_submodule(NULL,"QuasiNewton","opt_module"); + * Matrix data; + * index_t dim_param_space; + * + * ... + * QuasiNewton opt; + * opt.Init(obj_function, data, dim_param_space, opt_module); + * ... + * opt.Eval(init_pt); + * // init_pt contains the optimal point found + * @endcode + * + */ + +class QuasiNewton { + + private: + index_t dimension_; + Matrix data_; + long double (*func_ptr_)(Vector&, const Matrix&, Vector*); + datanode *opt_module_; + + public: + + QuasiNewton(){ + } + + ~QuasiNewton(){ + } + + void Init(long double (*fun)(Vector&, const Matrix&, Vector*), + Matrix& data, datanode *opt_module){ + + data_.Copy(data); + func_ptr_ = fun; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double *pt){ + + index_t n = dimension(), iters; + index_t i, its, MAXIMUM_ITERATIONS = fx_param_int(opt_module_,"MAX_ITERS",500); + long double temp_1, temp_2, temp_3, temp_4, f_previous, f_min, + maximum_step_length, sum = 0.0, sumdg, sumxi, temp, test; + Vector dgrad, grad, hdgrad, xi; + Vector pold, pnew; + Matrix hessian; + double EPSILON = fx_param_double(opt_module_, "EPSILON", 3.0e-8); + fx_format_param(opt_module_, "TOLERANCE", "%lf", 1.0e-5); + double TOLERANCE = fx_param_double_req(opt_module_, "TOLERANCE"); + double MAX_STEP_SIZE = fx_param_double(opt_module_, "MAX_STEP_SIZE", 100.0); + double g_tol = fx_param_double(opt_module_, "gtol", 1.0e-7); + + dgrad.Init(n); + grad.Init(n); + hdgrad.Init(n); + hessian.Init(n,n); + pnew.Init(n); + xi.Init(n); + pold.Copy(pt,n); + f_previous = (*func_ptr_)(pold, data(), &grad); + Vector tmp; + tmp.Init(n); + tmp.SetAll(1.0); + hessian.SetDiagonal(tmp); + la::ScaleOverwrite(-1.0, grad, &xi); + + sum = la::Dot(pold, pold); + double fmax; + if( sqrt(sum) > (float)n ) { + fmax = sqrt(sum); + } + else { + fmax = (float)n; + } + maximum_step_length = MAX_STEP_SIZE*fmax; + + for(its = 0; its < MAXIMUM_ITERATIONS; its++) { + + dgrad.CopyValues(grad); + LineSearch_(pold, f_previous, &grad, &xi, + &pnew, &f_min, maximum_step_length); + f_previous = f_min; + la::SubOverwrite(pold, pnew, &xi); + pold.CopyValues(pnew); + + for(i = 0; i < n; i++) { + pt[i] = pold.get(i); + } + + test = 0.0; + for(i = 0; i < n; i++){ + if(fabs(pold.get(i)) > 1.0) fmax = fabs(pold.get(i)); + else{ fmax = 1.0; } + temp = fabs(xi.get(i)) / fmax; + if(temp > test) test = temp; + } + if(test < TOLERANCE) { + iters = its; + fx_format_result(opt_module_, "iters", "%d", iters); + fx_format_result(opt_module_,"min_obtained","%Lf", f_previous); + return; + } + + test = 0.0; + if(f_min > 1.0) temp_1 = f_min; + else{ temp_1 = 1.0; } + + for(i = 0; i < n; i++) { + if(fabs(pold.get(i)) > 1.0) fmax = pold.get(i); + else{ fmax = 1.0; } + + temp = fabs(grad.get(i))*fmax / temp_1; + if(temp > test) test = temp; + } + if(test < g_tol) { + iters = its; + fx_format_result(opt_module_, "iters", "%d", iters); + fx_format_result(opt_module_,"min_obtained","%Lf", f_previous); + return; + } + + la::SubFrom(grad, &dgrad); + la::Scale(-1.0, &dgrad); + la::MulOverwrite(hessian,dgrad, &hdgrad); + + temp_2 = la::Dot(dgrad, xi); + temp_4 = la::Dot(dgrad, hdgrad); + sumdg = la::Dot(dgrad, dgrad); + sumxi = la::Dot(xi, xi); + + if (temp_2 > sqrt(EPSILON*sumdg*sumxi)) { + temp_2 = 1.0 / temp_2; + temp_3 = 1.0 / temp_4; + + la::ScaleOverwrite(temp_2, xi, &dgrad); + la::AddExpert((-1.0*temp_3), hdgrad, &dgrad); + + Matrix co, ro, tmp; + co.AliasColVector(xi); + ro.AliasRowVector(xi); + la::MulInit(co, ro, &tmp); + la::AddExpert(temp_2, tmp, &hessian); + + co.Destruct(); + ro.Destruct(); + tmp.Destruct(); + co.AliasColVector(hdgrad); + ro.AliasRowVector(hdgrad); + la::MulInit(co, ro, &tmp); + la::AddExpert((-1.0*temp_3), tmp, &hessian); + + co.Destruct(); + ro.Destruct(); + tmp.Destruct(); + co.AliasColVector(dgrad); + ro.AliasRowVector(dgrad); + la::MulInit(co, ro, &tmp); + la::AddExpert(temp_4, tmp, &hessian); + } + la::MulOverwrite(hessian, grad, &xi); + la::Scale((-1.0), &xi); + } + NOTIFY("Too many iterations in Quasi Newton\n"); + fx_format_result(opt_module_,"min_obtained","%Lf", f_previous); + } + + void LineSearch_(Vector pold, long double fold, Vector *grad, + Vector *xi, Vector *pnew, long double *f_min, + long double maximum_step_length){ + + index_t i, n = dimension(); + long double a, step_length, previous_step_length = 0.0, + minimum_step_length, b, disc, previous_f_value = 0.0, + rhs1, rhs2, slope, sum, temp, test, temp_step_length, + MIN_DECREASE = 1.0e-4, TOLERANCE = 1.0e-7; + + sum = la::Dot(*xi, *xi); + sum = sqrt(sum); + if(sum > maximum_step_length) { + la::Scale((maximum_step_length/sum), xi); + } + slope = la::Dot(*grad, *xi); + if(slope >= 0.0){ + return; + } + test = 0.0; + for(i = 0; i < n; i++) { + double fmax; + fmax = (fabs(pold.get(i)) > 1.0 ? fabs(pold.get(i)) : 1.0); + temp = fabs((*xi).get(i)) / fmax; + if(temp > test) test = temp; + } + + minimum_step_length = TOLERANCE/test; + step_length = 1.0; + for(;;) { + + pnew->CopyValues(pold); + la::AddExpert(step_length, *xi, pnew); + + *f_min = (*func_ptr_)((*pnew), data(), grad); + if(step_length < minimum_step_length) { + pnew->CopyValues(pold); + return; + } + else if( *f_min <= fold + MIN_DECREASE*step_length*slope) { + return; + } + else { + if (step_length == 1.0) { + temp_step_length = -slope/(2.0*(*f_min - fold - slope)); + } + else { + rhs1 = *f_min - fold - step_length*slope; + rhs2 = previous_f_value - fold - previous_step_length*slope; + a = (rhs1 / (step_length*step_length) + - rhs2/(previous_step_length*previous_step_length)) + / (step_length-previous_step_length); + b = (-previous_step_length*rhs1/(step_length*step_length) + +step_length*rhs2/(previous_step_length*previous_step_length)) + / (step_length - previous_step_length); + if(a == 0.0) { + temp_step_length = -slope / (2.0*b); + } + else { + disc = b*b - 3.0*a*slope; + if(disc < 0.0) { + temp_step_length = 0.5*step_length; + } + else if (b <= 0.0) { + temp_step_length = (-b+sqrt(disc))/(3.0*a); + } + else { + temp_step_length = -slope / (b+sqrt(disc)); + } + } + if(temp_step_length > 0.5*step_length) { + temp_step_length = 0.5*step_length; + } + } + } + previous_step_length = step_length; + previous_f_value = *f_min; + step_length = (temp_step_length > 0.1*step_length + ? temp_step_length : 0.1*step_length); + } + } +}; + +/** + * Normal Gradient Descent implemented here + * documentation later + * + */ + +class GradientDescent { + + private: + index_t dimension_; + Matrix data_; + long double (*func_ptr_)(Vector&, const Matrix&, Vector*); + datanode *opt_module_; + + public: + + GradientDescent(){ + } + + ~GradientDescent(){ + } + + void Init(long double (*fun)(Vector&, const Matrix&, Vector*), + Matrix& data, datanode *opt_module){ + + data_.Copy(data); + func_ptr_ = fun; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double *pt){ + + index_t iters; + index_t MAXIMUM_ITERATIONS = fx_param_int(opt_module_,"MAX_ITERS",100); + double EPSILON = fx_param_double(opt_module_, "EPSILON", 1.0e-5); + fx_format_param(opt_module_, "TOLERANCE", "%lf", 0.001); + double TOLERANCE = fx_param_double_req(opt_module_, "TOLERANCE"); + // double MAX_STEP_SIZE = fx_param_double(opt_module_, + // "MAX_STEP_SIZE", 100.0); + index_t dim = fx_param_int_req(opt_module_, "param_space_dim"); + Vector pold, pnew, grad; + long double f_old, f_new; + double scale, alpha = 0.1, gamma; + long double p_tol = 0.0, f_tol = 0.0; + + // have to decide what to assign alpha value as + // step lengths are crucial because this is + // ending up oscillating close to the optimal + // hence never actually reaching the optimal + + pold.Init(dim); + pnew.Init(dim); + pold.CopyValues(pt); + grad.Init(dim); + + f_old = (*func_ptr_)(pold, data(), &grad); + printf("first val: %Lf\n", f_old); + + // Here we are doing the normal gradient step + // scale = || - \nabla_\theta f(X, \theta_k) || + // \theta_{k+1} = \theta_k - + // alpha * \nabla_\theta f(X,\theta_k) / scale; + + for (iters = 0; iters < MAXIMUM_ITERATIONS; iters++) { + + scale = sqrt(la::Dot(grad, grad)); + gamma = - alpha / scale; + pnew.SetZero(); + la::AddTo(pold, &pnew); + la::AddExpert(gamma, grad, &pnew); + + Vector diff; + la::SubInit(pnew, pold, &diff); + p_tol = sqrt(la::Dot(diff, diff)); + + f_new = (*func_ptr_)(pnew, data(), &grad); + f_tol = fabs(f_new - f_old); + + if (((f_tol < EPSILON) && (p_tol < TOLERANCE)) || (scale < EPSILON)) { + fx_format_result(opt_module_, "iters", "%d", iters+1); + fx_format_result(opt_module_,"min_obtained","%Lf", f_old); + for (index_t i = 0; i < dim; i++) { + pt[i] = pold.get(i); + } + return; + } + + pold.CopyValues(pnew); + f_old = f_new; + } + + NOTIFY("Too many iterations in Gradient Descent\n"); + fx_format_result(opt_module_,"min_obtained","%Lf", f_old); + for(index_t i = 0; i < dim; i++) { + printf("%lf, ", pold.get(i)); + } + printf("\nfinal val: %Lf\n p_tol : %Lf, f_tol : %Lf, iters : %"LI"d\n", f_old, p_tol, f_tol, iters); + return; + } + +}; + +/** + * Stochastic Gradient Descent implemented here + * documentation later + * + */ + +class SGD { + + private: + index_t dimension_; + Matrix data_; + long double (*func_ptr_)(Vector&, const Matrix&, Vector*); + datanode *opt_module_; + + public: + + SGD(){ + } + + ~SGD(){ + } + + void Init(long double (*fun)(Vector&, const Matrix&, Vector*), + Matrix& data, datanode *opt_module){ + + data_.Copy(data); + func_ptr_ = fun; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double *pt){ + + index_t iters; + index_t MAXIMUM_ITERATIONS = fx_param_int(opt_module_,"MAX_ITERS",100); + double EPSILON = fx_param_double(opt_module_, "EPSILON", 1.0e-5); + fx_format_param(opt_module_, "TOLERANCE", "%lf", 0.001); + double TOLERANCE = fx_param_double_req(opt_module_, "TOLERANCE"); + // double MAX_STEP_SIZE = fx_param_double(opt_module_, + // "MAX_STEP_SIZE", 100.0); + index_t dim = fx_param_int_req(opt_module_, "param_space_dim"); + index_t num_batch = fx_param_int(opt_module_, "BATCHES",50); + Vector pold, pnew, grad; + long double f_old, f_new; + double scale, alpha = 0.1, gamma; + long double p_tol = 0.0, f_tol = 0.0; + Matrix data_batched; + index_t batch_size = data().n_cols() / num_batch; + + // have to decide what to assign alpha value as + // step lengths are crucial because this is + // ending up oscillating close to the optimal + // hence never actually reaching the optimal + + pold.Init(dim); + pnew.Init(dim); + pold.CopyValues(pt); + grad.Init(dim); + data_batched.Copy(data()); + + f_old = (*func_ptr_)(pold, data(), &grad); + printf("first val: %Lf\n", f_old); + + // Here we are doing the normal gradient step + // scale = || - \nabla_\theta f(X_t, \theta_t) || + // \theta_{t+1} = \theta_t - + // alpha * \nabla_\theta f(X_t,\theta_t) / scale; + + for (iters = 0; iters < MAXIMUM_ITERATIONS; iters++) { + + // Now going through the data batchwise + for (index_t in = 0; in < num_batch; in++) { + + scale = sqrt(la::Dot(grad, grad)); + gamma = - alpha / scale; + pnew.SetZero(); + la::AddTo(pold, &pnew); + la::AddExpert(gamma, grad, &pnew); + + Vector diff; + la::SubInit(pnew, pold, &diff); + p_tol = sqrt(la::Dot(diff, diff)); + + // using a batch + Matrix single_batch; + index_t st_pt = in * batch_size; + data_batched.MakeColumnSlice(st_pt, batch_size, &single_batch); + f_new = (*func_ptr_)(pnew, single_batch, &grad); + f_tol = fabs(f_new - f_old); + + if ((f_tol < EPSILON) && (p_tol < TOLERANCE)) { + fx_format_result(opt_module_, "iters", "%d", iters+1); + fx_format_result(opt_module_,"min_obtained","%Lf", f_old); + for (index_t i = 0; i < dim; i++) { + pt[i] = pold.get(i); + } + printf("iters: %"LI"d, min: %Lf\n", iters, f_old); + return; + } + + pold.CopyValues(pnew); + f_old = f_new; + + } + + // permuting the data matrix + data_batched.Destruct(); + PermuteMatrix_(data(), &data_batched); + //printf("data permuted\n"); + } + + NOTIFY("Too many iterations in Stochastic Gradient Descent\n"); + fx_format_result(opt_module_,"min_obtained","%Lf", f_old); + for(index_t i = 0; i < dim; i++) { + printf("%lf, ", pold.get(i)); + } + long double f_final = (*func_ptr_)(pold, data(), &grad); + printf("\nfinal val: %Lf\n p_tol : %Lf, f_tol : %Lf, iters : %"LI"d\n", + f_final, p_tol, f_tol, iters); + return; + } + + void PermuteMatrix_(const Matrix& input, Matrix *output) { + + ArrayList perm_array; + index_t size = input.n_cols(); + Matrix perm_mat; + + perm_mat.Init(size, size); + perm_mat.SetAll(0.0); + + math::MakeRandomPermutation(size, &perm_array); + for(index_t i = 0; i < size; i++) { + perm_mat.set(perm_array[i], i, 1.0); + } + + la::MulInit(input, perm_mat, output); + return; + } +}; + + +/** + * Stochastic Meta Descent with a + * Single step model implemented here + * documentation later + * + */ + +/* +class SMD_SingleStep { + + private: + index_t dimension_; + Matrix data_; + // This is the original way of calling the function + long double (*func_ptr_)(Vector&, const Matrix&, Vector*); + + // But since the L2 function is pretty awesome, we need + // to make the calling a little diff + // long double (*func_ptr_stoc_)(Vector&, const Matrix&, Vector*, index_t); + datanode *opt_module_; + + public: + + SMD_SingleStep(){ + } + + ~SMD_SingleStep(){ + } + + void Init(long double (*fun)(Vector&, const Matrix&, Vector*), + //long double (*fun_stoc)(Vector&, const Matrix&, Vector*, index_t), + Matrix& data, datanode *opt_module){ + + data_.Copy(data); + + // function pointer to the original function + func_ptr_ = fun; + + // pointer to the broken up function used + // for stochastic optimization + // func_ptr_stoc_ = fun_stoc; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double *pt){ + + index_t iters; + index_t MAXIMUM_ITERATIONS = fx_param_int(opt_module_,"MAX_ITERS",100); + // double EPSILON = fx_param_double(opt_module_, "EPSILON", 1.0e-2); + double TOLERANCE = fx_param_double(opt_module_, "TOLERANCE", 1.0e-5); + index_t dim = fx_param_int_req(opt_module_, "param_space_dim"); + index_t num_batch = fx_param_int(opt_module_, "BATCHES",50); + Vector pold, pnew, grad, prev_grad; + long double f_old, f_new; + double scale, scale_prev, eta = 0.1, gamma, mu = 0.05; + double p_tol = 0.0;//, f_tol = 0.0; + Matrix data_batched; + index_t batch_size = fx_param_int(opt_module_,"BATCH_SIZE", + data().n_cols() / num_batch); + + // fx_clear_param(opt_module_,"BATCHES"); + num_batch = fx_param_int(opt_module_,"BATCHES", + data().n_cols() / batch_size); + // to decide how to chose starting value + // of alpha (right now it is just 1). + // also have to decide the value for the + // meta parameter mu (right now it is + // arbitrarily chosen as 0.1) + + pold.Init(dim); + pnew.Init(dim); + pold.CopyValues(pt); + grad.Init(dim); + prev_grad.Init(dim); + data_batched.Copy(data()); + + // This is the original way of calling the function + // f_old = (*func_ptr_)(pold, data(), &grad); + + // But since the L2 function is pretty awesome, we need + // to make the calling a little diff + f_old = (*func_ptr_)(pold, data(), &grad); + // printf("first val: %Lf\n", f_old); + // scale = sqrt(la::Dot(grad, grad)); + + // Here we are doing the gradient step + // scale = || - \nabla_\theta f(X_t, \theta_t) || + // \theta_{t+1} = \theta_t - + // \eta_t * \nabla_\theta f(X_t,\theta_t) / scale; + + for (iters = 0; iters < MAXIMUM_ITERATIONS; iters++) { + + // Now going through the data batchwise + for (index_t in = 0; in < num_batch; in++) { + + // instead of scaling the gradient, how about using low values + // of the step sizes, because scaling the gradients result + // in the gradient being significant even when it is close + // to the optimal + // gamma = - eta / scale; + gamma = -eta; + pnew.SetZero(); + la::AddTo(pold, &pnew); + la::AddExpert(gamma, grad, &pnew); + + + // using a batch + Matrix single_batch; + index_t st_pt = in * batch_size; + data_batched.MakeColumnSlice(st_pt, batch_size, &single_batch); + prev_grad.CopyValues(grad); + // This is the original way of calling the function + f_new = (*func_ptr_)(pnew, single_batch, &grad); + + // But since the L2 function is pretty awesome, we need + // to make the calling a little diff + // f_new = (*func_ptr_stoc_)(pnew, single_batch, &grad, num_batch); + + // Terminating conditions + // |f_t+1 - f_t| < epsilon & ||\theta_t+1 - \theta_t|| < delta + // but since we are using stochastic method, it is + // better to only consider \theta instead of f + Vector diff; + la::SubInit(pnew, pold, &diff); + // f_tol = fabs(f_new - f_old); + p_tol = sqrt(la::Dot(diff, diff)); + + // but instead if we used the condition + // ||grad_t|| < epsilon' & ||\theta_t+1 - \theta_t|| < delta + // if ((f_tol < EPSILON) && (p_tol < TOLERANCE)) { + + // this doesn't work either, same problem + // if ((scale < EPSILON) && (p_tol < TOLERANCE)){ + + // using just the point in the param_space + // which refuses to move + if (p_tol < TOLERANCE) { + // rejected because stops too early, need the check + // the overall gradient is small + Vector temp_grad; + temp_grad.Init(dim); + long double f_final = (*func_ptr_)(pold, data(), &temp_grad); + // but maybe we can skip that now + // double temp_grad_val = sqrt(la::Dot(temp_grad, temp_grad)); + // if (temp_grad_val < EPSILON) { + + fx_format_result(opt_module_, "iters", "%d", iters+1); + fx_format_result(opt_module_,"min_obtained","%Lf", f_final); + for (index_t i = 0; i < dim; i++) { + pt[i] = pold.get(i); + } + // printf("iters: %"LI"d\n", iters); + // for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", pold.get(i)); + // } + printf("\nfinal val: %Lf\n", f_final); + printf("p_tol : %lf, iters : %"LI"d, batch_number : %"LI"d\n", + p_tol, iters, in); + + return; + // } + } + + pold.CopyValues(pnew); + f_old = f_new; + + // updating the step size as per the following + // \eta_i = \eta_{i-1} * max(0.5, 1 + mu * \eta_{i-1} * + // \nabla_\theta f_{i-1}'* + // \nabla_\theta f_i + // ) + scale_prev = scale; + scale = sqrt(la::Dot(grad, grad)); + // double temp_eta = 1 + mu * eta * (la::Dot(grad, prev_grad)) / + //(scale * scale_prev); + double temp_eta = 1 + mu * eta * (la::Dot(grad, prev_grad)); + eta = eta * ((0.5 > temp_eta)?0.5 : temp_eta); + } + + // permuting the data matrix + data_batched.Destruct(); + PermuteMatrix_(data(), &data_batched); + //printf("data permuted\n"); + } + + NOTIFY("Too many iterations in Stochastic Meta Descent\n"); + for(index_t i = 0; i < dim; i++) { + pt[i] = pold.get(i); + } + long double f_final = (*func_ptr_)(pold, data(), &grad); + // scale = sqrt(la::Dot(grad, grad)); + printf("\nfinal val: %Lf, p_tol : %lf, iters : %"LI"d\n", + f_final, p_tol, iters); + fx_format_result(opt_module_,"min_obtained","%Lf", f_final); + + return; + } + + void PermuteMatrix_(const Matrix& input, Matrix *output) { + + ArrayList perm_array; + index_t size = input.n_cols(); + Matrix perm_mat; + + perm_mat.Init(size, size); + perm_mat.SetAll(0.0); + + math::MakeRandomPermutation(size, &perm_array); + for(index_t i = 0; i < size; i++) { + perm_mat.set(perm_array[i], i, 1.0); + } + + la::MulInit(input, perm_mat, output); + return; + } +}; +*/ + +class SMD_SingleStep { + + private: + index_t dimension_; + Matrix data_; + // This is the original way of calling the function + long double (*func_ptr_)(Vector&, const Matrix&, Vector*); + + // But since the L2 function is pretty awesome, we need + // to make the calling a little diff + // long double (*func_ptr_stoc_)(Vector&, const Matrix&, Vector*, index_t); + datanode *opt_module_; + + public: + + SMD_SingleStep(){ + } + + ~SMD_SingleStep(){ + } + + void Init(long double (*fun)(Vector&, const Matrix&, Vector*), + //long double (*fun_stoc)(Vector&, const Matrix&, Vector*, index_t), + Matrix& data, datanode *opt_module){ + + data_.Copy(data); + + // function pointer to the original function + func_ptr_ = fun; + + // pointer to the broken up function used + // for stochastic optimization + // func_ptr_stoc_ = fun_stoc; + opt_module_ = opt_module; + dimension_ = fx_param_int_req(opt_module_, "param_space_dim"); + } + + const Matrix data() { + return data_; + } + + index_t dimension() { + return dimension_; + } + + void Eval(double *pt){ + + index_t iters; + index_t MAXIMUM_ITERATIONS = fx_param_int(opt_module_,"MAX_ITERS",100); + // double EPSILON = fx_param_double(opt_module_, "EPSILON", 1.0e-2); + double TOLERANCE = fx_param_double(opt_module_, "TOLERANCE", 1.0e-2); + index_t dim = fx_param_int_req(opt_module_, "param_space_dim"); + index_t num_batch = fx_param_int(opt_module_, "BATCHES",50); + Vector pold, pnew, eta, grad, prev_grad, one_vector; + long double f_old, f_new; + double scale, scale_prev, gamma; + double mu = fx_param_double(opt_module_, "MU", 0.1); + double p_tol = 0.0;//, f_tol = 0.0; + Matrix data_batched, matrix_for_had_prod; + index_t batch_size = fx_param_int(opt_module_,"BATCH_SIZE", + data().n_cols() / num_batch); + + // fx_clear_param(opt_module_,"BATCHES"); + num_batch = fx_param_int(opt_module_,"BATCHES", + data().n_cols() / batch_size); + // to decide how to chose starting value + // of alpha (right now it is just 1). + // also have to decide the value for the + // meta parameter mu (right now it is + // arbitrarily chosen as 0.1) + + pold.Init(dim); + pnew.Init(dim); + pold.CopyValues(pt); + grad.Init(dim); + prev_grad.Init(dim); + eta.Init(dim); + eta.SetAll(0.1); + one_vector.Init(dim); + one_vector.SetAll(1.0); + matrix_for_had_prod.Init(dim,dim); + data_batched.Copy(data()); + + // This is the original way of calling the function + // f_old = (*func_ptr_)(pold, data(), &grad); + + // But since the L2 function is pretty awesome, we need + // to make the calling a little diff + f_old = (*func_ptr_)(pold, data(), &grad); + // printf("first val: %Lf\n", f_old); + // scale = sqrt(la::Dot(grad, grad)); + + // Here we are doing the gradient step + // scale = || - \nabla_\theta f(X_t, \theta_t) || + // \theta_{t+1} = \theta_t - + // \eta_t * \nabla_\theta f(X_t,\theta_t) / scale; + + for (iters = 0; iters < MAXIMUM_ITERATIONS; iters++) { + + // Now going through the data batchwise + for (index_t in = 0; in < num_batch; in++) { + + // instead of scaling the gradient, how about using low values + // of the step sizes, because scaling the gradients result + // in the gradient being significant even when it is close + // to the optimal + // gamma = - eta / scale; + gamma = -1.0; + matrix_for_had_prod.SetDiagonal(eta); + pnew.SetZero(); + la::AddTo(pold, &pnew); + la::MulExpert(gamma, grad, matrix_for_had_prod, 1.0, &pnew); + //printf("\n pold:"); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", pold.get(i)); + //} + //printf("\n prev grad:"); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", grad.get(i)); + //} + //printf("\n pnew:"); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", pnew.get(i)); + //} + + + // using a batch + Matrix single_batch; + index_t st_pt = in * batch_size; + data_batched.MakeColumnSlice(st_pt, batch_size, &single_batch); + prev_grad.CopyValues(grad); + // This is the original way of calling the function + f_new = (*func_ptr_)(pnew, single_batch, &grad); + + // But since the L2 function is pretty awesome, we need + // to make the calling a little diff + // f_new = (*func_ptr_stoc_)(pnew, single_batch, &grad, num_batch); + + // Terminating conditions + // |f_t+1 - f_t| < epsilon & ||\theta_t+1 - \theta_t|| < delta + // but since we are using stochastic method, it is + // better to only consider \theta instead of f + Vector diff; + la::SubInit(pnew, pold, &diff); + // f_tol = fabs(f_new - f_old); + p_tol = sqrt(la::Dot(diff, diff)); + + // using just the point in the param_space + // which refuses to move + if (p_tol < TOLERANCE) { + // rejected because stops too early, need the check + // the overall gradient is small + Vector temp_grad; + temp_grad.Init(dim); + long double f_final = (*func_ptr_)(pold, data(), &temp_grad); + // but maybe we can skip that now + // double temp_grad_val = sqrt(la::Dot(temp_grad, temp_grad)); + // if (temp_grad_val < EPSILON) { + + fx_format_result(opt_module_, "iters", "%d", iters+1); + fx_format_result(opt_module_,"min_obtained","%Lf", f_final); + for (index_t i = 0; i < dim; i++) { + pt[i] = pold.get(i); + } + // printf("iters: %"LI"d\n", iters); + // for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", pold.get(i)); + // } + printf("\nfinal val: %Lf\n", f_final); + printf("p_tol : %lf, iters : %"LI"d, batch_number : %"LI"d\n", + p_tol, iters, in); + + return; + // } + } + + pold.CopyValues(pnew); + f_old = f_new; + + // updating the step size as per the following + // \eta_i = \eta_{i-1} * max(0.5, 1 + mu * \eta_{i-1} * + // \nabla_\theta f_{i-1}'* + // \nabla_\theta f_i + // ) + scale_prev = scale; + scale = sqrt(la::Dot(grad, grad)); + // double temp_eta = 1 + mu * eta * (la::Dot(grad, prev_grad)) / + //(scale * scale_prev); + //printf("\n grad:"); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", grad.get(i)); + //} + //printf("\n eta_old:"); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", eta.get(i)); + //} + Vector prev_grad_prod_eta; + la::MulInit(prev_grad, matrix_for_had_prod, &prev_grad_prod_eta); + matrix_for_had_prod.SetZero(); + matrix_for_had_prod.SetDiagonal(prev_grad_prod_eta); + Vector grad_prod_prev_grad_prod_eta; + la::MulInit(grad, matrix_for_had_prod, + &grad_prod_prev_grad_prod_eta); + Vector temp_eta; + temp_eta.Copy(one_vector); + la::AddExpert(mu, grad_prod_prev_grad_prod_eta, &temp_eta); + //printf("\n 1+mu*.. : "); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf ",temp_eta.get(i)); + //} + //printf("\n max(0.5, 1+mu*..) : "); + for (index_t i = 0; i < dim; i++) { + if ((temp_eta.ptr())[i] < 0.5) { + (temp_eta.ptr())[i] = 0.5; + } + } + //for(index_t i = 0; i < dim; i++) { + // printf("%lf ", temp_eta.get(i)); + //} + matrix_for_had_prod.SetZero(); + matrix_for_had_prod.SetDiagonal(temp_eta); + la::MulOverwrite(matrix_for_had_prod, eta, &temp_eta); + eta.CopyValues(temp_eta); + //printf("\n eta_new:"); + //for(index_t i = 0; i < dim; i++) { + // printf("%lf, ", eta.get(i)); + //} + //printf("\n"); + } + + // permuting the data matrix + data_batched.Destruct(); + PermuteMatrix_(data(), &data_batched); + //printf("data permuted\n"); + } + fflush(NULL); + NOTIFY("Too many iterations in Stochastic Meta Descent"); + for(index_t i = 0; i < dim; i++) { + pt[i] = pold.get(i); + } + long double f_final = (*func_ptr_)(pold, data(), &grad); + // scale = sqrt(la::Dot(grad, grad)); + printf("\nfinal val: %Lf, p_tol : %lf, iters : %"LI"d\n", + f_final, p_tol, iters); + fx_format_result(opt_module_,"min_obtained","%Lf", f_final); + + return; + } + + void PermuteMatrix_(const Matrix& input, Matrix *output) { + + ArrayList perm_array; + index_t size = input.n_cols(); + Matrix perm_mat; + + perm_mat.Init(size, size); + perm_mat.SetAll(0.0); + + math::MakeRandomPermutation(size, &perm_array); + for(index_t i = 0; i < size; i++) { + perm_mat.set(perm_array[i], i, 1.0); + } + + la::MulInit(input, perm_mat, output); + return; + } +}; + + +#endif diff --git a/fastlib2/contrib/pram/opt/phi.h b/fastlib2/contrib/pram/opt/phi.h new file mode 100644 index 0000000000..bdddb8720f --- /dev/null +++ b/fastlib2/contrib/pram/opt/phi.h @@ -0,0 +1,156 @@ +/** + * @author Parikshit Ram (pram@cc.gatech.edu) + * @file phi.h + * + * This file computes the Gaussian probability + * density function + */ +#include "fastlib/fastlib.h" +#include "fastlib/fastlib_int.h" +#include + +/** + * Calculates the multivariate Gaussian probability density function + * + * Example use: + * @code + * Vector x, mean; + * Matrix cov; + * .... + * long double f = phi(x, mean, cov); + * @endcode + */ + +long double phi(Vector& x , Vector& mean , Matrix& cov) { + + long double det, f; + double exponent; + index_t dim; + Matrix inv; + Vector diff, tmp; + + dim = x.length(); + la::InverseInit(cov, &inv); + det = la::Determinant(cov); + + if( det < 0){ + det = -det; + } + la::SubInit(mean,x,&diff); + la::MulInit(inv, diff, &tmp); + exponent = la::Dot(diff, tmp); + long double tmp1, tmp2, tmp3; + tmp1 = 1; + tmp2 = dim; + tmp2 = tmp2/2; + tmp2 = pow((2*(math::PI)),tmp2); + tmp1 = tmp1/tmp2; + tmp3 = 1; + tmp2 = sqrt(det); + tmp3 = tmp3/tmp2; + tmp2 = -exponent; + tmp2 = tmp2 / 2; + f = (tmp1*tmp3*exp(tmp2)); + + return f; +} + +/** + * Calculates the univariate Gaussian probability density function + * + * Example use: + * @code + * double x, mean, var; + * .... + * long double f = phi(x, mean, var); + * @endcode + */ + +long double phi(double x, double mean, double var) { + + long double f; + + f = exp(-1.0*((x-mean)*(x-mean)/(2*var)))/sqrt(2*math::PI*var); + return f; +} + +/** + * Calculates the multivariate Gaussian probability density function + * and also the gradients with respect to the mean and the variance + * + * Example use: + * @code + * Vector x, mean, g_mean, g_cov; + * ArrayList d_cov; // the dSigma + * .... + * long double f = phi(x, mean, cov, d_cov, &g_mean, &g_cov); + * @endcode + */ + +long double phi(Vector& x, Vector& mean, Matrix& cov, ArrayList& d_cov, Vector *g_mean, Vector *g_cov){ + + long double det, f; + double exponent; + index_t dim; + Matrix inv; + Vector diff, tmp; + + dim = x.length(); + la::InverseInit(cov, &inv); + det = la::Determinant(cov); + + if( det < 0){ + det = -det; + } + la::SubInit(mean,x,&diff); + la::MulInit(inv, diff, &tmp); + exponent = la::Dot(diff, tmp); + long double tmp1, tmp2, tmp3; + tmp1 = 1; + tmp2 = dim; + tmp2 = tmp2/2; + tmp2 = pow((2*(math::PI)),tmp2); + tmp1 = tmp1/tmp2; + tmp3 = 1; + tmp2 = sqrt(det); + tmp3 = tmp3/tmp2; + tmp2 = -exponent; + tmp2 = tmp2 / 2; + f = (tmp1*tmp3*exp(tmp2)); + + // Calculating the g_mean values which would be a (1 X dim) vector + la::ScaleInit(f,tmp,g_mean); + + // Calculating the g_cov values which would be a (1 X (dim*(dim+1)/2)) vector + double *g_cov_tmp; + g_cov_tmp = (double*)malloc(d_cov.size()*sizeof(double)); + for(index_t i = 0; i < d_cov.size(); i++){ + Vector tmp_d; + // Matrix inv_d; + Matrix inv_d, tmp_mat_1, tmp_mat_2; + double tmp_d_cov_d_r; + + la::MulInit(d_cov[i],inv,&tmp_mat_1); + la::MulInit(inv, tmp_mat_1, &tmp_mat_2); + la::MulInit(tmp_mat_2, diff, &tmp_d); + tmp_d_cov_d_r = la::Dot(diff, tmp_d); + + // la::MulInit(d_cov[i], tmp, &tmp_d); + // tmp_d_cov_d_r = la::Dot(tmp_d,tmp); + + la::MulInit(inv,d_cov[i],&inv_d); + + double trace = 0; + for(index_t j = 0; j < dim; j++) { + trace += inv_d.get(j,j); + } + + tmp_d_cov_d_r -= trace; + //printf("trace = %lf\n", trace); + + g_cov_tmp[i] = f*tmp_d_cov_d_r/2; + } + g_cov->Copy(g_cov_tmp,d_cov.size()); + + return f; +}