//
//-----------------------------------------------------------------------
/// @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	SolarEdgeModel.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 "SolarEdgeModel.h"
#include "ConfigFile.h"
#include "Event.h"
#include "XAxisFunctor.h"
#include "EdgeLocation.h"
#include "GATS_Utilities.hpp"
#include "lmRegression.h"

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



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

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


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

  cf_      = cf;
  section_ = section;

  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);
  }
}


void SolarEdgeModel::calcImpactAngleGrid(double &topedge, double &bottomedge)
{
  double m,b;

  m = (bottomedge - topedge);
  b = topedge;
  ImpactAngleAxis_.set(b,m);
  ImpactAngleGrid_.setValueZero();
  ImpactAngleGrid_ = ExtentNormalizedGrid_ * m + b;
}


void SolarEdgeModel::calcFPACoordinateGrid(double &topedge, double &bottomedge)
{
  double m,b;

  m = (bottomedge - topedge);
  b = topedge;
  FPACoordinateAxis_.set(b,m);
  FPACoordinateGrid_.setValueZero();
  FPACoordinateGrid_ = ExtentNormalizedGrid_ * m + b;
}


void SolarEdgeModel::attenuateEdgeModel(EventVar<double> &ImpactAngle, EventVar<double> &Transmission)
{
  int    nelem = ImpactAngleGrid_.size();
  int    i,ext;
  double ia,trans;
  
  ext = 1;
  EdgeModel_.resize(nelem);
  for(i=0; i<nelem; i++)
  {
    ia    = ImpactAngleGrid_[i];
    if(ia < 0.38)
      trans = 1.0;
    else
      trans = Transmission.Interpol(ImpactAngle,ia);
    EdgeModel_[i] = trans * AverageIntensity_[i];
  }
//  EdgeModel_.Plot("EdgeModel","");
//  EdgeModel_.Plot("EdgeModel","",ImpactAngleGrid_);
}

void SolarEdgeModel::plotEdgeModel(double &ts)
{
  std::string str = "time = " + GATS_Utilities::ConvertToString(ts);
  FPACoordinateGrid_.Plot("FPACoordinateGrid",str);
  EdgeModel_.Plot("EdgeModel",str);
  EdgeModel_.Plot("EdgeModel",str,FPACoordinateGrid_);
}



EventVarVect<double> SolarEdgeModel::findEdges(EventVar<double> &Xgrid, EventVar<double> &sample, double &lx, double &lxval, 
					   double &hx, double &hxval, double &te, double &be)
{
  EventVarVect<double> evv;
  EventVar<double> grid;
  EventVar<double> tmp(3);
  double edge;
  double r2;
  int i;
  int len = ExtentNormalizedGrid_.size();

  lmInitialize(EdgeModel_,ExtentNormalizedGrid_);
  FPACoordinateGrid_ = lmFitGrid(sample,lx,lxval,hx,hxval,r2);
//  std::cerr << "R2 = " << r2 << "\n";
  tmp[2] = r2;
//  i = 0;
//  while((i < len) && ExtentNormalizedGrid_[i] < 0)
//    i++;
//  tmp[0] = findLowXEdge(te,EdgeModel_,FPACoordinateGrid_);
  tmp[0] = FPACoordinateGrid_.Interpol(ExtentNormalizedGrid_,0.0);
//  std::cerr << "Edge:: " << ExtentNormalizedGrid_[i];
//  std::cerr << ", "     << FPACoordinateGrid_[i];
//  std::cerr << ", "     << edge;
//  std::cerr << ", "     << tmp[0] << "\n";
//  tmp[0] = edge;
//  i = len-1;
//  while((i > 0) && ExtentNormalizedGrid_[i] > 1)
//    i--;
//  tmp[1] = findHighXEdge(be,EdgeModel_,FPACoordinateGrid_);
  tmp[1] = FPACoordinateGrid_.Interpol(ExtentNormalizedGrid_,1.0);
//  std::cerr << "Edge:: " << ExtentNormalizedGrid_[i];
//  std::cerr << ", "     << FPACoordinateGrid_[i];
//  std::cerr << ", "     << edge;
//  std::cerr << ", "     << tmp[1] << "\n";
//  tmp[1] = edge;
  grid = EdgeModel_.VInterpol(FPACoordinateGrid_,Xgrid);
  evv.addEventVar(tmp);
  evv.addEventVar(grid);

  return evv;
//  std::cerr << "findEdges: testing lm_minimize \n";
//  lm_minimize(m_dat,n_par,lmEvaluate,lmPrint,&data,&control);
//  std::cerr << "findEdges:lm_minimize: status " << lm_shortmsg[control.info];
//  std::cerr << " after " << control.nfev << " evaluations\n";

}


