//   Read the documentation to learn more about C++ code generator
//   versioning.
//	  %X% %Q% %Z% %W%

// Fit
#include <XSFit/Fit/Fit.h>
// Step
#include <XSFit/Fit/Step.h>

#include <XSFit/Fit/StatMethod.h>
#include <XSModel/GlobalContainer/ModelContainer.h>
#include <XSUtil/Utils/XSstream.h>
#include <XSUtil/Parse/XSparse.h>
#include "XSContainer.h"
#include "XSsymbol.h"
#include "XSstreams.h"  
#include <iomanip>


// Class Step 

Step::Step (Fit* fit)
  : m_best(false),
    m_bestFit(.0),
    m_fit(fit),
    m_increasing()
{
}


Step::~Step()
{
}


void Step::doGrid ()
{
  using namespace std;      

  // suppress all fit messages - this isn't quite what we want

  SpecContainer& parameters = parameter();
  int N(parameters.size());
  m_increasing.clear();
  m_increasing.resize(N, true);
  const bool saveFitStatus = m_fit->isStillValid();
  SpecContainer::iterator p    = parameters.begin();
  SpecContainer::iterator pEnd = parameters.end();
  try
  {
     while ( p != pEnd )
     {
        ParameterSpec*& spec = *p; 

        // freeze each parameter by removing from list of fit's variable
        // parameters. They will be added on command exit. 
        m_fit->freezeByIndex(spec->index);
        ModParam* mp = spec->address;
        // Freeze call really only matters for norm param, since that will
        // be accessed even if it's not part of fit->variableParameters map
        // (see StatMethod::renorm function).  
        mp->freeze();
        mp->setValue(spec->lowRange,'v');
        ++p;
     } 
     m_fit->calculateModel();   
     std::vector<Real>& highestDimensionArray = parameters.back()->parameterValues;
     int Npt ( highestDimensionArray.size());
     IntegerArray coords(N,0);

     bool TITLE(true);
     bool NOTITLE(false);

     report(TITLE);

     if ( N > 1 )
     {
           for ( int pt = 0; pt < Npt; ++pt)
           {
                   parameters[N-1]->value = highestDimensionArray[pt];
                   coords[N-1] = pt;
                   stepThruDimension(N-1,pt,coords);
           }
     }
     else
     {
           for ( int pt = 0; pt < Npt; ++pt)
           {
                   parameters[0]->value = highestDimensionArray[pt];
                   doPoint();
                   Real stat(m_fit->statistic());

		   ios_base::fmtflags saveFmt(tpout.flags());
                   size_t savePrecision(tpout.precision());
		   tpout.setf(ios_base::showpoint|ios_base::right,ios_base::scientific);
		   tpout.precision(5);		

		   tpout << right << setw(16) << stat
		         << right << setw(12) << stat - m_bestFit;

		   tpout << right << setw(5) << pt
		         << right << setw(12) << parameters[0]->parameterValues[pt] << std::endl;

                   tpout.flags(saveFmt);
                   tpout.precision(savePrecision);
                   grid()[pt] = stat;
           }

     }
     restoreFit();
     m_fit->isStillValid(saveFitStatus);
  }
  catch (...)
  {
     restoreFit();
     m_fit->isStillValid(saveFitStatus);
     throw;
  }
}

void Step::retrieveParameter (Fit* fit, string& s, ModParam*& parameter, int& paramIndex, string& paramName)
{
  using namespace XSContainer;
  string modelName("");
  size_t index(0);      
  XSparse::stringIntPair(s,modelName,index);
  paramIndex = models->fullIndex(index,modelName);
  Parameter*  p(models->lookupParameter(paramIndex));

  if ( p )
  {
       if ( modelName.length())
       {
                paramName = modelName + string(":");       
       }
       paramName += p->name();
       parameter = fit->variableParameters(paramIndex);
       if ( !parameter )
       {
                string msg (" Step: Parameter ");
                msg += s;
                msg += " is not variable";
                throw InvalidParameter(msg);                       
       }               
  }
  else
  {
       string msg (" Step: Parameter ");
       msg += s;
       msg += " is not defined";
       throw InvalidParameter(msg);        
  }

}

