//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// G4ChordFinderDelegate inline methods implementation
//
// Author: Dmitry Sorokin (CERN, Google Summer of Code 2017), 12.09.2018
// --------------------------------------------------------------------

template <class Driver>
G4ChordFinderDelegate<Driver>::~G4ChordFinderDelegate()
{
#ifdef G4VERBOSE
    if (GetDriver().GetVerboseLevel() > 0)
    {
        PrintStatistics();
    }
#endif
}

template <class Driver>
void G4ChordFinderDelegate<Driver>::ResetStepEstimate()
{
    fLastStepEstimate_Unconstrained = DBL_MAX;  
}

template <class Driver>
Driver& G4ChordFinderDelegate<Driver>::GetDriver()
{
    return static_cast<Driver&>(*this);
}

template <class Driver>
G4double G4ChordFinderDelegate<Driver>::
AdvanceChordLimitedImpl(G4FieldTrack& yCurrent, G4double stepMax, 
                        G4double epsStep, G4double chordDistance)
{
    G4double dyErr;
    G4FieldTrack yEnd = yCurrent;
    G4double nextStep;

    const G4double stepPossible = FindNextChord(yCurrent, stepMax,
                                                epsStep, chordDistance,
                                                yEnd, dyErr, nextStep);
    if (dyErr < epsStep * stepPossible)
    {
         // Accept this accuracy.
         //
         yCurrent = yEnd;
         return stepPossible;
    }

    // Advance more accurately to "end of chord"
    //
    const G4double startCurveLen = yCurrent.GetCurveLength();
    const G4bool goodAdvance =
          GetDriver().AccurateAdvance(yCurrent,stepPossible,epsStep,nextStep);

    return goodAdvance ? stepPossible
                       : yCurrent.GetCurveLength() - startCurveLen;
}

// Returns Length of Step taken
//
template <class T>
G4double G4ChordFinderDelegate<T>::
FindNextChord(const G4FieldTrack& yStart,
              G4double stepMax,
              G4double epsStep,
              G4double chordDistance,
              G4FieldTrack& yEnd, // Endpoint
              G4double& dyErrPos, // Error of endpoint
              G4double& stepForAccuracy)
{
    //  1.)  Try to "leap" to end of interval
    //  2.)  Evaluate if resulting chord gives d_chord that is good enough.
    // 2a.)  If d_chord is not good enough, find one that is.

    G4double dydx[G4FieldTrack::ncompSVEC];
  
    G4bool validEndPoint = false;
    G4double dChordStep, lastStepLength;

    GetDriver().GetDerivatives(yStart, dydx);

    const G4double safetyFactor = fFirstFraction; //  0.975 or 0.99 ? was 0.999

    G4double stepTrial = std::min(stepMax,
                                  safetyFactor*fLastStepEstimate_Unconstrained);

    G4double newStepEst_Uncons = 0.0; 
    G4double stepForChord;

    G4int noTrials = 1;
    constexpr G4int maxTrials = 75; // Avoid endless loop for bad convergence 
    for (; noTrials < maxTrials; ++noTrials)
    {
        yEnd = yStart; // Always start from initial point  
        GetDriver().QuickAdvance(yEnd, dydx, stepTrial, dChordStep, dyErrPos);
        lastStepLength = stepTrial; 

        validEndPoint = dChordStep < chordDistance;
        stepForChord = NewStep(stepTrial, dChordStep,
                               chordDistance, newStepEst_Uncons);
        if (validEndPoint)
        {
            break;
        }

        if (stepTrial <= 0.0)
        {
            stepTrial = stepForChord;
        }
        else if (stepForChord <= stepTrial)
        {
            // Reduce by a fraction, possibly up to 20% 
            stepTrial = std::min( stepForChord, fFractionLast * stepTrial);
        }
        else
        {
            stepTrial *= 0.1;
        }
    }

    if (noTrials >= maxTrials)
    {
        std::ostringstream message;
        message << "Exceeded maximum number of trials= " << maxTrials << G4endl
                << "Current sagita dist= " << dChordStep << G4endl
                << "Max sagita dist= " << chordDistance << G4endl
                << "Step sizes (actual and proposed): " << G4endl
                << "Last trial =         " << lastStepLength  << G4endl
                << "Next trial =         " << stepTrial  << G4endl
                << "Proposed for chord = " << stepForChord  << G4endl;
        G4Exception("G4ChordFinder::FindNextChord()", "GeomField0003",
                    JustWarning, message);
    }

    if (newStepEst_Uncons > 0.0)
    {
        fLastStepEstimate_Unconstrained = newStepEst_Uncons;
    }

    AccumulateStatistics(noTrials);


    // Calculate the step size required for accuracy, if it is needed
    G4double dyErr_relative = dyErrPos / (epsStep * lastStepLength);
    stepForAccuracy = dyErr_relative > 1 ? 
        GetDriver().ComputeNewStepSize(dyErr_relative, lastStepLength) : 0;

  return stepTrial; 
}

