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

// ChainManager
#include <XSFit/MCMC/ChainManager.h>
// FitMethod
#include <XSFit/Fit/FitMethod.h>
// StatMethod
#include <XSFit/Fit/StatMethod.h>
// ModParam
#include <XSModel/Parameter/ModParam.h>
// SpectralData
#include <XSModel/Data/SpectralData.h>
// RandomizerBase
#include <XSFit/Randomizer/RandomizerBase.h>
// ModelContainer
#include <XSModel/GlobalContainer/ModelContainer.h>
// DataContainer
#include <XSModel/GlobalContainer/DataContainer.h>
// Fit
#include <XSFit/Fit/Fit.h>

#include <XSFit/Fit/FitErrorCalc.h>
#include <XSModel/Data/SpectralData.h>
#include <XSModel/Data/BackCorr/Background.h>
#include <XSModel/GlobalContainer/ResponseContainer.h>
#include <XSModel/Parameter/ResponseParam.h>
#include <XSUtil/Numerics/Numerics.h>
#include <XSContainer.h>
#include <XSstreams.h>
#include <iomanip>
#include <algorithm>

Fit*   XSContainer::fit = 0;
using namespace XSContainer;


// Class Fit::NoSuchFitMethod 

Fit::NoSuchFitMethod::NoSuchFitMethod (const string& name)
  : YellowAlert(" Fitting Algorithm not loaded: ")
{
  tcerr << name << '\n';
}


// Class Fit::NoSuchStatMethod 

Fit::NoSuchStatMethod::NoSuchStatMethod (const string& name)
  : YellowAlert(" Statistic not recognized: ")
{
  tcerr << name << std::endl;;
}


// Class Fit::FitOverDetermined 

Fit::FitOverDetermined::FitOverDetermined()
  : YellowAlert() 
{
  tcerr << "***Warning: Ill-formed Fit problem - number of variable parameters exceeds number of bins"
     << std::endl;
}


// Class Fit::ParameterLookupError 

Fit::ParameterLookupError::ParameterLookupError()
  : RedAlert(" parameter indexing error ")               
{
}


// Class Fit::NotVariable 

Fit::NotVariable::NotVariable()
  : YellowAlert(" No variable parameters for fit ")
{
}


// Class Fit::CantInitialize 

Fit::CantInitialize::CantInitialize (const string& diag)
  : YellowAlert(" initializing statistic: ")
{
  tcerr << diag << '\n';
}


// Class Fit::HoldSpecVar 

Fit::HoldSpecVar::HoldSpecVar()
  : m_spectrum(), m_rawVariance(), m_variance(), m_bckSpectrum(), 
    m_bckRawVariance(), m_bckVariance()
{
}


Fit::HoldSpecVar::~HoldSpecVar()
{
}


void Fit::HoldSpecVar::saveSpecVar (const SpectralData& data)
{
  size_t nChans = data.channels();
  m_spectrum.resize(nChans);
  m_rawVariance.resize(nChans);
  m_variance.resize(nChans);
  m_spectrum = data.spectrum();
  m_rawVariance = data.rawVariance();
  m_variance = data.variance();
  const Background* bck = data.background();
  if (bck)
  {
     m_bckSpectrum.resize(nChans);
     m_bckRawVariance.resize(nChans);
     m_bckVariance.resize(nChans);
     m_bckSpectrum = bck->spectrum();
     m_bckRawVariance = bck->data()->rawVariance();
     m_bckVariance = bck->variance();
  }
}

void Fit::HoldSpecVar::restoreSpecVar (SpectralData& data) const
{
  data.setSpectrum(m_spectrum);
  data.setRawVariance(m_rawVariance);
  data.setVariance(m_variance);
  Background* bck = data.getBackground();
  if (bck)
  {
     data.setBackgroundSpectrum(m_bckSpectrum);
     bck->setRawVariance(m_bckRawVariance);
     bck->setVariance(m_bckVariance);
  } 
}

// Additional Declarations

// Class Fit::FitInterrupt 

Fit::FitInterrupt::FitInterrupt ()
  : YellowAlert()
{
  tcerr << "\n*** Warning: User interrupted fit, fit not valid."
        << std::endl;
}


// Class Fit 
Fit* Fit::s_instance = 0;
bool Fit::s_errorCalc = false;
Real Fit::s_deltaStat = 2.706;
int Fit::s_errorTry = 20;
Real Fit::s_tolerance = 0.01;
Real Fit::s_chiMax = 2.0;
int Fit::s_MCrealizations = 100;
Real Fit::s_lastFtest = .0;
const int Fit::s_RESPAR_INDEX = (1 << XSContainer::ModelContainer::SHIFT()-1) - 1;
Real Fit::s_lastGoodness = .0;

Fit::Fit (const string& fitMethodName, const string& statMethodName)
  : m_isStillValid(false),
    m_chainManager(ChainManager::Instance()),
    m_nativeRandomizerNames(),
    m_useNumericalDifferentiation(false),
    m_fitMethod(0),
    m_statMethod(0),
    m_oldParameterValues(),
    m_variableParameters(),
    m_spectraToFit(),
    m_queryMode(ON),
    m_stepGrid(0),
    m_activeModels(),
    m_randomizingStrategies(),
    m_renormType(PREFIT)
{
  m_fitMethod.reset(FitMethod::get(fitMethodName));
  m_statMethod.reset(StatMethod::get(statMethodName));

  // Register Fit as an observer of ModelContainer.  If
  // ModelContainer doesn't already exist, it will be created here.
  ModelContainer::Instance()->Attach(this);  
}


Fit::~Fit()
{
  // let's deallocate the resources anyway, but actually
  // Fit's dtor is only called by program exit.
  // Fit's fitMethod and statMethod are auto_ptrs thus will die automatically.
  deleteSavedParams();
  delete m_chainManager; 
  clearRandomizingStrategies();     
}


void Fit::renormalize ()
{

  m_statMethod->renormalize(this);
}

void Fit::reinitialize (bool saveParameters)
{
  cleanup();

  initializeFitModels();

  if (m_activeModels.empty()) 
  {
          throw CantInitialize(" no models defined");
  }

  initializeFitParameters();      

  if ( saveParameters)     
  {      
        // now create an (owned) copy of the parameters for resetting
        // purposes.
        std::map<int,ModParam*>::const_iterator vp    = m_variableParameters.begin();
        std::map<int,ModParam*>::const_iterator vpEnd = m_variableParameters.end();
        for ( ; vp != vpEnd; ++vp )
        {
                ModParam* savePar = vp->second->clone();
                m_oldParameterValues.insert(std::map<int,ModParam*>::value_type(vp->first,savePar));

        }
  }     

  std::vector<size_t> missingModels;   
  initializeFitSpectra(missingModels);
  if (missingModels.size())
  {
     std::ostringstream oss;
     oss << "Spectra with responses but no models: ";
     for (size_t i=0; i<missingModels.size(); ++i)
        oss << missingModels[i] << " ";
     throw CantInitialize(oss.str());
  }

  if ( m_spectraToFit.empty() )  throw CantInitialize(" no data defined ");

  // initialize the statistical calculation, allocating for the arrays
  // for intermediate results.
  m_statMethod->initialize(this);
}

FitMethod* Fit::fitMethod ()
{
  return m_fitMethod.get();
}

StatMethod* Fit::statMethod ()
{
  return m_statMethod.get();
}

void Fit::fitMethod (FitMethod* method)
{
  bool wasStillValid = m_isStillValid;
  m_fitMethod.reset(method);
  Update();
  // Update will automatically mark fit as invalid, but if
  // it was previously valid and all we did is change the
  // fit method, let's still call it valid.
  m_isStillValid = wasStillValid;    
}

void Fit::statMethod (StatMethod* statistic)
{
  bool wasStillValid = m_isStillValid;
  m_statMethod.reset(statistic);  
  Update();
  // Update will automatically mark fit as invalid, but if
  // it was previously valid and all we did is change the
  // stat method, let's still call it valid.
  m_isStillValid = wasStillValid;    
}

