//
//-----------------------------------------------------------------------
/// @copyright
///            (c) Copyright 2008 by GATS, Inc.,
///     11864 Canon Blvd, Suite 101, Newport News VA 23606
///
///  All Rights Reserved. No part of this software or publication may be
///  reproduced, stored in a retrieval system, or transmitted, in any form
///  or by any means, electronic, mechanical, photocopying, recording, or
///  otherwise without the prior written permission of GATS, Inc.
///
//-----------------------------------------------------------------------
///
/// @file	EdgeModel.cpp
///
/// @author	John Burton	<j.c.burton@gats-inc.com>
///
/// @date  	Thu May  8 14:07:43 2008
///
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//

#include <string>
#include <valarray>
#include "EdgeModel.h"
#include "ConfigFile.h"
#include "Event.h"
#include "XAxisFunctor.h"
#include "EdgeLocation.h"
#include "GATS_Utilities.hpp"
#include "lmmin.h"

//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//

typedef struct {
  double *y;
  double *x;
  double R2;
  double (*func) (double *par, double user_x);
} lm_data_type;

static int ModelType = 0;
static int AttenuationMethod = 0;
EventVar<double> AverageIntensity_;
EventVar<double> ExtentNormalizedGrid_;
EventVar<double> ImpactAngle_;
EventVar<double> Transmission_;
double TopEdgeImpactAngle_;
double BottomEdgeImpactAngle_;


//
//-----------------------------------------------------------------------
// Utility Functions:
//-----------------------------------------------------------------------
//


static double Boltzmann(double &x, double &a0, double &a1, double &a2, double &a3)
{
  double y;

  y = a0 + a1 / (1 + exp((x - a2)/a3));

  return y;
}


static double TopEdge(double &x)
{
  double y;
  double a0 = 0.464446;
  double a1 = -0.4144044;
//  double a2 = 0.00534758;
  double a2 = 0.0;         // shifted so max positive slope is at 0.0
  double a3 = 0.00404353;
 
  y = Boltzmann(x,a0,a1,a2,a3);

  return y;
}

static double BottomEdge(double &x)
{
  double y;
  double a0 = -0.000650605;
  double a1 = 0.596465605;
//  double a2 = 0.991427;
  double a2 = 1.0;         // shifted so max negative slope is at 1.0;
  double a3 = 0.00716971;

  y = Boltzmann(x,a0,a1,a2,a3);

  return y;
}

static double Center(double &x)
{
  double y;
  double a0 = 0.532345;
  double a1 = 3.50656;
  double a2 = -10.4222;
  double a3 = 14.2721;
  double a4 = -7.51987;

  y = a0 + x * (a1 + x * (a2 + x * (a3 + x * (a4))));

  return y;
}


static double TopCenter(double &x)
{
  double y;
  double a0 = 0.197303;
  double a1 = 17.8989;
  double a2 = -161.485;

  y = a0 + x * (a1 + x * (a2));

  return y;
}



static double BottomCenter(double &x)
{
  double y;
  double a0 = 136.795;
  double a1 = -504.052;
  double a2 = 620.074;
  double a3 = -252.696;

  y = a0 + x * (a1 + x * (a2 + x * (a3)));

  return y;
}


static double Model(double &x)
{
  double y;
  
  if(ModelType == 1)
    y =  AverageIntensity_.Interpol(ExtentNormalizedGrid_,x);
  else
  {
    if(x <= 0.01)
      y = TopEdge(x);
    else if(x < 0.04)
      y = TopCenter(x);
    else if(x < 0.75)
      y = Center(x);
    else if(x < 0.98)
      y = BottomCenter(x);
    else
      y = BottomEdge(x);
  }

  return y;
}


static double lmFitEdge(double *par, double loc)
{
  double x,y;
  double a0,a1,a2,a3;
  double attn_factor;
  
  attn_factor = -1;

  x = par[0] +  loc * par[1];

  y = Model(x);

  if(AttenuationMethod == 1)
  {
    double m = (BottomEdgeImpactAngle_ - TopEdgeImpactAngle_);
    double b = TopEdgeImpactAngle_;
//    double ia = loc * m + b;
    double ia = x * m + b;
    attn_factor = Transmission_.Interpol(ImpactAngle_,ia);
    y = attn_factor * y;
  }
  else
  {
    if(x <= 0.01)
      y = par[2] * y;
    else if(x < 0.04)
      y = par[3] * y;
    else if(x < 0.75)
      y = par[4] * y;
    else if(x < 0.98)
      y = par[5] * y;
    else
      y = par[6] * y;
  }
  return y;
}

static void lmEvaluate(double *par, int m_dat, double *fvec,
		       void *data, int *info)
{
  int i;
  lm_data_type *mydata;

  mydata = (lm_data_type *)data;

  for(i=0; i<m_dat; i++)
    fvec[i] = mydata->y[i] - mydata->func(par,mydata->x[i]);
}



