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

// Weight
#include <XSModel/GlobalContainer/Weight.h>
// StatMethod
#include <XSFit/Fit/StatMethod.h>

#include <XSFit/Fit/Fit.h>
#include <XSFit/Fit/FitMethod.h>
#include <XSModel/Model/Model.h>
#include <XSModel/Data/SpectralData.h>
#include <XSModel/Data/DataUtility.h>
#include <XSModel/Data/Detector/Response.h>
#include <XSModel/GlobalContainer/DataContainer.h>
#include <XSModel/GlobalContainer/ModelContainer.h>
#include <XSModel/Parameter/ModParam.h>
#include <XSModel/Parameter/ResponseParam.h>
#include <XSUtil/Utils/XSstream.h>
#include <XSstreams.h>
#include <XSContainer.h>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <set>
using namespace XSContainer;


// Class StatMethod::InvalidSpectrumNumber 

StatMethod::InvalidSpectrumNumber::InvalidSpectrumNumber (size_t num)
  : RedAlert(" initializing spectra for fit:\n ")  
{
  std::ostringstream message;
  message << " Statistic tried to access spectrum number " << num << " not loaded ";
  reportAndExit(message.str(),-5);
}


// Class StatMethod::InvalidParameterError 

StatMethod::InvalidParameterError::InvalidParameterError ()
  : RedAlert(" retrieving model derivative array:\n")              
{
  std::ostringstream message;
  message << " Statistic calculation attempted to retrieve parameter derivative not allocated ";
  reportAndExit(message.str(),-5);
}


// Class StatMethod::FitNotDefined 

StatMethod::FitNotDefined::FitNotDefined (const string& diag)
  : YellowAlert (" no model defined or no data read ")
{
  tcerr << diag << '\n';
}


// Class StatMethod::NoSuchSpectrum 

StatMethod::NoSuchSpectrum::NoSuchSpectrum (const string& diag)
  : YellowAlert(" spectrum not loaded ")
{
}


// Class StatMethod 
Bayes StatMethod::s_bayes = Bayes();

std::map<std::string, StatMethod*> StatMethod::s_statMethods;
std::map<string, XSContainer::Weight*> StatMethod::s_weightScheme;

StatMethod::StatMethod(const StatMethod &right)
  : m_statistic(right.m_statistic),
    m_isDerivAnalytic(right.m_isDerivAnalytic),
    m_alpha(right.m_alpha),
    m_name(right.m_name),
    m_weight(0),
    m_beta(right.m_beta),
    m_dModel(right.m_dModel),
    m_protoWorkspace(right.m_protoWorkspace),
    m_dRM(right.m_dRM)               
{
  if (right.m_weight) m_weight = right.m_weight->clone();
}

StatMethod::StatMethod (const std::string& name)
  : m_statistic(0.),
    m_isDerivAnalytic(false),
    m_alpha(),
    m_name(name),
    m_weight(0),
    m_beta(),
    m_dModel(),
    m_protoWorkspace(),
    m_dRM()
{
}


StatMethod::~StatMethod()
{
  delete m_weight;
}


void StatMethod::differentiate (Fit* fit, const IntegerArray& which)
{
   if (m_isDerivAnalytic)
      analyticDifferentiate(fit, which);
   else
      numericalDifferentiate(fit, which);
}

StatMethod* StatMethod::get (const std::string& name)
{
  std::map<std::string,StatMethod*>::iterator f(s_statMethods.find(name));
  if (f != s_statMethods.end()) return (*f).second->clone();
  else return 0;
}

void StatMethod::registerMethod (const string& name, StatMethod* method)
{
  s_statMethods.insert(std::map<std::string,StatMethod*>::value_type(name,method));
}

void StatMethod::registerWeightingMethod (const std::string& name, XSContainer::Weight* weightMethod)
{
  s_weightScheme.insert(std::map<std::string, XSContainer::Weight*>::value_type(name,weightMethod));
}

XSContainer::Weight* StatMethod::weightScheme (const std::string& name)
{
  string lcName = XSutility::lowerCase(name);
  XSContainer::Weight* requested = 0;
  std::map<std::string,XSContainer::Weight*>::const_iterator iw = s_weightScheme.lower_bound(lcName);
  if (iw != s_weightScheme.end()) 
  {
     if (iw->first.find(lcName) == 0)
     {
        requested = iw->second;
     }
  }
  return requested;
}