Fit* Fit::Instance (const string& fitMethodName, const string& statMethodName)
{
  if (s_instance == 0) {s_instance = new Fit(fitMethodName,statMethodName);} return s_instance;
}

void Fit::Update (Subject* changed)
{
  m_isStillValid = false;     
  cleanup();
  initializeFitModels();
  initializeFitParameters();
  m_chainManager->isSynchedWithFitParams(checkChainsForSynch());
  // In this context, we don't care if some spectra have responses
  // but no models since they aren't trying to do a fit at this
  // point.
  std::vector<size_t> dummyMissingMods;
  initializeFitSpectra(dummyMissingMods);
  calculateModel();

  if (!m_spectraToFit.empty() && !m_activeModels.empty()) 
  {
          m_statMethod->initialize(this);
          try
          {
             initializeStatistic(true);
  	     if (!s_errorCalc) m_statMethod->report(this);
          }
          catch (FitOverDetermined&)
          {
             // initializeStatistic can throw this.  We'll catch
             // it when called from Update to allow models to be
             // loaded which may have more free params than data
             // bins at the moment.  
          }
  }
  else
  {
     m_statMethod->statistic(0.0);
  }
}

void Fit::deleteSavedParams ()
{
  if (!m_oldParameterValues.empty())      
  {      
        std::map<int,ModParam*>::const_iterator op    = m_oldParameterValues.begin();
        std::map<int,ModParam*>::const_iterator opEnd = m_oldParameterValues.end();
        for ( ; op != opEnd; ++op )
        {
                delete  op->second;
        }    
        m_oldParameterValues.clear();
  }    
}

void Fit::cleanup ()
{
  deleteSavedParams();

  m_spectraToFit.clear();

  m_variableParameters.clear();
}

void Fit::resetParameters ()
{
  std::map<int,ModParam*>::const_iterator p (m_oldParameterValues.begin());
  std::map<int,ModParam*>::const_iterator pEnd (m_oldParameterValues.end());
  std::map<int,ModParam*>::iterator varEnd(m_variableParameters.end());
  while (p != pEnd )
  {
     // Note: m_variableParameters may not have all the ModParam
     // objects that were originally stored in m_oldParameterValues
     // (see for example its usage in FitErrorCalc::calcUncertainty).
     // In this case, simply do nothing.

     std::map<int,ModParam*>::iterator itVar = 
                m_variableParameters.find(p->first);
     if (itVar != varEnd)
     {
        itVar->second->setValue(p->second->value(), 'v');
     }
     ++p;
  }
}

void Fit::reportParameters (std::ostream& s, bool header)
{
  if (s_errorCalc) return;
  using namespace std;
  ios_base::fmtflags saveFormat(s.flags());
  const int savePrecision(s.precision());
  int nvpar = m_variableParameters.size();
  std::map<int,ModParam*>::const_iterator vp = m_variableParameters.begin();
  if (header)
  {

        // write header for output
        s << "\n" << setw(12) << left << m_statMethod->fullName() << setw(20)
                << " Lvl Par #";



        int j (0);
        // apes fortran code: print out five at a time except for the first
        // row which prints the statistical method name and the value.                
        while ( j < std::min(nvpar,4) )
        {
                s  << setw(14) << left << vp->second->getParameterLabel();
                ++j, ++vp;
        }
        s << std::endl;

        while ( j < nvpar )
        {
                int jfirst = j;
                for (int k = jfirst; k < std::min(jfirst + 5,nvpar); ++k)
                {  
                        s  << setw(14) << left <<  vp->second->getParameterLabel();
                        ++j,++vp;
                }   
                s << std::endl;    
        }          
  }      
  else
  {
        m_fitMethod->reportProgress(s,this); 
        s.precision(6);
        int j (0);  

          while ( j < std::min(nvpar,4) )
          {
                s << left << setw(14) << vp->second->value('v');
                ++j, ++vp;
          }
          s << std::endl;

          while ( j < nvpar )
          {
                int jfirst = j;
                for (int k = jfirst; k < std::min(jfirst + 5,nvpar); ++k)
                {  
                        s << left << setw(14) << vp->second->value('v');
                        ++j,++vp;
                }   
                s << std::endl;    
          }        

  }
  s.flags(saveFormat);
  s.precision(savePrecision);      
}

void Fit::report ()
{
  if (s_errorCalc) return;
  m_fitMethod->reportCovariance();
  // Only want to report active/on models.
  std::map<string,bool>::const_iterator itModNames = models->activeModelNames().begin();
  std::map<string,bool>::const_iterator itNamesEnd = models->activeModelNames().end();
  while (itModNames != itNamesEnd)
  {
     if (itModNames->second)
     {
        // It's active, is it on?
        std::vector<Model*> mods = models->lookupModelGroup(itModNames->first);
        if (mods.size() && mods[0]->isActive())
        {
           // It's on.
           mods[0]->printHeading();
           for (size_t i=0; i<mods.size(); ++i)
           {
              tcout << *mods[i];
           }
           mods[0]->printMixComp();           
           tcout << string(72,'_') << '\n' << std::endl;
        }
     }
     ++itModNames;
  }

  if (responses->totalResponseParams())
  {
     tcout << "\nResponse Parameters:" <<std::endl;
     responses->reportResponseParams();
  }                
  m_statMethod->report(this);

  if (m_useNumericalDifferentiation)
  {
     tcout << "\nNote that fit is using the full (slower) numerical differentiation method."
        <<"\nThis setting may be changed in the user's ~/.xspec/Xspec.init start-up file.\n"
        << std::endl;
  }      
}

void Fit::checkPegged ()
{

  // this is much simpler than the XSPEC11 routine which is supposed
  // to be sensitive as to whether the parameter is a norm param or not.
  // ... not clear that that code actually works, however...      
  std::map<int,ModParam*>::const_iterator vp = m_variableParameters.begin();
  std::map<int,ModParam*>::const_iterator vpEnd = m_variableParameters.end();

  while ( vp != vpEnd )
  {
     ModParam* par = vp->second;
     if (!par->isPeriodic())
     {
        if ( par->value() <= par->value('l') ||
             par->value() >= par->value('h') )
        {
           int savConVerbose = tpout.conVerbose();
           int savLogVerbose = tpout.logVerbose();
           XSstream::verbose(tpout, 15, 15);
           tcout << " Fit Parameter: "  << vp->second->getParameterLabel()
                 << " has pegged due to running into its hard limit." << std::endl;
           XSstream::verbose(tpout, savConVerbose, savLogVerbose);
           par->peg();
           break;       
        }
     }
     ++vp;       
  }
}

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

  m_statMethod->peakResidual(this,max,min,spectrum,range);
}