static void lmDegreeOfFit(int n_par, double *par, int m_dat, double *fvec,
			  void *data, int iflag, int iter, int nfev)
{
  double f, y, t;
  int i;
  double ss_reg = 0.0;
  double ss_tot = 0.0;
  double tmp_reg;
  double tmp_tot;
  double sum = 0.0;
  double ave = 0.0;
  lm_data_type *mydata;
  mydata = (lm_data_type *) data;
  
  if (iflag == -1) 
  {
    sum = 0.0;
    for(i=0; i<m_dat; i++)
      sum = sum + (mydata->y)[i];
    ave = sum / (double)m_dat;
    ss_tot = 0.0;
    ss_reg = 0.0;
    for (i = 0; i < m_dat; ++i) 
    {
      t = (mydata->x)[i];
      y = (mydata->y)[i];
      f = mydata->func(par,t);
      tmp_tot = ave - y;
      tmp_reg = f   - y;
      ss_tot = ss_tot + (tmp_tot * tmp_tot);
      ss_reg = ss_reg + (tmp_reg * tmp_reg);
    }
    mydata->R2 = 1.0 - (ss_reg/ss_tot);
  }
}

//
//-----------------------------------------------------------------------
// Class Methods:
//-----------------------------------------------------------------------
//


EdgeModel::EdgeModel(ConfigFile &cf, std::string section)
{
  std::string PixName33  = "PixName33";
  std::string IntensityName = "EdgeModelInt";
  std::string PositionName  = "EdgeModelPos";
  std::string fname;

  cf_      = cf;
  section_ = section;

  ModelType = 1;

  if(!(cf_.ValidEntry(section_,IntensityName)))
  {
    std::string errstr = "SolarEdgeModel::SolarEdgeModel: ConfigFile Entry not found: ";
    errstr = errstr + IntensityName;
    throw std::runtime_error(errstr);
  }
  fname  = cf_.GetStr(section_,IntensityName);
  AverageIntensity_.readVar(fname);

  if(!(cf_.ValidEntry(section_,PositionName)))
  {
    std::string errstr = "SolarEdgeModel::SolarEdgeModel: ConfigFile Entry not found: ";
    errstr = errstr + PositionName;
    throw std::runtime_error(errstr);
  }
  fname  = cf_.GetStr(section_,PositionName);
  ExtentNormalizedGrid_.readVar(fname);
  if(AverageIntensity_.size() != ExtentNormalizedGrid_.size())
  {
    std::string errstr = "SolarEdgeModel::SolarEdgeModel: Size of EventVars does not match";
    throw std::runtime_error(errstr);
  }
}




EdgeModel::EdgeModel(EventVar<double> &sample, EventVar<double> &location, 
                     double &topedge, double &bottomedge)
{
  int len;
  double extent;

  Sample_   = sample;      // normalized pixel intensity
  Location_ = location;    // FPA row location of pixel/sum

  ModelType = 0;

  len = Sample_.size();

  EdgeModel_.resize(len);
  Grid_.resize(len);

  extent = (bottomedge - topedge);
  Grid_  = (Location_ - topedge) / extent;

  if(Grid_.Mean() < 0.5)  // top edge model
  {
//    a0_ = 0.0500416;
//    a1_ = 0.464446;
    a0_ = 0.464446;
    a1_ = -0.4144044;
    a2_ = 0.00534758;
    a3_ = 0.00404353;
  }
  else                    // bottom edge model
  {
//    a0_ = 0.595815;
//    a1_ = -0.000650605;
    a0_ = -0.000650605;
    a1_ = 0.596465605;
    a2_ = 0.991427;
    a3_ = 0.00716971;
  }
}


void EdgeModel::setTransmissionModel(EventVar<double> &ImpactAngle, EventVar<double> &Transmission)
{
  ImpactAngle_  = ImpactAngle;
  Transmission_ = Transmission;
}


void EdgeModel::attenuateEdgeModel(void)
{
  AttenuationMethod = 0;
}


void EdgeModel::attenuateEdgeModel(double &topedge_ia, double &bottomedge_ia)
{

  AttenuationMethod      = 1;
  TopEdgeImpactAngle_    = topedge_ia;
  BottomEdgeImpactAngle_ = bottomedge_ia;

}