// Is called to estimate the next step size, even for successful steps,
// in order to predict an accurate 'chord-sensitive' first step
// which is likely to assist in more performant 'stepping'.
//
template <class T>
G4double G4ChordFinderDelegate<T>::
NewStep(G4double stepTrialOld, 
        G4double dChordStep, // Curr. dchord achieved
        G4double fDeltaChord,
        G4double& stepEstimate_Unconstrained)  
{
    G4double stepTrial;

    if (dChordStep > 0.0)
    {
        stepEstimate_Unconstrained =
            stepTrialOld * std::sqrt(fDeltaChord / dChordStep);
        stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
    }
    else
    {
        // Should not update the Unconstrained Step estimate: incorrect!
        stepTrial =  stepTrialOld * 2.; 
    }

    if (stepTrial <= 0.001 * stepTrialOld)
    {
        if (dChordStep > 1000.0 * fDeltaChord)
        {
            stepTrial = stepTrialOld * 0.03;   
        }
        else
        {
            if (dChordStep > 100. * fDeltaChord)
            {
                stepTrial = stepTrialOld * 0.1;   
            }
            else   // Try halving the length until dChordStep OK
            {
                stepTrial = stepTrialOld * 0.5;   
            }
        }
    }
    else if (stepTrial > 1000.0 * stepTrialOld)
    {
        stepTrial = 1000.0 * stepTrialOld;
    }

    if (stepTrial == 0.0)
    {
        stepTrial= 0.000001;
    }

    // A more sophisticated chord-finder could figure out a better
    // stepTrial, from dChordStep and the required d_geometry
    //   e.g.
    //      Calculate R, r_helix (eg at orig point)
    //      if( stepTrial < 2 pi  R )
    //          stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
    //      else    
    //          ??

    return stepTrial;
}

template <class T>
void G4ChordFinderDelegate<T>::AccumulateStatistics(G4int noTrials) 
{
    fTotalNoTrials += noTrials; 
    ++fNoCalls; 
      
    if (noTrials > fmaxTrials) 
    { 
        fmaxTrials = noTrials; 
    }
}

template <class T>
void G4ChordFinderDelegate<T>::PrintStatistics()
{
    // Print Statistics
    G4cout << "G4ChordFinder statistics report: \n"
           << "  No trials: " << fTotalNoTrials
           << "  No Calls: "  << fNoCalls
           << "  Max-trial: " <<  fmaxTrials << "\n"
           << "  Parameters: " 
           << "  fFirstFraction "  << fFirstFraction
           << "  fFractionLast "   << fFractionLast
           << "  fFractionNextEstimate " << fFractionNextEstimate
           << G4endl; 
}

template <class T>
G4int G4ChordFinderDelegate<T>::GetNoCalls()
{
    return fNoCalls;
}

template <class T>
G4int G4ChordFinderDelegate<T>::GetNoTrials()
{
    return fTotalNoTrials;
}

template <class T>
G4int G4ChordFinderDelegate<T>::GetNoMaxTrials()
{
    return fmaxTrials;
}

