diff --git a/fastlib/u/pram/mog_l2e/README.txt b/fastlib/u/pram/mog_l2e/README.txt
index 5867d1e629..ed3c1e933a 100644
--- a/fastlib/u/pram/mog_l2e/README.txt
+++ b/fastlib/u/pram/mog_l2e/README.txt
@@ -14,12 +14,9 @@ The files are the following:
 
 4. phi.h - this contains the functions that calculates the value of the multivariate Gaussian PDF, and also the gradients of the gaussian PDF with respect to the mean and variance when d(sigma) is provided.
 
-5. optimizers.h - this file contains the implementation of two optimizer
-classes
+5. build.py
 
-6. build.py
-
-7. The .arff and .csv file - the data files on which you can run the program.
+6. The .arff and .csv file - the data files on which you can run the program.
 
 -> An example run would be:
 fl-build mog_l2e_main
diff --git a/fastlib/u/pram/mog_l2e/build.py b/fastlib/u/pram/mog_l2e/build.py
index a46a43abc1..a0bc847cdf 100644
--- a/fastlib/u/pram/mog_l2e/build.py
+++ b/fastlib/u/pram/mog_l2e/build.py
@@ -10,7 +10,7 @@ librule(
 binrule(
     name = "mog_l2e_main",                  # the executable name
     sources = ["mog_l2e_main.cc"],          # compile main.cc
-    headers = ["optimizers.h"], # no extra headers
+    headers = ["../opt/optimizers.h"], # no extra headers
     deplibs = [":mog_l2e","fastlib:fastlib"]           #
     )
 
diff --git a/fastlib/u/pram/mog_l2e/mog_l2e_main.cc b/fastlib/u/pram/mog_l2e/mog_l2e_main.cc
index 7599fc69c4..466f4767c9 100644
--- a/fastlib/u/pram/mog_l2e/mog_l2e_main.cc
+++ b/fastlib/u/pram/mog_l2e/mog_l2e_main.cc
@@ -22,7 +22,7 @@
  */
 
 #include "mog.h"
-#include "optimizers.h"
+#include "../opt/optimizers.h"
 
 
 int main(int argc, char* argv[]) {
diff --git a/fastlib/u/pram/mog_l2e/optimizers.h b/fastlib/u/pram/mog_l2e/optimizers.h
deleted file mode 100644
index 4ec569be8d..0000000000
--- a/fastlib/u/pram/mog_l2e/optimizers.h
+++ /dev/null
@@ -1,474 +0,0 @@
-/**
- * @author Parikshit Ram (pram@cc.gatech.edu)
- * @file optimizers.h
- *
- * Implements classes for two types of optimizer
- *
- */
-
-#ifndef OPTIMIZERS_H
-#define OPTIMIZERS_H
-
-#include <fastlib/fastlib.h>
-
-/**
- * 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, &number_of_function_evaluations);
- * // 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-5);
-    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;
-	}
-	break;
-      }
-      if(num_func_eval > NMAX){
-	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, &number_of_function_evaluations);
- * // 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",200);
-    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);
-	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);
-	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");
-  }
-
-  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);
-    }
-  }
-};
-
-#endif