void StatMethod::setWeight (const std::string& weightName)
{
  static const string STANDARD ("standard");
  // Default setWeight only accepts standard weighting.
  const string lcName = XSutility::lowerCase(weightName);
  if (STANDARD.find(lcName) != 0)
  {
        tcout << "*** Warning: statistic " << m_name << " does not recognize weighting \n"
                << " scheme " << weightName << ": using standard weighting" << std::endl; 
  }
  else
  {
     tcout << " " << m_name << " statistic: variance weighted using standard weighting"
          << std::endl;
  }
  weight(weightScheme(std::string("standard")));
}

void StatMethod::initialize (Fit* fit)
{

  doInitialize(fit);
}

void StatMethod::doInitialize (Fit* fit)
{

  // clear all and start again.
  deallocate();
  // allocate memory for derivative arrays. If derivatives are not
  // required we don't need to do anything.

  // create proto workspace

  std::map<size_t,SpectralData*>::const_iterator sp (fit->spectraToFit().begin());
  std::map<size_t,SpectralData*>::const_iterator spEnd (fit->spectraToFit().end());

  while (sp != spEnd)
  {
        size_t N = sp->second->indirectNotice().size();
        m_protoWorkspace.insert(ArrayContainer::value_type(sp->first,RealArray(0.,N)));  
        ++sp;       
  }


  if (fit->fitMethod()->firstDerivativeRequired())
  {
        size_t numFreeParam = fit->variableParameters().size();
        m_beta.resize(numFreeParam,0.);


        std::map<int,ModParam*>::const_iterator vp = fit->variableParameters().begin();
        std::map<int,ModParam*>::const_iterator vpEnd = fit->variableParameters().end();

        // one derivative arrayContainer per variable parameter per source
        // and one real array per spectrum is required. 
        // Because of linked parameters, more than one model
        // can contribute to the derivative with respect to an individual spectrum.
        while (vp != vpEnd)
        {
            m_dModel.insert(std::map<size_t,ArrayContainer>::value_type(vp->first,protoWorkspace()));                                
            ++vp;        
        }
        //the original plan was to size the derivative arrays here, but that
        // would be rather complicated because in practice it would mean extracting
        // data group information and comparing it with parameter indexing to
        // see which parameters belonged to each data group. In any case many of the
        // derivative arrays will be zero, so it's probably more efficient to 
        // resize on the first iteration and then to check the size at subsequent
        // iterations. But we can determine how large the derivative array container
        // will be here.


        if (fit->fitMethod()->secondDerivativeRequired())
        {
                m_alpha.resize(numFreeParam*numFreeParam,0.0);       
        }       

  }
}

void StatMethod::deallocate ()
{

        doDeallocate();
}

void StatMethod::doDeallocate ()
{
        m_alpha.resize(0);
        m_beta.resize(0);
        m_dModel.clear(); 
        m_protoWorkspace.clear(); 
        m_dRM.clear();  
}

