{-# LANGUAGE OverloadedStrings #-}
import Control.Applicative
import Control.Monad
import Numeric.SpecFunctions
import Numeric.MathFunctions.Constants
import CPython.Sugar
import CPython.MPMath
import qualified CPython as Py

import HEP.ROOT.Plot


----------------------------------------------------------------


viewBetaDelta = runPy $ do
  addToPythonPath "."
  m  <- loadMPMath
  mpmSetDps m 100
  xs <- forM pqBeta $ \(p,q) -> do x <- fromMPNum =<< mpmLog m =<< mpmBeta m (MPDouble p) (MPDouble q)
                                   return (p,q, relErr x (logBeta p q))
  draws $ do
    -- let xs = [ (p,q, logBeta p q `relErr` (logGammaL p + logGammaL q - logGammaL (q+p)))
    --          | (p,q) <- pqBeta
    --          ]
    add $ Graph2D xs


pqBeta = [ (p,q)
         | p <- logRange 50 0.3 0.6
         , q <- logRange 50 5 6
         ]
  where




viewIBeta x = runPy $ do
  addToPythonPath "."
  m <- loadMPMath
  mpmSetDps m 30
  --
  let n  = 40
  let pq =  (,)
        <$> logRange n 100 1000
        <*> logRange n 100 1000
  --
  xs <- forM pq $ \(p,q) -> do
          i <- fromMPNum =<< mpmIncompleteBeta m (MPDouble p) (MPDouble q) (MPDouble x)
          return (p,q, incompleteBeta p q x `relErr` i)
  --
  draws $ do
    add $ Graph2D xs


go = runPy $ do
  addToPythonPath "."
  m <- loadMPMath
  mpmSetDps m 16
  --
  print =<< fromMPNum =<< mpmIncompleteBeta m (MPDouble 10) (MPDouble 10) (MPDouble 0.4)
  print $ incompleteBeta 10 10 0.4




viewLancrox = runPy $ do
  addToPythonPath "."
  m <- loadMPMath
  mpmSetDps m 50
  --
  let xs = logRange 10000 (1e-8) (1e-1)
  pl <- forM xs $ \x -> do y0 <- fromMPNum =<< mpmLog m =<< mpmGamma m (MPDouble x)
                           return (x, y0)
  draws $ do
    add $ Graph $ [ (x, abs $ y `relErr` logGammaL x) | (x,y) <- pl ]
    set $ lineColor RED
    --
    add $ Graph $ [ (x, abs $ y `relErr` logGamma x) | (x,y) <- pl ]
    set $ lineColor BLUE
    --
    set $ xaxis $ logScale ON
    -- set $ yaxis $ logScale ON
    --
    add $ HLine m_epsilon
    add $ HLine $ negate m_epsilon


----------------------------------------------------------------

relErr :: Double -> Double -> Double
relErr 0 0 = 0
relErr x y = (x - y) / max (abs x) (abs y)



logRange :: Int -> Double -> Double -> [Double]
logRange n a b
  = [ a * r^i | i <- [0 .. n] ]
  where
    r = (b / a) ** (1 / fromIntegral n)