template <class T>
void G4ChordFinderDelegate<T>::SetFractions_Last_Next(G4double fractLast, 
                                                      G4double fractNext)
{
    // Use -1.0 as request for Default.
    if (fractLast == -1.0) { fractLast = 1.0; }  // 0.9;
    if (fractNext == -1.0) { fractNext = 0.98; } // 0.9; 

    // fFirstFraction  = 0.999; // Safe value, range: ~ 0.95 - 0.999
    if (GetDriver().GetVerboseLevel() > 0)
    { 
        G4cout << " ChordFnd> Trying to set fractions: "
               << " first " << fFirstFraction
               << " last " <<  fractLast
               << " next " <<  fractNext
               << G4endl;
    } 

    if (fractLast > 0 && fractLast <= 1) 
    {
        fFractionLast = fractLast;
    } else
    {
        std::ostringstream message;
        message << "Invalid fraction Last = " << fractLast
                << "; must be  0 <  fractionLast <= 1 ";
        G4Exception("G4ChordFinderDelegate::SetFractions_Last_Next()",
                    "GeomField1001", JustWarning, message);
    }
    if (fractNext > 0. && fractNext < 1)
    {
        fFractionNextEstimate = fractNext;
    } else
    {
        std::ostringstream message;
        message << "Invalid fraction Next = " << fractNext
                << "; must be  0 <  fractionNext < 1 ";
        G4Exception("G4ChordFinderDelegate::SetFractions_Last_Next()",
                    "GeomField1001", JustWarning, message);
    }
}

template <class T>
void G4ChordFinderDelegate<T>::SetFirstFraction(G4double fractFirst)
{
    fFirstFraction = fractFirst;
}

template <class T>
G4double G4ChordFinderDelegate<T>::GetFirstFraction()
{
    return fFirstFraction;
}

template <class T>
G4double G4ChordFinderDelegate<T>::GetFractionLast()
{
    return fFractionLast;
}

template <class T>
G4double G4ChordFinderDelegate<T>::GetFractionNextEstimate()
{
    return fFractionNextEstimate;
}

template <class T>
G4double G4ChordFinderDelegate<T>::GetLastStepEstimateUnc()
{
    return fLastStepEstimate_Unconstrained;   
}

template <class T>
void G4ChordFinderDelegate<T>::SetLastStepEstimateUnc(G4double stepEst)
{
    fLastStepEstimate_Unconstrained = stepEst;
}

template <class T>
void G4ChordFinderDelegate<T>::TestChordPrint(G4int noTrials, 
                                              G4int lastStepTrial, 
                                              G4double dChordStep,
                                              G4double fDeltaChord,
                                              G4double nextStepTrial)
{
     G4int oldprec = G4cout.precision(5);
     G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials 
            << " this_step= "       << std::setw(10) << lastStepTrial;
     if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
     {
       G4cout.precision(8);
     }
     else
     {
       G4cout.precision(6);
     }
     G4cout << " dChordStep=  " << std::setw(12) << dChordStep;
     if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
     else                           { G4cout << " d-"; }
     G4cout.precision(5);
     G4cout <<  " new_step= "       << std::setw(10)
            << fLastStepEstimate_Unconstrained
            << " new_step_constr= " << std::setw(10)
            << lastStepTrial << G4endl;
     G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
     G4cout.precision(oldprec);
}

template <class T>
void G4ChordFinderDelegate<T>::StreamDelegateInfo( std::ostream& os ) const
{
// Write out the parameters / state of the driver
  os << "State of G4ChordFinderDelegate: " << std::endl;
  os << "--Parameters: " << std::endl;
  os << "    First Fraction = " << fFirstFraction << std::endl;
  os << "    Last  Fraction = " << fFractionLast << std::endl;
  os << "    Fract Next est = " << fFractionNextEstimate << std::endl;

  os << "--State (fungible): " << std::endl;
  os << "    Maximum No Trials (seen)         = " << fmaxTrials  << std::endl;  
  os << "    LastStepEstimate (Unconstrained) = " << fLastStepEstimate_Unconstrained
     << std::endl;
  // os << "    Statistics NOT printed. " << std::endl;
  os << "--Statistics: trials= " << fTotalNoTrials
     << "  calls= " << fNoCalls << std::endl;
}