void StatMethod::renormalize (Fit* fit)
{

        if (  models->activeModelNames().empty() || datasets->numberOfSpectra() == 0)
        {
                throw Fit::CantInitialize(" renorm: either no model defined or no data read ");      

        }
        else
        {
           if (m_name == "cstat" || m_name == "lstat")
           {
              // These stat methods aren't set up for renorming with
              // background data.  
              for (DataArrayConstIt itDs = datasets->dataArray().begin(); 
                    itDs != datasets->dataArray().end(); ++itDs)
              {
                 const DataSet* ds = itDs->second;
                 if (ds->isMultiple())
                 {
                    SpectralDataMapConstIt itSd = ds->multiSpectralData().begin();
                    while (itSd != ds->multiSpectralData().end())
                    {
                       const SpectralData* sd = itSd->second;
                       if (sd && sd->background())  return;
                       ++itSd;
                    }
                 }
                 else if (const SpectralData* sd = ds->spectralData())
                 {
                    if (sd->background())  return;
                 }
              }
           }
           // perform the renormalization...
           ModelMapConstIter mod (models->modelSet().begin());
           ModelMapConstIter modEnd (models->modelSet().end());

           // first, compute and store current model.
           std::map<Model*,bool> frozenComponents;

           size_t numberOfNorms (0);
           size_t numberOfFrozenNorms(0);
           while ( mod != modEnd )
           {
              Model* currentModel = mod->second;
              frozenComponents[currentModel] = false;
              if ( currentModel->isActive() ) 
              {                   
                 std::list<ModParam*> norms (currentModel->normParams());
                 numberOfNorms += norms.size();
                 for (std::list<ModParam*>::const_iterator np = norms.begin();
                      np != norms.end(); ++np)
                 {
                     if ((*np)->isFrozen() || Parameter::isLinkedToFrozen(*np))
                     {
	                frozenComponents[currentModel] = true; 
                        ++numberOfFrozenNorms;
                     }            
                 }  
              }
              ++mod;
           }
           if (numberOfFrozenNorms == numberOfNorms )
           {
               if (!fit->errorCalc())
               {        
                   tcout << " Warning: renorm - no variable model to allow  renormalization"
                         << std::endl; 
               }   
           }
           else
           {
              mod    = models->modelSet().begin();

              std::pair<Real,Real> fSums;
              Real ofSum(0.);    // =>     fSums[0]
              Real mfSum(0.);    // =>     fSums[1]
              doReset(fit);

              // Want to analyze spectra grouped by data group number.
              // DataContainer currently provides no ready-made way to
              // do this, so in order to avoid O(N^2) behavior let's
              // collate things here.
              const size_t nDgs = datasets->numberOfGroups();
              const size_t nSpecs = datasets->numberOfSpectra();
              std::vector<std::vector<const SpectralData*> > specsByGroup(nDgs);
              for (size_t iSp=1; iSp<=nSpecs; ++iSp)
              {
                 const SpectralData* spec = datasets->lookup(iSp);
                 const size_t dgNum = spec->parent()->dataGroup();
                 specsByGroup[dgNum-1].push_back(spec);
              }

              for (size_t iDg=0; iDg<nDgs; ++iDg)
              {
                 std::vector<Model*> modsForGroup(datasets->numSourcesForSpectra(),0);
                 std::map<Model*,ArrayContainer > saveFrozenModVals;
                 const std::set<size_t>& sourcesForGroup =  
                        datasets->dgToSources().find(iDg+1)->second;
                 std::set<size_t>::const_iterator itSource =
                        sourcesForGroup.begin();
                 std::set<size_t>::const_iterator itSourceEnd =
                        sourcesForGroup.end();
                 while (itSource != itSourceEnd)
                 {
                    const string modName = models->lookupModelForSource(*itSource);
                    if (modName.length())
                    {
                       // Model must be active/on if in here.
                       Model* mod = models->lookup(modName, iDg+1);
                       modsForGroup[*itSource-1] = mod;
                       mod->fold();
                       if (frozenComponents[mod])
                       {
                          ArrayContainer& savModVal = saveFrozenModVals[mod];
                          ArrayContainer::const_iterator it = mod->foldedModel().begin();
                          ArrayContainer::const_iterator itEnd = mod->foldedModel().end();
                          while (it != itEnd)
                          {
                             savModVal.insert(ArrayContainer::value_type(it->first,it->second));
                             ++it;
                          }
                          mod->calculate(false, true);
                          mod->fold();
                       }
                    }
                    ++itSource;
                 }

                 const std::vector<const SpectralData*>& specsForGroup = 
                                specsByGroup[iDg];
                 for (size_t iSp=0; iSp<specsForGroup.size(); ++iSp)
                 {
                    const SpectralData* spec = specsForGroup[iSp];
                    const std::vector<Response*>& detectors = spec->detector();
                    std::vector<Model*> modsForSpec;
                    for (size_t iDet=0; iDet<detectors.size(); ++iDet)
                    {
                       if (detectors[iDet])
                       {
                          if (modsForGroup[iDet])
                             modsForSpec.push_back(modsForGroup[iDet]);
                       }
                    }
                    // adjustForFrozen functions ASSUME at least one 
                    // model is applied to spec, so filter out no model 
                    // case here.
                    if (modsForSpec.size())
                    {
                       fSums = adjustForFrozen(spec, modsForSpec, saveFrozenModVals);
                       ofSum += fSums.first;
                       mfSum += fSums.second;
                    }
                 } 
              } 


              Real renorm(1.);
              if(std::abs(ofSum - mfSum) < SMALL)
              {
                  if (!fit->errorCalc())
                  {
                      tcout << " renorm: no renormalization necessary" <<std::endl;   
                  }
              }
              else if(std::abs(mfSum) > SMALL)
              { 
	          renorm = ofSum / mfSum;
                  mod = models->modelSet().begin();
                  while (mod != modEnd)
                  {
	              Model* currentModel (mod->second);
                      if ( currentModel->isActive()) 
                      {
	                 std::list<ModParam*> norms (currentModel->normParams());
	                 std::list<ModParam*>::iterator np = norms.begin();
	                 std::list<ModParam*>::iterator npEnd = norms.end();
	                 while (np != npEnd)
	                 {
		             if (!(*np)->isFrozen() && !(*np)->isLinked())
		             {
		                 Real newValue = (*np)->value()*renorm;
		                 newValue = std::max( (*np)->value('l'),
					      std::min(newValue,(*np)->value('h')));
		                 (*np)->setValue(newValue,'v');
		             }
		             ++np;       
	                 }
                      }
	              ++mod;
                  }
              }
           } // end frozenNorms != numNorms

        }
}