void Fit::simulate (std::vector<Real>& trialValues, bool simulateParams, int statMode)
{

  ArrayContainer modelsSim;
  Real stat ( m_statMethod->statistic() );  
  size_t higher(0), lower(0);

  size_t nPars (m_variableParameters.size());
  RealArray saveParameterValues(nPars);
  if (!simulateParams)
  {
        int dontRandomizeParameters(0);
        simulateModel(modelsSim,dontRandomizeParameters);       
  }
  else
  {
     reportErrorSimulationMethod();
     std::map<int,ModParam*>::const_iterator itVp = m_variableParameters.begin();
     std::map<int,ModParam*>::const_iterator itVpEnd = m_variableParameters.end();
     size_t k=0;
     while (itVp != itVpEnd)
     {
        saveParameterValues[k] = itVp->second->value('a');
        ++itVp, ++k;    
     }
  } 

  // Simulated spectrum is initialized for each trial by modelsSim.
  // Background spectrum though must be restored each time.
  std::map<size_t,RealArray> saveBackSpecs;
  std::map<size_t,SpectralData*>::const_iterator itCSpec = m_spectraToFit.begin();
  std::map<size_t,SpectralData*>::const_iterator itCSpecEnd = m_spectraToFit.end();
  while (itCSpec != itCSpecEnd)
  {
     const Background* bck = itCSpec->second->background();
     if (bck)
     {
        saveBackSpecs[itCSpec->first].resize(bck->spectrum().size());
        saveBackSpecs[itCSpec->first] = bck->spectrum();
     }
     ++itCSpec;
  }  

  const size_t nRealize = trialValues.size();
  int savConVerbose = tpout.conVerbose();
  int savLogVerbose = tpout.logVerbose();  
  try
  {
     // just to stress that the operation of randomizing the model parameters
     // and randomizing the background are quite different.
     for ( size_t j = 0; j < nRealize; ++j)
     {          
        if (simulateParams)
        {
           // First time through, flag = 2 tells model randomizer 
           // to initialize for run.
           int randomizeFlag = j ? 1 : 2;
           simulateModel(modelsSim,randomizeFlag);
           std::map<int,ModParam*>::const_iterator itVp = m_variableParameters.begin();
           for ( size_t k = 0; k < nPars; ++k, ++itVp)
           {
              itVp->second->setValue(saveParameterValues[k],'a');       
           }
        }

        std::map<size_t,SpectralData*>::iterator itSpec = m_spectraToFit.begin();
        std::map<size_t,SpectralData*>::iterator itSpecEnd = m_spectraToFit.end();
        XSstream::verbose(tpout,15,15);         
        while (itSpec != itSpecEnd )
        {
           size_t spectrumNumber = itSpec->first;
           SpectralData* data = itSpec->second;

           // Background must be restored each time, see note above.
           Background* bck = data->getBackground();
           if (bck)
           {
              bck->setSpectrum(saveBackSpecs[spectrumNumber]);
           }
           // These two functions will perform the randomization and store
           // the new spectra (including variances).
           data->simulateSpectrum(modelsSim[spectrumNumber]);
           data->simulateBackground();


           if ( tpout.maxChatter() >= 15 )
           {
              const RealArray& model = modelsSim[spectrumNumber];
              const RealArray& simSpec = data->spectrum();
              const RealArray& simVar = data->variance();
              const size_t nChans = simSpec.size();
              tcout << "Simulation: " << j << std::endl; 
              if (data->background())
              {
                 const RealArray& simBack = data->background()->spectrum();
                 const RealArray& bckVar = data->background()->variance();
                 tcout << " simModel " << " randSpec " << " randVar " 
                    << " rand B " << " randBVar " << std::endl; 
                 for (size_t k = 0; k < nChans ; ++k)
                 {
                    tcout << k << " " << model[k] << " " << simSpec[k] << " " 
                       << simVar[k] << "  " << simBack[k] << " " 
                       << bckVar[k] << std::endl;                                        
                 }
              }
              else
              {
                 tcout << " simModel " << " randSpec " << " randVar " << std::endl;
                 for (size_t k = 0; k < nChans ; ++k)
                 {
                    tcout << k << " " << model[k] << " " << simSpec[k] << " " 
                       << simVar[k] << std::endl;
                 }
              }
              tcout << " simvarTot " << simVar.sum() << std::endl;
           }
           ++itSpec;       
        }     
        XSstream::verbose(tpout,savConVerbose,savLogVerbose);
        m_statMethod->perform(this);
        trialValues[j] = m_statMethod->statistic();
        // tcout  << " ChiSquare " << trialValues[j] ;
        if (trialValues[j] > stat) ++higher;
        if (trialValues[j] <= stat) ++lower;
        // tcout << " H " << higher << " L " << lower;
        // tcout << std::endl;

        // reset the statistic after the fact.
        m_statMethod->statistic(stat);          
     }
  }
  catch (...)
  {
     XSstream::verbose(tpout,savConVerbose,savLogVerbose);
     // Only restore changes to the models.  The restoring of the 
     // original spectra and fit statistic will be handled
     // higher up.
     if (simulateParams)
     {
        std::map<int,ModParam*>::const_iterator itVp = m_variableParameters.begin();
        std::map<int,ModParam*>::const_iterator itEnd = m_variableParameters.end();
        size_t k=0;
        while(itVp != itEnd)
        {
           itVp->second->setValue(saveParameterValues[k],'a'); 
           ++itVp, ++k;      
        }
     }
     throw;
  }      
}

void Fit::goodness (int realizations, bool simulateParams, int statMode)
{
  using namespace std; 
  map<size_t,HoldSpecVar*> saveData;

  map<size_t,SpectralData*>::iterator f = m_spectraToFit.begin();
  map<size_t,SpectralData*>::iterator fEnd = m_spectraToFit.end();

  // Save the spectrum and variance (raw and processed) of every 
  // spectral data object that is fit.
  // We are going to need to replace them temporarily with simulated
  // realizations. 

  Real stat = m_statMethod->statistic();

  while ( f != fEnd)
  {
     HoldSpecVar *specVar = new HoldSpecVar();
     specVar->saveSpecVar(*f->second);
     saveData.insert(map<size_t,HoldSpecVar*>::value_type(f->first, specVar));
     ++f;
  }

  std::vector<Real> simResults(realizations,0.);     

  try
  {  
     simulate(simResults,simulateParams,statMode);
  }
  catch (...)
  {
     // Assume any changes to models have been restored at 
     // lower levels.  Restore the spectra here.
     f = m_spectraToFit.begin();
     while ( f != fEnd )
     {
        HoldSpecVar *specVar = saveData[f->first];
        specVar->restoreSpecVar(*f->second);
        delete specVar;
        ++f; 
     }
     m_statMethod->statistic(stat);
     throw;
  }
#ifndef STD_COUNT_DEFECT
        int n (count_if(simResults.begin(),simResults.end(),bind2nd(less<Real>(),stat)));
#else

        int n (0);
        count_if(simResults.begin(),simResults.end(),bind2nd(less<Real>(),stat),n);
#endif
  std::ios_base::fmtflags currentSetting(tcout.flags());
  size_t savePrecision = tcout.precision();
  Real pct = Real(n)/realizations*100.;
  tcout.precision(2);
  tcout <<  fixed << pct << "% of realizations are < best fit statistic " << stat;
  tcout.precision(7);
  simulateParams ? tcout << "  (sim)" : tcout << "  (nosim)";
  tcout << std::endl;
  s_lastGoodness = pct;
  tcout.flags(currentSetting);
  tcout.precision(savePrecision);
  // replace the original data. Got to ensure exception safety here.

  f = m_spectraToFit.begin();
  while ( f != fEnd )
  {
     HoldSpecVar *specVar = saveData[f->first];
     specVar->restoreSpecVar(*f->second);
     delete specVar;
     ++f; 
  }

}

void Fit::simulateModel (ArrayContainer& simSpec, int randomizeIndicator)
{
  // prepare one model spectrum for each real spectrum. Sum all the relevant
  // models, scale by exposure time and area, and add in the background
  // and any correction factor.

  switch (randomizeIndicator)
  {
     case 0:
        break;
     case 1:
        randomizeModelParameters(false);
        break;
     case 2:
     default:
        randomizeModelParameters(true);
        break;
  }

  calculateModel();  

  std::map<size_t,SpectralData*>::const_iterator f = m_spectraToFit.begin();
  std::map<size_t,SpectralData*>::const_iterator fEnd = m_spectraToFit.end();
  while ( f != fEnd )
  {
        SpectralData* data = f->second;
        const size_t NC (data->spectrum().size());
        RealArray sim(0.,NC);
        ModelMapConstIter mm (models->modelSet().begin());
        ModelMapConstIter mmEnd (models->modelSet().end());
        // sum all the folded models that are fit to a single spectrum
        // ( the normal case is one model per spectrum, but the following
        //   code supports the multisource case).
        while ( mm != mmEnd )
        {
                Model* m (mm->second);

                if (!m->isActive()) 
                {
                   ++mm;
                   continue;
                }
                ArrayContainer::const_iterator s = m->foldedModel().find(f->first);

                if ( s != m->foldedModel().end())
                {
                        sim += s->second;       
                }
                ++mm;

        }

        simSpec.insert(ArrayContainer::value_type(f->first,sim));
        ++f;      
  }  
}