//EventVar<double> SolarEdgeModel::doEdgeModel(EventVar<double> &sample, double &lx, double &lxval, 
//					   double &hx, double &hxval, double &te, double &be)
EventVar<double> SolarEdgeModel::doEdgeModel(EventVar<double> &Xgrid, EventVar<double> &sample,
					     double &lx, double &lxval, double &hx, double &hxval)
{
  int i;
  int len = ExtentNormalizedGrid_.size();
  int nrows = sample.size();
  double r2;
  EventVar<double> tmp(nrows);
//  EventVar<double> Xgrid(0.0,(double)(nrows-1),1.0);

  lmInitialize(EdgeModel_,ExtentNormalizedGrid_);
  FPACoordinateGrid_ = lmFitGrid(sample,lx,lxval,hx,hxval,r2);

//  for(i=0;i<nrows;i++)
//  {
//    tmp[i] = EdgeModel21_.Interpol(FPACoordinateGrid_,(double)i);
//  }
  tmp = EdgeModel_.VInterpol(FPACoordinateGrid_,Xgrid);
  return tmp;
}




EventVar<double> SolarEdgeModel::findIntensityCoeff(EventVar<double> &vals, EventVar<int> &locs)
{
  EventVar<double> tmp(2);
  int i;
//  int len = ExtentNormalizedGrid_.size();
  int len = locs.size(); 
  double m,b;
  double x,y;

  std::cerr << "findIntensityCoeff::len = " << len << "\n";
  lmInitialize(EdgeModel_,FPACoordinateGrid_);
//  tmp = lmFitIntensities(vals,locs);
  tmp = lmFitNoiseFloorAndBlur(vals,locs);

//  x = (double)locs[0];
//  y = EdgeModel21_.Interpol(FPACoordinateGrid_,x);
//  b = (y - vals[0]);
//  std::cerr << "findIntensityCoeff::x,y,vals[x],b = " << x << ", " << y << ", " << vals[0] << ", " << b << "\n";

//  x = (double)locs[len-1];
//  y = EdgeModel21_.Interpol(FPACoordinateGrid_,x);
//  m = (y / vals[len-1]);
//  std::cerr << "findIntensityCoeff::x,y,vals[x],m = " << x << ", " << y << ", " << vals[len-1] << ", " << m << "\n";

//  tmp[0] = m;
//  tmp[1] = b;
  return tmp;
}


EventVar<double> SolarEdgeModel::findNoiseFloorAndBlur(EventVar<double> &vals, EventVar<int> &locs)
{
  EventVar<double> tmp;
  int i;
//  int len = ExtentNormalizedGrid_.size();
  int len = locs.size(); 
  double m,b;
  double x,y;

  std::cerr << "findNoiseFloorAndBlur::len = " << len << "\n";
  lmInitialize(EdgeModel_,FPACoordinateGrid_);
  tmp = lmFitNoiseFloorAndBlur(vals,locs);

//  x = (double)locs[0];
//  y = EdgeModel_.Interpol(FPACoordinateGrid_,x);
//  b = (y - vals[0]);
//  std::cerr << "findIntensityCoeff::x,y,vals[x],b = " << x << ", " << y << ", " << vals[0] << ", " << b << "\n";

//  x = (double)locs[len-1];
//  y = EdgeModel_.Interpol(FPACoordinateGrid_,x);
//  m = (y / vals[len-1]);
//  std::cerr << "findIntensityCoeff::x,y,vals[x],m = " << x << ", " << y << ", " << vals[len-1] << ", " << m << "\n";

//  tmp[0] = m;
//  tmp[1] = b;
  return tmp;
}


void SolarEdgeModel::applyNoiseFloorAndBlur(EventVar<double> &coeffs)
{
  double blur = coeffs[0];
  double nf   = coeffs[1];
  double mx   = coeffs[2];

  int maf = (int)(1000 * blur);

  std::cerr << "applyNoiseFloorAndBlur::nf   = " << nf   << "\n";
  std::cerr << "applyNoiseFloorAndBlur::blur = " << blur << "\n";
  std::cerr << "applyNoiseFloorAndBlur::maf  = " << maf  << "\n";
  std::cerr << "applyNoiseFloorAndBlur::mx   = " << mx   << "\n";
  

  AverageIntensity_ = AverageIntensity_.Smooth(maf);
//  AverageIntensity21Pixel_ = AverageIntensity21Pixel_ * mx + nf;
  AverageIntensity_ = (AverageIntensity_ + nf)* mx;
//  EdgeModel21_ = EdgeModel21_.Smooth(maf);
//  EdgeModel21_ = EdgeModel21_  + nf;
}