string StatMethod::weightNames ()
{
  std::map<string,XSContainer::Weight*>::const_iterator ws = s_weightScheme.begin();
  std::map<string,XSContainer::Weight*>::const_iterator wsEnd = s_weightScheme.end();
  string names("");
  while ( ws != wsEnd)
  {
          names += ws->first;
          ++ws;       
          if (ws != wsEnd) names += " | ";
          else names += " ";
  }

  return names;      
}

string StatMethod::statNames ()
{
  std::map<string,StatMethod*>::const_iterator sm = s_statMethods.begin();
  std::map<string,StatMethod*>::const_iterator smEnd = s_statMethods.end();
  string names("");
  while ( sm != smEnd)
  {
          names += sm->first;
          ++sm;       
          if (sm != smEnd) names += " | ";
          else names += " ";
  }

  return names;      
}

void StatMethod::degreesOfFreedom (Fit* fit, int& degrees, size_t& bins) const
{
        // compute the number of degrees of freedom for a fit.
        std::map<size_t,SpectralData*>::const_iterator   sm (fit->spectraToFit().begin());
        std::map<size_t,SpectralData*>::const_iterator   smEnd (fit->spectraToFit().end());

        while (  sm != smEnd )
        {
                bins += sm->second->indirectNotice().size();
                ++sm;
        }      

        degrees = bins - fit->variableParameters().size();
}

void StatMethod::report (Fit* fit) const
{

  int dof = 0;
  size_t bins = 0;
  degreesOfFreedom(fit, dof, bins);
  if (dof < 0)
  {
     tcout << "Ill-formed Fit problem - number of variable parameters exceeds number of bins"
        << std::endl;
  }
  else
  { 
     if (!fit->spectraToFit().empty() && !fit->activeModels().empty())
        doReport(fit);
     if (!fit->isStillValid())
        tcout << " Current data and model not fit yet." << std::endl;
  }
}

void StatMethod::resetDerivativeArrays ()
{
  // is it quicker to trash the arrays and reallocate arrays initialized to
  // zero? I doubt it. Anyway, all of the arrays here are valarrays and
  // so the code is very simple.
  m_alpha = 0.;
  m_beta  = 0.;

  std::map<size_t,ArrayContainer >::iterator  dm = m_dModel.begin();
  std::map<size_t,ArrayContainer >::iterator  dmEnd = m_dModel.end();   
  while ( dm != dmEnd )
  {
        ArrayContainer::iterator dmsp    = dm->second.begin();
        ArrayContainer::iterator dmspEnd = dm->second.end();
        while ( dmsp != dmspEnd)
        {
                dmsp->second = 0.;
                ++dmsp;  
        }
        ++dm;       
  } 
}

void StatMethod::peakResidual (const Fit* fit, std::pair<Real,Real>& max, std::pair<Real,Real>& min, size_t spectrum, const std::pair<Real,Real> range) const
{

  ArrayContainer difference (getDifferences(fit));

  doPeakResidual(fit,difference,max,min,spectrum,range);
}