void Fit::initializeFitModels ()
{
  // are there any models that have attached responses? 

  ModelMapConstIter mi (models->modelSet().begin());
  ModelMapConstIter miEnd (models->modelSet().end());
  m_activeModels.clear();
  for ( ; mi != miEnd; ++mi)
  {
        if (mi->second->isActive()) 
        {
                // active models is a set, will ignore duplicates.
                m_activeModels.insert(std::set<string>::value_type(mi->first));
                mi->second->prepareForFit();  
        }    
  }  
}

void Fit::initializeFitSpectra (std::vector<size_t>& missingModels)
{
  missingModels.clear();
  size_t N = datasets->numberOfSpectra();
  size_t nSources = datasets->numSourcesForSpectra();
  bool allHaveResp = true;
  bool thisHasResp = false;
  bool respHasMod = false;
  for ( size_t j = 1; j <= N; ++j)
  {
        thisHasResp = false;
        respHasMod = false;
        SpectralData* sp (datasets->lookup(j));
        // Requirement for spectrum to be added to fit:  at
        // least 1 detector slot must be filled with an actual
        // response suitable for fitting, and that slot must have
        // a model associated with it.  

        // If a spectrum has no responses, we'll let that go with just
        // a warning and won't add it to the missingModels container.
        // If it has 1 or more responses but no models, it gets added
        // to missingModels.
	for (size_t i=0; i<nSources; ++i)
	{
           if (sp->responseLoaded(i) && !sp->isDummyrspMode2(i)) 
           {
	      thisHasResp = true;
              if (models->lookupModelForSource(i+1).size())
              {
                 respHasMod = true;
                 break;
              }
           }
	}
	if (thisHasResp)
	{
           if (respHasMod)
           {
              if (!sp->spectrumIsZeroed())
              {
                 m_spectraToFit.insert(std::map<size_t,SpectralData*>::
		   				              value_type(j,sp));
              }
              sp->prepareForFit(); 
           }
           else
              missingModels.push_back(j);
	}
        else
           allHaveResp = false;      
  }
  if (!allHaveResp)
  {
     tcout << "***Warning!  One or more spectra are missing responses,\n";
     tcout << "               and are not suitable for fit.\n";
     tcout.flush();
  }
}

void Fit::initializeFitParameters ()
{
  // second, get an array of variable parameters, if saveParameters
  // is true (and therefore we are fitting rather than renormalizing).
  const ParamMap& pars = models->parameterList();
  ParamMapConstIter gp = pars.begin();
  ParamMapConstIter gpEnd = pars.end();

  for ( ; gp != gpEnd ; ++gp)
  {
     const string modelName = gp->first.substr(0,gp->first.find_first_of(':'));
     const std::vector<Model*> mods = models->lookupModelGroup(modelName);
     const size_t nParsInMod = mods[0]->numberOfParameters();
     const size_t groupNum = (gp->second->index()-1)/nParsInMod;     
     const Model* mod = mods[groupNum];

     if (mod && mod->isActive())
     {
          ModParam* mp = dynamic_cast<ModParam*>(gp->second);
          if (mp && !(mp->isFrozen() || mp->isLinked()) )
          {
             // create a map entry with the full index as
             // key and ModParam* value
             int mpIndex = models->keyToIndex(gp->first);
             // shallow copy
             variableParameters(mpIndex,mp); 
             //  tcout << " Parameter: " << (gp->first) << " index " << mpIndex 
             //        << " val: " << mp->value() << " address " << mp << std::endl;
          }
      }
  }

  // Now add any response parameters, give index a high enough
  // value so that it is extremely unlikely to ever conflict with
  // a model index.
  // This scheme should allow a max of 2^16 response params.
  if (responses->totalResponseParams())
  {
     const size_t MAX = 1 << (31 - ModelContainer::SHIFT() + 1);
     const size_t nSources = responses->respParContainer().size();
     if (responses->totalResponseParams() >= MAX)
     {
        std::ostringstream msg;
        msg << "Number of response parameters has exceeded the maximum of " << MAX;
        throw RedAlert(msg.str());
     }
     int k=1;
     for (size_t i=0; i<nSources; ++i)
     {
        const std::vector<ResponseParam*>& parsForSource = 
                responses->respParContainer()[i];
        for (size_t j=0; j<parsForSource.size(); ++j)
        {
           ResponseParam* rp = parsForSource[j];
           if (rp && !(rp->isFrozen() || rp->isLinked()))
           {
              // This puts the respFloor at (2^15 - 1)*2^16 = 2^31 - 2^16
              int respIndex = (s_RESPAR_INDEX << ModelContainer::SHIFT()) + k;
              // shallow copy
              variableParameters(respIndex, rp);
           }
           ++k;
        }
     }

  }
}

void Fit::calculateModel ()
{
  m_isStillValid = false;
  std::set<string>::iterator ss (m_activeModels.begin());
  std::set<string>::iterator ssEnd (m_activeModels.end());

  while ( ss != ssEnd)
  {
        models->calculate(*ss);
        models->fold(*ss);
	++ss;
  }         
}

void Fit::randomizeModelParameters (bool callInitialize, Real fSigma)
{
  RealArray parVals(.0, m_variableParameters.size());
  std::map<int,ModParam*>::iterator itVp = m_variableParameters.begin();
  std::map<int,ModParam*>::iterator itVpEnd = m_variableParameters.end();
  size_t i=0;
  while (itVp != itVpEnd)
  {
     parVals[i] = itVp->second->value('a');
     ++itVp, ++i;
  }  

  if (m_chainManager->isSynchedWithFitParams())
  {
     m_chainManager->getRandomPoint(parVals);
  }
  else
  {
     RandomizerBase* rand = 
                m_randomizingStrategies.find("gaussian fit")->second;
     if (callInitialize)
        rand->initializeRun(this);

     rand->randomize(parVals, this);
  }

  itVp = m_variableParameters.begin();
  i=0;
  while (itVp != itVpEnd)
  {
     itVp->second->setValue(parVals[i],'a');
     ++itVp, ++i;
  }

  // If fit was valid before, it sure isn't now.  Leave it up
  // to client to decide if it wants to return fit to valid state.
  m_isStillValid = false;      
}