void Step::stepThruDimension (int dimension, int point, IntegerArray& coordinates)
{
    SpecContainer& parameters = parameter();
    const size_t N (parameters.size());
    int nextDimension (dimension - 1);
    std::vector<Real>& currentParamArray = parameters[dimension - 1]->parameterValues; 
    const size_t NCD ( currentParamArray.size() );        

    using namespace std;

    if ( dimension != 1)
    {
	// iterate forwards
	if  ( m_increasing[dimension - 1])
	{
	    for ( size_t j = 0; j < NCD; ++j )
	    {
		parameters[dimension-1]->value = currentParamArray[j];
		coordinates[dimension-1] = j;
		stepThruDimension(nextDimension,j,coordinates);
	    }
	    m_increasing[dimension-1] = false;
	}
	else
	{
	    for ( int j = NCD - 1; j >= 0; --j)
	    {
		parameters[dimension-1]->value = currentParamArray[j];
		coordinates[dimension-1] = j;
		stepThruDimension(nextDimension,j,coordinates);
	    }
	    m_increasing[dimension-1] = true;
	}
    }      
    else
    {
	int j = 0, limit = NCD - 1, incr = 1;

	if  ( !m_increasing[0] )
	{
	    j = NCD - 1;
	    incr = -1;
	    limit = 0;
	}

	ios_base::fmtflags saveFmt(tpout.flags());
        size_t savePrecision(tpout.precision());
	tpout.setf(ios_base::showpoint|ios_base::right,ios_base::scientific);
	tpout.precision(5);

	for(; (incr > 0 ? j <= limit : j >= limit); j += incr)
	{
	    parameters[0]->value = currentParamArray[j];
	    coordinates[0] = j;
	    // here and below we should have a correct set of
	    // parameter values and be able to replace this 
	    // output with the function that actually fits.
	    doPoint();
	    //XSutility::starPrint("Calculating...",j,tcout);
	    Real stat(m_fit->statistic());

	    tpout << right << setw(16) << stat
		  << right << setw(12) << stat - m_bestFit;

	    for(int i = 0; i < N; ++i)
		tpout << right << setw(5) << coordinates[i]
		      << right << setw(12) << parameters[i]->parameterValues[coordinates[i]];

	    tpout << endl;

	    int arrayPoint(coordinates[0]);
	    int block(1);
	    for (size_t k = 1; k < N; ++k)
	    {
		block *= (parameters[k - 1]->intervals + 1);
		arrayPoint += block*coordinates[k];       
	    }
	    grid()[arrayPoint] = stat;
	}
        tpout.flags(saveFmt);
        tpout.precision(savePrecision);
	m_increasing[0] = !m_increasing[0];
    }        

}

void Step::restoreFit ()
{
  SpecContainer& parameters = parameter();     
  SpecContainer::iterator sp = parameters.begin();
  SpecContainer::iterator spEnd = parameters.end();

  // put back parameters into fit's variable parameters set. 
  while ( sp != spEnd )
  {
        ParameterSpec*& p = *sp;
        p->address->thaw();  
        m_fit->variableParameters(p->index,p->address);         
        ++sp;       
  }

  std::map<int,ModParam*>::const_iterator v (m_fit->variableParameters().begin());       
  std::map<int,ModParam*>::const_iterator vEnd (m_fit->variableParameters().end());
  while ( v != vEnd )
  {
          int index ( v->first );
          Real oldVal = m_fit->oldParameterValues(index)->value('v');
          m_fit->variableParameters(index)->setValue(oldVal,'v'); 
          ++v;     
  }       

  // go through again and set the value fields of the parameter specs to the 
  // best fit value.
  sp = parameters.begin();
  for ( ; sp != spEnd; ++sp ) (*sp)->value = (m_fit->variableParameters((*sp)->index)->value('v'));

  m_fit->calculateModel();

  // We want the parameters set back to their original values prior to steppar,
  // so no renormalization.
  Fit::RenormSetting saveRenorm = m_fit->renormType();
  m_fit->renormType(Fit::NONE);
  try
  {
     m_fit->initializeStatistic();
  }
  catch (...)
  {
     m_fit->renormType(saveRenorm);
     throw;
  }
  m_fit->renormType(saveRenorm);
}

