#include <LinearInterp.h>
#include <XSUtil/Error/Error.h>

namespace Numerics {

    namespace Rebin {


    void initializeBins(const RealArray& interpolant, const RealArray& inputArray,
             const Real FUZZY, size_t& inputBin, size_t& outputBin,
             IntegerArray& startBin, IntegerArray& endBin, 
             RealArray& startWeight, RealArray& endWeight)
    {
       // The assumption for original starting values is that 
       // inputBin and outputBin are the lowest values which
       // satisfy:
       //    interpolant[outputBin]*(1.0+FUZZY) >= inputArray[inputBin]
       // Therefore if there's any overlap between the two arrays,
       // at least one of inputBin,outputBin must be 0.
       size_t Nout (interpolant.size() - 1);
       size_t Nin  (inputArray.size() - 1);
       startWeight.resize(Nout,0);
       endWeight.resize(Nout,0);
       startBin.clear();
       endBin.clear();
       startBin.resize(Nout,0);
       endBin.resize(Nout,0);

       while ( outputBin < Nout )
       {
           if ( interpolant[outputBin+1] <= inputArray[inputBin+1] )
           {
              endBin[outputBin] = 0;
              startBin[outputBin] = inputBin;
              startWeight[outputBin] 
                   = (interpolant[outputBin+1]- interpolant[outputBin])
                     / ( inputArray[inputBin+1] - inputArray[inputBin]);
              if (interpolant[outputBin+1]*(1.0+FUZZY) >= inputArray[inputBin+1]) 
              {
                 ++inputBin;
              }
              ++outputBin;
           }
           else
           {
              startBin[outputBin] = inputBin;
              startWeight[outputBin] 
                   = ( inputArray[inputBin+1] - interpolant[outputBin] )
                     / (inputArray[inputBin+1] - inputArray[inputBin]);
              while ( inputBin < Nin 
                      && interpolant[outputBin+1] > inputArray[inputBin+1] ) 
              {
                 ++inputBin;
              }
              if ( inputBin >= Nin )
              {
                 endBin[outputBin] = Nin - 1;
                 endBin[outputBin] != startBin[outputBin] 
                                ?  endWeight[outputBin] = 1. 
                                : endBin[outputBin] = 0;
                 ++outputBin;       
              }
              else
              {
                 endBin[outputBin] = inputBin;
                 endWeight[outputBin] 
                     = (interpolant[outputBin+1] - inputArray[inputBin])
                        / ( inputArray[inputBin+1] - inputArray[inputBin]);
                 if (interpolant[outputBin+1]*(1.0+FUZZY) >= inputArray[inputBin+1]) 
                 {
                    ++inputBin;
                 }
                 ++outputBin;
              }
           }

           if ( inputBin >= Nin ) outputBin = Nout + 1;
        }
    }

        void rebin(const RealArray& inputArray, const IntegerArray& startBin, 
              const IntegerArray& endBin, const RealArray& startWeight, const RealArray& endWeight,
              RealArray& outputArray)
        {
              size_t Nout(outputArray.size());
              for (size_t j = 0; j < Nout; ++j)
              {
                    if ( startBin[j] >= 0 ) 
                    {
                            outputArray[j] = startWeight[j]*inputArray[startBin[j]];
                            if ( endBin[j] > startBin[j] )
                            {
                                    for ( int k = startBin[j] + 1; k < endBin[j]; ++k)
                                    {
                                            outputArray[j] += inputArray[k];
                                    } 
                                    outputArray[j] += endWeight[j]*inputArray[endBin[j]];
                            }      
                    }
                    else outputArray[j] = 0;          
              }      

        }

   void interpolate(const RealArray& inputArray,const IntegerArray& startBin, 
         const IntegerArray& endBin, const RealArray& startWeight, const RealArray& endWeight,
         RealArray& outputArray, bool exponential)
   {
      size_t Nout (outputArray.size());
      for ( size_t j = 0; j < Nout; ++j)
      {
         Real denominator (0.);
         outputArray[j] = startWeight[j]*inputArray[startBin[j]];
         denominator += startWeight[j];
         if ( endBin[j] > startBin[j] )
         {
            for ( int k = startBin[j] + 1; k < endBin[j]; ++k)
            {
               outputArray[j] += inputArray[k];
               denominator += 1;
            }
            outputArray[j] 
                    += endWeight[j]*inputArray[endBin[j]];
            denominator += endWeight[j];
         }
         if ( denominator > 0 )
         {
            outputArray[j] /= denominator;
         }
         else
         {
            outputArray[j] = ( exponential ? 0 :  1.) ;
         }
      }
      if ( exponential )
      {
         for (size_t j = 0; j < Nout; ++j) 
         { 
            outputArray[j] > 100. 
              ?  outputArray[j] = 0.
              : outputArray[j] = exp(-outputArray[j]);
         }           
      }

   }

        bool findFirstBins(const RealArray& inputArray, const RealArray& outputArray,
                    const Real FUZZY, size_t& inputStart, size_t& outputStart)
        {
           // Find first bin in outputArray which which is contained
           // in range of inputArray, and the bin inputArray to
           // which it corresponds.  To conform with inibin.f function,
           // an output bin only partially contained on LOW end of 
           // inputArray is NOT accepted.  It IS accepted if it is
           // partially overlapping the HIGH end.  If no output bin
           // is found, function returns false.
           bool isOverlap = false;
           inputStart = outputStart = 0;
           if (inputArray.size() && outputArray.size())
           {
              size_t nInBins = inputArray.size() - 1;
              size_t nOutBins = outputArray.size() - 1;
              while (outputStart < nOutBins && outputArray[outputStart]*(1.+FUZZY)
                        < inputArray[0])
              {
                 ++outputStart;
              }
              if (outputStart < nOutBins)
              {
                 while (inputStart < nInBins && inputArray[inputStart+1] 
                        < outputArray[outputStart]*(1.+FUZZY))
                 {
                    ++inputStart;
                 }
                 if (inputStart < nInBins)
                 {
                    isOverlap = true;
                 }  
              }
           }
           return isOverlap;
        }

       void gainRebin(const RealArray& inputArray, const IntegerArray &startBin,
               const IntegerArray& endBin, const RealArray& startWeight,
               const RealArray& endWeight, RealArray& outputArray)
       {
          // ASSUMPTION: All arrays passed here are the SAME SIZE
          size_t sz = outputArray.size();
          // Partial check.
          if (sz != inputArray.size() || sz != startBin.size())
          {
             throw RedAlert("Array size mismatch occured in gainRebin function.");
          }
          for (size_t i=0; i<sz; ++i)
          {
             Real total = .0;
             Real& output_i = outputArray[i];
             const int startBin_i = startBin[i];
             const int endBin_i = endBin[i];
             output_i = startWeight[i]*inputArray[startBin_i];
             total += startWeight[i];
             if (endBin_i > startBin_i)
             {
                for (int j=startBin_i+1; j<endBin_i; ++j)
                {
                   // If we're here, output bin spreads over more than
                   // just 2 input bins.
                   output_i += inputArray[j];
                   total += 1.0;
                }
                output_i += endWeight[i] *inputArray[endBin_i];
                total += endWeight[i];
             }
             if (total > .0)
             {
                output_i /= total;
             }
             else
             {
                output_i = 1.0;
             }
          }
       }


    }  // namespace Rebin

}  // namespace Numerics