void Fit::fluxes (bool lumin, bool error, int nTrials, Real level, Real eMin, Real eMax, Real redshift)
{
    using std::setw;
    static const Real LUMCON (1.07057E13);
    static std::pair<Real,Real> ZERO(0,0);

    // for each spectrum, and for each model that fits that 
    // spectrum, integrate.


    size_t nSpectra = datasets->numberOfSpectra();
    size_t nSources = datasets->numSourcesForSpectra();
    // The idea here is 2 modes of operation for flux/lumin command.  If ANY spectra
    // are found with ANY kind of response attached, flux will be calculated for those
    // spectra, using only those models that correspond to the attached responses.
    // Results are then stored with the SpectralData object.
    // If NO spectra are loaded, or no loaded spectra have any responses, flux/lumin
    // is calculated for active/off models using singleton dummy response, and result
    // stored with the model object.  
    std::vector<SpectralData*> specs;
    for (size_t i=0; i<nSpectra; ++i)
    {
       SpectralData* sd = datasets->lookup(i+1);
       for (size_t j=0; j<nSources; ++j)
       {
          if (sd->responseLoaded(j))
          {
             specs.push_back(sd);
             break;
          }
       }
    }
    if (specs.empty())  specs.push_back(0);


    std::ios_base::fmtflags currentSetting(tcout.flags());
    size_t savePrecision (tcout.precision());

    Numerics::FZSQ fzsq;
    Real H0 (models->cosmo().H0/100.);
    Real q0 (models->cosmo().q0);
    Real lambda0 = (models->cosmo().lambda0);
    if  ( H0 <= 0 ) 
    {
	H0 = 0.50;
	tcout << "  Warning: stored H0 value is zero, taken as 50 km/s/Mpc " << std::endl;
    }

    std::map<size_t,RealArray> prevEngsForSpec;
    size_t prevDataGroupNum = 0;

    std::vector<std::vector<Model*> > modsForSpectra;
    std::vector<SpectralData*>::const_iterator
	begSpecs = specs.begin(),
	endSpecs = specs.end();

    while ( begSpecs != endSpecs)
    {
       SpectralData* sd = *begSpecs; 
       size_t spectrumNumber = 0;
       size_t currDataGroupNum = 0;
       std::vector<SpectralData::FluxCalc> modelFluxCalc;
       bool diffFound = false;
       std::map<size_t,const RealArray*> engsForSpec;
       modsForSpectra.push_back(getModsForFlux(sd));
       const std::vector<Model*>& modsForSpec = modsForSpectra.back();       
       if (sd)
       {
          spectrumNumber = sd->spectrumNumber();
          currDataGroupNum = sd->parent()->dataGroup();
       }
       for (size_t i=0; i<modsForSpec.size(); ++i)
       {
          Model* mod = modsForSpec[i];
          ArrayContainer::const_iterator itEng = mod->energy().find(spectrumNumber);
          if (itEng != mod->energy().end())
          {
             // Should always get in here.
             engsForSpec[mod->sourceNumber()] = &itEng->second;
          }
       }

       size_t sz = engsForSpec.size();
       if (sz == prevEngsForSpec.size())
       {
          std::map<size_t,RealArray>::const_iterator itPrev = prevEngsForSpec.begin();
          std::map<size_t,RealArray>::const_iterator itPrevEnd = prevEngsForSpec.end();
          std::map<size_t,const RealArray*>::const_iterator itEngsEnd = engsForSpec.end();
          while (itPrev != itPrevEnd && !diffFound)
          {
             size_t iSource = itPrev->first;
             std::map<size_t,const RealArray*>::const_iterator itEngs = 
                        engsForSpec.find(iSource);
             if (itEngs != itEngsEnd)
             {
                const RealArray& eng_iSource = *itEngs->second;
                const RealArray& prevEng_iSource = itPrev->second;
                const size_t engSize = eng_iSource.size();
                if (engSize != prevEng_iSource.size())
                {
                   diffFound = true;
                }
                else
                {
                   for (size_t i=0; i<engSize; ++i)
                   {
                      if (eng_iSource[i] != prevEng_iSource[i])
                      {
                         diffFound = true;
                         break;
                      }
                   }
                }
             }
             else
             {
                diffFound = true;
             }
             ++itPrev;
          }
       }
       else
       {
	   diffFound = true;
       }
       if (diffFound)
       {
	   prevEngsForSpec.clear();
           std::map<size_t,const RealArray*>::const_iterator itEngs =
                engsForSpec.begin();
           std::map<size_t,const RealArray*>::const_iterator itEngsEnd =
                engsForSpec.end();
           while (itEngs != itEngsEnd)
           {
              prevEngsForSpec.insert(std::map<size_t,RealArray>::value_type
                        (itEngs->first, *itEngs->second));
              ++itEngs;
           }
       }

       const bool showOutput = diffFound || (currDataGroupNum != prevDataGroupNum);
       prevDataGroupNum = currDataGroupNum;
       bool isRangeError = false;
       for (size_t i=0; i<sz && !isRangeError; ++i)
       {
	   Model* m = modsForSpec[i];
	   Real energyLow (eMin);
	   Real energyHigh (eMax);
	   const RealArray& energyArray = m->energy().find(spectrumNumber)->second;
	   const size_t N (energyArray.size());
	   const Real& arrayMax = energyArray[N-1];
	   const Real& arrayMin = energyArray[0];

	   if ( energyLow > arrayMax || energyHigh < arrayMin )
	   {
	       if (showOutput)
	       {
		   tcout << " No overlap between matrix range  ( "
			 << setw(12) << arrayMin << ',' << setw(12)
			 << arrayMax << " ) \nand the requested range "
			 << "( " << setw(12) << energyLow << ',' 
			 << setw(12) << energyHigh << " )"<<std::endl;
	       }
	       isRangeError = true;
               break;
	   }

	   if ( energyLow < arrayMin ) 
	   {
	       if (showOutput)
	       {
		   tcout << " Lower range bound " << setw(12) << energyLow
			 << " reset by matrix bound to " << setw(12)
			 << arrayMin << std::endl;
	       }
	       energyLow = arrayMin;
	   } 

	   if ( energyHigh > arrayMax ) 
	   {
	       if (showOutput)
	       {
		   tcout << " Upper range bound " << setw(12) << energyHigh
			 << " reset by matrix bound to " << setw(12)
			 << arrayMax << std::endl;
	       }
	       energyHigh = arrayMax;
	   } 

	   Real kFlux (0);
	   Real eFlux (0);

	   m->integrateFlux(spectrumNumber,energyLow, energyHigh, kFlux, eFlux);
	   m->keVFlux(spectrumNumber,kFlux);
	   if (lumin)
	   {
	       Real lumNorm  = LUMCON*fzsq(redshift,q0,lambda0)/H0/H0;
	       eFlux *= lumNorm;
	   }
	   m->ergFlux(spectrumNumber,eFlux);

	   modelFluxCalc.push_back(SpectralData::FluxCalc(eFlux,.0,.0,kFlux));
           // If fit is not valid, don't go into the error block later on
           // regardless if error flag is set to true.
	   if (!error || !m_isStillValid)
	   {
	       if ( datasets->numberOfSpectra() > 1 && showOutput)
	       {
		   tcout << " Spectrum Number: " << spectrumNumber << std::endl;
                   tcout << " Data Group Number: " << currDataGroupNum << std::endl;
	       }
	       m->keVFluxRange(spectrumNumber,ZERO);
	       m->ergFluxRange(spectrumNumber,ZERO);
	       if (showOutput)
	       {
		   m->reportFluxes(spectrumNumber,1.+redshift,
				   lumin,energyLow,energyHigh);
	       }
	   }                               
       } // end mods for spectrum loop

       if (!isRangeError)
       {
          if(spectrumNumber)
          {
	          SpectralData* spec = datasets->lookup(spectrumNumber);

	          if (lumin)
		      spec->setLastModelLuminCalc(modelFluxCalc);
	          else
		      spec->setLastModelFluxCalc(modelFluxCalc);
          }
          else
          {
	      for(size_t i = 0; i < sz; ++i)
	      {
	          Model* m = modsForSpec[i];

	          if(lumin)
		      m->lastModelLuminCalc(modelFluxCalc[i]);
	          else
		      m->lastModelFluxCalc(modelFluxCalc[i]);
	      }
          }
       }

       ++begSpecs;
    } // end spectra loop

    // the error loop has to be broken out ( cf. xspec11) because
    // the randomizeModelParameters acts on all of the models.

    if ( error && m_isStillValid)
    {
	// one vector per trial containing a vector of
	// size #spectra in fit containing a vector of # sources
	// being simultaneously fit.
	std::vector<std::vector< std::vector<Real> > > evFF(nTrials);
	std::vector<std::vector< std::vector<Real> > > ergFF(nTrials);    
	size_t N (m_variableParameters.size());
	size_t nFluxSpecs = specs.size();
	std::vector<Real> saveParameterValues(N);
	std::map<int,ModParam*>::const_iterator vp 
	    = m_variableParameters.begin();
	for ( size_t k = 0; k < N; ++k, ++vp)
	{
	    saveParameterValues[k] =    vp->second->value('a');    
	}

        reportErrorSimulationMethod();
        if (m_chainManager->isSynchedWithFitParams())
        {
           const std::vector<size_t>& chainLengths = m_chainManager->accumulatedLengths();
           size_t totLength = chainLengths[chainLengths.size()-1];
           if (static_cast<size_t>(nTrials) > totLength)
           {
              tcout<<"***Warning: The number of error samples: "<< nTrials 
                 <<"\n    is larger than the total lengths of loaded chains: "<< totLength
                 << std::endl;
           }
        }
	for (int j = 0; j < nTrials; ++j)
	{
	    std::vector<std::vector<Real> >& evf = evFF[j];  
	    std::vector<std::vector<Real> >& ergf = ergFF[j];  

	    evf.resize(nFluxSpecs);  
	    ergf.resize(nFluxSpecs);

            bool firstTime = !j;               
	    randomizeModelParameters(firstTime);

	    calculateModel();

	    // iterate through datasets/models  again.
	    for (size_t iSp=0; iSp<nFluxSpecs; ++iSp)
	    {
	       size_t spectrumNumber = specs[iSp] ?
                           specs[iSp]->spectrumNumber() : 0;
               const std::vector<Model*> modsForSpec = modsForSpectra[iSp];
               const size_t nMods = modsForSpec.size();
	       std::vector<Real>& evff = evf[iSp];
	       std::vector<Real>& ergff = ergf[iSp];
	       evff.resize(nMods,0);
	       ergff.resize(nMods,0);
               for (size_t iMod=0; iMod<nMods; ++iMod)
	       {
		   Model* m = modsForSpec[iMod];

		   Real& evt = evff[iMod];
		   Real& ergt = ergff[iMod];
		   ArrayContainer::const_iterator 
		       s (m->energy().find(spectrumNumber));
		   ArrayContainer::const_iterator sEnd (m->energy().end());
		  if ( s != sEnd)
		  {
		     const RealArray& energyArray = s->second;
		     const size_t N (energyArray.size());
		     const Real& arrayMax = energyArray[N-1];
		     const Real& arrayMin = energyArray[0];
		     Real energyLow (eMin);
		     Real energyHigh (eMax);
		     if ( energyLow > arrayMax ||
		          energyHigh < arrayMin ) break;
		     energyLow = std::max(energyLow,arrayMin);
		     energyHigh = std::min(energyHigh,arrayMax);
		     m->integrateFlux(spectrumNumber,energyLow, 
				      energyHigh, evt, ergt);
		     if (lumin)
		     {
		        Real lumNorm  
			    = LUMCON*fzsq(redshift,q0,lambda0)/H0/H0;
		        ergt *= lumNorm;
		     }
		  }
	       }  // foreach model fitting the spectrum
	    } // foreach spectrum

	    vp = m_variableParameters.begin();
	    for ( size_t k = 0; k < N; ++k, ++vp)
	    {
		vp->second->setValue
		    (saveParameterValues[k],'a');       
	    }
	} // foreach trial

	// collect results and find confidence ranges
        try
        {
	   for (size_t iSp=0; iSp<nFluxSpecs; ++iSp)
	   {
	       size_t spectrumNumber = specs[iSp] ?
                           specs[iSp]->spectrumNumber() : 0;
               const std::vector<Model*>& modsForSpec = modsForSpectra[iSp];
               const size_t nMods = modsForSpec.size();
	       for (size_t iMod=0; iMod<nMods; ++iMod)
	       {
	          Model* m = modsForSpec[iMod];
	          ArrayContainer::const_iterator s 
		      (m->energy().find(spectrumNumber));
	          ArrayContainer::const_iterator sEnd (m->energy().end());
	          if ( s != sEnd )
	          {
		      const RealArray& energyArray = s->second;
		      std::vector<Real> teV(nTrials);
		      std::vector<Real> terg(nTrials);
		      for (int j = 0; j < nTrials; ++j)
		      {
		          teV[j] = evFF[j][iSp][iMod];
		          terg[j] = ergFF[j][iSp][iMod];
		      } 
                      // This can throw.
		      std::pair<Real,Real> ceVRange = 
		          XSutility::confidenceRange(level,teV);
		      std::pair<Real,Real> cergRange  =
		          XSutility::confidenceRange(level,terg);
		      m->keVFluxRange(spectrumNumber,ceVRange);
		      m->ergFluxRange(spectrumNumber,cergRange);

		      SpectralData::FluxCalc tmp; 
		      tmp.errLow = cergRange.first;
		      tmp.errHigh = cergRange.second;
                      tmp.photonLow = ceVRange.first;
                      tmp.photonHigh = ceVRange.second;

		      if(spectrumNumber)
		      {
		          SpectralData* spec = specs[iSp];

		          if (lumin)
		          {
			      tmp.value = spec->lastModelLuminCalc(iMod).value;
                              tmp.photonValue = spec->lastModelLuminCalc(iMod).photonValue;
			      spec->lastModelLuminCalc(iMod,tmp);
		          }
		          else
		          {
			      tmp.value = spec->lastModelFluxCalc(iMod).value;
                              tmp.photonValue = spec->lastModelFluxCalc(iMod).photonValue;
			      spec->lastModelFluxCalc(iMod,tmp);
		          }
		      }
		      else
		      {
		          if(lumin)
		          {
			      tmp.value = m->lastModelLuminCalc().value;
			      tmp.photonValue = m->lastModelLuminCalc().photonValue;
			      m->lastModelLuminCalc(tmp);
		          }
		          else
		          {
			      tmp.value = m->lastModelFluxCalc().value;
			      tmp.photonValue = m->lastModelFluxCalc().photonValue;
			      m->lastModelFluxCalc(tmp);
		          }
		      }

		      if ( spectrumNumber && datasets->numberOfSpectra() > 1)
		      {
		          tcout << " Spectrum Number: " 
			        << spectrumNumber << std::endl;
		      }
		      const Real& arrayMax = energyArray[energyArray.size()-1];
		      const Real& arrayMin = energyArray[0];
		      Real energyLow (eMin);
		      Real energyHigh (eMax);
		      if ( energyLow > arrayMax ||
		           energyHigh < arrayMin ) break;
		      energyLow = std::max(energyLow,arrayMin);
		      energyHigh = std::min(energyHigh,arrayMax);
		      m->reportFluxes(spectrumNumber, 1.+redshift, lumin,
				      energyLow, energyHigh, level);
		      tcout << std::flush;  
	          }
	       }       
	   } // end spectra loop
        }
        catch (YellowAlert&)
        {
           calculateModel();
           m_isStillValid = true;
           throw;
        }
        // restore model
        calculateModel();
        m_isStillValid = true;
    } // end if error

    if ( lumin )
    {
	tcout.precision(3);
	tcout << "     (z = " << std::showpoint << setw(5) << redshift ;
	tcout << " H0 = "  << setw(5) << H0*100. << " q0 = "  << setw(5) << q0;
	if ( lambda0 > 0 ) tcout << " Lambda0 = " << setw(5) << lambda0 << ")" << std::endl;
    }

    tcout.flags(currentSetting);
    tcout.precision(savePrecision);
}