void StatMethod::doPeakResidual (const Fit* fit, const ArrayContainer& difference, std::pair<Real,Real>& max, std::pair<Real,Real>& min, size_t spectrum, const std::pair<Real,Real> range) const
{

  if ( fit->activeModels().empty() || fit->spectraToFit().empty()) 
  {

          throw StatMethod::FitNotDefined(" Peak Residual");
  }


  std::map<size_t,SpectralData*>::const_iterator sm = fit->spectraToFit().find(spectrum);
  if (sm == fit->spectraToFit().end()) 
  {
        std::ostringstream s;
        s << " Peak residual " << spectrum;
        throw StatMethod::NoSuchSpectrum(s.str());      
  }

  SpectralData* sp = sm->second;
  ArrayContainer::const_iterator f = difference.find(spectrum);
  const RealArray& residual  = f->second;
  const std::valarray<size_t>& IN = sp->indirectNotice();
  int M (IN.size());
  int lower(0);
  int upper(M-1);
  if ( range.first != 0 && range.second != 0)
  {
     sp->energyChannel(range.first,lower);
     sp->energyChannel(range.second,upper);
  }
  Real maxResid (-LARGE);
  Real maxEnergy (0);


  Real minResid (LARGE);
  Real minEnergy (0);
  int jmax (0),jmin (0);
  const RealArray& effectiveArea = sp->detector(0)->efficiency();
  for (int j = lower; j <= upper; ++j)
  {
        if  (effectiveArea[j] > SMALL)
        {        
                Real resid = residual[j];
                if (resid > maxResid) 
                {
                        maxResid = resid;
                        maxEnergy = sp->channelEnergy(j);
                        jmax = j;

                }         
                if (resid < minResid) 
                {
                        minResid = resid;
                        minEnergy = sp->channelEnergy(j);
                        jmin = j;

                } 
        }  
  }

  int idx;
  const Response* rsp = sp->detector(0);
  XSutility::find(rsp->energies(), maxEnergy, idx);
  if (idx < 0 || idx >= effectiveArea.size())
  {
     // This should never happen
     throw YellowAlert("Peak resid energy is outside effective area array");
  }
  max.first = maxResid/effectiveArea[idx];
  max.second = maxEnergy;
  XSutility::find(rsp->energies(), minEnergy, idx);
  if (idx < 0 || idx >= effectiveArea.size())
  {
     // This should never happen
     throw YellowAlert("Peak resid energy is outside effective area array");
  }
  min.first = minResid/effectiveArea[idx];
  min.second = minEnergy;

  // now, ...        
}

Real StatMethod::sumDerivs (int parameterIndex, bool isRespar)
{
   return .0;
}

void StatMethod::getWeightDerivs (Fit* fit)
{
}

Real StatMethod::plotChi (Real obs, Real mod, Real error, Real areaTime, const Real& NODATA) const
{

  return 0;
}

Real StatMethod::plotDeltaChi (Real obs, Real mod, Real error, const Real& NODATA) const
{

  return 0;
}

Real StatMethod::sumSecondDerivs (int pi, int pj, bool isPiResp, bool isPjResp)
{
   return .0;
}

void StatMethod::clearMethods ()
{
   std::map<std::string, StatMethod*>::iterator itSm = s_statMethods.begin();
   std::map<std::string, StatMethod*>::iterator itEnd = s_statMethods.end();
   while (itSm != itEnd)
   {
      delete itSm->second;
      ++itSm;
   }
   s_statMethods.clear();
}

void StatMethod::clearWeightingMethods ()
{
   std::map<std::string, XSContainer::Weight*>::iterator itWs = s_weightScheme.begin();
   std::map<std::string, XSContainer::Weight*>::iterator itEnd = s_weightScheme.end();
   while (itWs != itEnd)
   {
      delete itWs->second;
      ++itWs;
   }
   s_weightScheme.clear();
}

