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

Copyright (C) 2018

A. Ronald Gallant
Post Office Box 659
Raleigh 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.

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

#ifndef __FILE_LPEQNS_H_SEEN__
#define __FILE_LPEQNS_H_SEEN__

#include "libscl.h"
#include "lprior.h"

class lpeqns : public scl::nleqns_base {
private:
  INTEGER M;
  INTEGER L;
  INTEGER mrs_pos;
  INTEGER cf_pos;
  INTEGER years;
  REAL lambda;
  scl::realmat target;
  scl::intvec idx;
  scl::realmat stats;
public:
  lpeqns(INTEGER M_, INTEGER L_, INTEGER mp, INTEGER cp, INTEGER yr, 
  REAL lam, const scl::realmat& tar, const scl::intvec& dx) 
  : M(M_),L(L_),mrs_pos(mp),cf_pos(cp),years(yr),lambda(lam), 
    target(tar), idx(dx) {}
  bool get_f(const scl::realmat& x, scl::realmat& f)
  {
    REAL ss = ssd(M,L,mrs_pos,cf_pos,years,x,lambda,target,idx,stats);
    f.resize(1,1);
    f[1] = ss;
    if (!scl::IsFinite(ss)) return false;
    return true;
  }
  bool get_F(const scl::realmat& x, scl::realmat& f, scl::realmat& F)
  {
    if (get_f(x,f)) {
      return df(x,F);
    }
    else {
      return false;
    }
  }
  bool df(scl::realmat x, scl::realmat& F)
  {
    INTEGER d = x.nrow();
    scl::realmat f0;
    scl::realmat f1;
    if (!get_f(x,f0)) return false;
    INTEGER r = f0.nrow();
    F.resize(r,d);
    REAL eps = std::pow(double(REAL_EPSILON),0.33333333);
    for (INTEGER j=1; j<=d; j++) {
      REAL tmp = x[j];
      REAL h = eps*std::fabs(tmp);
      if (h == 0) h = eps; 
      REAL hi = tmp + h;
      REAL lo = tmp - h;
      x[j] = hi;
      if (!get_f(x,f1) ) return false;
      x[j] = lo;
      if (!get_f(x,f0) ) return false;
      REAL difference = hi - lo;
      x[j] = tmp;
      for (INTEGER i=1; i<=r; i++) {
        F(i,j) = (f1[i] - f0[i])/difference;
      }
    }
    return true;
  }
};
  
  
#endif