EventVar<double> EdgeModel::findEdge(EventVar<double> &sample, EventVar<double> &location, 
				     double &topedge, double &bottomedge)
{
  int len, m_dat, n_par;
  int i,j;
  double extent;

  Sample_   = sample;      // normalized pixel intensity
  Location_ = location;    // FPA row location of pixel/sum

  m_dat = Sample_.size();

  EdgeModel_.resize(m_dat);
  Grid_.resize(m_dat);

  extent = (bottomedge - topedge);
  Grid_  = (Location_ - topedge) / extent;

  n_par = 7;
  len = m_dat + n_par + 3;

  EventVar<double> tmp(len);

  double edge0, edge1;
  double loc0,loc1;
  double par[n_par];
  double val[m_dat];
  double loc[m_dat];

  par[0] = 0.0;
  par[1] = 1.0;
  par[2] = 1.0;
  par[3] = 1.0;
  par[4] = 1.0;
  par[5] = 1.0;
  par[6] = 1.0;

  lm_control_type control;
  lm_data_type    data;

  lm_initialize_control(&control);

  for(i=0; i<m_dat; i++)
  {
    val[i] = Sample_[i];
    loc[i] = Grid_[i];
  }

  data.func = lmFitEdge;
  data.y    = val;
  data.x    = loc;
  lm_minimize(m_dat,n_par,par,lmEvaluate,lmDegreeOfFit,&data,&control);
//  edge0 = par[0];
//  edge1 = par[0] + par[1];
  edge0 = -(par[0]/par[1]);
  edge1 = (1-par[0])/par[1];

  j = 0;
  loc0 = Location_.Interpol(Grid_,edge0);  // Top edge location in FPA space
  loc1 = Location_.Interpol(Grid_,edge1);  // Bottom edge location in FPA space
//  std::cerr << "EdgeModel: " << edge0 << ", " << edge1 ;
//  std::cerr << ", " << loc0 << ", " << loc1 << "\n";
  tmp[j++] = loc0;
  tmp[j++] = loc1;
  tmp[j++] = data.R2;                           // degreee of fit
  for(i=0; i<n_par; i++)
    tmp[j++] = par[i];                            // final a0 value;
  for(i=0; i<m_dat; i++)
    tmp[j++] = lmFitEdge(par,loc[i]);           // final a1 value;
  return tmp;
}



EventVar<double> EdgeModel::findEdge(void)
{
  int m_dat = Grid_.size();
  int n_par = 7;
  int len = m_dat + n_par + 3;
  int i,j;

  EventVar<double> tmp(len);

  double edge0, edge1;
  double loc0,loc1;
  double par[n_par];
  double val[m_dat];
  double loc[m_dat];

//  par[0] =(Grid_.Mean() < 0.5) ? 0.00534758 : 0.991427;
//  par[1] = 1.0;
//  par[2] = 0.0;
//  par[3] = 0.0;
  par[0] = 0.0;
  par[1] = 1.0;
  par[2] = 1.0;
  par[3] = 1.0;
  par[4] = 1.0;
  par[5] = 1.0;
  par[6] = 1.0;

  lm_control_type control;
  lm_data_type    data;

  lm_initialize_control(&control);

  for(i=0; i<m_dat; i++)
  {
    val[i] = Sample_[i];
    loc[i] = Grid_[i];
  }

  data.func = lmFitEdge;
  data.y    = val;
  data.x    = loc;
  lm_minimize(m_dat,n_par,par,lmEvaluate,lmDegreeOfFit,&data,&control);
  edge0 = par[0];
  edge1 = par[0] + par[1];

  j = 0;
  loc0 = Location_.Interpol(Grid_,edge0);  // Top edge location in FPA space
  loc1 = Location_.Interpol(Grid_,edge1);  // Bottom edge location in FPA space
  tmp[j++] = loc0;
  tmp[j++] = loc1;
  tmp[j++] = data.R2;                           // degreee of fit
  for(i=0; i<n_par; i++)
    tmp[j++] = par[i];                            // final a0 value;
  for(i=0; i<m_dat; i++)
    tmp[j++] = lmFitEdge(par,loc[i]);           // final a1 value;
  return tmp;
}





EventVar<double> EdgeModel::findEdge(double &a0, double &a1, double &a2, double &a3)
{
  int m_dat = Grid_.size();
  int n_par = 4;
  int len = m_dat + n_par + 2;
  int i,j;

  EventVar<double> tmp(len);

  double edge;
  double par[n_par];
  double val[m_dat];
  double loc[m_dat];

  par[0] =(Grid_.Mean() < 0.5) ? 0.00534758 : 0.991427;
  par[1] = 0.0;
  par[2] = 1.0;
  par[3] = 0.0;

  lm_control_type control;
  lm_data_type    data;

  lm_initialize_control(&control);

  for(i=0; i<m_dat; i++)
  {
    val[i] = Sample_[i];
    loc[i] = Grid_[i];
  }

  data.func = lmFitEdge;
  data.y    = val;
  data.x    = loc;
  lm_minimize(m_dat,n_par,par,lmEvaluate,lmDegreeOfFit,&data,&control);
  j = 0;
  tmp[j++] = Location_.Interpol(Grid_,par[0]);  // edge location in FPA space
  tmp[j++] = data.R2;                           // degreee of fit
  for(i=0; i<n_par; i++)
    tmp[j++] = par[i];                            // final a0 value;
  for(i=0; i<m_dat; i++)
    tmp[j++] = lmFitEdge(par,loc[i]);           // final a1 value;
  return tmp;
}