module Main where

import qualified Numeric.Interpolation.NodeList as Nodes
import qualified Numeric.Interpolation.Piecewise as Piecewise
import qualified Numeric.Interpolation.Type as Type

import qualified Data.Packed.Matrix as Matrix
import qualified Data.Packed.Vector as Vector
import qualified Numeric.Container as Container
import Numeric.Container ((<\>))

import qualified Graphics.Gnuplot.Advanced as GP
import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D

import System.Random (randomRs, mkStdGen)
import Control.Monad.HT (void)
import Data.Monoid ((<>))


noisy :: [(Double, Double)]
noisy =
   take 100 $
   zipWith
      (\x d -> (x, sin x + d))
      (randomRs (0,7) (mkStdGen 23))
      (randomRs (-0.2,0.2) (mkStdGen 42))

fit ::
   Type.T Double Double ny ->
   [Double] -> [(Double, Double)] -> Nodes.T Double ny
fit typ xs target =
   let txs = Vector.fromList $ map fst target
       tys = Vector.fromList $ map snd target
       matrix =
          Matrix.fromColumns $
          map (flip Container.cmap txs . Piecewise.interpolateConstantExt typ) $
          Type.basisFunctions typ xs
   in  Type.coefficientsToInterpolator typ xs $
       Vector.toList $ matrix <\> tys

plotBasisFunctions ::
   Type.T Double Double ny -> [Double] -> Plot2D.T Double Double
plotBasisFunctions nodeType xs =
   let abscissa = Plot2D.linearScale 1000 (minimum xs, maximum xs)
   in  Plot2D.functions Graph2D.lines abscissa $
       map (Piecewise.interpolateConstantExt nodeType) $
       Type.basisFunctions nodeType xs


main :: IO ()
main = do
   let xs = [0, 1, 3, 4, 6, 7]
       exs = (-1) : xs ++ [8]
   void $ GP.plotDefault $ plotBasisFunctions Type.linear xs
   void $ GP.plotDefault $ plotBasisFunctions Type.cubic xs
   void $ GP.plotDefault $ plotBasisFunctions Type.cubicLinear exs
   void $ GP.plotDefault $ plotBasisFunctions Type.cubicParabola exs
   let linearNodes = fit Type.linear xs noisy
       hermite1Nodes = fit Type.cubic xs noisy
       cubicLinearNodes = fit Type.cubicLinear exs noisy
       cubicParabolaNodes = fit Type.cubicParabola exs noisy
   void $ GP.plotDefault $
      Plot2D.list Graph2D.points noisy
      <>
      (Plot2D.functions Graph2D.lines (Plot2D.linearScale 1000 (-2,10)) $
       Piecewise.interpolateConstantExt Type.linear linearNodes :
       Piecewise.interpolateConstantExt Type.cubic hermite1Nodes :
       Piecewise.interpolateConstantExt Type.cubicLinear cubicLinearNodes :
       Piecewise.interpolateConstantExt Type.cubicParabola cubicParabolaNodes :
       [])