/*-----------------------------------------------------------------------------

Copyright (C) 2005, 2006, 2009, 2010.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA   

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/
#include "pathname.h"

#include "mpi.h"

#include "snp.h"
#include "snpusr.h"

using namespace scl;
using namespace libsnp;
using namespace snp;
using namespace std;

namespace {

  int my_rank;
  int no_procs;

  void mpi_error(string msg) 
  {
    cerr << msg << endl; MPI_Abort(MPI_COMM_WORLD, my_rank);
  }
  void mpi_warn(string msg) 
  {
    cerr << msg << endl;
  }

  const INTEGER bufsz = 512;

  snpll nlopt_ll;

  class objective_function : public nleqns_base {
  private:
    INTEGER count;
  public:
    objective_function() : count(0) { }
    void reset_count() { count = 0; }
    INTEGER get_count() { return count; }
    bool get_f(const realmat& rho, realmat& f)
    {
      if (f.nrow()!=1) f.resize(1,1);
      nlopt_ll.set_rho(rho);
      REAL log_like = nlopt_ll.log_likehood();
      f[1] = -log_like/nlopt_ll.get_datparms().n;
      ++count;
      return true;
    }
    bool get_F(const realmat& rho, realmat& f, realmat& F)
    {
      if (f.nrow()!=1) f.resize(1,1);
      nlopt_ll.set_rho(rho);
      realmat dllwrho;
      REAL log_like = nlopt_ll.log_likehood(dllwrho);
      f[1] = -log_like/nlopt_ll.get_datparms().n;
      F = -dllwrho/nlopt_ll.get_datparms().n;
      ++count;
      return true;
    }
  };
}

int main(int argc, char** argp, char** envp)
{
  MPI_Init(&argc, &argp);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &no_procs);

  LIB_ERROR_HANDLER_PTR previous_error = set_lib_error_handler(&mpi_error);
  LIB_WARN_HANDLER_PTR previous_warn = set_lib_warn_handler(&mpi_warn);

  string pathname = "/tmp/";

  if (my_rank == 0) {
    pathname = string(PATHNAME) + string("/");
  }

  ostream* detail_ptr = &cout;
  if (my_rank == 0) {
    string filename = pathname + "detail.dat";
    detail_ptr = new(nothrow) ofstream(filename.c_str());
    if ( (detail_ptr==0) || (!*detail_ptr) )
      error("Error, snp, " + filename + " open failed");
  }  
  ostream& detail = *detail_ptr;

  int  dim[4];          // dim[0] size, dim[1] rows, dim[2] cols, dim[3] extra

  dim[0] = dim[1] = dim[2] = dim[3] = 0;

  int dest;
  int src;
  int tag;
  int dum = 0;

  MPI_Status status;
  
  if (my_rank == 0) {

    string filename = pathname + "control.dat";
    ifstream ctrl_ifs(filename.c_str());
    if(!ctrl_ifs) error("Error, snp, control.dat open failed");
    control ctrl(ctrl_ifs);
  
    filename = pathname + "summary.dat";
    #if defined COMPILER_HAS_BASE_IOS 
      ofstream summary_ofs(filename.c_str(), ios_base::app);
    #else
      ofstream summary_ofs(filename.c_str(), ios::app);
    #endif
    if(!summary_ofs) error("Error, snp, summary.dat open failed");
    sumry summary(summary_ofs);
    summary.write_header();
  
    INTEGER active_procs = no_procs;
    while (active_procs > 1) { 

      MPI_Recv(&dum,1,MPI_INT,MPI_ANY_SOURCE,MPI_ANY_TAG,
        MPI_COMM_WORLD,&status);
      dest = src = status.MPI_SOURCE; tag = status.MPI_TAG;
  
      if (tag == 20) { //print optimization results

        char control_line[bufsz];
        MPI_Recv(&control_line,bufsz,MPI_CHAR,src,tag,MPI_COMM_WORLD,&status);
 
        string line = control_line;
        cutstr(line);
        string output_file = cutstr(line);

        vector<string> pfvec; // Parmfile as a vector of strings
        vecstrbuf pfbuf;      // Parmfile as a char array
        parmfile pf;          // Parmfile as an object

        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        pfbuf.resize(dim[1],dim[2]);
        MPI_Recv(pfbuf.get_ptr(),dim[0],MPI_CHAR,src,tag,
          MPI_COMM_WORLD,&status);
        pfvec = pfbuf.get_vec();

        if (!pf.set_parms(pfvec, detail)) {
           detail.flush();
           error("Error, snp, cannot interpret parmfile");
        }
  
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        INTEGER iter_count = dim[0];
        INTEGER termination_code = dim[1];
  
        realmat rho_stop;
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        rho_stop.resize(dim[1],dim[2]);
        MPI_Recv (rho_stop.get_x(),dim[0],MPI_DOUBLE,src,tag,
          MPI_COMM_WORLD,&status); 
  
        vector<string> detail_vec;
        vecstrbuf detail_buf;
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        detail_buf.resize(dim[1],dim[2]);
        MPI_Recv(detail_buf.get_ptr(),dim[0],MPI_CHAR,src,tag,
          MPI_COMM_WORLD,&status);
        detail_vec = detail_buf.get_vec();

        realmat X,Y;
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        X.resize(dim[1],dim[2]);
        MPI_Recv(X.get_x(),dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        Y.resize(dim[1],dim[2]);
        MPI_Recv(Y.get_x(),dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
  
        if (pf.get_optparms().print) {
          vector<string>::const_iterator itr;
          for (itr=detail_vec.begin(); itr!=detail_vec.end(); ++itr) {
            detail << *itr;
          }
        }
  
        snpll ll(pf.get_datparms(),
          pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(),
          pf.get_rfunc(), pf.get_afunc_mask(), pf.get_ufunc_mask(),
          pf.get_rfunc_mask());
 
        ll.set_XY(&X,&Y);
        summary.set_ll_ptr(&ll);  //Allows writing a full summary line

        nlopt_ll = ll;
        objective_function obj;
  
        realmat sn;
        obj.get_f(rho_stop,sn);

        realmat rho = rho_stop;

        ll.set_rho(rho);
        pf.set_afunc(ll.get_afunc());
        pf.set_ufunc(ll.get_ufunc());
        pf.set_rfunc(ll.get_rfunc());

        summary.set_filename(output_file.c_str());
        summary.set_pname(pf.get_optparms().pname);
        summary.set_n(pf.get_datparms().n);
    
        summary.set_sn(sn[1]);
        summary.set_iter(iter_count);
        summary.set_termination_code(termination_code);
        summary.write_line();

        if (pf.get_optparms().print) {
          detail << starbox("/Estimated value of rho//");
          detail << rho;
          detail << starbox("/Stats and diagnostics//");
          detail << '\n';
          ll.write_stats(detail);
        }
                                  
        string filename = pathname + output_file;
        if (!pf.write_parms(filename.c_str(),string(control_line),ll)){
          detail.flush();
          error("Error, snp, cannot write parmfile"); 
        }
      }

      else if (ctrl.read_line()) {

        vector<string> pfvec; // Parmfile as a vector of strings
        vecstrbuf pfbuf;      // Parmfile as a char array
        parmfile pf;          // Parmfile as an object

        string filename = pathname + ctrl.get_input_file();
  
        ifstream pf_ifs(filename.c_str());
        if (!pf_ifs) error("Error, snp, cannot open " + filename);
  
        string line;
  
        while (getline(pf_ifs, line)) pfvec.push_back(line);
  
        if (!pf.set_parms(pfvec, detail)) {
           detail.flush();
           error("Error, snp, cannot interpret parmfile");
        }
  
        if (pf.get_optparms().print) {
          detail << starbox("/Control line for above parmfile//");
          detail << '\n' << ctrl.get_line() << endl;
        }
    
        summary.set_filename(ctrl.get_output_file());
        summary.set_pname(pf.get_optparms().pname);
        summary.set_n(pf.get_datparms().n);
    
        datparms dp = pf.get_datparms();
        if (dp.dsn[0] != '/') dp.dsn = pathname + dp.dsn;
        pf.set_datparms(dp);
        
        realmat Y;
        datread_type dr; dr.initialize(pf.get_datparms());
        if (!dr.read_data(Y)) error("Error, snp, read_data failed");
    
        if (pf.get_optparms().print) {
          detail << starbox("/First 12 observations//");
          detail << Y("",seq(1,12));
          detail << starbox("/Last 12 observations//");
          detail << Y("",seq(pf.get_datparms().n-11, pf.get_datparms().n));
        }
    
        trnfrm tr(pf.get_datparms(), pf.get_tranparms(), Y);
        pf.set_tranparms(tr.get_tranparms());
        
        if (pf.get_optparms().print) {
          detail << starbox("/Mean and variance of data//");
          if (pf.get_tranparms().diag) {
            detail << "(Variance has been diagonalized.)\n";
          }
          detail << pf.get_tranparms().mean << pf.get_tranparms().variance;
        }

        if (!pf.get_tranparms().useold) {
          keyword kw;
          string header;
          vector<string>::iterator kw_ptr;
          header=kw.set_keyword("TRANSFORM START VALUES FOR mean");
          kw_ptr = find_if(pfvec.begin(), pfvec.end(), kw);
          if (kw_ptr == pfvec.end()) {
            pfvec.push_back("TRANSFORM START VALUES FOR mean (optional)");
            realmat mean = pf.get_tranparms().mean;
            for (INTEGER i=1; i<=mean.size(); ++i) {
              pfvec.push_back(fmt('e',26,17,mean[i]).get_ostr());
            }
          }
          else {
            ++kw_ptr;
            realmat mean = pf.get_tranparms().mean;
            for (INTEGER i=1; i<=mean.size(); ++i) {
              *kw_ptr++ = fmt('e',26,17,mean[i]).get_ostr(); 
            }
          }
          header=kw.set_keyword("TRANSFORM START VALUES FOR variance");
          kw_ptr = find_if(pfvec.begin(), pfvec.end(), kw);
          if (kw_ptr == pfvec.end()) {
            pfvec.push_back("TRANSFORM START VALUES FOR variance (optional)");
            realmat variance = pf.get_tranparms().variance;
            for (INTEGER i=1; i<=variance.size(); ++i) {
              pfvec.push_back(fmt('e',26,17,variance[i]).get_ostr());
            }
          }
          else {
            ++kw_ptr;
            realmat variance = pf.get_tranparms().variance;
            for (INTEGER i=1; i<=variance.size(); ++i) {
              *kw_ptr++ = fmt('e',26,17,variance[i]).get_ostr(); 
            }
          }
        }

        vecstrbuf send_buf(pfvec);
        pfbuf = send_buf;      

        realmat X = Y;
        tr.normalize(Y);
        tr.normalize(X);
        
        if (pf.get_tranparms().squash == 1) {
          tr.spline(X);
        } else if (pf.get_tranparms().squash == 2) {
          tr.logistic(X);
        }
    
        if (pf.get_optparms().print) {
          detail << starbox("/First 12 normalized observations//");
          detail << Y("",seq(1,12));
          detail << starbox("/Last 12 normalized observations//");
          detail << Y("",seq(pf.get_datparms().n-11, pf.get_datparms().n));
          if (pf.get_tranparms().squash > 0) {
            detail << starbox("/First 12 transformed observations//");
            detail << X("",seq(1,12));
            detail << starbox("/Last 12 transformed observations//");
            detail << X("",seq(pf.get_datparms().n-11, pf.get_datparms().n));
          }
        }
    
        if (pf.get_optparms().task == 0) { //fit

           tag = 10;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           char control_line[bufsz];
           string cls = ctrl.get_line();
           INTEGER limit = cls.size() < bufsz-1 ? cls.size() : bufsz-1;
           for (INTEGER i=0; i<limit; ++i) control_line[i] = cls[i];
           control_line[limit] = '\0';
           MPI_Send(&control_line,limit,MPI_CHAR,dest,tag,MPI_COMM_WORLD);

           dim[0] = pfbuf.size();
           dim[1] = pfbuf.get_rows();
           dim[2] = pfbuf.get_cols();
           MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
           MPI_Send(pfbuf.get_ptr(),dim[0],MPI_CHAR,dest,tag,MPI_COMM_WORLD);
  
           dim[0] = X.size();
           dim[1] = X.get_rows();
           dim[2] = X.get_cols();
           MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
           MPI_Send(X.get_x(),dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
           MPI_Send(Y.get_x(),dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
  
         }
         else if (pf.get_optparms().task == 1) { //res
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           residual_type rt
           (pf.get_optparms(),pf.get_datparms(),pf.get_tranparms(),
             pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(), pf.get_rfunc(), 
             pf.get_afunc_mask(), pf.get_ufunc_mask(), pf.get_rfunc_mask(),
             detail, &tr);
    
           string filename = pathname + ctrl.get_output_file();
           ofstream os(filename.c_str());
           if (!os) error("Error, snp, can't open " + filename);
    
           rt.initialize(&os);
           rt.initialize(ctrl.get_output_file());
           rt.set_XY(&X,&Y);
    
           if (!rt.calculate()) error("Error, snp, residual write failed");
    
           summary.write_partial_line();  // Information in ll not set.
         }
         else if (pf.get_optparms().task == 2) { //mu
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           mean_type mt(pf.get_optparms(),pf.get_datparms(),pf.get_tranparms(),
             pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(), pf.get_rfunc(), 
             pf.get_afunc_mask(), pf.get_ufunc_mask(), pf.get_rfunc_mask(),
             detail, &tr);
    
           string filename = pathname + ctrl.get_output_file();
           ofstream os(filename.c_str());
           if (!os) error("Error, snp, can't open " + filename);
    
           mt.initialize(&os);
           mt.initialize(ctrl.get_output_file());
           mt.set_XY(&X,&Y);
    
           if (!mt.calculate()) error("Error, snp, mean write failed");
    
           summary.write_partial_line();  // Information in ll not set.
         }
         else if (pf.get_optparms().task == 3) { //sig
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           variance_type vt
           (pf.get_optparms(),pf.get_datparms(),pf.get_tranparms(),
             pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(), pf.get_rfunc(), 
             pf.get_afunc_mask(), pf.get_ufunc_mask(), pf.get_rfunc_mask(),
             detail, &tr);
    
           string filename = pathname + ctrl.get_output_file();
           ofstream os(filename.c_str());
           if (!os) error("Error, snp, can't open " + filename);
    
           vt.initialize(&os);
           vt.initialize(ctrl.get_output_file());
           vt.set_XY(&X,&Y);
    
           if (!vt.calculate()) error("Error, snp, variance write failed");
    
           summary.write_partial_line();  // Information in ll not set.
         }
         else if (pf.get_optparms().task == 4) { //plt
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           plot_type pt(pf.get_optparms(),pf.get_datparms(),pf.get_tranparms(),
             pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(), pf.get_rfunc(), 
             pf.get_afunc_mask(), pf.get_ufunc_mask(), pf.get_rfunc_mask(),
             detail, &tr);
    
           string filename = pathname + ctrl.get_output_file();
           ofstream os(filename.c_str());
           if (!os) error("Error, snp, can't open " + filename);
    
           pt.initialize(&os);
           pt.initialize(ctrl.get_output_file());
           pt.set_XY(&X,&Y);
    
           if (!pt.calculate()) error("Error, snp, plot write failed");
    
           summary.write_partial_line();  // Information in ll not set.
         }
         else if (pf.get_optparms().task == 5) { //sim
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           simulate_type st
           (pf.get_optparms(),pf.get_datparms(),pf.get_tranparms(),
             pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(), pf.get_rfunc(), 
             pf.get_afunc_mask(), pf.get_ufunc_mask(), pf.get_rfunc_mask(),
             detail, &tr);
    
           string filename = pathname + ctrl.get_output_file();
           ofstream os(filename.c_str());
           if (!os) error("Error, snp, can't open " + filename);
    
           st.initialize(&os);
           st.initialize(ctrl.get_output_file());
           st.set_XY(&X,&Y);
    
           if (!st.calculate()) error("Error, snp, simulation write failed");
    
           summary.write_partial_line();  // Information in ll not set.
         }
         else if (pf.get_optparms().task == 6) { //usr
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           
           user_type ut(pf.get_optparms(),pf.get_datparms(),pf.get_tranparms(),
             pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(), pf.get_rfunc(), 
             pf.get_afunc_mask(), pf.get_ufunc_mask(), pf.get_rfunc_mask(),
             detail, &tr);
    
           string filename = pathname + ctrl.get_output_file();
           ofstream os(filename.c_str());
           if (!os) error("Error, snp, can't open " + filename);
    
           ut.initialize(&os);
           ut.initialize(ctrl.get_output_file());
           ut.set_XY(&X,&Y);
    
           if (!ut.calculate()) error("Error, snp, user write failed");
    
           summary.write_partial_line();  // Information in ll not set.
         }
         else { //error
           tag = 30;
           MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
           error("Error, snp, no such task");
         }
      }
      else { // no more control lines to read
        tag = 90;
        MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
        --active_procs;
      }
    }
  }  // end my_rank == 0

  else {  // my_rank != 0

    vector<string> pfvec; // Parmfile as a vector of strings
    vecstrbuf pfbuf;      // Parmfile as an char array
    parmfile pf;          // Parmfile as an object

    dest = 0; tag = 10; 
    MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
    MPI_Recv(&dum,1,MPI_INT,dest,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
    src = status.MPI_SOURCE; tag = status.MPI_TAG;

    while (tag != 90) {
    
      if (tag == 10) {
        
        char control_line[bufsz];
        MPI_Recv(&control_line,bufsz,MPI_CHAR,src,tag,MPI_COMM_WORLD,&status);

        string line = control_line;
        string input_file = cutstr(line);
        string output_file = cutstr(line);
        string word = cutstr(line);
        REAL fnew = atof(word.c_str());
        word = cutstr(line);
        REAL fold = atof(word.c_str());
        word = cutstr(line);
        INTEGER nstart = atoi(word.c_str());
        word = cutstr(line);
        INT_32BIT jseed = atoi(word.c_str());
 
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        pfbuf.resize(dim[1],dim[2]);
        MPI_Recv(pfbuf.get_ptr(),dim[0],MPI_CHAR,src,tag,
          MPI_COMM_WORLD,&status);
        pfvec = pfbuf.get_vec();

        keyword kw;  // Must suppress print if my_rank != 0
        vector<string>::iterator kw_ptr;
        string header=kw.set_keyword("OPTIMIZATION DESCRIPTION");
        kw_ptr = find_if(pfvec.begin(), pfvec.end(), kw);
        if (kw_ptr == pfvec.end() || kw_ptr + 6 > pfvec.end())
        error("Error, snp_mpi, " + header);
        kw_ptr += 6;
        bool print = atoi(kw_ptr->substr(0,10).c_str());
        *kw_ptr = "0 ";

        if (!pf.set_parms(pfvec, detail)) {
           detail.flush();
           error("Error, snp, cannot interpret parmfile");
        }
  
        realmat X,Y;
        MPI_Recv(&dim,4,MPI_INT,src,tag,MPI_COMM_WORLD,&status);
        X.resize(dim[1],dim[2]);
        MPI_Recv(X.get_x(),dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);
        Y.resize(dim[1],dim[2]);
        MPI_Recv(Y.get_x(),dim[0],MPI_DOUBLE,src,tag,MPI_COMM_WORLD,&status);

        snpll ll(pf.get_datparms(),
          pf.get_snpden(), pf.get_afunc(), pf.get_ufunc(),
          pf.get_rfunc(), pf.get_afunc_mask(), pf.get_ufunc_mask(),
          pf.get_rfunc_mask());

        ll.set_XY(&X,&Y);

        realmat rho = ll.get_rho();

        nlopt_ll = ll;
     
        objective_function obj;
     
        realmat sn;
        obj.get_f(rho,sn);
        
        nlopt minimizer(obj);
        minimizer.set_output(false);
        minimizer.set_warning_messages(false);
        minimizer.set_iter_limit(pf.get_optparms().itmax0);
        minimizer.set_solution_tolerance(pf.get_optparms().toler);
     
        realmat best_rho = rho;
     
        INTEGER best_iter = 0; 
        INTEGER evaluations = 0;
        
        vector<string> detail_vec;
     
        if (print) {
          string header_line = "\nFor output file ";
          header_line += output_file;
          header_line += "\n";
          detail_vec.push_back(header_line);
        }

        if (nstart > 0) {
          if (print) {
            detail_vec.push_back("\nStarting value of rho\n");
            for (INTEGER i=1; i<=rho.size(); ++i) {
              detail_vec.push_back(fmt('e',27,17,rho[i]).get_ostr()+ "\n");
            }
            detail_vec.push_back("\nTrials\n");
          }
          best_iter = pf.get_optparms().itmax0;
          REAL best_f = REAL_MAX;  
          for (INTEGER trial=1; trial<=nstart; ++trial) {
            ll.set_rho(rho, fold, fnew, jseed);
            realmat rho_start = ll.get_rho();
            realmat rho_stop;
            realmat sn;
            obj.reset_count();
            minimizer.minimize(rho_start, rho_stop);
            obj.get_f(rho_stop,sn);
     
            evaluations += obj.get_count();
     
            if (print) {
              detail_vec.push_back
              ("     trial =" + fmt('d',4,trial).get_ostr() 
              +",  obj =" + fmt('e',23,16,sn[1]).get_ostr() 
              +",  eval = " + fmt('d',4,obj.get_count()).get_ostr()
              +",  iter = " + fmt('d',4,minimizer.get_iter_count()).get_ostr()
              +"\n");
            }
            if (sn[1] < best_f) { 
              best_iter = minimizer.get_iter_count();
              best_f = sn[1];
              best_rho = rho_stop;
            }
          }
        }
     
        realmat rho_start = best_rho;
        realmat rho_stop = rho_start;
        
        minimizer.set_iter_limit(pf.get_optparms().itmax1);
     
        obj.reset_count();
     
        if (pf.get_optparms().itmax1 > 0) {
          minimizer.minimize(rho_start, rho_stop);
          if (print) {
            detail_vec.push_back("\nFinal\n");
            realmat sn;
            obj.get_f(rho_stop,sn);
            detail_vec.push_back
            (string("                ")
            +"   obj =" + fmt('e',23,16,sn[1]).get_ostr() 
            +",  eval = " + fmt('d',4,obj.get_count()).get_ostr()
            +",  iter = " + fmt('d',4,minimizer.get_iter_count()).get_ostr()
            +"\n");
          }
        }
  
        dest = 0; tag = 20;
        MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);

        MPI_Send(&control_line,bufsz,MPI_CHAR,dest,tag,MPI_COMM_WORLD);

        dim[0] = pfbuf.size();
        dim[1] = pfbuf.get_rows();
        dim[2] = pfbuf.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        MPI_Send(pfbuf.get_ptr(),dim[0],MPI_CHAR,dest,tag,MPI_COMM_WORLD);
  
        dim[0] = best_iter + minimizer.get_iter_count();
        dim[1] = minimizer.get_termination_code();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
  
        dim[0] = rho_stop.size();
        dim[1] = rho_stop.get_rows();
        dim[2] = rho_stop.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        MPI_Send(rho_stop.get_x(),dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
        
        vecstrbuf send_buf(detail_vec);
        dim[0] = send_buf.size();
        dim[1] = send_buf.get_rows();
        dim[2] = send_buf.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        MPI_Send(send_buf.get_ptr(),dim[0],MPI_CHAR,dest,tag,MPI_COMM_WORLD); 

        dim[0] = X.size();
        dim[1] = X.get_rows();
        dim[2] = X.get_cols();
        MPI_Send(&dim,4,MPI_INT,dest,tag,MPI_COMM_WORLD);
        MPI_Send(X.get_x(),dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
        MPI_Send(Y.get_x(),dim[0],MPI_DOUBLE,dest,tag,MPI_COMM_WORLD);
      }
      
      dest = 0; tag = 10; 
      MPI_Send(&dum,1,MPI_INT,dest,tag,MPI_COMM_WORLD);
      MPI_Recv(&dum,1,MPI_INT,dest,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
      src = status.MPI_SOURCE; tag = status.MPI_TAG;
    }
  } //end my_rank != 0   

  if (my_rank == 0) delete detail_ptr;

  MPI_Finalize();

  previous_error = set_lib_error_handler(previous_error);
  previous_warn = set_lib_warn_handler(previous_warn);
    
  return 0;
}