Real StatMethod::setOutputFormat (Fit* fit, int degrees) const
{
   using namespace std;
   Real deltaCrit = fit->fitMethod()->deltaCrit();
   if (deltaCrit <= 0.0)  
      throw YellowAlert("Zero or neg critical delta detected in StatMethod::setOutputFormat.\n");
   Real truncatedStat = .0;
   Real absStatistic = m_statistic;
   if (m_statistic == 0.0)
   {
      // just print 0.0
      tcout << fixed << setprecision(1);
   }
   else
   {
      if (m_statistic < .0)
         absStatistic *= -1.0;
      int expDelta = static_cast<int>(floor(log10(deltaCrit)));
      int expStat = static_cast<int>(floor(log10(absStatistic)));
      int nDigits = expStat - expDelta + 1;
      if (nDigits <= 0)
      {
         // This case shouldn't really ever happen, but on the outside
         // chance, adjust expDelta to give 1 digit.
         expDelta = expStat;
         nDigits = 1;
      }
      Real stat = absStatistic;
      int lowestExp = expDelta;
      if (degrees)
      {
         // We're dealing with a reduced-chisq calc.  
         // nDigits should remain the same.
         stat /= degrees;
         expStat = static_cast<int>(floor(log10(stat)));
         lowestExp = expStat - nDigits + 1;
      }
      double truncator = pow(10.0, lowestExp);
      truncatedStat = floor(stat/truncator + .5)*truncator;
      if (expStat > 5 || expStat < -3)
      {
         if (nDigits == 1)
         {
            // use general format
            tcout.setf(ios_base::fmtflags(0), ios_base::floatfield);
            tcout.precision(1);
         }
         else
         {
            int prec = min(nDigits-1, 6);
            tcout << scientific << setprecision(prec);
         }
      }
      else
      {
         int prec = (lowestExp < 0) ? -lowestExp : 0;
         prec = min(prec, 6-expStat);
         tcout << fixed << setprecision(prec);
      }
      // Put back the '-' if there is one.
      if (m_statistic < .0)
         truncatedStat *= -1.0;
   }
   return truncatedStat;
}

void StatMethod::numericalDifferentiate (Fit* fit, const IntegerArray& parNums)
{
   using namespace XSContainer;
   using namespace std;

   const size_t nPars = parNums.size();
   RealArray dPars(.0,nPars);
   for (size_t i=0; i<nPars; ++i)
   {
      ModParam* par = fit->variableParameters(parNums[i]);
      Real origVal = par->value('a');
      Real delta = par->value('d');      
      par->setValue(origVal+delta,'a');      
      perform(fit);
      Real plusDeltaStat = m_statistic;
      par->setValue(origVal-delta,'a');
      perform(fit);
      Real minusDeltaStat = m_statistic;
      dPars[i] = delta;
      par->setValue(origVal,'a');
      beta(i,-.25*(plusDeltaStat-minusDeltaStat)/delta);
   }
   if (fit->fitMethod()->secondDerivativeRequired())
   {
      size_t i=0;
      RealArray alphaElems(0.0,nPars*nPars);
      // Calculate:
      // [(S(pi+di,pj+dj)-S(pi-di,pj+dj))-(S(pi+di,pj-dj)-S(pi-di,pj-dj))]/(4di*dj)
      while (i < nPars)
      {
         ModParam* par_i = fit->variableParameters(parNums[i]);
         const Real orig_i = par_i->value('a');
         const size_t iPos = i*nPars;
         const Real dPars_i = dPars[i];
         for (int iDel=1; iDel>=-1; iDel-=2)
         {
            par_i->setValue(orig_i+iDel*dPars_i,'a');
            size_t j=0;
            while (j <= i)
            {
               ModParam* par_j = fit->variableParameters(parNums[j]);
               Real orig_j = par_j->value('a');
               for (int jDel=1; jDel>=-1; jDel-=2)
               {
                  par_j->setValue(orig_j+jDel*dPars[j],'a');
                  perform(fit);
                  const size_t jElem = j+iPos;
                  alphaElems[jElem] += iDel*jDel*m_statistic/(dPars_i*dPars[j]);
                  if (i != j)
                     alphaElems[i + nPars*j] = alphaElems[jElem];
               }
               par_j->setValue(orig_j,'a');
               ++j;
            }
         }
         par_i->setValue(orig_i,'a');
         ++i;
      }
      alphaElems /= 8.0;
      setAlpha(alphaElems);
   }
   // This will reset all model fluxes as well as statistic.
   perform(fit);
}

