// Created on: 2005-12-08
// Created by: Sergey KHROMOV
// Copyright (c) 2005-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <math.hxx>
#include <math_Function.hxx>
#include <math_KronrodSingleIntegration.hxx>
#include <math_Vector.hxx>
#include <TColStd_SequenceOfReal.hxx>

//==========================================================================
// function : An empty constructor.
//==========================================================================
math_KronrodSingleIntegration::math_KronrodSingleIntegration()
    : myIsDone(Standard_False),
      myValue(0.),
      myErrorReached(0.),
      myAbsolutError(0.),
      myNbPntsReached(0),
      myNbIterReached(0)
{
}

//==========================================================================
// function : Constructor
//
//==========================================================================

math_KronrodSingleIntegration::math_KronrodSingleIntegration(math_Function&         theFunction,
                                                             const Standard_Real    theLower,
                                                             const Standard_Real    theUpper,
                                                             const Standard_Integer theNbPnts)
    : myIsDone(Standard_False),
      myValue(0.),
      myErrorReached(0.),
      myAbsolutError(0.),
      myNbPntsReached(0),
      myNbIterReached(0)
{
  Perform(theFunction, theLower, theUpper, theNbPnts);
}

//==========================================================================
// function : Constructor
//
//==========================================================================

math_KronrodSingleIntegration::math_KronrodSingleIntegration(math_Function&         theFunction,
                                                             const Standard_Real    theLower,
                                                             const Standard_Real    theUpper,
                                                             const Standard_Integer theNbPnts,
                                                             const Standard_Real    theTolerance,
                                                             const Standard_Integer theMaxNbIter)
    : myIsDone(Standard_False),
      myValue(0.),
      myErrorReached(0.),
      myAbsolutError(0.),
      myNbPntsReached(0),
      myNbIterReached(0)
{
  Perform(theFunction, theLower, theUpper, theNbPnts, theTolerance, theMaxNbIter);
}

//==========================================================================
// function : Perform
//           Computation of the integral.
//==========================================================================

void math_KronrodSingleIntegration::Perform(math_Function&         theFunction,
                                            const Standard_Real    theLower,
                                            const Standard_Real    theUpper,
                                            const Standard_Integer theNbPnts)
{
  // const Standard_Real aMinVol = Epsilon(1.);
  const Standard_Real aPtol = 1.e-9;
  myNbIterReached           = 0;

  if (theNbPnts < 3)
  {
    myIsDone = Standard_False;
    return;
  }

  if (theUpper - theLower < aPtol)
  {
    myIsDone = Standard_False;
    return;
  }

  // Get an odd value of number of initial points.
  myNbPntsReached = (theNbPnts % 2 == 0) ? theNbPnts + 1 : theNbPnts;
  myErrorReached  = RealLast();

  Standard_Integer aNGauss = myNbPntsReached / 2;
  math_Vector      aKronrodP(1, myNbPntsReached);
  math_Vector      aKronrodW(1, myNbPntsReached);
  math_Vector      aGaussP(1, aNGauss);
  math_Vector      aGaussW(1, aNGauss);

  if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW)
      || !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW))
  {
    myIsDone = Standard_False;
    return;
  }

  myIsDone = GKRule(theFunction,
                    theLower,
                    theUpper,
                    aGaussP,
                    aGaussW,
                    aKronrodP,
                    aKronrodW,
                    myValue,
                    myErrorReached);

  if (!myIsDone)
    return;

  // Standard_Real anAbsVal = Abs(myValue);

  myAbsolutError = myErrorReached;

  // if (anAbsVal > aMinVol)
  // myErrorReached /= anAbsVal;

  myNbIterReached++;
}

//=================================================================================================