bool Fit::getErrors (const IntegerArray& paramNums, const string& modelName)
{
  FitErrorCalc::setAllErrParamStrings(this, paramNums, FitErrorCalc::OK);
  if (!m_fitMethod->getErrors(this, modelName, paramNums))
  {
     bool doAnotherCalc = true;
     while (doAnotherCalc)
     {
        FitErrorCalc  errCalcObj(this, modelName, paramNums);
        try
        {
           errCalcObj.perform();
           doAnotherCalc = false;
        }
        catch (FitErrorCalc::NewMinFound)
        {
           if (m_queryMode == ON)
           {
              if (XSutility::yesToQuestion("Re-start error calc with new best fit parameters? ",
                        0, tcin) <= 0)
              {
                 doAnotherCalc = false;
              }
           }
           else if (m_queryMode == NO)
           {
              doAnotherCalc = false;
           }
           // Need to replace m_oldParameters with new best fit values:
           reinitialize(true);
        }
        catch (...)
        {
           throw;
        }
     }
  }
  return true;
}

Step* Fit::stepGrid ()
{
  return m_stepGrid.get();
}

void Fit::stepGrid (Step* value)
{
  m_stepGrid.reset(value);
}

void Fit::initializeStatistic (bool isUpdate)
{

    int degrees(0);
    size_t bins(0);
    m_statMethod->degreesOfFreedom(fit,degrees,bins);
    if ( degrees < 0) 
    {
            throw FitOverDetermined();       
    }


    if (m_renormType == AUTO || (m_renormType == PREFIT && !isUpdate)) 
    {
            // can throw "can't initialize".
            m_statMethod->renormalize(fit);       
    }
    // this has the byproduct of computing the initial
    // model difference array, which we need.
    m_statMethod->perform(fit);
}