void Step::doPoint ()
{
  SpecContainer& parameters = parameter();     
  SpecContainer::iterator sp = parameters.begin();
  SpecContainer::iterator spEnd = parameters.end();
  // the parameters held in SpecContainer have been taken out of Fit's variable
  // parameter set and are therefore frozen. Now we set those parameter values to
  // the grid array values and fit .
  while ( sp != spEnd )
  {
        ModParam* mp ( (*sp)->address );
        mp->setValue((*sp)->value ,'v' ); 
        ++sp;       
  }
  if ( m_best )
  {
        // the set of indices of variableParameters is a subset of that of
        // oldParameterValues
        std::map<int,ModParam*>::const_iterator v (m_fit->variableParameters().begin());       
        std::map<int,ModParam*>::const_iterator vEnd (m_fit->variableParameters().end());
        while ( v != vEnd )
        {
                int index ( v->first );
                Real oldVal = m_fit->oldParameterValues(index)->value('v');
                m_fit->variableParameters(index)->setValue(oldVal,'v'); 
                ++v;     
        }       
  }

  // 
  m_fit->calculateModel();

  if (m_fit->variableParameters().size())  
  {
     // renorm, DOF check, recompute current statistic
     m_fit->initializeStatistic();

     m_fit->perform();
  }
  else
  {
     m_fit->statMethod()->perform(m_fit);
  }
}

void Step::report (bool title) const
{
    using namespace std;
    // preserve current chatter level settings and force output.
    int con (tpout.consoleChatterLevel());
    int log (tpout.logChatterLevel());
    tpout.consoleChatterLevel( tpout.conVerbose() );
    tpout.logChatterLevel( tpout.logVerbose() );

    StatMethod* stat = m_fit->statMethod();

    const SpecContainer& parameters = getParameter();
    int Npar  (parameters.size());
    if ( title )
    {
            tpout << "\n" << setw(16) <<  stat->fullName()
                  << setw(12)   <<   "  Delta   "; 
            for (int  j = 0; j < Npar; ++j )
            {
                tpout << setw(14) << parameters[j]->name;       
            }

            tpout << "\n";
            tpout << setw(16) <<    " "
                  << setw(12)   <<   stat->fullName(); 
            for (int  j = 0; j < Npar; ++j )
            {
                tpout << setw(14) << parameters[j]->address->index();       
            }
            tpout << endl << endl;

    }
    else
    {
            ios_base::fmtflags save(tpout.flags());
            tpout.setf(ios_base::showpoint|ios_base::right,ios_base::scientific);
            tpout.precision(5);
            int Ngrid (getGrid().size());
            IntegerArray coordinates(Npar,0);
            IntegerArray offsets (Npar,1);
            for (int k = 1; k < Npar; ++k) 
            {
                    offsets[k] = parameters[k]->parameterValues.size()*offsets[k-1];
            }
            // arrayPoint, iteration.

            for (int j = 0; j < Ngrid; ++j)
            { 
                    int k (Npar-1);
                    int remainder = j;
                    while ( k >=  0)
                    {

                            coordinates[k] = remainder / ( offsets[k] );
                            remainder -= coordinates[k] * offsets[k];
                            --k;                        
                    }
                    Real stat = getGrid()[j];
                    tpout << right << setw(16)  << stat  << right << setw(12)  << stat - m_bestFit; 
                    for (k = 0 ; k < Npar; ++k)
                    {
                            tpout << right << setw(4)  << coordinates[k]  
                                  << right << setw(8)  
                                  << parameters[k]->parameterValues[coordinates[k]]; 
                    }
                    tpout << '\n';
            }
            tpout.flags(save);                               

    }
    tpout << flush;


    tpout.consoleChatterLevel(con);
    tpout.logChatterLevel(log);
}

// Additional Declarations