void math_KronrodSingleIntegration::Perform(math_Function&         theFunction,
                                            const Standard_Real    theLower,
                                            const Standard_Real    theUpper,
                                            const Standard_Integer theNbPnts,
                                            const Standard_Real    theTolerance,
                                            const Standard_Integer theMaxNbIter)
{
  Standard_Real aMinVol = Epsilon(1.);
  myNbIterReached       = 0;

  // Check prerequisites.
  if (theNbPnts < 3 || theTolerance <= 0.)
  {
    myIsDone = Standard_False;
    return;
  }
  // Get an odd value of number of initial points.
  myNbPntsReached = (theNbPnts % 2 == 0) ? theNbPnts + 1 : theNbPnts;

  Standard_Integer aNGauss = myNbPntsReached / 2;
  math_Vector      aKronrodP(1, myNbPntsReached);
  math_Vector      aKronrodW(1, myNbPntsReached);
  math_Vector      aGaussP(1, aNGauss);
  math_Vector      aGaussW(1, aNGauss);

  if (!math::KronrodPointsAndWeights(myNbPntsReached, aKronrodP, aKronrodW)
      || !math::OrderedGaussPointsAndWeights(aNGauss, aGaussP, aGaussW))
  {
    myIsDone = Standard_False;
    return;
  }

  // First iteration
  myIsDone = GKRule(theFunction,
                    theLower,
                    theUpper,
                    aGaussP,
                    aGaussW,
                    aKronrodP,
                    aKronrodW,
                    myValue,
                    myErrorReached);

  if (!myIsDone)
    return;

  Standard_Real anAbsVal = Abs(myValue);

  myAbsolutError = myErrorReached;
  if (anAbsVal > aMinVol)
    myErrorReached /= anAbsVal;

  myNbIterReached++;

  if (myErrorReached <= theTolerance)
    return;
  if (myNbIterReached >= theMaxNbIter)
    return;

  TColStd_SequenceOfReal anIntervals;
  TColStd_SequenceOfReal anErrors;
  TColStd_SequenceOfReal aValues;

  anIntervals.Append(theLower);
  anIntervals.Append(theUpper);

  anErrors.Append(myAbsolutError);
  aValues.Append(myValue);

  Standard_Integer i, nint, nbints;

  Standard_Real    maxerr;
  Standard_Integer count = 0;

  while (myErrorReached > theTolerance && myNbIterReached < theMaxNbIter)
  {
    // Searching interval with max error
    nbints = anIntervals.Length() - 1;
    nint   = 0;
    maxerr = 0.;
    for (i = 1; i <= nbints; ++i)
    {
      if (anErrors(i) > maxerr)
      {
        maxerr = anErrors(i);
        nint   = i;
      }
    }

    Standard_Real a = anIntervals(nint);
    Standard_Real b = anIntervals(nint + 1);
    Standard_Real c = 0.5 * (a + b);

    Standard_Real v1, v2, e1, e2;

    myIsDone = GKRule(theFunction, a, c, aGaussP, aGaussW, aKronrodP, aKronrodW, v1, e1);

    if (!myIsDone)
      return;

    myIsDone = GKRule(theFunction, c, b, aGaussP, aGaussW, aKronrodP, aKronrodW, v2, e2);

    if (!myIsDone)
      return;

    myNbIterReached++;

    Standard_Real deltav = v1 + v2 - aValues(nint);
    myValue += deltav;
    if (Abs(deltav) <= Epsilon(Abs(myValue)))
      ++count;

    Standard_Real deltae = e1 + e2 - anErrors(nint);
    myAbsolutError += deltae;
    if (myAbsolutError <= Epsilon(Abs(myValue)))
      ++count;

    if (Abs(myValue) > aMinVol)
      myErrorReached = myAbsolutError / Abs(myValue);
    else
      myErrorReached = myAbsolutError;

    if (count > 50)
      return;

    // Inserting new interval

    anIntervals.InsertAfter(nint, c);
    anErrors(nint) = e1;
    anErrors.InsertAfter(nint, e2);
    aValues(nint) = v1;
    aValues.InsertAfter(nint, v2);
  }
}

//=================================================================================================

Standard_Boolean math_KronrodSingleIntegration::GKRule(math_Function&      theFunction,
                                                       const Standard_Real theLower,
                                                       const Standard_Real theUpper,
                                                       const math_Vector& /*theGaussP*/,
                                                       const math_Vector& theGaussW,
                                                       const math_Vector& theKronrodP,
                                                       const math_Vector& theKronrodW,
                                                       Standard_Real&     theValue,
                                                       Standard_Real&     theError)
{

  Standard_Boolean IsDone;

  Standard_Integer aNKronrod = theKronrodP.Length();

  Standard_Real    aGaussVal;
  Standard_Integer aNPnt2 = (aNKronrod + 1) / 2;
  Standard_Integer i;
  Standard_Real    aDx;
  Standard_Real    aVal1;
  Standard_Real    aVal2;

  math_Vector f1(1, aNPnt2 - 1);
  math_Vector f2(1, aNPnt2 - 1);

  Standard_Real aXm = 0.5 * (theUpper + theLower);
  Standard_Real aXr = 0.5 * (theUpper - theLower);

  // Compute Gauss quadrature
  aGaussVal = 0.;
  theValue  = 0.;

  for (i = 2; i < aNPnt2; i += 2)
  {
    aDx = aXr * theKronrodP.Value(i);

    if (!theFunction.Value(aXm + aDx, aVal1) || !theFunction.Value(aXm - aDx, aVal2))
    {
      IsDone = Standard_False;
      return IsDone;
    }

    f1(i) = aVal1;
    f2(i) = aVal2;
    aGaussVal += (aVal1 + aVal2) * theGaussW.Value(i / 2);
    theValue += (aVal1 + aVal2) * theKronrodW.Value(i);
  }

  // Compute value in the middle point.
  if (!theFunction.Value(aXm, aVal1))
  {
    IsDone = Standard_False;
    return IsDone;
  }

  Standard_Real fc = aVal1;
  theValue += aVal1 * theKronrodW.Value(aNPnt2);

  if (i == aNPnt2)
    aGaussVal += aVal1 * theGaussW.Value(aNPnt2 / 2);

  // Compute Kronrod quadrature
  for (i = 1; i < aNPnt2; i += 2)
  {
    aDx = aXr * theKronrodP.Value(i);

    if (!theFunction.Value(aXm + aDx, aVal1) || !theFunction.Value(aXm - aDx, aVal2))
    {
      IsDone = Standard_False;
      return IsDone;
    }

    f1(i) = aVal1;
    f2(i) = aVal2;
    theValue += (aVal1 + aVal2) * theKronrodW.Value(i);
  }

  Standard_Real mean = 0.5 * theValue;

  Standard_Real asc = Abs(fc - mean) * theKronrodW.Value(aNPnt2);
  for (i = 1; i < aNPnt2; ++i)
  {
    asc += theKronrodW.Value(i) * (Abs(f1(i) - mean) + Abs(f2(i) - mean));
  }

  asc *= aXr;
  theValue *= aXr;
  aGaussVal *= aXr;

  // Compute the error and the new number of Kronrod points.

  theError = Abs(theValue - aGaussVal);

  Standard_Real scale = 1.;
  if (asc != 0. && theError != 0.)
    scale = Pow((200. * theError / asc), 1.5);
  if (scale < 1.)
    theError = Min(theError, asc * scale);

  // theFunction.GetStateNumber();

  return Standard_True;
}