include/Geant4/G4DNARevertProbability.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /*
0027  * G4DNARevertFunction.hh
0028  *
0029  *  Created on: Aug 25, 2015
0030  *      Author: mkaramit
0031  */
0032
0033 #ifndef SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_UTILS_INCLUDE_G4DNAREVERTPROBABILITY_HH_
0034 #define SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_UTILS_INCLUDE_G4DNAREVERTPROBABILITY_HH_
0035
0036 #include <functional>
0037
0038 //------------------------------------------------------------------------------
0039 // *****************************************************
0040 // Revert a probability by binary search
0041 // *****************************************************
0042
0043 class G4DNARevertProbability
0044 {
0045 public:
0046   double fXmax;
0047   double fXmin;
0048   double fXmaxDef;
0049   double fXminDef;
0050   double fToleranceProba;
0051   double fIncreasingCumulativeFunction;
0052
0053   G4DNARevertProbability(double toleranceY) :
0054       fToleranceProba(toleranceY)
0055   {
0056     fXmax = 0.;
0057     fXmin = 0.;
0058     fXminDef = 0.;
0059     fXmaxDef = 0.;
0060     fIncreasingCumulativeFunction = 0.;
0061   }
0062
0063   G4DNARevertProbability(double xmin, double xmax, double toleranceY) :
0064       fXmax(xmax), fXmin(xmin), fToleranceProba(toleranceY)
0065   {
0066     if(fXmax < fXmin)
0067     {
0068       double tmp = fXmin;
0069       fXmin = fXmax;
0070       fXmax = tmp;
0071     }
0072
0073     fXminDef = fXmin;
0074     fXmaxDef = fXmax;
0075     fIncreasingCumulativeFunction = 0;
0076   }
0077
0078   void Reset()
0079   {
0080     fXmin = fXminDef;
0081     fXmax = fXmaxDef;
0082   }
0083
0084   void SetBoundaries(double xmin, double xmax)
0085   {
0086     fXmax = xmax;
0087     fXmin = xmin;
0088     if(fXmax < fXmin)
0089     {
0090       double tmp = fXmin;
0091       fXmin = fXmax;
0092       fXmax = tmp;
0093     }
0094
0095     fXminDef = fXmin;
0096     fXmaxDef = fXmax;
0097   }
0098
0099   void SetTolerance(double toleranceY)
0100   {
0101     fToleranceProba = toleranceY;
0102   }
0103
0104   bool Propose(double proposedXValue,
0105                double proposedProba,
0106                double nextProba,
0107                double& returnedValue)
0108   {
0109     bool returnFlag = false;
0110
0111     if(proposedProba < nextProba - fToleranceProba) // proba trop petite ==> augmente
0112     {
0113       if(fIncreasingCumulativeFunction > 0) // croissant
0114       {
0115         if(proposedXValue > fXmin) fXmin = proposedXValue;
0116       }
0117       else if(fIncreasingCumulativeFunction < 0) // decroissant
0118       {
0119         if(proposedXValue < fXmax) fXmax = proposedXValue;
0120       }
0121
0122       returnedValue = (fXmax + fXmin) / 2.;
0123       returnFlag = false;
0124     }
0125     else if(proposedProba > nextProba + fToleranceProba) // proba trop grande
0126     {
0127       if(fIncreasingCumulativeFunction > 0)
0128       {
0129         if(proposedXValue < fXmax)
0130         {
0131           fXmax = proposedXValue;
0132         }
0133       }
0134       else if(fIncreasingCumulativeFunction < 0)
0135       {
0136         if(proposedXValue > fXmin)
0137         {
0138           fXmin = proposedXValue;
0139         }
0140       }
0141
0142       returnedValue = (fXmax + fXmin) / 2.;
0143       returnFlag = false;
0144     }
0145     else
0146     {
0147       // Assuming search for next proba is increasing
0148       if(fIncreasingCumulativeFunction < 0)
0149       {
0150         fXmin = fXminDef;
0151         fXmax = proposedXValue;
0152       }
0153       else if(fIncreasingCumulativeFunction > 0)
0154       {
0155         fXmin = proposedXValue;
0156         fXmax = fXmaxDef;
0157       }
0158       returnFlag = true;
0159     }
0160
0161     return returnFlag;
0162   }
0163
0164   //**********************************
0165   // *** REVERT THE FUNCTION "FUNCT"
0166   // *********************************
0167
0168   double Revert(double probaForSearchedTime,
0169                 std::function<double(double)>& funct)
0170   {
0171     Reset();
0172     bool notFound = true;
0173     double proposedX;
0174     double x = fXmax;
0175
0176     fIncreasingCumulativeFunction = (funct(fXmax) - funct(fXmin))
0177         / (fXmax - fXmin);
0178
0179     while(notFound)
0180     {
0181       double newProba = funct(x);
0182
0183       if(Propose(x, newProba, probaForSearchedTime, proposedX))
0184       {
0185         notFound = false;
0186       }
0187       else
0188       {
0189         if(x == proposedX)
0190         {
0191           return x;
0192 //          G4cout << "BREAK : x= " << x << G4endl;
0193 //          G4cout << "nextProba= " << probaForSearchedTime << G4endl;
0194 //          G4cout << "newProba= " << newProba << G4endl;
0195 //          abort();
0196         }
0197         x = proposedX;
0198       }
0199     }
0200     return x;
0201   }
0202 };
0203 #endif /* SOURCE_PROCESSES_ELECTROMAGNETIC_DNA_UTILS_INCLUDE_G4DNAREVERTPROBABILITY_HH_ */