void StatMethod::analyticDifferentiate (Fit* fit, const IntegerArray& parNums)
{

    // differentiate the statistic with respect to parameters with
    // indices specified by parNums.

    // The m_difference array will contain the quantity
    // d = (D-F), the m_variance array will contain the variance var, 
    // corrected for any model systematic terms. Calculation of m_difference
    // is done in Model::difference(ArrayContainer&) which is invoked by
    // ChiSquare::doPerform

    // NOTE finally that since each model's parameters are independent of the other
    // models, the result of all this  will be one set of derivative arrays 
    // per model per spectrum. This is different from  the case of chisquare itself,
    // in which the difference can depend on two models simultaneously.
    using namespace XSContainer;
    using std::vector;
    using std::set;      
    getWeightDerivs(fit);

    const size_t N = parNums.size();
    size_t np (0);
    bool saveComponentFlux(true);
    static const int respFloor = Fit::RESPAR_INDEX() << ModelContainer::SHIFT();
    m_dRM.clear();
    for (size_t k = 0; k < N; ++k)
    {
       bool isRespar = false;
       int parNum = parNums[k];
       ModParam* pj (fit->variableParameters(parNum));
       Real p (pj->value('a'));
       Real delta (pj->value('d'));

        // results for each spectrum at parameter indexed by k.

        ArrayContainer& dm = m_dModel[parNum];


        set<string>::const_iterator ss    (fit->activeModels().begin());
        set<string>::const_iterator ssEnd (fit->activeModels().end());

        while ( ss != ssEnd )
        {
           // forward difference  
           const string& modelName = *ss;

           pj->setValue(p + delta,'a');
           // save the model and component fluxes to avoid unnecessary recalculation.
           GroupFluxContainer saveFlux;
           models->storeFluxes(modelName, saveFlux);
           models->calculate(modelName, saveComponentFlux);

           // save results in 'difference' container
           GroupFluxContainer forwardDifference;
           models->storeFluxes(modelName, forwardDifference);

           // backward difference
           pj->setValue(p - delta,'a');
           models->calculate(modelName);

           // difference contains the forward difference and the models
           // contain the backward difference. Combine these into a derivative,
           // fold, and return the value in dm

           models->makeDerivatives(modelName, forwardDifference, delta, dm);

           // ... and put the model back straight afterwards.

           models->restoreFluxes (modelName, saveFlux);
           models->resetComponentFluxes(modelName);

           ++ss;       
        }
        pj->setValue(p,'z');
        // Rather than do a dynamic_cast to tell if this is a
        // ResponseParam, simply rely on the fact that ResponseParams
        // are given an index above a set value.
        if (parNum >= respFloor)
        {
           ResponseParam* rpj = dynamic_cast<ResponseParam*>(pj);
           if (!rpj)
           {
              throw RedAlert("Expecting but did not find response parameter in StatMethod.cxx");
           }
           isRespar = true;
           Response* resp = rpj->responseParent();
           size_t specNum = resp->spectrumNumber();
           // calc (dR/dp)*M
           RealArray dr;
           resp->makeDerivative(delta, dr, rpj->parType());
           // calc dF/dp = d(R*M)/dp = (dR/dp)*M + R*(dM/dp)
           m_dRM[parNum] = ArrayContainer();
           m_dRM[parNum][specNum] = RealArray();
           RealArray& drm = m_dRM[parNum][specNum];
           drm.resize(dr.size());
           drm = dr;
           if (drm.size() != dm[specNum].size())
           {
              throw RedAlert("Response and folded model derivative array size mismatch");
           }
           drm += dm[specNum];

        }

	Real sum = sumDerivs(parNum, isRespar);
        // beta_i is defined as -.5*d(stat)/d(par_i)
        beta(np,-0.5*sum - s_bayes.dlPrior(fit, parNum));        
        ++np;
    }


    if ( fit->fitMethod()->secondDerivativeRequired() ) 
    {
        size_t ni (0);
        const size_t N(parNums.size());
        while ( ni <  N )
        {
	    int i = parNums[ni];
	    size_t nj (0);
	    while ( nj <= ni )
	    {
		int j = parNums[nj];
                Real alpha_ij = sumSecondDerivs(i, j, i>=respFloor, j>=respFloor);
		alpha(nj + N*ni,0.5*alpha_ij - s_bayes.d2lPrior(fit, i, j)) ;
		if (ni != nj ) alpha(ni + N*nj,0.5*alpha_ij);
		++nj; 
	    }  
	    ++ni; 
        }   
    }
}

void StatMethod::weight (XSContainer::Weight* value)
{
  delete m_weight;
  m_weight = value->clone();
  DataUtility::weightScheme(m_weight);
}

// Additional Declarations