Poisson.hs

{-# LANGUAGE UndecidableInstances #-}

{-# OPTIONS_GHC -Wno-type-defaults #-}

module Poisson (main) where

import Prelude hiding (id, (.))
import qualified Prelude
import Data.Kind
import Math.Polynomial hiding (x)
import Math.Polynomial.Interpolation
import Math.Polynomial.Legendre
import qualified Math.Polynomial

default (Int)



-- | A constrained category, so that we can compose functions
class Category k where
  type Ok k :: Type -> Constraint
  id :: Ok k a => k a a
  (.) :: (Ok k a, Ok k b, Ok k c) => k b c -> k a b -> k a c
  eval :: (Ok k a, Ok k b) => k a b -> a -> b

class Unconstrained a
instance Unconstrained a

-- | Hask is a category
instance Category (->) where
  type Ok (->) = Unconstrained
  id = Prelude.id
  (.) = (Prelude..)
  eval = id

data Fun a b where
  FunHask :: (a -> b) -> Fun a b
  FunPoly :: Poly a -> Fun a a
  FunId :: Fun a a
  FunComp :: FunOk b => Fun b c -> Fun a b -> Fun a c

class (Eq a, Num a) => FunOk a
instance (Eq a, Num a) => FunOk a

evalFun :: (FunOk a, FunOk b) => Fun a b -> a -> b
evalFun (FunHask f) = f
evalFun (FunPoly p) = evalPoly p
evalFun FunId = id
evalFun (FunComp g f) = evalFun g . evalFun f

-- | Fun is a category
instance Category Fun where
  type Ok Fun = FunOk
  id = FunId
  (.) = FunComp
  eval = evalFun



-- | Logistic function
-- > 1/2 + 1/2 tanh (x / 2)
logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

-- | Gaussian
gaussian :: Floating a => a -> a
gaussian x = exp (-1/2 * x^2)

-- | Sigmoid function
sigmoid :: Floating a => a -> a
sigmoid = logistic



xpoly :: (Eq a, Num a) => Poly a
xpoly = Math.Polynomial.x

approxPoly :: (Fractional a, Ord a) => Int -> a -> (a -> a) -> Poly a
approxPoly n eps f =
  let xs = [-1] ++ legendreRoots (n - 2) eps ++ [1]
      ys = map f xs
  in lagrangePolyFit (zip xs ys)



initial :: (Floating a, Ord a) => Poly a
initial = approxPoly 10 1.0e-12 gaussian



main :: IO ()
main = let f0 = FunPoly initial
           f1 = FunHask sigmoid
           f = f1 . f0
           xs = [-1 + i / 10 | i <- [0..20]] :: [Double]
           ys = map (eval f) xs
       in mapM_ putStrLn $ show <$> zip [0..] (zip xs ys)