void Fit::freezeByIndex (int index)
{

    // std::map<K,V>::erase returns 1 if erasure was successful.
    if ( !m_variableParameters.erase(index))
    {
            tcerr << "*** Warning: attempt to freeze non-fit parameter "<<std::endl;       
    }      
}

Real Fit::statistic () const
{

  return m_statMethod->statistic();
}

void Fit::perform ()
{
  m_fitMethod->perform(this);
}

void Fit::calcEqErrors (const XSContainer::EqWidthRecord& currComp, Model* mod, size_t nIter, Real confLevel)
{
      using XSContainer::datasets;
      // Calculate errors based on gaussian distribution of parameters.
      if (!m_isStillValid)
      {
         tcout <<"\nCannot determine error of equiv width without a valid fit."<<std::endl;
         return;
      }
      if (!nIter)
      {
         return;
      }
      RealArray eqWidths(.0, nIter);
      size_t nParams = m_variableParameters.size();
      std::vector<Real> saveParameterValues(nParams);
      std::map<int,ModParam*>::const_iterator vp = m_variableParameters.begin();
      for (size_t i=0; i<nParams; ++i, ++vp)
      {
         saveParameterValues[i] = vp->second->value();
      }
      int savConVerbose = tpout.conVerbose();
      int savLogVerbose = tpout.logVerbose();
      Real savLastEqWidth = .0;
      // Find the first spectrum in the first Dataset that matches mod's
      // data group.  All lastEqWidths for this group should be the same.
      DataArrayIt itDs = datasets->dataArray().begin();
      DataArrayIt itDsEnd = datasets->dataArray().end();
      while (itDs != itDsEnd)
      {
         DataSet* ds = itDs->second;
         if (ds->dataGroup() == mod->dataGroupNumber())
         {
            if (ds->isMultiple())
            {
               savLastEqWidth = ds->multiSpectralData().begin()->second->
                        lastEqWidthCalc().value;
            }
            else
            {
               savLastEqWidth = ds->spectralData()->lastEqWidthCalc().value;
            }
            break;
        } 
         ++itDs;
      }

      if (m_chainManager->isSynchedWithFitParams())
      {
         const std::vector<size_t>& chainLengths = m_chainManager->accumulatedLengths();
         size_t totLength = chainLengths[chainLengths.size()-1];
         if (nIter > totLength)
         {
            tcout<<"***Warning: The number of error samples: "<< nIter 
               <<"\n    is larger than the total lengths of loaded chains: "<< totLength
               << std::endl;
         }
      }

      try
      {
         for (size_t i=0; i<nIter; ++i)
         {
            bool firstTime = !i;
            randomizeModelParameters(firstTime);
            calculateModel();

            // debug info
      /*      std::map<int,ModParam*>::const_iterator ip = m_variableParameters.begin();
            while (ip != fit->variableParameters().end())
            {
               tcout << ip->second->name() << "  " << ip->second->value() << std::endl;
               ++ip;
            }
            m_statMethod->perform(fit);
            tcout <<"Chisq: " <<statistic() <<std::endl;
     */
            // Do not print the various eqwidth output messages
            // while merely building up an errors array.
            XSstream::verbose(tpout, 14, 14);
            eqWidths[i] = XSContainer::models->calcEqWidths(currComp, mod);
            vp = m_variableParameters.begin();
            for (size_t i=0; i<nParams; ++i, ++vp)
            {
               vp->second->setValue(saveParameterValues[i]);
            }
            // reset chatter levels
            XSstream::verbose(tpout, savConVerbose, savLogVerbose);
         }
      }
      catch (YellowAlert&)
      {
         // calcEqWidths changes what is stored in lastEqWidthCalc, so
         // it needs to be restored.
         itDs = datasets->dataArray().begin();
         SpectralData::FluxCalc eqWidthStruct(savLastEqWidth, 0., 0.);
         while (itDs != itDsEnd)
         {
            DataSet* ds = itDs->second;
            if (ds->dataGroup() == mod->dataGroupNumber())
            {
               if (ds->isMultiple())
               {
                  SpectralDataMapConstIt itSd = ds->multiSpectralData().begin();
                  SpectralDataMapConstIt itSdEnd = ds->multiSpectralData().end();
                  while (itSd != itSdEnd)
                  {
                     itSd->second->lastEqWidthCalc(eqWidthStruct);
                     ++itSd;
                  }
               }
               else
               {
                  SpectralData* sd = ds->spectralData();
                  if (sd)
                  {
                     sd->lastEqWidthCalc(eqWidthStruct);
                  }
               }
            }
            ++itDs;
         }

         vp = m_variableParameters.begin();
         for (size_t i=0; i<nParams; ++i, ++vp)
         {
            vp->second->setValue(saveParameterValues[i]);
         }
         calculateModel();
         m_isStillValid = true;
         XSstream::verbose(tpout);
         throw;
      }
      // Note: fit statistic hasn't been changed during the errors
      // loop, so it doesn't need to be recalculated.
      calculateModel();
      m_isStillValid = true;

      std::sort(&eqWidths[0], &eqWidths[0] + nIter);
      size_t lowErr = static_cast<size_t>(.5*(nIter-1)*(1.0-confLevel/100.0)+.5);
      size_t highErr = static_cast<size_t>(.5*(nIter-1)*(1.0+confLevel/100.0)+.5);
      if (nIter < static_cast<size_t>(s_MCrealizations))
      {
         tcout << "\n***Warning: The number of simulations for error determination is low."
              <<std::endl;
      }
      tcout << "Equiv width error range:  " << eqWidths[lowErr] << " - "
            << eqWidths[highErr] << " keV" << std::endl;

      itDs = datasets->dataArray().begin();
      SpectralData::FluxCalc eqWidthStruct(savLastEqWidth, eqWidths[lowErr],
                                eqWidths[highErr]);
      while (itDs != itDsEnd)
      {
         DataSet* ds = itDs->second;
         if (ds->dataGroup() == mod->dataGroupNumber())
         {
            if (ds->isMultiple())
            {
               SpectralDataMapConstIt itSd = ds->multiSpectralData().begin();
               SpectralDataMapConstIt itSdEnd = ds->multiSpectralData().end();
               while (itSd != itSdEnd)
               {
                  itSd->second->lastEqWidthCalc(eqWidthStruct);
                  ++itSd;
               }
            }
            else
            {
               SpectralData* sd = ds->spectralData();
               if (sd)
               {
                  sd->lastEqWidthCalc(eqWidthStruct);
               }
            }
         }
         ++itDs;
      }
}

std::vector<Model*> Fit::getModsForFlux (const SpectralData* sd)
{
   // If sd exists, models vector will return with mods in order
   // sorted by sourceNum.  Else, models vector will contain
   // all active/off mods sorted alphabetically.
   std::vector<Model*> mods;
   if (sd)
   {
      const size_t dataGroup = sd->parent()->dataGroup();
      const std::vector<Response*>& detectors = sd->detector();
      for (size_t i=0; i<detectors.size(); ++i)
      {
         if (detectors[i])
         {
            const string modName(models->lookupModelForSource(i+1));
            if (modName.length())
            {
               Model* mod = models->lookup(modName, dataGroup);
               mods.push_back(mod);
            }
         }
      }
   }
   else
   {
      std::map<string,bool>::const_iterator itModName =
                models->activeModelNames().begin();
      std::map<string,bool>::const_iterator itModNameEnd =
                models->activeModelNames().end();
      while (itModName != itModNameEnd)
      {
         if (itModName->second)
         {
            Model* mod = models->lookup(itModName->first);
            mods.push_back(mod);
         }
         ++itModName;
      }
   }
   return mods;
}

const RealArray& Fit::getSVDevalue () const
{
   const FitMethod* cfitMethod = 
        const_cast<const FitMethod*>(m_fitMethod.get());
   return cfitMethod->evalue();
}

const RealArray& Fit::getSVDevector () const
{
   const FitMethod* cfitMethod = 
        const_cast<const FitMethod*>(m_fitMethod.get());
   return cfitMethod->evector();
}

