//
//  XSPEC12  November 2003
//
//

#include <XSUser/Global/Global.h>  
#include <XSUser/Global/XSGlobal.h>  
#include <XSFit/Fit/Fit.h>
#include <XSFit/Fit/FitErrorCalc.h>
#include <XSFit/Fit/StatMethod.h>
#include <XSFit/Fit/FitMethod.h>
#include <XSFit/MCMC/ChainManager.h>
#include <XSModel/GlobalContainer/ModelContainer.h>
#include <XSUtil/Parse/XSparse.h>   
#include <XSUtil/Signals/SignalHandler.h>   

#include <XSContainer.h>
#include <XSstreams.h>
#include <sstream>
#include <iomanip>

namespace
{
   void reportParameter(const ModParam* par, bool isFitValid);
}

int
XSGlobal::xsError(ClientData cdata,Tcl_Interp* tclInterp,int objc, Tcl_Obj* CONST  objv[])
{
   using XSContainer::fit;
   if (objc == 2 && string(Tcl_GetString(objv[1])) == "?") 
   {
           printDocs(static_cast<char*>(cdata),"?");
           return globalData->autoSave(TCL_OK);
   }

   const string ss = "sS";
   const string mm = "mM";
   const string ir = ".Ee";

   int status (globalData->autoSave(TCL_OK));
   // store i/o flags.

   StringArray inArgs;
   StringArray params;
   IntegerArray iParams;
   for (int i=1; i<objc; ++i)
   {
      inArgs.push_back(Tcl_GetString(objv[i]));
   }
   
   static IntegerArray savErrorParams(1,1);
   static string savModName("");
   int tmpNtrial = fit->errorTry();
   Real tmpTolerance = fit->tolerance();
   Real tmpChiMax = fit->chiMax();
   Real tmpDeltaStat = fit->deltaStat();
   string tmpModName = savModName;
   IntegerArray tmpErrorParams(savErrorParams);
      
   XSparse::collectParams(inArgs, iParams, params);
   size_t nPars = iParams.size();
   bool parseError = false;
   bool stopatProcessed = false;
   bool deltaProcessed = false;
   bool isFirstParRange = true;
   for (size_t i=0; i<nPars && !parseError; ++i)
   {
      const string& parString = params[i];
      // qualified range here means a parameter range of form  
      // modelName:{range string}. We need to do this to allow
      // model names beginning with 's' or 'S'.
      bool isQualified = (parString.find_first_of(':') != string::npos);
      bool isStopat = (!isQualified && ss.find(parString[0]) != string::npos);
      bool isMaximum = (!isQualified && mm.find(parString[0]) != string::npos);
      // Foremost rule: if stopat appears, it better be the first param.
      if (isStopat)
      {
         if (iParams[i] != 0)
         {
            tcout << "stopat option must be specified as the first parameter."
              << std::endl;
            parseError = true;
         }
         else
         {
            // Process stopat w/ possible pars 2 and 3
            size_t itmp = i+1;
            while (itmp < nPars && iParams[itmp] < 3 && !parseError)
            {
               std::istringstream testIss(params[itmp]);
               if (iParams[itmp] == 1)
               {
                  // numTrials
                  if (!(testIss >> tmpNtrial) || !testIss.eof() || tmpNtrial < 1)
                  {
                     tcout << "Improper value for nTrial parameter." << std::endl;
                     parseError = true;
                  }
               }
               else if (iParams[itmp] == 2)
               {
                  // tolerance
                  if (!(testIss >> tmpTolerance) || !testIss.eof() || tmpTolerance <= .0)
                  {
                     tcout << "Improper value for tolerance parameter." << std::endl;
                     parseError = true;
                  }
               }
               // i will get incremented again at end of outer for loop, 
               // so it should be lagging itmp by 1 here.
               ++i;
               ++itmp;
            }
            stopatProcessed = true;
         }
      } // end if isStopat
      else if (isMaximum)
      {
         if ((stopatProcessed && iParams[i] != 3) || (!stopatProcessed &&
               iParams[i] != 0))
         {
            tcout << "Improper location of \"maximum\" option specifier."
               << std::endl;
            parseError = true;
         }
         else
         {
            if (i+1 < nPars && iParams[i+1] == iParams[i]+1)
            {
               // chi max
               ++i;
               std::istringstream testIss(params[i]);
               if (!(testIss >> tmpChiMax) || !testIss.eof() || tmpChiMax <=.0)
               {
                  tcout << "Improper maximum chi-square value." << std::endl;
                  parseError = true;
               }
            }
         }
      } // end if isMaximum
      else
      {
         // Only constraint here is that delta stat (identified
         // as a Real) must appear just once, and it must 
         // appear before ANY parameter range indicators.
         size_t testInt = XSutility::isInteger(parString);
         Real testReal=.0;
         if (testInt != string::npos)
         {
            // testInt must be positive if it makes it here.
            // Assume it must be trying to give a parameter range.
            if (isFirstParRange)
            {
               tmpErrorParams.clear();
               isFirstParRange = false;
               tmpModName.erase();
            }
            if (testInt == 0)
            {
               tcout << "Improper parameter specifier: 0" << std::endl;
               parseError = true;
            }
            else
            {
               tmpErrorParams.push_back(static_cast<int>(testInt));
            }
         }
         else if (XSutility::isReal(parString, testReal) && isFirstParRange)
         {
            if (deltaProcessed)
            {
               tcout << "Improper parameter range specifier: " 
                  << testReal << std::endl;
               parseError = true;
            }
            else if (testReal <= .0)
            {
               tcout << "Improper delta fit stat value: " << testReal
                 << std::endl;
               parseError = true;
            }
            else
            {
               tmpDeltaStat = testReal;
               deltaProcessed = true;
            }
         }
         else
         {
            try
            {
               IntegerArray paramRange;
               // Current implementation only allows 1 setting
               // of modName.  Use from first range specifier. 
               string testModName;
               // This can throw
               XSparse::stringRangePair(parString, testModName, paramRange);
               if (isFirstParRange)
               {
                  tmpErrorParams.clear();
                  tmpModName = testModName;
                  isFirstParRange = false;
               }
               for (size_t j=0; j<paramRange.size(); ++j)
               {
                  if (paramRange[j] <= 0)
                  {
                     tcout << "Improper parameter specifier: " 
                        << paramRange[j] << std::endl;
                     parseError = true;
                  }
                  tmpErrorParams.push_back(paramRange[j]);
               }
            }
            catch (YellowAlert&)
            {
               parseError = true;
            }
            
         }
      }
   } // end parameters for loop
   if (parseError)
   {
      printDocs("error","?");
      return globalData->autoSave(TCL_ERROR);           
   }
   else
   {
      fit->errorTry() = tmpNtrial;
      fit->tolerance() = tmpTolerance;
      fit->chiMax() = tmpChiMax;
      fit->deltaStat() = tmpDeltaStat;
      savModName = tmpModName;
      savErrorParams = tmpErrorParams;   
   }
   
   // Parsing syntax is validated by this point, though still not
   // home free.  Check if model parameters are valid.
   const size_t nSpecifiedPar = savErrorParams.size();
   // These two arrays will contain a subset of savErrorParams.
   // Errors from chains need just the short par index nums (which is
   // what the user sees), while fit errors can perform easier 
   // lookup with the full parameter index.
   IntegerArray fullErrorParams;
   IntegerArray shortErrorParams;
   
   const bool useChains = fit->chainManager()->isSynchedWithFitParams();
   
   for (size_t i=0; i<nSpecifiedPar; ++i)
   {
      int parNum = savErrorParams[i];
      string modErrString = savModName.length() ? 
           savModName + string(":") : string("");
      Parameter* testPar = 
           XSContainer::models->lookupParameter(parNum,savModName);
      if (testPar)
      {
         ModParam* modPar = dynamic_cast<ModParam*>(testPar);
         if (modPar)
         {
            std::map<int,ModParam*>::const_iterator vp =
                 fit->variableParameters().begin();
            std::map<int,ModParam*>::const_iterator vpEnd =
                 fit->variableParameters().end(); 
            while (vp != vpEnd && modPar != vp->second)
            {
               ++vp;
            }
            if (vp == vpEnd)
            {
               tcout << "*** Parameter " << modErrString << parNum
                  << " is not a variable model parameter and no confidence range will be calculated." 
                  << std::endl;
               modPar->lastErrorStatus(FitErrorCalc::errCodeToString(FitErrorCalc::FROZEPAR));
            }
            else
            {
               fullErrorParams.push_back(vp->first);
               shortErrorParams.push_back(parNum);
            }           
         }
         else
         {
            // Not a ModParam - Switch or Scale
            tcout << "*** Parameter " << modErrString << parNum
               << " is not a variable model parameter and no confidence range will be calculated."
               << std::endl;
         }
      }
      else
      {
         tcout << "***Error:  Unable to locate parameter " << modErrString
              << parNum << "\nIf parameter belongs to a named model,"
              << " make sure model name is specified." << std::endl;
         return globalData->autoSave(TCL_ERROR);
      }
   }
   
   // Now checking for fit validity much later.  This is done after
   // pars are verified so we can easily set their status strings.
   if (!useChains && !fit->isStillValid())
   {
      tcout << "A valid fit is first required in order to run error command."
           <<std::endl;
      FitErrorCalc::setAllErrParamStrings(fit, fullErrorParams, 
                FitErrorCalc::GENPROB);
      return globalData->autoSave(TCL_ERROR);           
   }


   SIGINT_Handler intHandler;
   EventHandler* oldHandler = 0;

   std::ios_base::fmtflags save (tcout.flags());
   try 
   {
      if (useChains)
      {
         oldHandler = SignalHandler::instance()->registerHandler(
                         SIGINT, &intHandler);
         tcout <<"Errors calculated from chains\n"     
           << " Parameter   Confidence Range (" << fit->deltaStat() << ")"
           << std::endl;
         for (size_t i=0; i<shortErrorParams.size(); ++i)
         {
            // This can throw
            std::pair<Real,Real> range = 
                fit->chainManager()->getParErrorRange(fit->deltaStat(),
                        shortErrorParams[i], savModName);
            // This has already been verified in loop above:
            ModParam* modPar = 
               fit->variableParameters().find(fullErrorParams[i])->second;
            modPar->setValue(range.first,'m');
            modPar->setValue(range.second,'p');
            modPar->lastErrorStatus(FitErrorCalc::errCodeToString(FitErrorCalc::OK));
            reportParameter(modPar, fit->isStillValid());
         }
         SignalHandler::instance()->registerHandler(SIGINT,oldHandler);
      }
      else
      {
         fit->reinitialize(true);

         // check fit is not overdetermined (more variable
         // parameters than data bins.
         int degrees(0);
         size_t bins(0);

         StatMethod* stat = fit->statMethod();


         stat->degreesOfFreedom(fit,degrees,bins);

         if ( degrees <= 0 ) 
         {
            FitErrorCalc::setAllErrParamStrings(fit, fullErrorParams,
                      FitErrorCalc::GENPROB);
            throw Fit::FitOverDetermined();       
         }
         Real redChi (stat->statistic()/degrees);
         if ( stat->name() == "chi" && redChi > fit->chiMax())
         {
            tcout << "Cannot do error calc: Reduced Chi^2 (= " << redChi << ") > maximum (" 
                  << fit->chiMax() << ")"<< std::endl; 
            FitErrorCalc::setAllErrParamStrings(fit, fullErrorParams,
                      FitErrorCalc::TOOLARGE);
            return globalData->autoSave(TCL_ERROR);
         }
         fit->errorCalc(true);

         tcout << " Parameter   Confidence Range (" << fit->deltaStat() << ")\n";
         tcout <<std::flush;

         oldHandler = SignalHandler::instance()->registerHandler(
                         SIGINT, &intHandler);
         fit->getErrors(fullErrorParams, savModName);
         SignalHandler::instance()->registerHandler(SIGINT,oldHandler);
      }
   }
   catch (YellowAlert&)
   {
      if (useChains)
      {
         // If something threw by this point using chain errors,
         // can assume the problem exists for all parameters.
         FitErrorCalc::setAllErrParamStrings(fit, fullErrorParams, 
                FitErrorCalc::GENPROB);
      }
      if (oldHandler)
      {
         SignalHandler::instance()->registerHandler(SIGINT,oldHandler);
      }
      tcout.flags(save); 
      fit->errorCalc(false);
      return globalData->autoSave(TCL_ERROR);
   }
   tcout << std::flush;
   tcout.flags(save);        
   fit->errorCalc(false);
   return globalData->autoSave(TCL_OK);       
}

namespace
{
   void reportParameter(const ModParam* par, bool isFitValid)
   {
      // This is very similar to FitErrorCalc::reportParameter,
      // which is not accessible when getting errors from chains.
      using std::setw;
      
      std::streamsize savePrec(tcout.precision());
      std::ios_base::fmtflags save(tcout.flags());
      tcout.precision(6);
      tcout << setw(6) << par->index() << setw(13) << par->value('m')
            << setw(13) << par->value('p');
      if (isFitValid)
      {      
         tcout << "    (" << par->value('m') - par->value('a') << ","
	    << par->value('p') - par->value('a') << ")";
      }
      tcout << std::endl;
      tcout.flags(save);
      tcout.precision(savePrec);
   }
}