const RealArray& Fit::getSVDcovariance () const
{
   const FitMethod* cfitMethod = 
        const_cast<const FitMethod*>(m_fitMethod.get());
   return cfitMethod->covariance();
}

bool Fit::checkChainsForSynch () const
{
  // Naturally this should be called only after variableParameters
  // map has been updated.
  // The case of no chains loaded is defined as synch = false.  This
  // makes ChainManager's isSynchedWithFitParams flag the only thing
  // that needs to be checked to determine whether to get distribution
  // from covariance or from chains.
   bool areAllParsMatched = false;
   if (m_chainManager->chains().size())
   {
      // All loaded chains must have the same pars, so just
      // grap the info from the first chain.
      const std::vector<Chain::ParamID>& chainPars = 
                m_chainManager->chains().begin()->second->paramIDs();
      const size_t nPars = m_variableParameters.size();
      if (chainPars.size() == nPars)
      {
         areAllParsMatched = true;
         std::vector<Chain::ParamID>::const_iterator itChainPar = 
                   chainPars.begin();
         std::map<int,ModParam*>::const_iterator itVar = 
                   m_variableParameters.begin();
         std::map<int,ModParam*>::const_iterator itVarEnd = 
                   m_variableParameters.end();
         while (itVar != itVarEnd)
         {
            const Chain::ParamID& chainPar = *itChainPar;
            const ModParam* fitPar = itVar->second;
            // Note that passing the test below does NOT guarantee that 
            // parameters in the chain file were actually generated from 
            // those currently stored in the variableParameters map.  
            if (chainPar.modName != fitPar->modelName() ||
                chainPar.parName != fitPar->name() ||
                chainPar.index != fitPar->index() ||
                chainPar.units != fitPar->unit())
            {
               areAllParsMatched = false;
               break;
            }
            ++itVar, ++itChainPar;
         }
      }
   }
   return areAllParsMatched;
}

void Fit::reportErrorSimulationMethod () const
{
   tcout << "Parameter distribution is derived from ";
   if (m_chainManager->isSynchedWithFitParams())
      tcout << "loaded chain files." << std::endl;
   else
      tcout <<"fit covariance matrix." << std::endl;
}

void Fit::randomizeForChain (bool onlyInit)
{
  RandomizerBase* rand = m_chainManager->chainProposal();
  if (onlyInit)
  {
     rand->initializeRun(this);
  }
  else
  {
     RealArray parVals(.0, m_variableParameters.size());
     std::map<int,ModParam*>::iterator itVp = m_variableParameters.begin();
     std::map<int,ModParam*>::iterator itVpEnd = m_variableParameters.end();
     size_t i=0;
     while (itVp != itVpEnd)
     {
        parVals[i] = itVp->second->value('a');
        ++itVp, ++i;
     }  

     rand->randomize(parVals, this);
     itVp = m_variableParameters.begin();
     i=0;
     while (itVp != itVpEnd)
     {
        itVp->second->setValue(parVals[i],'a');
        ++itVp, ++i;
     }

     // If fit was valid before, it sure isn't now.  Leave it up
     // to client to decide if it wants to return fit to valid state.
     m_isStillValid = false; 
  }     
}

string Fit::getChainProposalNames () const
{
   string names;
   std::map<string,RandomizerBase*>::const_iterator itRand = m_randomizingStrategies.begin();
   std::map<string,RandomizerBase*>::const_iterator itEnd = m_randomizingStrategies.end();
   while (itRand != itEnd)
   {
      RandomizerBase* rand = itRand->second;
      if (!names.empty())
         names += " | ";
      string::size_type cmdLoc = rand->name().find("<cmdline>");
      if (cmdLoc != string::npos)
      {
         // Give a more precise description.
         // ie. "<dist> <cmdline>" becomes  
         // "<dist> diagonal <values> | <dist> matrix <values>"
         const string distName(rand->name().substr(0, cmdLoc));
         names += distName + "diagonal <values> | " + distName + "matrix <values>";
      }
      else
        names += rand->name();
      ++itRand;
   }
   return names;
}

void Fit::registerRandomizingStrategy (const string& name, RandomizerBase* strategy)
{
   const string ws(" \t\n");
   if (name.find_first_not_of(ws) == string::npos)
   {
      throw YellowAlert("No name provided for add-on randomizer class.\n");
   }
   // Assuming that the list of reserved names in m_nativeRandomizerNames
   // will be filled AFTER all native strategies have been entered.  
   // Therefore this check should really only ever apply to user-entered 
   // strategies.
   if (m_nativeRandomizerNames.size())
   {
      // User classes should only have one word names, no whitespace.
      if (name.find_first_of(ws) != string::npos)
      {
         throw YellowAlert("Whitespace is not allowed in name member of add-on randomizer classes\n");
      }
      if (m_nativeRandomizerNames.find(name) != m_nativeRandomizerNames.end())
      {
         string errMsg("Cannot add randomizing strategy named ");
         errMsg += name;
         errMsg += ".\nThe name is reserved for built-in randomizing options.\n";
         throw YellowAlert(errMsg);
      }
   }

   std::map<string,RandomizerBase*>::iterator it = m_randomizingStrategies.find(name);
   if (it != m_randomizingStrategies.end())
   {
      tcout <<"***Warning: An add-on randomizing strategy named "<< name 
        <<" already exists.  It will be replaced." <<std::endl;
      delete it->second;
      m_randomizingStrategies.erase(it);
   }
   m_randomizingStrategies[name] = strategy;
}

RandomizerBase* Fit::getRandomizingStrategy (const string& name)
{
   RandomizerBase* pStrategy = 0;
   // If strategy doesn't exist, simply return a null pointer and
   // let calling function decide if it wants to issue a message.
   // This allows for abbreviated names to match.
   std::map<string,RandomizerBase*>::iterator it =
                m_randomizingStrategies.lower_bound(name);
   if (it != m_randomizingStrategies.end() && it->first.find(name) == 0)
      pStrategy = it->second;
   return pStrategy;
}

void Fit::clearRandomizingStrategies ()
{
   std::map<string,RandomizerBase*>::iterator itRand = m_randomizingStrategies.begin();
   std::map<string,RandomizerBase*>::iterator itRandEnd = m_randomizingStrategies.end();
   while (itRand != itRandEnd)
   {
      delete itRand->second;
      ++itRand;
   }
   m_randomizingStrategies.clear();
}

void Fit::initNativeRandomizerNames ()
{
   // This should only be called once, during start-up.  If it detects
   // existing entries it will do nothing.  Only the FIRST part of
   // compound names are stored.
   if (m_nativeRandomizerNames.empty())
   {
      std::map<string,RandomizerBase*>::const_iterator itRand = 
                m_randomizingStrategies.begin();
      std::map<string,RandomizerBase*>::const_iterator itEnd =
                m_randomizingStrategies.end();
      while (itRand != itEnd)
      {
         string::size_type firstWS = itRand->first.find(' ');
         string primaryName(itRand->first.substr(0, firstWS));
         m_nativeRandomizerNames.insert(primaryName);
         ++itRand;
      }
   }
}

void Fit::reportChainAcceptance (bool isAccepted) const
{
   RealArray parVals(.0, m_variableParameters.size());
   std::map<int,ModParam*>::const_iterator itVp = m_variableParameters.begin();
   std::map<int,ModParam*>::const_iterator itVpEnd = m_variableParameters.end();
   size_t i=0;
   while (itVp != itVpEnd)
   {
      parVals[i] = itVp->second->value();
      ++itVp, ++i;
   }  
   m_chainManager->chainProposal()->acceptedRejected(parVals, isAccepted);
}

ModParam* Fit::oldParameterValues (int index) const
{
  std::map<int,ModParam*>::const_iterator p (m_oldParameterValues.find(index));
  if (p != m_oldParameterValues.end()) return p->second;
  else return 0;
}

ModParam* Fit::variableParameters (int index) const
{
  std::map<int,ModParam*>::const_iterator p (m_variableParameters.find(index));
  if (p != m_variableParameters.end()) return p->second;
  else return 0;
}

// Additional Declarations