text:
{-# LANGUAGE TypeOperators, DeriveDataTypeable, FlexibleInstances, ScopedTypeVariables, ImplicitParams, PackageImports, ViewPatterns #-}
module Target.Heston where

import Target.Prelude
import qualified Data.Record as R
import Prelude hiding (scanl,repeat,replicate)
import Data.Record.Combinators ((!!!))
import Data.Kind
import Data.List (transpose)
import Data.TypeFun
import Numeric.LinearAlgebra hiding (diag, linspace, svd, )
import Math.Probably.Sampler hiding (uniform,primOneOf,logNormal,invGamma,binomial,gamma,oneOf,bernoulli, normal, unit, unormal,multiNormal)
import qualified Data.Text as T
import Foreign.Storable (Storable)
import Data.STRef
import Control.Monad.ST
import Control.Monad (forM, forM_)
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.Vector.Storable as VS
import Foreign.Storable.Tuple

(.==.) :: (Eq a) => (a -> ((a -> Double)))
(.==.) = \((x::a)) -> \((y::a)) -> if (x==y) then 1.000 else 0.000

scanl :: (((a -> ((b -> a)))) -> ((a -> ((([b]) -> ([a]))))))
scanl = \(((_arg0)::(a -> ((b -> a))))) -> \((_arg1)) -> \(((_arg2)::[b])) -> case (_arg0,((_arg1,_arg2))) of {((op::(a -> ((b -> a)))),(acc,[])) -> acc:[]; ((op::(a -> ((b -> a)))),(acc,(x:(xs::[b])))) -> acc:((scanl op (op acc x)) xs)}

fromTo :: (BayNum a) => (a -> ((a -> ([a]))))
fromTo = \((n::a)) -> \((m::a)) -> if (n<m) then (n:(fromTo (n+1) m)) else (n:[])

bigSum :: (Int -> ((Int -> ((((Int -> Double)) -> Double)))))
bigSum = \((lo::Int)) -> \((hi::Int)) -> \((f::(Int -> Double))) -> sum (map f (fromTo lo hi))

replicate :: (Int -> ((a -> ([a]))))
replicate = \(((_arg0)::Int)) -> \((_arg1)) -> case (_arg0,_arg1) of {(0,x) -> []; ((n::Int),x) -> x:(replicate (n-1) x)}

neg :: (BayNum a) => (a -> a)
neg = \((x::a)) -> 0-x

ix :: (Int -> ((([a]) -> a)))
ix = \(((_arg0)::Int)) -> \(((_arg1)::[a])) -> case (_arg0,_arg1) of {(0,(x:_)) -> x; ((n::Int),(_:(xs::[a]))) -> ix (n-1) xs}

linspace :: (Double -> ((Double -> ((Double -> ([Double]))))))
linspace = \((from::Double)) -> \((to::Double)) -> \((num::Double)) -> let {(dt::Double) = (to-from)/num;
 (is::[Double]) = fromTo 0 (num-1);
 } in map (\((i::Double)) -> ((unround i)*dt)+from) is

countSamples :: ((Prob a) -> (Maybe Int))
countSamples = \(((_arg0)::Prob a)) -> case _arg0 of {Samples (xs::[a]) -> Just (length xs); Sampler _ -> Nothing}

unit :: Prob Double
unit = Sampler primUnit

primOneOf :: (([a]) -> ((Seed -> ((a,Seed)))))
primOneOf = \((xs::[a])) -> \((seed::Seed)) -> let {((u::Double),(nextSeed::Seed)) = primUnit seed;
 (idx::Int) = floor (u*(unround (length xs)));
 } in ((ix idx xs),nextSeed)

append :: (([a]) -> ((([a]) -> ([a]))))
append = \(((_arg0)::[a])) -> \(((_arg1)::[a])) -> case (_arg0,_arg1) of {([],(ys::[a])) -> ys; ((x:(xs::[a])),(ys::[a])) -> x:(append xs ys)}

invlogit :: (Double -> Double)
invlogit = \((x::Double)) -> 1/(1+(exp (0.000-x)))

logit :: (Double -> Double)
logit = \((x::Double)) -> log (x/(1-x))

boolToReal :: (Bool -> Double)
boolToReal = \(((_arg0)::Bool)) -> case _arg0 of {True  -> 1.000; False  -> 0.000}

for :: (Int -> ((Int -> ((((Int -> (Prob a))) -> (Prob ([a])))))))
for = \((n::Int)) -> \((m::Int)) -> \((s::(Int -> (Prob a)))) -> if (n<m) then ((s n)>>=(\(x) -> ((for (n+1) m) s)>>=(\((xs::[a])) -> return (x:xs)))) else ((s n)>>=(\(v) -> return (v:[])))

repeat :: (Int -> (((Prob a) -> (Prob ([a])))))
repeat = \((n::Int)) -> \((sam::Prob a)) -> (for 1 n) (\((i::Int)) -> sam)

square :: (BayNum a) => (a -> a)
square = \((x::a)) -> x*x

step :: (BayNum a,BayNum b) => (a -> b)
step = \((x::a)) -> if (x<0) then 0 else 1

fac :: (Int -> Int)
fac = \(((_arg0)::Int)) -> case _arg0 of {1 -> 1; (n::Int) -> n*(fac (n-1))}

zipWithNats :: (BayNum b) => (([a]) -> ((b -> ([(b,a)]))))
zipWithNats = \(((_arg0)::[a])) -> \(((_arg1)::b)) -> case (_arg0,_arg1) of {([],_) -> []; ((x:(xs::[a])),(n::b)) -> ((n,x)):(zipWithNats xs (n+1))}

unSamples :: ((Prob a) -> ([a]))
unSamples = \(Samples (xs::[a])) -> xs

chainPlot :: ((Prob Double) -> Plot)
chainPlot = \(Samples (xs::[Double])) -> Plot [] (return ((Points (zipWithNats xs 0)):[]))

style :: (([(T.Text,T.Text)]) -> ((Plot -> Plot)))
style = \((opts::[(T.Text,T.Text)])) -> \(Plot (pos::[(T.Text,T.Text)]) (plr::Prob ([Radian]))) -> Plot pos (fmap (\((lrs::[Radian])) -> (Options opts lrs):[]) plr)

distPlot0 :: ((Prob Double) -> Plot)
distPlot0 = \(((_arg0)::Prob Double)) -> case _arg0 of {Samples (xs::[Double]) -> Plot ((((T.pack "range-y"),(T.pack "0"))):[]) (return ((Histogram xs):[])); (sampler::Prob Double) -> Plot ((((T.pack "range-y"),(T.pack "0"))):[]) ((repeat 2000 sampler)>>=(\((xs::[Double])) -> return ((Histogram xs):[])))}

distPlot :: ((Prob Double) -> Plot)
distPlot = \((p::Prob Double)) -> style ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "fill-opacity"),(T.pack "0.3"))):((((T.pack "stroke"),(T.pack "none"))):((((T.pack "bar-width"),(T.pack "0.8"))):[])))) (distPlot0 p)

histogram :: (BayNum a) => (([a]) -> Plot)
histogram = \((xs::[a])) -> Plot [] (return ((Options ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "fill-opacity"),(T.pack "0.3"))):((((T.pack "stroke"),(T.pack "none"))):((((T.pack "bar-width"),(T.pack "0.8"))):[])))) ((Histogram (map unround xs)):[])):[]))

unPlot :: (Plot -> (Prob ([Radian])))
unPlot = \(Plot _ (x::Prob ([Radian]))) -> x

over :: (([Plot]) -> Plot)
over = \((plots::[Plot])) -> Plot [] ((mapM unPlot plots)>>=(\((items::[[Radian]])) -> return (map (\((Plot (os::[(T.Text,T.Text)]) _,(rdns::[Radian]))) -> Options os rdns) (zip plots items))))

scatterPlot :: (BayNum a,BayNum b) => (([(a,b)]) -> Plot)
scatterPlot = \((xys::[(a,b)])) -> style ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "marker-size"),(T.pack "30"))):((((T.pack "marker"),(T.pack "circle"))):((((T.pack "stroke"),(T.pack "none"))):[])))) (Plot [] (return ((Points (map (\(((x::a),(y::b))) -> ((unround x),(unround y))) xys)):[])))

linePlot :: (BayNum a,BayNum b) => (([(a,b)]) -> Plot)
linePlot = \((xys::[(a,b)])) -> style ((((T.pack "stroke"),(T.pack "pop_colour"))):((((T.pack "stroke-width"),(T.pack "2"))):[])) (Plot [] (return ((Lines (map (\(((x::a),(y::b))) -> ((unround x),(unround y))) xys)):[])))

plines :: (BayNum a,BayNum b) => ((Prob ([(a,b)])) -> Plot)
plines = \((plns::Prob ([(a,b)]))) -> style ((((T.pack "stroke-opacity"),(T.pack "0.2"))):((((T.pack "stroke-width"),(T.pack "2"))):[])) (Plot [] (repeat 50 (fmap (\((xys::[(a,b)])) -> Lines (map (\(((x::a),(y::b))) -> ((unround x),(unround y))) xys)) plns)))

ppoints :: ((Prob ([(Double,Double)])) -> Plot)
ppoints = \((ppts::Prob ([(Double,Double)]))) -> Plot ((((T.pack "fill"),(T.pack "pop_colour"))):[]) (repeat 50 (fmap Points ppts))

sigPlot :: (((Double -> Double)) -> Plot)
sigPlot = \((sig::(Double -> Double))) -> Plot ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "stroke-width"),(T.pack "2"))):[])) (return ((Timeseries sig):[]))

thin :: (Int -> ((([a]) -> ([a]))))
thin = \((skip::Int)) -> \((xs::[a])) -> map snd (filter (\(((i::Int),x)) -> (mod i (skip+1))==0) (zip (fromTo 0 ((length xs)-1)) xs))

thinTo :: (Int -> ((([a]) -> ([a]))))
thinTo = \((n::Int)) -> \((xs::[a])) -> let {(nxs::Int) = length xs;
 (ratio::Int) = round ((unround nxs)/(unround n));
 } in thin ratio xs

psigPlot :: ((Prob ((Double -> Double))) -> Plot)
psigPlot = \(((_arg0)::Prob ((Double -> Double)))) -> case _arg0 of {Samples (sigs::[(Double -> Double)]) -> Plot ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "stroke-opacity"),(T.pack "0.2"))):((((T.pack "stroke-width"),(T.pack "2"))):[]))) (return (map Timeseries (thinTo 20 sigs))); Sampler (f::(Seed -> ((((Double -> Double)),Seed)))) -> Plot ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "stroke-opacity"),(T.pack "0.2"))):((((T.pack "stroke-width"),(T.pack "2"))):[]))) ((repeat 20 (Sampler f))>>=(\((sigs::[(Double -> Double)])) -> return (map Timeseries sigs)))}

psigNPlot :: (Int -> (((Prob ((Double -> Double))) -> Plot)))
psigNPlot = \(((_arg0)::Int)) -> \(((_arg1)::Prob ((Double -> Double)))) -> case (_arg0,_arg1) of {((n::Int),Samples (sigs::[(Double -> Double)])) -> Plot ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "stroke-opacity"),(T.pack "0.2"))):((((T.pack "stroke-width"),(T.pack "2"))):[]))) (return (map Timeseries (thinTo n sigs))); ((n::Int),Sampler (f::(Seed -> ((((Double -> Double)),Seed))))) -> Plot ((((T.pack "fill"),(T.pack "pop_colour"))):((((T.pack "stroke-opacity"),(T.pack "0.1"))):((((T.pack "stroke-width"),(T.pack "2"))):[]))) ((repeat n (Sampler f))>>=(\((sigs::[(Double -> Double)])) -> return (map Timeseries sigs)))}

probBoolToP :: ((Prob Bool) -> (Prob Double))
probBoolToP = \(((_arg0)::Prob Bool)) -> case _arg0 of {Sampler (f::(Seed -> ((Bool,Seed)))) -> (repeat 200 (Sampler f))>>=(\((bs::[Bool])) -> probBoolToP (Samples bs)); Samples (bs::[Bool]) -> let {(yeas::[Bool]) = filter id bs;
 } in return ((unround (length yeas))/(unround (length bs)))}

pcurve :: (([(Double,(Prob Bool))]) -> Plot)
pcurve = \((xps::[(Double,(Prob Bool))])) -> Plot [] (let {(xs::[Double]) = map fst xps;
 } in (mapM (probBoolToP.snd) xps)>>=(\((ps::[Double])) -> return ((Lines (zip xs ps)):[])))

plotStyle :: (([(T.Text,T.Text)]) -> ((Plot -> Plot)))
plotStyle = \(((_arg0)::[(T.Text,T.Text)])) -> \(((_arg1)::Plot)) -> case (_arg0,_arg1) of {((opts::[(T.Text,T.Text)]),Plot (pos::[(T.Text,T.Text)]) (plr::Prob ([Radian]))) -> Plot (append opts pos) plr; ((opts::[(T.Text,T.Text)]),PlotRow (pos::[(T.Text,T.Text)]) (plr::[Plot])) -> PlotRow (append opts pos) plr; ((opts::[(T.Text,T.Text)]),PlotColumn (pos::[(T.Text,T.Text)]) (plr::[Plot])) -> PlotColumn (append opts pos) plr; ((opts::[(T.Text,T.Text)]),PlotStack (pos::[(T.Text,T.Text)]) (plr::[(T.Text,Plot)])) -> PlotStack (append opts pos) plr; ((opts::[(T.Text,T.Text)]),PlotGrid (pos::[(T.Text,T.Text)]) (x::Int) (y::Int) (plr::[(T.Text,Plot)])) -> ((PlotGrid (append opts pos) x) y) plr}

wide :: (Plot -> Plot)
wide = plotStyle ((((T.pack "width"),(T.pack "750"))):[])

aspect :: (Double -> ((Plot -> Plot)))
aspect = \((x::Double)) -> plotStyle ((((T.pack "aspect"),(showReal x))):[])

besides :: (([Plot]) -> Plot)
besides = PlotRow []

above :: (([Plot]) -> Plot)
above = PlotColumn []

axisLabels :: (T.Text -> ((T.Text -> ((Plot -> Plot)))))
axisLabels = \((xlab::T.Text)) -> \((ylab::T.Text)) -> \(Plot (popts::[(T.Text,T.Text)]) (plns::Prob ([Radian]))) -> Plot ((((T.pack "axis-x-label"),xlab)):((((T.pack "axis-y-label"),ylab)):popts)) plns

unPlotOpts :: (Plot -> ([(T.Text,T.Text)]))
unPlotOpts = \(Plot (os::[(T.Text,T.Text)]) (x::Prob ([Radian]))) -> os

sigLast :: (((Double -> Double)) -> Double)
sigLast = \((sig::(Double -> Double))) -> case observeSig sig of {ObservedSignal (dt::Double) (t0::Double) (pts::Vector Double) -> pts@>((dim pts)-1); ObservedXYSignal (pts::Vector ((Double,Double))) -> snd (pts@>((dim pts)-1))}

sigTail :: (((Double -> Double)) -> ([Double]))
sigTail = \((sig::(Double -> Double))) -> case observeSig sig of {ObservedSignal (dt::Double) (t0::Double) (pts::Vector Double) -> vecToList pts}

sigNPts :: (((Double -> Double)) -> Int)
sigNPts = \((sig::(Double -> Double))) -> case observeSig sig of {ObservedSignal (dt::Double) (t0::Double) (pts::Vector Double) -> dim pts; ObservedXYSignal (pts::Vector ((Double,Double))) -> dim pts}

between :: (BayNum a) => (a -> ((a -> ((a -> Bool)))))
between = \((lo::a)) -> \((hi::a)) -> \((x::a)) -> (x>lo)&&(x<hi)

runP :: (a -> ((((a -> ((c,d)))) -> c)))
runP = \(xs) -> \((my::(a -> ((c,d))))) -> fst (my xs)

runP1 :: (((a -> ((b,c)))) -> ((a -> b)))
runP1 = \((my::(a -> ((b,c))))) -> \(xs) -> fst (my xs)

returnP :: (a -> ((b -> ((a,b)))))
returnP = \(x) -> \(xs) -> (x,xs)

bindP :: (((a -> ((b,c)))) -> ((((b -> ((c -> f)))) -> ((a -> f)))))
bindP = \((f::(a -> ((b,c))))) -> \((g::(b -> ((c -> f))))) -> \(xs) -> let {(x,xs') = f xs;
 } in g x xs'

headP :: (([a]) -> ((a,([a]))))
headP = \(((_arg0)::[a])) -> case _arg0 of {(x:(xs::[a])) -> (x,xs); [] -> bayError (T.pack "headP: empty list")}

takeP :: (Int -> ((([a]) -> ((([a]),([a]))))))
takeP = \((n::Int)) -> \((xs::[a])) -> ((take n xs),(drop n xs))

forP :: (([a]) -> ((((a -> ((c -> ((d,c)))))) -> ((c -> ((([d]),c)))))))
forP = \(((_arg0)::[a])) -> \(((_arg1)::(a -> ((c -> ((d,c))))))) -> \((_arg2)) -> case (_arg0,((_arg1,_arg2))) of {([],((f::(a -> ((c -> ((d,c)))))),xs)) -> ([],xs); ((a:(as::[a])),((f::(a -> ((c -> ((d,c)))))),xs)) -> (bindP (f a) (\(y) -> bindP (forP as f) (\((ys::[d])) -> returnP (y:ys)))) xs}

noP :: (a -> (((),a)))
noP = \(xs) -> ((),xs)

fmapP :: (((a -> b)) -> ((((c -> ((a,e)))) -> ((c -> ((b,e)))))))
fmapP = \((f::(a -> b))) -> \((mx::(c -> ((a,e))))) -> bindP mx (\(x) -> returnP (f x))

fixP :: (Int -> (((Prob a) -> (Prob (Prob a)))))
fixP = \(((_arg0)::Int)) -> \(((_arg1)::Prob a)) -> case (_arg0,_arg1) of {((n::Int),Sampler (f::(Seed -> ((a,Seed))))) -> (repeat n (Sampler f))>>=(\((xs::[a])) -> return (Samples xs)); ((n::Int),Samples (xs::[a])) -> return (Samples (thinTo n xs))}

diag :: ((Vector Double) -> (Matrix Double))
diag = \((v::Vector Double)) -> fillM (((dim v),(dim v))) (\(((i::Int),(j::Int))) -> if (i==j) then (v@>i) else 0.000)

diagL :: (([Double]) -> (Matrix Double))
diagL = diag.listToVec

mmap :: (BayNum a,BayNum b) => (((a -> b)) -> (((Matrix a) -> (Matrix b))))
mmap = \((f::(a -> b))) -> \((m::Matrix a)) -> fillM (mdims m) (\(((i::Int),(j::Int))) -> f (m@@>((i,j))))

transM :: (BayNum a) => ((Matrix a) -> (Matrix a))
transM = \((m::Matrix a)) -> fillM (mdims m) (\(((i::Int),(j::Int))) -> m@@>((j,i)))

ident :: (BayNum a) => (a -> (Matrix Double))
ident = \((n::a)) -> diagL (map (const 1) (fromTo 1 n))

meanL :: (([Double]) -> Double)
meanL = \((xs::[Double])) -> (sum xs)/(unround (length xs))

varL :: (([Double]) -> Double)
varL = \((xs::[Double])) -> let {(mu::Double) = meanL xs;
 } in (sum (map (\((x::Double)) -> (x-mu)*(x-mu)) xs))/(unround ((length xs)-1))

popvarL :: (([Double]) -> Double)
popvarL = \((xs::[Double])) -> let {(mu::Double) = meanL xs;
 } in (sum (map (\((x::Double)) -> (x-mu)*(x-mu)) xs))/(unround (length xs))

minL :: (([Double]) -> Double)
minL = \(((x::Double):(xs::[Double]))) -> (foldl min x) xs

maxL :: (([Double]) -> Double)
maxL = \(((x::Double):(xs::[Double]))) -> (foldl max x) xs

getT :: ((a) -> (((b) -> ((a,b)))))
getT = \(dt) -> \(tmax) -> (dt,tmax)

pearson :: (([(Double,Double)]) -> Double)
pearson = \((xys::[(Double,Double)])) -> let {(xs::[Double]) = map fst xys;
 (ys::[Double]) = map snd xys;
 (mx::Double) = meanL xs;
 (my::Double) = meanL ys;
 (sx::Double) = sqrt (varL xs);
 (sy::Double) = sqrt (varL ys);
 (invN::Double) = 1.000/(unround ((length xys)-1));
 } in invN*(sum (map (\(((x::Double),(y::Double))) -> ((x-mx)/sx)*((y-my)/sy)) xys))

data Strategy a =
   GStrategy (((((Vector Double) -> ((Double,(Vector Double))))) -> (((Vector Double) -> ((a -> ((((Double,(Vector Double))) -> (Prob (((((Vector Double),a)),((Double,(Vector Double)))))))))))))) (((Vector Double) -> a))
   |VStrategy (((((Vector Double) -> Double)) -> (((Vector Double) -> ((a -> ((Double -> (Prob (((((Vector Double),a)),Double))))))))))) (((Vector Double) -> a))
intBetweenLogPdf :: (BayNum d) => (a -> ((b -> ((c -> d)))))
intBetweenLogPdf = \(lo) -> \(hi) -> \(x) -> 1

anyLogPdf :: (a -> Double)
anyLogPdf = \(x) -> 1.000

d :: ((Double) -> ((((Double -> Double)) -> ((Double -> Double)))))
d = \((dt::Double)) -> \((w::(Double -> Double))) -> \((t::Double)) -> ((w t)-(w (t-dt)))/dt

unormal :: Prob Double
unormal = unit>>=(\((u1::Double)) -> unit>>=(\((u2::Double)) -> return ((sqrt ((0.000-2.000)*(log u1)))*(cos ((2.000*pi)*u2)))))

gammaAux :: (Double -> ((Double -> (Prob Double))))
gammaAux = \((a::Double)) -> \((b::Double)) -> let {(d::Double) = a-(1.000/3.000);
 (c::Double) = 1.000/(3.000*(sqrt d));
 } in unormal>>=(\((x::Double)) -> let {(cx::Double) = c*x;
 (v::Double) = (1.000+cx)**3.000;
 (x_2::Double) = x*x;
 (x_4::Double) = x_2*x_2;
 } in if (cx<(-1.000)) then (gammaAux a b) else (unit>>=(\((u::Double)) -> if ((u<(1.000-(3.310e-2*x_4)))||((log u)<((0.500*x_2)+(d*((1.000-v)+(log v)))))) then (return ((b*d)*v)) else (gammaAux a b))))

gamma :: (Double -> ((Double -> (Prob Double))))
gamma = \((k::Double)) -> \((theta::Double)) -> if (k<1.000) then (unit>>=(\((u::Double)) -> (gamma (1.000+k) theta)>>=(\((x::Double)) -> return (x*(u**(1.000/k)))))) else (gammaAux k theta)

improper_uniform :: Prob Double
improper_uniform = gamma 1 0.100

improper_uniformLogPdf :: (a -> Double)
improper_uniformLogPdf = \(_) -> 1.000

improper_uniform_positive :: Prob Double
improper_uniform_positive = gamma 1 1

unfoldN :: (Int -> ((Int -> ((a -> ((((Int -> ((a -> (Prob a))))) -> (Prob ([a])))))))))
unfoldN = \((n::Int)) -> \((m::Int)) -> \(lastx) -> \((s::(Int -> ((a -> (Prob a)))))) -> if (n<m) then ((s n lastx)>>=(\(x) -> (((unfoldN (n+1) m) x) s)>>=(\((xs::[a])) -> return (x:xs)))) else ((s n lastx)>>=(\(v) -> return (v:[])))

unfold :: (Int -> ((a -> ((((a -> (Prob a))) -> (Prob ([a])))))))
unfold = \((n::Int)) -> \(x0) -> \((s::(a -> (Prob a)))) -> ((unfoldN 1 n) x0) (\((i::Int)) -> s)

improper_uniformInit :: Double
improper_uniformInit = 1.000

ser :: Double
ser = 1.000

normal :: (Double -> ((Double -> (Prob Double))))
normal = \((mean::Double)) -> \((variance::Double)) -> unormal>>=(\((u::Double)) -> return ((u*(sqrt variance))+mean))

rwmTrans :: (BayNum a) => ((((Vector Double) -> Double)) -> (((Vector Double) -> ((((Double,((Double,a)))) -> ((Double -> (Prob (((((Vector Double),((Double,((Double,a)))))),Double))))))))))
rwmTrans = \((posterior::((Vector Double) -> Double))) -> \((xi::Vector Double)) -> \(((sigma::Double),((i::Double),(iaccept::a)))) -> \((pi::Double)) -> (fmap listToVec (mapM (\((x::Double)) -> normal x sigma) (vecToList xi)))>>=(\((xstar::Vector Double)) -> let {(pstar::Double) = posterior xstar;
 (ratio::Double) = exp (pstar-pi);
 } in unit>>=(\((u::Double)) -> let {(accept::Bool) = u<ratio;
 (sigmaNext::Double) = if accept then (sigma*((min 1.400 (1+(5.000/i)))**3)) else (sigma*(max 0.714 (1-(5.000/i))));
 } in return (if accept then ((((xstar,((sigmaNext,(((i+1),(iaccept+1))))))),pstar)) else ((((xi,((sigmaNext,(((i+1),iaccept)))))),pi)))))

rwmIni :: (a -> ((Double,((Double,Double)))))
rwmIni = \(_) -> (0.100,((1.000,0.000)))

rwm :: Strategy ((Double,((Double,Double))))
rwm = VStrategy rwmTrans rwmIni

malaTrans :: ((((Vector Double) -> ((Double,a)))) -> (((Vector Double) -> ((Double -> ((((Double,a)) -> (Prob (((((Vector Double),Double)),((Double,a))))))))))))
malaTrans = \((postgrad::((Vector Double) -> ((Double,a))))) -> \((xi::Vector Double)) -> \((sigma::Double)) -> \(((pi::Double),gradienti)) -> let {(xstarMean::Vector Double) = xi;
 } in (fmap listToVec (mapM (\((x::Double)) -> normal x sigma) (vecToList xi)))>>=(\((xstar::Vector Double)) -> let {((pstar::Double),gradientStar) = postgrad xstar;
 (ratio::Double) = exp (pstar-pi);
 } in unit>>=(\((u::Double)) -> let {(accept::Bool) = u<ratio;
 (sigmaNext::Double) = sigma;
 } in return (if accept then ((((xstar,sigmaNext)),((pstar,gradientStar)))) else ((((xi,sigmaNext)),((pi,gradienti)))))))

malaIni :: (BayNum b) => (a -> b)
malaIni = \(vini) -> 1

mala :: Strategy Double
mala = GStrategy malaTrans malaIni

uniform :: (Double -> ((Double -> (Prob Double))))
uniform = \((lo::Double)) -> \((hi::Double)) -> unit>>=(\((x::Double)) -> return ((x*(hi-lo))+lo))

oneOf :: (([a]) -> (Prob a))
oneOf = \((xs::[a])) -> (fmap floor (uniform 0.000 (unround (length xs))))>>=(\((idx::Int)) -> return (ix idx xs))

oneOfLogPdf :: (([a]) -> ((b -> Double)))
oneOfLogPdf = \((xs::[a])) -> \(_) -> 1.000/(unround (length xs))

uniformLogPdf :: (Double -> ((Double -> ((Double -> Double)))))
uniformLogPdf = \((lo::Double)) -> \((hi::Double)) -> \((x::Double)) -> if ((x<hi)&&(x>lo)) then ((log 1)-(log (hi-lo))) else (0.000-1.000e20)

intBetween :: (BayNum a) => (a -> ((a -> (Prob Int))))
intBetween = \((lo::a)) -> \((hi::a)) -> unit>>=(\((x::Double)) -> return (floor ((x*(unround ((hi+1)-lo)))+(unround lo))))

any :: Prob a
any = undefined

improper_uniform_positiveLogPdf :: (BayNum a) => (a -> Double)
improper_uniform_positiveLogPdf = \((x::a)) -> if (x>0) then 1.000 else (0.000-1.000e20)

oneTo :: (BayNum a) => (a -> (Prob Int))
oneTo = \((n::a)) -> (uniform 0.500 ((unround n)+0.500))>>=(\((x::Double)) -> return (round x))

oneToLogPdf :: (BayNum a) => (a -> ((a -> Double)))
oneToLogPdf = \((hi::a)) -> \((x::a)) -> if ((x<(hi+1))&&(x>0)) then 1.000 else (0.000-1.000e10)

normalLogPdf :: (Double -> ((Double -> ((Double -> Double)))))
normalLogPdf = \((mean::Double)) -> \((variance::Double)) -> \((x::Double)) -> ((log 1)-(0.500*(log ((2.000*pi)*variance))))-(((x-mean)**2)/(2*variance))

logNormal :: (Double -> ((Double -> (Prob Double))))
logNormal = \((mean::Double)) -> \((variance::Double)) -> fmap exp (normal mean variance)

logNormalLogPdf :: (Double -> ((Double -> ((Double -> Double)))))
logNormalLogPdf = \((mean::Double)) -> \((variance::Double)) -> \((x::Double)) -> (log (1/(sqrt ((2.000*pi)*variance))))+(0.000-((((log x)-mean)*((log x)-mean))/(2*variance)))

normalLines :: (Double -> ((Double -> ([(Double,Double)]))))
normalLines = \((mean::Double)) -> \((v::Double)) -> let {(f::(Double -> ((Double,Double)))) = \((x::Double)) -> (x,(exp ((normalLogPdf mean v) x)));
 } in map f ((linspace (mean-(3*(sqrt v))) (mean+(3*(sqrt v)))) 50)

binomialProb :: (Int -> ((Double -> ((Int -> Double)))))
binomialProb = \((n::Int)) -> \((p::Double)) -> \((k::Int)) -> ((choose n k)*(p^k))*((1.000-p)^(n-k))

binomialLogProb :: (Int -> ((Double -> ((Int -> Double)))))
binomialLogProb = \((n::Int)) -> \((p::Double)) -> \((k::Int)) -> ((log (choose n k))+((unround k)*(log p)))+((unround (n-k))*(log (1.000-p)))

bernoulli :: (Double -> (Prob Bool))
bernoulli = \((p::Double)) -> unit>>=(\((u::Double)) -> return (u<p))

bernoulli01 :: (Double -> (Prob Int))
bernoulli01 = \((p::Double)) -> unit>>=(\((u::Double)) -> if (u<p) then (return 1) else (return 0))

bernoulliLogPdf :: (Double -> ((Bool -> Double)))
bernoulliLogPdf = \((p::Double)) -> \((b::Bool)) -> log (((p*((2*(boolToReal b))-1))+1)-(boolToReal b))

bernoulli01LogPdf :: (Double -> ((Double -> Double)))
bernoulli01LogPdf = \((p::Double)) -> \((b::Double)) -> log (((p*((2*b)-1))+1)-b)

countTrue :: (BayNum a) => (([Bool]) -> a)
countTrue = \(((_arg0)::[Bool])) -> case _arg0 of {[] -> 0; (True :(bs::[Bool])) -> 1+(countTrue bs); (False :(bs::[Bool])) -> countTrue bs}

binomial :: (BayNum a) => (Int -> ((Double -> (Prob a))))
binomial = \((n::Int)) -> \((p::Double)) -> (repeat n (unit>>=(\((u::Double)) -> return (u<p))))>>=(\((bools::[Bool])) -> return (countTrue bools))

exponential :: (Double -> (Prob Double))
exponential = \((lam::Double)) -> unit>>=(\((u::Double)) -> return (neg ((log u)/lam)))

exponentialLogPdf :: (Double -> ((Double -> Double)))
exponentialLogPdf = \((lam::Double)) -> \((x::Double)) -> lam*(exp ((0.000-lam)*x))

poissonAux :: (BayNum a) => (Double -> ((a -> ((Double -> (Prob a))))))
poissonAux = \((bigl::Double)) -> \((k::a)) -> \((p::Double)) -> if (p>bigl) then (unit>>=(\((u::Double)) -> (poissonAux bigl (k+1)) (p*u))) else (return (k-1))

poisson :: (Double -> (Prob Int))
poisson = \((lam::Double)) -> (poissonAux (exp (0.000-lam)) 0) 1

poissonLogPdf :: (Double -> ((Int -> Double)))
poissonLogPdf = \((lam::Double)) -> \((x::Int)) -> ((lam**(unround x))*(exp (0.000-lam)))/(unround (fac x))

betaAux :: (Int -> (Prob Double))
betaAux = \((n::Int)) -> (repeat n unit)>>=(\((us::[Double])) -> return (log (product us)))

beta :: (Int -> ((Int -> (Prob Double))))
beta = \((a::Int)) -> \((b::Int)) -> (betaAux a)>>=(\((g1::Double)) -> (betaAux b)>>=(\((g2::Double)) -> return (g1/(g1+g2))))

cof :: [Double]
cof = 76.180:((-86.505):(24.014:((-1.232):(1.209e-3:((-5.395e-6):[])))))

gammaln :: (Double -> Double)
gammaln = \((xx::Double)) -> let {(tmp'::Double) = (xx+5.500)-((xx+0.500)*(log (xx+5.500)));
 (ser'::Double) = ser+(sum (map (\(((y::Double),(c::Double))) -> c/(xx+y)) (zip (fromTo 1 7) cof)));
 } in (0.000-tmp')+(log ((2.507*ser')/xx))

betaf :: (Double -> ((Double -> Double)))
betaf = \((z::Double)) -> \((w::Double)) -> exp (((gammaln z)+(gammaln w))-(gammaln (z+w)))

betaLogPdf :: (Int -> ((Int -> ((Double -> Double)))))
betaLogPdf = \((a::Int)) -> \((b::Int)) -> \((x::Double)) -> log (((1.000/(betaf (unround a) (unround b)))*(x^(a-1)))*((1.000-x)^(b-1)))

invGamma :: (Double -> ((Double -> (Prob Double))))
invGamma = \((a::Double)) -> \((b::Double)) -> (gamma a (1.000/b))>>=(\((g::Double)) -> return (1.000/g))

wiener :: ((Double) -> (((Double) -> (Prob ((Double -> Double))))))
wiener = \((dt::Double)) -> \((tmax::Double)) -> let {(n::Int) = (round (tmax/dt))+1;
 } in (repeat n unormal)>>=(\((ns::[Double])) -> let {(etas::[Double]) = map (\((u::Double)) -> u*(sqrt dt)) ns;
 } in return ((pack dt (0.000-dt)) (listToVec ((scanl (\((x::Double)) -> \((y::Double)) -> x+y) 0.000) etas))))

multiNormal :: ((Vector Double) -> (((Matrix Double) -> (Prob (Vector Double)))))
multiNormal = \((vmean::Vector Double)) -> \((cov::Matrix Double)) -> (repeat (dim vmean) (normal 0 1))>>=(\((ns::[Double])) -> let {((u::Matrix Double):((s::Matrix Double):((v::Matrix Double):[]))) = svd cov;
 (j::Matrix Double) = mXm (mXm v (mmap sqrt s)) (transM v);
 } in return (vmean+(mXv j (listToVec ns))))

wieners :: ((Double) -> (((Double) -> (((Matrix Double) -> (Prob ([(Double -> Double)])))))))
wieners = \((dt::Double)) -> \((tmax::Double)) -> \((cov::Matrix Double)) -> let {(n_time_pts::Int) = (round (tmax/dt))+2;
 (n_dims::Int) = fst (mdims cov);
 (zeroV::Vector Double) = fillV n_dims (const 0.000);
 (scale_cov::Matrix Double) = dt*%cov;
 } in ((unfold n_time_pts zeroV) (\((vlast::Vector Double)) -> (multiNormal zeroV scale_cov)>>=(\((etas::Vector Double)) -> return (vlast+etas))))>>=(\((vs::[Vector Double])) -> let {(ll::[[Double]]) = transL (map vecToList vs);
 } in return (map (\((ys::[Double])) -> (pack dt (0.000-dt)) (listToVec ys)) ll))

diff :: ((Double) -> ((((Double -> Double)) -> ((Double -> Double)))))
diff = \((dt::Double)) -> \((w::(Double -> Double))) -> \((t::Double)) -> ((w t)-(w (t-dt)))/dt

decide :: (Double -> (((Vector Double) -> (((Prob a) -> (((((Vector Double) -> ((a -> Double)))) -> (Vector Double))))))))
decide = \((tol::Double)) -> \((ini::Vector Double)) -> \((dist::Prob a)) -> \((util::((Vector Double) -> ((a -> Double))))) -> (optimise tol ini) (\((vaction::Vector Double)) -> expect (fmap (util vaction) dist))

nsig :: (((Double -> Double)) -> ((Double -> (Prob ((Double -> Double))))))
nsig = \((sig::(Double -> Double))) -> \((v::Double)) -> case observeSig sig of {ObservedSignal (dt::Double) (t0::Double) (vpts::Vector Double) -> (repeat (dim vpts) (normal 0 v))>>=(\((ns::[Double])) -> return ((pack dt t0) (listToVec (map (\(((x::Double),(y::Double))) -> x+y) (zip ns (vecToList vpts))))))}

quantile :: (Double -> (((Prob Double) -> Double)))
quantile = \((x::Double)) -> \(Samples (xs::[Double])) -> let {(total::Int) = length xs;
 (under::Int) = length (filter (\((y::Double)) -> x<y) xs);
 } in (unround under)/(unround total)

gammaLogPdf :: (Double -> ((Double -> ((Double -> Double)))))
gammaLogPdf = \((k::Double)) -> \((theta::Double)) -> \((x::Double)) -> ((((k-1)*(log x))+((0.000-x)/theta))-(k*(log theta)))-(gammaln k)

invGammaLogPdf :: (Double -> ((Double -> ((Double -> Double)))))
invGammaLogPdf = \((a::Double)) -> \((b::Double)) -> \((x::Double)) -> log ((((b**a)/(exp (gammaln a)))*(x**((0.000-a)-1)))*(exp ((0.000-b)/x)))

cookAssert :: ((Prob Double) -> SamplerDensity)
cookAssert = \((s::Prob Double)) -> SamplerDensity s (uniformLogPdf 0 1)

multiNormalLogPdf :: (a -> ((b -> c)))
multiNormalLogPdf = \(vmean) -> \(cov) -> undefined

invWishart :: (Double -> (((Matrix Double) -> (Prob (Matrix Double)))))
invWishart = \((nu::Double)) -> \((s::Matrix Double)) -> undefined

invWishartLogPdf :: (a -> ((b -> ((c -> d)))))
invWishartLogPdf = \(nu) -> \(s) -> \(w) -> undefined

scandyn' :: (BayNum a) => (((Int -> ((odes -> ((sdes -> ((sdes -> ((a,odes)))))))))) -> ((a -> ((Int -> ((sdes -> ((odes -> ((([odes]) -> ((([sdes]) -> ((a,([odes]))))))))))))))))
scandyn' = \(((_arg0)::(Int -> ((odes -> ((sdes -> ((sdes -> ((a,odes))))))))))) -> \(((_arg1)::a)) -> \(((_arg2)::Int)) -> \((_arg3)) -> \((_arg4)) -> \(((_arg5)::[odes])) -> \(((_arg6)::[sdes])) -> case (_arg0,((_arg1,((_arg2,((_arg3,((_arg4,((_arg5,_arg6))))))))))) of {((f::(Int -> ((odes -> ((sdes -> ((sdes -> ((a,odes)))))))))),((p::a),(_,(_,(_,((odeacc::[odes]),[])))))) -> (p,(reverse odeacc)); ((f::(Int -> ((odes -> ((sdes -> ((sdes -> ((a,odes)))))))))),((p0::a),((i::Int),(sdelast,(odecurr,((odeacc::[odes]),(sde:(sdes::[sdes])))))))) -> let {((p1::a),odenext) = ((f i odecurr) sdelast) sde;
 } in (((((scandyn' f (p0+p1)) (i+1)) sde) odenext) (odenext:odeacc)) sdes}

scandyn :: (BayNum a) => (((Int -> ((odes -> ((sdes -> ((sdes -> ((a,odes)))))))))) -> ((odes -> ((sdes -> ((([sdes]) -> ((a,([odes]))))))))))
scandyn = \((f::(Int -> ((odes -> ((sdes -> ((sdes -> ((a,odes))))))))))) -> \(initode) -> \(initsde) -> \((sdes::[sdes])) -> (((((scandyn' f 0) 0) initsde) initode) []) sdes

tmax :: Double
tmax = 10.000

dt :: Double
dt = 1.000e-2

samples :: (BayNum a) => a
samples = 2000

heston :: Prob ((R.X R.:& V R.::: (Double -> Double) R.:& S R.::: (Double -> Double)) (Id KindStar))
heston = (gamma 1 1)>>=(\((k::Double)) -> (gamma 1 0.100)>>=(\((th::Double)) -> (gamma 1 0.100)>>=(\((eta::Double)) -> (normal 0 0.100)>>=(\((mu::Double)) -> (wiener dt tmax)>>=(\((w1::(Double -> Double))) -> (wiener dt tmax)>>=(\((w2::(Double -> Double))) -> (gamma 1 0.100)>>=(\((v_0::Double)) -> (uniform 0.000 2.000)>>=(\((s_0::Double)) -> let {v = solveODE (\v-> \((t::Double)) -> (k*(th-v))+((eta*(sqrt v))*((d dt w1) t))) tmax dt v_0;
s = solveODE (\s-> \((t::Double)) -> (mu*s)+(((sqrt (v t))*s)*((d dt w2) t))) tmax dt s_0;
} in return (R.X R.:& V R.:= (v) R.:& S R.:= (s))))))))))

heston1 :: Prob ((R.X R.:& V R.::: (Double -> Double) R.:& S R.::: (Double -> Double)) (Id KindStar))
heston1 = (return (R.X R.:& S_0 R.:= (1.000) R.:& V_0 R.:= (2.000e-2) R.:& Mu R.:= (5.000e-2) R.:& Eta R.:= (0.100) R.:& Th R.:= (2.000e-2) R.:& K R.:= (1.000))::Prob ((R.X R.:& S_0 R.::: Double R.:& V_0 R.::: Double R.:& Mu R.::: Double R.:& Eta R.::: Double R.:& Th R.::: Double R.:& K R.::: Double) (Id KindStar)))>>=(\(((_pars)::((R.X R.:& S_0 R.::: Double R.:& V_0 R.::: Double R.:& Mu R.::: Double R.:& Eta R.::: Double R.:& Th R.::: Double R.:& K R.::: Double) (Id KindStar)))) -> let {(k::Double) = _pars!!!K;
 } in let {(th::Double) = _pars!!!Th;
 } in let {(eta::Double) = _pars!!!Eta;
 } in let {(mu::Double) = _pars!!!Mu;
 } in let {(v_0::Double) = _pars!!!V_0;
 } in let {(s_0::Double) = _pars!!!S_0;
 } in (wiener dt tmax)>>=(\((_w_v)) -> (wiener dt tmax)>>=(\((_w_s)) -> let {v = solveODE (\v-> \((t::Double)) -> (k*(th-v))+((eta*(sqrt v))*((d dt _w_v) t))) tmax dt v_0;
s = solveODE (\s-> \((t::Double)) -> (mu*s)+(((sqrt (v t))*s)*((d dt _w_s) t))) tmax dt s_0;
} in return (R.X R.:& V R.:= (v) R.:& S R.:= (s)))))

data S = S deriving Show
instance R.Name S where
   name = S
data V = V deriving Show
instance R.Name V where
   name = V
data K = K deriving Show
instance R.Name K where
   name = K
data Th = Th deriving Show
instance R.Name Th where
   name = Th
data Eta = Eta deriving Show
instance R.Name Eta where
   name = Eta
data Mu = Mu deriving Show
instance R.Name Mu where
   name = Mu
data V_0 = V_0 deriving Show
instance R.Name V_0 where
   name = V_0
data S_0 = S_0 deriving Show
instance R.Name S_0 where
   name = S_0
data Posterior = Posterior deriving Show
instance R.Name Posterior where
   name = Posterior
data Postgrad = Postgrad deriving Show
instance R.Name Postgrad where
   name = Postgrad
data VToRec = VToRec deriving Show
instance R.Name VToRec where
   name = VToRec
data Inisam = Inisam deriving Show
instance R.Name Inisam where
   name = Inisam

target = do
  ((fakedata::((R.X R.:& V R.::: (Double -> Double) R.:& S R.::: (Double -> Double)) (Id KindStar))))::((R.X R.:& V R.::: (Double -> Double) R.:& S R.::: (Double -> Double)) (Id KindStar)) <- sample (heston1::Prob ((R.X R.:& V R.::: (Double -> Double) R.:& S R.::: (Double -> Double)) (Id KindStar)))
  --fail "done"
  let prims :: ((R.X R.:& Inisam R.::: Prob ([Double]) R.:& VToRec R.::: ((Vector Double) -> ((R.X R.:& V R.::: (Double -> Double) R.:& V_0 R.::: Double R.:& Mu R.::: Double R.:& Eta R.::: Double R.:& Th R.::: Double R.:& K R.::: Double) (Id KindStar))) R.:& Postgrad R.::: ((Vector Double) -> ((Double,(Vector Double)))) R.:& Posterior R.::: (([Double]) -> Double)) (Id KindStar))
      prims = let {s::(Double -> Double) = (fakedata!!!S);
 posterior = runP1 (bindP (fmapP (\((_x)) -> exp _x) headP) (\(k) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(th) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(eta) -> bindP headP (\(mu) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(v_0) -> bindP (fmapP (\(s) -> v_0:(map (\(_x ) -> exp _x) s)) (takeP$((round (tmax/dt))-1))) (\(v__obsSig) -> bindP (returnP ((packL dt 0) v__obsSig)) (\(v) -> returnP (((gammaLogPdf 1 1) k)+(((gammaLogPdf 1 0.100) th)+(((gammaLogPdf 1 0.100) eta)+(((normalLogPdf 0 0.100) mu)+(((gammaLogPdf 1 0.100) v_0)+(((uniformLogPdf 0.000 2.000) (s 0))+(let {scanfun = \(i) -> \([]) -> \((v_last:(s_last:[]))) -> \((v_next:(s_next:[]))) -> let {t = (unround i)*dt;
 v = const v_last;
 s = const s_last;
 pthis = (0+((normalLogPdf (v_last+(dt*(k*(th-(v t))))) (dt*((eta*(sqrt (v t)))*(eta*(sqrt (v t)))))) v_next))+((normalLogPdf (s_last+(dt*(mu*(s t)))) (dt*(((sqrt (v t))*(s t))*((sqrt (v t))*(s t))))) s_next);
 } in (pthis,[]);
 (pdynsys,dynsysres) = ((scandyn scanfun []) ((v 0):((s 0):[]))) (transL ((tail v__obsSig):((sigTail s):[])));
 } in pdynsys+0)))))))))))))));
 postgrad = \(_v ) -> runST$((newSTRef 0)>>=(\(postref) -> (VSM.replicate (dim _v) 0)>>=(\(gradref) -> (newSTRef (0::Int))>>=(\(countref) -> let {post_incr = addSTRef postref;
 grad_incr = addMV gradref;
 count_incr = incrSTRef countref;
 } in (count_incr 1)>>=(\((_k_pos)) -> let {k = exp (_v@>_k_pos);
 } in (((post_incr ((\((x::Double)) -> ((((1-1)*(log x))+((0.000-x)/1))-(1*(log 1)))-(gammaln 1)) k))>>(grad_incr _k_pos (((\((_x)) -> (exp _x)*1) (_v@>_k_pos))*((((0+((1-1)*((1/k)*1)))+((0+(1*(0-1)))/(1*1)))-0)-0))))>>(return ()))>>((count_incr 1)>>=(\((_th_pos)) -> let {th = exp (_v@>_th_pos);
 } in (((post_incr ((\((x::Double)) -> ((((1-1)*(log x))+((0.000-x)/0.100))-(1*(log 0.100)))-(gammaln 1)) th))>>(grad_incr _th_pos (((\((_x)) -> (exp _x)*1) (_v@>_th_pos))*((((0+((1-1)*((1/th)*1)))+((0+(0.100*(0-1)))/(0.100*0.100)))-0)-0))))>>(return ()))>>((count_incr 1)>>=(\((_eta_pos)) -> let {eta = exp (_v@>_eta_pos);
 } in (((post_incr ((\((x::Double)) -> ((((1-1)*(log x))+((0.000-x)/0.100))-(1*(log 0.100)))-(gammaln 1)) eta))>>(grad_incr _eta_pos (((\((_x)) -> (exp _x)*1) (_v@>_eta_pos))*((((0+((1-1)*((1/eta)*1)))+((0+(0.100*(0-1)))/(0.100*0.100)))-0)-0))))>>(return ()))>>((count_incr 1)>>=(\((_mu_pos)) -> let {mu = _v@>_mu_pos;
 } in (((post_incr ((\((x::Double)) -> ((log 1)-(0.500*(log ((2.000*pi)*0.100))))-(((x-0)**2)/(2*0.100))) mu))>>(grad_incr _mu_pos (0-((0+((2*0.100)*((2*(mu-0))*(1-0))))/((2*0.100)*(2*0.100))))))>>(return ()))>>((count_incr 1)>>=(\((_v_0_pos)) -> let {v_0 = exp (_v@>_v_0_pos);
 } in (((post_incr ((\((x::Double)) -> ((((1-1)*(log x))+((0.000-x)/0.100))-(1*(log 0.100)))-(gammaln 1)) v_0))>>(grad_incr _v_0_pos (((\((_x)) -> (exp _x)*1) (_v@>_v_0_pos))*((((0+((1-1)*((1/v_0)*1)))+((0+(0.100*(0-1)))/(0.100*0.100)))-0)-0))))>>(return ()))>>(((post_incr ((\(_) -> (log 1)-(log (2.000-0.000))) (s 0.000)))>>(return ()))>>((count_incr ((round (tmax/dt))-1))>>=(\(v__startpos) -> let {v__obsSig = vcons v_0 (vmap (\((_x)) -> exp _x) ((slice v__startpos ((round (tmax/dt))-1)) _v));
 v = (pack dt 0) v__obsSig;
 } in (forM (fromTo 1 ((round (tmax/dt))-1)) (\((_timeix)) -> let {t = (unround _timeix)*dt;
 t__last = (unround (_timeix-1))*dt;
 v__last = v t__last;
 s__last = s t__last;
 } in ((return ())>>(let {v__next = v t;
 v__mean = v__last+(dt*(k*(th-v__last)));
 v__var = dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)));
 } in (post_incr (((log 1)-(0.500*(log ((2*pi)*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last))))))))-(((v__next-(v__last+(dt*(k*(th-v__last)))))**2)/(2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last))))))))>>((((((return ())>>(let {(_v__next_pos) = v__startpos+(_timeix-1);
 } in grad_incr _v__next_pos (((\((_x)) -> (exp _x)*1) (_v@>_v__next_pos))*(0-((0+((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*((2*(v__next-(v__last+(dt*(k*(th-v__last))))))*(1-0))))/((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*(2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))))))))>>(grad_incr _eta_pos (((\((_x)) -> (exp _x)*1) (_v@>_eta_pos))*((0-(0+(0.500*((1/((2*pi)*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last))))))*(0+((2*pi)*(0+(dt*((0+((eta*(sqrt v__last))*(0+((sqrt v__last)*1))))+((eta*(sqrt v__last))*(0+((sqrt v__last)*1))))))))))))-((0+(0-(((v__next-(v__last+(dt*(k*(th-v__last)))))**2)*(0+(2*(0+(dt*((0+((eta*(sqrt v__last))*(0+((sqrt v__last)*1))))+((eta*(sqrt v__last))*(0+((sqrt v__last)*1)))))))))))/((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*(2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))))))))>>(grad_incr _th_pos (((\((_x)) -> (exp _x)*1) (_v@>_th_pos))*(0-((0+((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*((2*(v__next-(v__last+(dt*(k*(th-v__last))))))*(0-(0+(0+(dt*(0+(k*(1-0))))))))))/((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*(2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))))))))>>(grad_incr _k_pos (((\((_x)) -> (exp _x)*1) (_v@>_k_pos))*(0-((0+((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*((2*(v__next-(v__last+(dt*(k*(th-v__last))))))*(0-(0+(0+(dt*(0+((th-v__last)*1)))))))))/((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*(2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))))))))>>(let {(_v__last_pos) = if (_timeix==1) then _v_0_pos else (v__startpos+(_timeix-2));
 } in grad_incr _v__last_pos (((\((_x)) -> (exp _x)*1) (_v@>_v__last_pos))*((0-(0+(0.500*((1/((2*pi)*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last))))))*(0+((2*pi)*(0+(dt*((0+((eta*(sqrt v__last))*(0+(eta*((0.500*(1/(sqrt v__last)))*1)))))+((eta*(sqrt v__last))*(0+(eta*((0.500*(1/(sqrt v__last)))*1)))))))))))))-(((0+((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*((2*(v__next-(v__last+(dt*(k*(th-v__last))))))*(0-(1+(0+(dt*(0+(k*(0-1))))))))))+(0-(((v__next-(v__last+(dt*(k*(th-v__last)))))**2)*(0+(2*(0+(dt*((0+((eta*(sqrt v__last))*(0+(eta*((0.500*(1/(sqrt v__last)))*1)))))+((eta*(sqrt v__last))*(0+(eta*((0.500*(1/(sqrt v__last)))*1))))))))))))/((2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))*(2*(dt*((eta*(sqrt v__last))*(eta*(sqrt v__last)))))))))))))>>(let {s__next = s t;
 s__mean = s__last+(dt*(mu*s__last));
 s__var = dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last));
 } in (post_incr (((log 1)-(0.500*(log ((2*pi)*s__var))))-(((s__next-s__mean)**2)/(2*s__var))))>>(((return ())>>(grad_incr _mu_pos (1*(0-((0+((2*(dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last))))*((2*(s__next-(s__last+(dt*(mu*s__last)))))*(0-(0+(0+(dt*(0+(s__last*1)))))))))/((2*(dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last))))*(2*(dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last))))))))))>>(let {(_v__last_pos) = if (_timeix==1) then _v_0_pos else (v__startpos+(_timeix-2));
 } in grad_incr _v__last_pos (((\((_x)) -> (exp _x)*1) (_v@>_v__last_pos))*((0-(0+(0.500*((1/((2*pi)*(dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last)))))*(0+((2*pi)*(0+(dt*((0+(((sqrt v__last)*s__last)*(0+(s__last*((0.500*(1/(sqrt v__last)))*1)))))+(((sqrt v__last)*s__last)*(0+(s__last*((0.500*(1/(sqrt v__last)))*1)))))))))))))-((0+(0-(((s__next-(s__last+(dt*(mu*s__last))))**2)*(0+(2*(0+(dt*((0+(((sqrt v__last)*s__last)*(0+(s__last*((0.500*(1/(sqrt v__last)))*1)))))+(((sqrt v__last)*s__last)*(0+(s__last*((0.500*(1/(sqrt v__last)))*1))))))))))))/((2*(dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last))))*(2*(dt*(((sqrt v__last)*s__last)*((sqrt v__last)*s__last)))))))))))))>>((readSTRef postref)>>=(\(postval) -> (unsafeFreeze gradref)>>=(\(gradval) -> return ((postval,gradval)))))))))))))))))))));
 vToRec = \(theta) -> runP (toList theta) (bindP (fmapP (\((_x)) -> exp _x) headP) (\(k) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(th) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(eta) -> bindP headP (\(mu) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(v_0) -> bindP (fmapP (\(s) -> v_0:(map (\(_x ) -> exp _x) s)) (takeP$((round (tmax/dt))-1))) (\(v__obsSig) -> bindP (returnP ((packL dt 0) v__obsSig)) (\(v) -> returnP (R.X R.:& V R.:= (v) R.:& V_0 R.:= (v_0) R.:& Mu R.:= (mu) R.:& Eta R.:= (eta) R.:& Th R.:= (th) R.:& K R.:= (k))))))))));
 inisam = (gamma 1 1)>>=(\(k) -> (gamma 1 0.100)>>=(\(th) -> (gamma 1 0.100)>>=(\(eta) -> (normal 0 0.100)>>=(\(mu) -> (gamma 1 0.100)>>=(\(v_0) -> (uniform 0.000 2.000)>>=(\(s_0) -> (wiener dt tmax)>>=(\(w2) -> (wiener dt tmax)>>=(\(w1) -> let {v = solveODE (\v-> \((t::Double)) -> (k*(th-v))+((eta*(sqrt v))*((d dt w1) t))) tmax dt v_0;
s = solveODE (\s-> \((t::Double)) -> (mu*s)+(((sqrt (v t))*s)*((d dt w2) t))) tmax dt s_0;
} in return (concat (((log k):[]):(((log th):[]):(((log eta):[]):((mu:[]):(((log v_0):[]):((map (\((_x)) -> log _x) (sigTail v)):[])))))))))))))));

 } in R.X R.:& Inisam R.:= (inisam) R.:& VToRec R.:= (vToRec) R.:& Postgrad R.:= (postgrad) R.:& Posterior R.:= (posterior)
  let post :: ((Vector Double) -> Double)
      post = (prims!!!Posterior) . VS.toList
  let postgrad :: ((Vector Double) -> ((Double,(Vector Double))))
      postgrad = prims!!!Postgrad
  let vtorec :: ((Vector Double) -> ((R.X R.:& V R.::: (Double -> Double) R.:& V_0 R.::: Double R.:& Mu R.::: Double R.:& Eta R.::: Double R.:& Th R.::: Double R.:& K R.::: Double) (Id KindStar)))
      vtorec = prims!!!VToRec
  let inisam :: Prob (Vector Double)
      inisam = fmap (VS.fromList) $ prims!!!Inisam
  return (post, postgrad,vtorec,inisam)

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
import System.Exit

import Quipper
import Quipper.Internal

import Quantum.Synthesis.Matrix
import Quantum.Synthesis.Ring
import Quipper.Libraries.Synthesis

import Quipper.Libraries.Decompose.GateBase
import Quipper.Libraries.Decompose

import Data.Complex
import Data.Ratio
import Data.Tuple

import Quipper.Utils.RandomSource
import System.Random

-- declare sample_oracle's data type
data Oracle = Oracle {
   qubit_num :: Int,
   function :: ([Qubit], Qubit) -> Circ ([Qubit], Qubit)
}

-- declare circuit function
circuit_function :: Oracle -> Circ ([Bit], Bit)
circuit_function oracle = do
     -- initialize string of qubits
     top_qubits <- qinit (replicate (qubit_num oracle) False)
     bottom_qubit <- qinit False
     label (top_qubits, bottom_qubit) ("|0>","|0>")

     --prepare the initial states
     --mapUnary hadamard top_qubits
     --mapUnary hadamard bottom_qubit

     comment "before oracle"
     -- call oracle
     function oracle (top_qubits, bottom_qubit)
     comment "after oracle"

     -- measure qubits
     (top_qubits, bottom_qubit) <- measure (top_qubits, bottom_qubit)
     -- discard unnecessary output and return result
     return (top_qubits, bottom_qubit)

-- * From a matrix
type Sixteen = Ten_and Six

mymatrix :: Matrix Sixteen Sixteen (Integer)
mymatrix = matrix [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([],0)
                   [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([],1)
                   [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([0],0)
                   [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([0],1)
                   [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([0,0],0)
                   [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([0,0],1)
                   [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], -- ([1,0],0)
                   [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], -- ([1,0],1)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], -- ([1],0)
                   [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], -- ([1],1)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], -- ([0,1],0)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], -- ([0,1],1)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], -- ([1,1],0)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], -- ([1,1],1)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], -- ([0,0,0],0)
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]] -- ([0,0,0],1)

synthesized = exact_synthesis mymatrix

circuit ::([Qubit], Qubit)-> Circ ([Qubit], Qubit)
circuit ([a, b, c], d) = do
  synthesized [a,b,c,d]
  return ([a,b,c],d)

-- Standard decompose in X, Y , Z, H , S, S*, T, T* and CNOT
-- needs : precision and randomSource
rand = RandomSource(fst(split(mkStdGen 10)))
prec = 20 * bits
circuit_decompose = decompose_generic (Standard prec rand) circuit

--
original_pdf :: IO ()
original_pdf = do
  print_generic Preview (circuit_function my_oracle)
    where
    -- declare empty_oracle's data type
    my_oracle :: Oracle
    my_oracle = Oracle {
    -- set the length of qubit string
        qubit_num = 3,
        function = circuit
    }

--
original_ASCII :: IO ()
original_ASCII = do
  print_generic ASCII (circuit_function my_oracle)
    where
     -- declare empty_oracle's data type
     my_oracle :: Oracle
     my_oracle = Oracle {
      -- set the length of qubit string
       qubit_num = 3,
       function = circuit
     }

decompose_pdf :: IO ()
decompose_pdf = do
  print_generic Preview (circuit_function my_oracle)
    where
   -- declare empty_oracle's data type
   my_oracle :: Oracle
   my_oracle = Oracle {
     -- set the length of qubit string
      qubit_num = 3,
      function = circuit_decompose
   }
 -- print mymatrix

decompose_ASCII :: IO ()
decompose_ASCII = do
  print_generic ASCII (circuit_function my_oracle)
    where
   -- declare empty_oracle's data type
   my_oracle :: Oracle
   my_oracle = Oracle {
     -- set the length of qubit string
      qubit_num = 3,
      function = circuit_decompose
   }

main_menu = do
  putStrLn "\n \n \n \n \n \n \n \n choose an option:"
  putStrLn " 1 - PDF original circuit\n"
  putStrLn " 2 - text description of original circuit\n"
  putStrLn " 3 - PDF decomposed circuit\n"
  putStrLn " 4 - text description of decomposed circuit\n"
  putStrLn " 5 - exit\n"
  line <- getLine
  case line of
    "1" -> do original_pdf
              main
    "2" -> do original_ASCII
              main
    "3" -> do decompose_pdf
              main
    "4" -> do decompose_ASCII
              main
    "5" -> do exitSuccess
    _ -> do main

main_qfold_cnot_7x2 = do
  decompose_ASCII

main = do
  --main_menu

  -- or
  -- if you want the circuit description a .txt file
  -- comment the main_menu above
  -- replacing it by:

  main_qfold_cnot_7x2

  -- and run:
  -- $ ./qfold_cnot_7x2_quipper > circuit_cnot_7x2_quipper.txt

ext: hs

max_stars_repo_name: AnaNeri/quantamorphismsGuide

max_issues_repo_name: AnaNeri/quantamorphismsGuide

max_forks_repo_name: AnaNeri/quantamorphismsGuide


text:
{-# LANGUAGE FlexibleContexts #-}
module Legacy.QuantumCircuit where

import Data.Complex
import Data.Either
import Data.List
import Data.Function
import System.Random
import Control.Monad.IO.Class

import Helper

import Tensor

-- applyGate :: (QuantumRegister q Num RealFloat, RealFloat a) => q (Complex a) -- ^ Gate tensor
--   -> Int -- ^ Index
--   -> q (Complex a) -- ^ State vector
--   -> Either String (q (Complex a))
-- applyGate x i y
--     | not (isSquare x) = Left "Invalid Gate Tensor"
--     | not (isSquare y) = Left "Invalid Statevector"
--     | i < 0 || i >= rank y = Left "Index out of bounds"
--     | otherwise = Right $ norm $ build (rank y) [(i, x)] pauliId * y

-- applyGateAll :: (QuantumRegister q Num RealFloat, RealFloat a) => q (Complex a) -> q (Complex a) -> Either String (q (Complex a))
-- applyGateAll x y
--     | not (isSquare x) = Left "Invalid Gate Tensor"
--     | not (isSquare y) = Left "Invalid Statevector"
--     | otherwise = Right $ norm $ build (rank y) [] x * y

-- applyControl :: (QuantumRegister q, RealFloat a) => q (Complex a) -- ^ Gate tensor
--   -> Int -- ^ Control qubit index
--   -> Int -- ^ Apply qubit index
--   -> q (Complex a) -- ^ State vector
--   -> Either String (q (Complex a))
-- applyControl x ctl i y
--     | not (isSquare x) = Left "Invalid gate tensor"
--     | not (isSquare y) = Left "Invalid state vector"
--     | i < 0 || i >= rank y = Left "Application qubit index out of bounds"
--     | ctl < 0 || ctl >= rank y = Left "Control qubit index out of bounds"
--     | ctl == i = Left "Control and application qubit cannot be the same"
--     | otherwise = Right $ norm $ (build (rank y) [(ctl, mask0)] pauliId
--         + build (rank y) (sortBy (compare `on` fst) [(ctl, mask1), (i, x)]) pauliId) * y

-- collapse :: (RealFloat a, QuantumRegister q) => Int -- ^ Index of the qubit
--         -> Bool -- ^ Collapsed State
--         -> q (Complex a) -- ^ Qubits to collapse
--         -> Either String (q (Complex a)) -- ^ Resulting Tensor
-- collapse i state q
--     | i < 0 || i >= rank q = Left "Index out of bounds"
--     | not (isSquare q) = Left "Invalid state vector"
--     | state = norm <$> applyGate mask1 i q
--     | otherwise = norm <$> applyGate mask0 i q

-- measure i v = do
--     a <- randomRIO (0.0, 1.0)
--     if a < fst p then
--         return $ fmap (\x -> (False, x)) (collapse i False =<< v)
--     else
--         return $ fmap (\x -> (True, x)) (collapse i True =<< v)
--     where
--         prob = getProb i =<< v
--         p = fromRight (0.0, 0.0) prob


text:
{-|
 Module      : Data.Boltzmann.System.Tuner
 Description : Interface utilities with the Paganini tuner.
 Copyright   : (c) Maciej Bendkowski, 2017-2020

 License     : BSD3
 Maintainer  : maciej.bendkowski@tcs.uj.edu.pl
 Stability   : experimental

 General utilities managing the IO interface between Boltzmann Brain
 and the Paganini tuner script.
 -}
module Data.Boltzmann.System.Tuner
  ( readPaganini
  , runPaganini
  )
where

import           Data.Void
import qualified Data.Map.Strict               as M

import           Numeric.LinearAlgebra   hiding ( size
                                                , double
                                                )

import           Text.Megaparsec

import           Data.Boltzmann.Internal.Logging
import           Data.Boltzmann.Internal.Parser
import           Data.Boltzmann.Internal.Tuner
import           Data.Boltzmann.System
import           Data.Boltzmann.System.Paganini
import qualified Data.Boltzmann.System.Tuner.Algebraic
                                               as A
import           Data.Boltzmann.Internal.Annotations
                                                ( withString )

type Parser = Parsec Void String

-- | Communicates with Paganini and collects the respective
--   tuning vector for the given system. If communication is not possible,
--   for instance due to the missing Paganini script, the current process
--   is terminated with an error message on the standard error output.
runPaganini
  :: System Int
  -> Parametrisation
  -> IO (Either (ParseErrorBundle String Void) (PSystem Double))

runPaganini sys paramT = do

  -- FIXME use arg
  info "Running paganini..."
  let gen = withString (annotations sys) "generate" (initType sys)
  result <- toPaganini gen sys
  case result of
    Left  err           -> fail' (show err)
    Right (rho, ts, us) -> do
      info "Computing branching probabilities..."
      let ts'  = fromList ts
      let sys' = parametrise sys paramT rho ts' us
      return $ Right sys'

-- | Parses the given input string as a Paganini tuning vector.
readPaganini
  :: System Int
  -> Parametrisation
  -> String
  -> IO (Either (ParseErrorBundle String Void) (PSystem Double))

readPaganini sys paramT f = do
  let spec = A.toPSpec sys
  pag <- parsePaganini spec f
  case pag of
    Left  err           -> return $ Left err
    Right (rho, us, ts) -> do
      let ts' = fromList ts
      return (Right $ parametrise sys paramT rho ts' us)

paganiniStmt :: PSpec -> Parser (Double, [Double], [Double])
paganiniStmt spec = do
  rho <- double
  us  <- parseN double $ numFreqs spec
  ts  <- parseN double $ numTypes spec
  return (rho, us, ts)

-- | Parses the given Paganini specification.
parsePaganini
  :: PSpec
  -> String
  -> IO
       ( Either
           (ParseErrorBundle String Void)
           (Double, [Double], [Double])
       )

parsePaganini spec = parseFromFile (paganiniStmt spec)

accumulateCons :: Double -> [Cons Double] -> [Cons Double]
accumulateCons _   []         = []
accumulateCons acc (con : xs) = con' : accumulateCons acc' xs
 where
  con' = con { weight = w + acc }
  acc' = acc + w
  w    = weight con

accumulateType :: (String, [Cons Double]) -> (String, [Cons Double])
accumulateType (t, cons) = (t, accumulateCons 0 cons)

accumulate :: PSystem Double -> PSystem Double
accumulate psys = psys { system = sys { defs = M.fromList types' } }
 where
  sys    = system psys
  typs   = M.toList (defs sys)
  types' = map accumulateType typs

-- | Compute the numerical branching probabilities for the given system.
parametrise
  :: System Int
  -> Parametrisation
  -> Double
  -> Vector Double
  -> [Double]
  -> PSystem Double

parametrise sys paramT rho ts us =
  let sys' = A.paramSystem sys rho ts us
  in  case paramT of
        Regular     -> sys'
        Cummulative -> accumulate sys'

max_stars_repo_name: timsears/boltzmann-brain

max_issues_repo_name: timsears/boltzmann-brain

max_forks_repo_name: timsears/boltzmann-brain


text:
-- Rummikub.hs
-- Alex Striff

module Rummikub where

import Prelude hiding
  (Enum, succ, pred, toEnum, fromEnum, enumFrom, enumFromThen, enumFromTo, enumFromThenTo)
import Prelude.SafeEnum
import Control.Applicative
import Control.Arrow
import Control.Monad
import Data.Function
import Data.Functor
import Data.Maybe
import Data.Function.Pointless
import Data.Traversable
import Data.Graph
import Data.List
import Data.List.Split
import Text.Read
import Numeric.LinearAlgebra hiding ((<#), (<>))

newtype TileNum = TileNum Int deriving (Eq, Ord)

instance Show TileNum where
  show (TileNum n) = show n

instance Bounded TileNum where
  minBound = TileNum 1
  maxBound = TileNum 12

instance DownwardEnum TileNum where
  pred t@(TileNum n)
    | minBound < t && t <= maxBound = TileNum <$> pred n
    | otherwise = Nothing
  precedes (TileNum n) (TileNum m) = minBound <= n && n < m && m <= maxBound

instance UpwardEnum TileNum where
  succ t@(TileNum n)
    | minBound <= t && t < maxBound = TileNum <$> succ n
    | otherwise = Nothing
  succeeds (TileNum n) (TileNum m) = maxBound >= n && n > m && m >= minBound

instance Enum TileNum where
  toEnum n
    | minBound <= tn && tn <= maxBound = Just tn
    | otherwise = Nothing
    where tn = TileNum n
  fromEnum (TileNum n) = Just n

tilenum :: Int -> Maybe TileNum
tilenum = toEnum

data Color = Black | Red | Orange | Blue | AnyColor

instance Show Color where
  show Black    = "B"
  show Red      = "r"
  show Orange   = "o"
  show Blue     = "b"
  show AnyColor = "*"

instance Eq Color where
  Black    /= Black    = False
  Red      /= Red      = False
  Orange   /= Orange   = False
  Blue     /= Blue     = False
  _ /= _ = True

instance Ord Color where
  compare Black    Red      = LT
  compare Red      Orange   = LT
  compare Orange   Blue     = LT
  compare Blue     Black    = LT
  compare Black    Black    = EQ
  compare Red      Red      = EQ
  compare Orange   Orange   = EQ
  compare Blue     Blue     = EQ
  compare AnyColor _        = EQ
  compare _        AnyColor = EQ
  compare _        _        = GT

color :: Char -> Maybe Color
color s
  | s == 'B'  = Just Black
  | s == 'r'  = Just Red
  | s == 'o'  = Just Orange
  | s == 'b'  = Just Blue
  | s == '*'  = Just AnyColor
  | otherwise = Nothing

data Tile = NumTile TileNum Color | Joker deriving (Ord)

instance Show Tile where
  show (NumTile n c) = show c ++ show n
  show Joker = "j"

instance Eq Tile where
  (NumTile n c) == (NumTile m d) = n == m && c == d
  _ == _ = True -- Comparision with Joker

(<#)   :: Tile -> Tile -> Bool
(NumTile n c) <# (NumTile m d) = succ n == pure m && c == d
Joker         <# (NumTile m _) = minBound < m
(NumTile n _) <# Joker         = n < maxBound
Joker         <# Joker         = True

(<#<)  :: Tile -> Tile -> Bool
(NumTile n c) <#< (NumTile m d) = (succ <=< succ) n == pure m && c == d
Joker         <#< (NumTile m _) = pure minBound < pred m
(NumTile n _) <#< Joker         = succ n < pure maxBound
Joker         <#< Joker         = True

(<##<)  :: Tile -> Tile -> Bool
(NumTile n c) <##< (NumTile m d) = (succ <=< succ <=< succ) n == pure m && c == d
Joker         <##< (NumTile m _) = pure minBound < (pred <=< pred) m
(NumTile n _) <##< Joker         = (succ <=< succ) n < pure maxBound
Joker         <##< Joker         = True

(<@)   :: Tile -> Tile -> Bool
(NumTile n c) <@ (NumTile m d) = c /= d && n == m
_ <@ _ = True -- Comparison with Joker

(<@<)  :: Tile -> Tile -> Bool
(<@<) = (<@)

tileNum :: Tile -> Maybe TileNum
tileNum (NumTile n _) = pure n
tileNum Joker = Nothing

tileColor (NumTile _ c) = c
tileColor _ = AnyColor

tile :: String -> Maybe Tile
tile [] = Nothing
tile s@(c:n)
  | s == "j"  = Just Joker
  | otherwise = liftA2 NumTile (tilenum =<< readMaybe n) (color c)

tiles :: String -> Maybe [Tile]
tiles = traverse tile . splitOn " "

allRel :: (b -> [(a, a)]) -> (a -> a -> Bool) -> b -> Bool
allRel p r = all (uncurry r) . p

dpairs :: [a] -> [(a, a)]
dpairs = sequence <=< ap zip (tail <$> tails)

neighbors :: [a] -> [(a, a)]
neighbors = ap zip tail

acquaintances :: [a] -> [(a, a)]
acquaintances [_] = []
acquaintances xs  = ap zip (tail . tail) xs

randos :: [a] -> [(a, a)]
randos xs
  | length xs < 4 = []
  | otherwise     = ap zip (tail . tail . tail) xs

validNumRun :: [Tile] -> Bool
validNumRun = (. (&)) . flip all $
  [ (3 <=) . length
  , allRel neighbors (<#)
  , allRel acquaintances (<#<)
  , allRel randos (<##<)]

validColorRun :: [Tile] -> Bool
validColorRun = (. (&)) . flip all $
  [ liftA2 (||) (3 ==) (4 ==) . length
  , allRel dpairs ((/=) `on` tileColor)
  , allRel neighbors (==) . mapMaybe tileNum]

validRun :: [Tile] -> Bool
validRun = liftA2 (||) validNumRun validColorRun

type Run    = [Tile]
type Board  = [Run]
type Frag   = Run
type PBoard = ([Frag], Board)

showBoard :: Board -> String
showBoard = join . intersperse "\n" . fmap (join . intersperse " " . fmap show)

rightFriends :: Tile -> Tile -> Bool
rightFriends = (liftA2 . liftA2) (||) (<#) (<@)

leftFriends :: Tile -> Tile -> Bool
leftFriends = flip rightFriends

splitBin :: (a -> b -> Bool) -> a -> [b] -> [([b], [b])]
splitBin = (.) $ liftA2 fmap (flip splitAt) . findIndices

leftFriendSplits :: Tile -> Frag -> [(Frag, Frag)]
leftFriendSplits  = splitBin leftFriends
rightFriendSplits = splitBin rightFriends

friendSplits :: Frag -> Run -> ([(Frag, Frag)], [(Frag, Frag)])
friendSplits = liftA2 (&&&) (leftFriendSplits . head) (rightFriendSplits . last)

leftSplitRuns :: Frag -> (Frag, Frag) -> [Frag]
leftSplitRuns fs (xs, (f:ys))
  | null ys   = [left]
  | otherwise = [left, ys]
  where left  = xs ++ [f] ++ fs

rightSplitRuns :: Frag -> (Frag, Frag) -> [Frag]
rightSplitRuns fs (xs, xs')
  | null xs   = [right]
  | otherwise = [xs, right]
  where right = fs ++ xs'

fragPlays' :: Run -> Frag -> [[Frag]]
fragPlays' ts fs = (leftSplitRuns fs <$> ls) ++ (rightSplitRuns fs <$> rs)
  where (ls, rs) = friendSplits fs ts

fragPlays :: Run -> Frag -> [PBoard]
fragPlays = fmap sortRuns .: fragPlays'

frags :: PBoard -> [Frag]
frags = fst

runs :: PBoard -> Board
runs = snd

pboard :: [Frag] -> Board -> PBoard
pboard = (,)

parsePBoard :: [String] -> [String] -> Maybe PBoard
parsePBoard = liftA2 pboard `on` traverse tiles

uniq :: Eq a => [a] -> [a]
uniq = nub -- lame O(n^2)

altTopBot :: Int -> [Int] -- altTopBot 4 == [1, 5, 2, 4, 3]
altTopBot n
  | n > 0 = take n . join . ap (zipWith $ \a b -> [a, b]) reverse $ [1..n]
  | otherwise = []

sumSet :: Int -> [[Int]]
sumSet n = sumSet' n [[]] where
  sumSet' n ss
    | n > 0 = [1..n] >>= \a -> sumSet' (n - a) ((a:) <$> ss)
    | otherwise = ss

fragments :: [a] -> [[[a]]]
fragments xs = reverse $ flip splitPlaces xs <$> sumSet (length xs)

sortRuns :: [Frag] -> PBoard
sortRuns = foldl (\(a, b) r -> if validRun r then (a, r:b) else (r:a, b)) mempty

addSplitRuns :: PBoard -> PBoard -> PBoard
addSplitRuns = mappend

playFrags :: PBoard -> PBoard
playFrags p = addSplitRuns (pboard mempty . runs $ p) . sortRuns . frags $ p

select :: Eq b => (a -> [b]) -> a -> [(b, [b])]
select = (ap (zipWith $ liftA2 (.) (,) delete) repeat .)

dropUnfinished :: [PBoard] -> [PBoard]
dropUnfinished = filter $ null . frags

liftFrag :: Functor f => ([Frag] -> f [Frag]) -> PBoard -> f PBoard
liftFrag f = uncurry (<&>) . (f *** flip pboard)

fragment :: PBoard -> [(Frag, [Frag])]
fragment p = do
  (frag, fs)   <- select frags p -- select a fragment
  ffs          <- fragments frag -- smash it
  (frag', fs') <- select id ffs  -- select one of the pieces
  return $ (frag', fs' ++ fs)

stepPlay :: ([PBoard] -> [PBoard]) -> PBoard -> [PBoard]
stepPlay filt p
  | null $ frags p = return p
  | otherwise = do
  (frag, fs)   <- fragment p     -- select a fragment (by smashing)
  (run,  rs)   <- select runs p  -- select a run
  -- fp           <- filt . uniq $ fragPlays run frag -- play the fragment in the run
  fp           <- filt $ fragPlays run frag -- play the fragment in the run
  return $ addSplitRuns fp $ pboard fs rs

stepPlay' :: Int -> ([PBoard] -> [PBoard]) -> PBoard -> [PBoard]
stepPlay' n = foldr1 (<=<) . take n . repeat . stepPlay

fuse :: Eq a => [a] -> [a] -> [[a]]
fuse x@(a:as) y@(b:bs)
  | last as == b = [x ++ bs]
  | a == last bs = [y ++ as]
  | otherwise    = [x, y]

-- TODO: solve problem of turning
-- [[b1], [b2], [b3]] into [[b1 b2 b3]] without duplicates or using uniq/nub

-- fuseFrags' :: [Frag] -> [[Frag]]
-- fuseFrags' fl = do
--   (f, f')   <- dpairs fl
--   return $ fuse f f' ++ (fl \\ [f, f'])

groupFrags' :: [Frag] -> [[Frag]]
groupFrags' fl = do
  (f, f')   <- dpairs fl
  fg        <- fragPlays' f f'
  return $ fg ++ (fl \\ [f, f'])

groupFrags :: PBoard -> [PBoard]
groupFrags = liftFrag groupFrags'

-- Graph stuff

type NGraph node key = (Graph, Vertex -> (node, key, [key]), key -> Maybe Vertex)

graph :: Ord key => NGraph node key -> Graph
graph (g, _, _) = g

nodeFromVertex :: Ord key => NGraph node key -> (Vertex -> (node, key, [key]))
nodeFromVertex (_, nfv, _) = nfv

vertexFromKey :: Ord key => NGraph node key -> (key -> Maybe Vertex)
vertexFromKey (_, _, vfk) = vfk

split2 :: Int -> Int -> [a] -> ([a], [a])
split2 i j xs = (ys, zs)
  where (ys, ys') = splitAt i xs
        zs = drop j ys'

constrain1 :: Int -> [Bool] -> (Int, Int) -> [[Bool]]
constrain1 m xs (i, h)
  | m' > l = [[]]
  | otherwise = fmap (\i -> sandwich (split2 i m' xs) ts) [a..b]
  where m' = m - h + 1
        l = length xs
        ts = replicate m' True
        a = max 0 $ i - (m' - 1)
        b = min i $ l - m'
        sandwich (x, y) z = x ++ z ++ y

-- True represents a constrained vertex that will be propagated
-- m is the minimum number of tiles required to form a run
-- TODO: Rewrite to always be minimal
constrain :: Int -> [Bool] -> [[Bool]]
constrain m xs
  | length xs < m  = [[]]
  | not (or xs) = [xs]
  | otherwise = nub
    . join -- . fmap (minimumBy howConstrained)
    . groupBy ((==) `on` parts)
    . foldr1 (liftA2 $ zipWith (||))
    . fmap (constrain1 m xs)
    $ boundaryTrues xs
    where trues = length . filter id
          howConstrained = compare `on` trues
          sigbool b = if b then 1 else -1
          parts bs = ((*) . sigbool . head $ bs)
            . length . filter (uncurry (/=)) . ap zip tail $ bs

boundaryTrues :: [Bool] -> [(Int, Int)]
boundaryTrues bs = reverse . (\(_, _, ts) -> ts)
  . foldl go (head bs, 0, []) $ groups
  where go (b, i, ts) h = (not b, i+h, if b then (i+h-1,h):(i,h):ts else ts)
        groups = length <$> group bs

-- Needs tests to check all cases
-- TODO: add more counting stuff to deal with cases like splitting
-- *oo...oo* (or maybe even *oo..oo*)
growthComponents :: Int -> [Bool] -> [[Bool]]
growthComponents m bs = flip splitPlaces bs . reverse $ maxBound:gs -- maxBound avoids calculating length
  where gc (c, n, i, gs, j) True  = (True, 1, i + 1, maybe gs (:gs) j, Nothing)
        gc (c, n, i, gs, j) False = if n == 2 * (m - 1)
                                    then (False, 1, m, gs, Just $ i - (m - 1))
                                    else (c, n', i + 1, gs, j)
                                    where n' = if c then n + 1 else 0
        (_, _, i, gs, _) = foldl gc (False, 1, 1, [], Nothing) bs

showConstraint :: [Bool] -> String
showConstraint = fmap $ \b -> if b then 'o' else '.'

parseConstraint :: String -> Maybe [Bool]
parseConstraint = traverse parseConstraint'
  where parseConstraint' c
          | c == '.'  = Just False
          | c == 'o'  = Just True
          | otherwise = Nothing

parseConstraint' = fromJust . parseConstraint

mpow :: Matrix Z -> Int -> Matrix Z
mpow = mconcat .: flip replicate

adjacencyMatrix :: (a -> a -> Bool) -> [a] -> Matrix Z
adjacencyMatrix r = ap (join (><) . length) (fmap edge . join (liftA2 r))
  where edge b = (if b then 1 else 0) :: Z

paths :: Int -> Matrix Z -> Int -> Z
paths i = liftA2 ((+) `on` sumElements . (!! i)) toColumns toRows .: mpow

minpaths :: Int -> Int -> Matrix Z -> Z
minpaths m i a = sum $ fmap (paths i a) [m..cols a]

tileGraph :: (Tile -> Tile -> Bool) -> [Tile]
          -> NGraph Tile Int
tileGraph r ts = graphFromEdges
  $ zipWith (\t i -> (t, i, fmap snd . filter fst
    $ zipWith (\t' j -> (r t t', j)) ts [0..])) ts [0..]

-- Here m is the number of tiles
partitionRuns m = join (***) (concat . fmap (foldMap (:[])))
  . partition ((m <=) . length) . components . graph

constraints :: Ord key => Int -> [NGraph node key] -> [[[(Vertex, Bool)]]]
constraints m gs = fmap ((\(a, b) -> (fmap . fmap) (id &&& (`elem` b)) a)
  . (fst *** outcasts . fmap snd)) . select id
  . fmap ((id *** concat) . partRuns . comps) $ gs
  where comps = fmap (foldMap (:[])) . components . graph
        partRuns = partition ((m <=) . length)
        outcasts = foldl1 intersect

max_stars_repo_name: jfjhh/rummikub

max_issues_repo_name: jfjhh/rummikub

max_forks_repo_name: jfjhh/rummikub


text:
{-# LANGUAGE DeriveDataTypeable, FlexibleContexts, Rank2Types #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Packed.Base
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Packed matrices.
--
module Numeric.LinearAlgebra.Packed.Base
    where

import Control.Monad( when )
import Control.Monad.ST( ST, RealWorld, runST, unsafeIOToST )
import Data.Typeable( Typeable )
import Foreign( Storable, Ptr )
import Text.Printf( printf )
import Unsafe.Coerce( unsafeCoerce )

import Numeric.LinearAlgebra.Types( Herm(..) )
import Numeric.LinearAlgebra.Vector( Vector, RVector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V
import Foreign.BLAS( BLAS2 )
import qualified Foreign.BLAS as BLAS

-- | Immutable packed matrices, stored in column-major order.
data Packed e = Packed !Int !(Vector e)
    deriving (Typeable)

-- | Mutable packed matrices in the 'ST' monad.
newtype STPacked s e = STPacked { unSTPacked :: Packed e }
    deriving (Typeable)

-- | Mutable packed matrices in the 'IO' monad.
type IOPacked = STPacked RealWorld

-- | The dimension of the packed matrix.
dim :: (Storable e) => Packed e -> Int
dim (Packed n _) = n
{-# INLINE dim #-}

-- | Allocate a mutable packed matrix of the given dimension.
new_ :: (Storable e) => Int -> ST s (STPacked s e)
new_ n
    | n < 0 = error $
        printf "new_ %d: negative dimension" n
    | otherwise = do
        mx <- V.new_ (n*(n+1) `div` 2)
        x <- V.unsafeFreeze mx
        return $ STPacked $ Packed n x
{-# INLINE new_ #-}

-- | Create a packed matrix view of a vector, ensurint that the
-- vector has dimension @n * (n+1)/2@, where @n@ is the desired dimension.
fromVector :: (Storable e) => Int -> Vector e -> Packed e
fromVector n x
    | not $ 2 * nx == n * (n+1) = error $
        printf ("fromVector %d <vector with dim %d>: dimension mismatch")
               n nx
    | otherwise =
        unsafeFromVector n x
  where
    nx = V.dim x
{-# INLINE fromVector #-}

-- | Create a packed matrix view of a vector, wihtout checking
-- the dimension of the vector.
unsafeFromVector :: (Storable e) => Int -> Vector e -> Packed e
unsafeFromVector = Packed
{-# INLINE unsafeFromVector #-}

-- | Returns the dimension and underlying vector storage of a
-- packed matrix.
toVector :: (Storable e) => Packed e -> (Int, Vector e)
toVector (Packed n v) = (n,v)
{-# INLINE toVector #-}

{-
-- | Create a packed matrix view of a vector, ensurint that the
-- vector has dimension @n * (n+1)/2@, where @n@ is the desired dimension.
fromSTVector :: (Storable e) => Int -> STVector s e -> STPacked s e
fromSTVector n x
    | not $ 2 * nx == n * (n+1) = error $
        printf ("fromVectorST %d <vector with dim %d>: dimension mismatch")
               n nx
    | otherwise =
        STPacked $ unsafeFromVector n x
  where
    nx = V.dim x
{-# INLINE fromSTVector #-}

-- | Create a packed matrix view of a vector, wihtout checking
-- the dimension of the vector.
unsafeFromSTVector :: (Storable e) => Int -> STVector s e -> STPacked s e
unsafeFromSTVector = STPacked
{-# INLINE unsafeFromSTVector #-}

-- | Returns the dimension and underlying vector storage of a
-- packed matrix.
toSTVector :: (Storable e) => STPacked s e -> (Int, STVector s e)
toSTVector (STPacked n v) = (n,v)
{-# INLINE toSTVector #-}
-}

-- | Read-only packed matrices.
class RPacked p where
    -- | Returns the dimension of the packed matrix.
    getDim :: (Storable e) => p e -> ST s Int

    -- | Perform an action with the underlying vector storage of
    -- the packed matrix.
    withVector :: (Storable e)
               => p e
               -> (forall v . RVector v => v e -> ST s a)
               -> ST s a

    -- | Perform an IO action with a pointer to the first element of
    -- the packed matrix.
    unsafeWith :: (Storable e) => p e -> (Ptr e -> IO a) -> IO a

    -- | Converts a read-only packed matrix into an immutable one. This simply
    -- casts the matrix from one type to the other without copying the matrix.
    -- Note that because the matrix is possibly not copied, any subsequent
    -- modifications made to the mutable version of the matrix may be shared
    -- with the immutable version. It is safe to use, therefore, if the mutable
    -- version is never modified after the freeze operation.
    unsafeFreeze :: (Storable e) => p e -> ST s (Packed e)

    unsafeThaw :: (Storable e) => p e -> ST s (STPacked s e)

-- | View a vector as a packed matrix and pass it to a function.
withFromVector :: (RVector v, Storable e)
               => Int
               -> v e
               -> (forall p . RPacked p => p e -> ST s a)
               -> ST s a
withFromVector n v f = do
    iv <- V.unsafeFreeze v
    f $ fromVector n iv
{-# INLINE withFromVector #-}

-- | View a mutable vector as a mutable packed matrix and pass it
-- to a function.
withFromVectorM :: (Storable e)
                => Int
                -> STVector s e
                -> (STPacked s e -> ST s a)
                -> ST s a
withFromVectorM n v f =
    withFromVector n v $ \p -> do
        mp <- unsafeThaw p
        f mp
{-# INLINE withFromVectorM #-}

-- | Perform an action with the underlying vector storage of
-- the mutable packed matrix.  See also 'withVectorView'.
withVectorM :: (Storable e)
            => STPacked s e
            -> (STVector s e -> ST s a)
            -> ST s a
withVectorM mp f =
    withVector mp $ \v -> do
        mv <- V.unsafeThaw v
        f mv
{-# INLINE withVectorM #-}

instance RPacked Packed where
    getDim = return . dim
    {-# INLINE getDim #-}
    withVector (Packed _ v) f = f v
    {-# INLINE withVector #-}
    unsafeWith (Packed _ v) = V.unsafeWith v
    {-# INLINE unsafeWith #-}
    unsafeFreeze = return
    {-# INLINE unsafeFreeze #-}
    unsafeThaw = return . STPacked
    {-# INLINE unsafeThaw #-}

instance RPacked (STPacked s) where
    getDim = getDim . unSTPacked
    {-# INLINE getDim #-}
    withVector = withVector . unSTPacked
    {-# INLINE withVector #-}
    unsafeWith = unsafeWith . unSTPacked
    {-# INLINE unsafeWith #-}
    unsafeFreeze = return . unSTPacked
    {-# INLINE unsafeFreeze #-}
    unsafeThaw p = return $ cast p
      where
        cast :: STPacked s e -> STPacked s' e
        cast = unsafeCoerce
    {-# INLINE unsafeThaw #-}

-- | Create a new copy of a packed matrix.
newCopy :: (RPacked p, Storable e)
        => p e -> ST s (STPacked s e)
newCopy p = do
    n <- getDim p
    p' <- new_ n
    withVector p $ \v ->
        withVectorM p' $ \v' ->
            V.unsafeCopyTo v' v
    return p'
{-# INLINE newCopy #-}

-- | Converts a mutable packed matrix to an immutable one by taking a complete
-- copy of it.
freeze :: (RPacked p, Storable e) => p e -> ST s (Packed e)
freeze p = do
    p' <- newCopy p
    unsafeFreeze p'

-- | A safe way to create and work with a mutable Packed before returning
-- an immutable one for later perusal.
create :: (Storable e)
       => (forall s. ST s ((STPacked s) e))
       -> Packed e
create stmp = runST $ do
    mp <- stmp
    unsafeFreeze mp

-- | A safe way to create and work with a mutable Herm Packed before returning
-- an immutable one for later perusal.
hermCreate :: (Storable e)
           => (forall s. ST s (Herm (STPacked s) e))
           -> Herm Packed e
hermCreate stmh = runST $ do
    (Herm u mp) <- stmh
    p <- unsafeFreeze mp
    return $ Herm u p

-- | @hermRank1Update alpha x a@ returns
-- @alpha * x * x^H + a@.
hermRank1Update :: (BLAS2 e)
                => Double -> Vector e -> Herm Packed e -> Herm Packed e
hermRank1Update alpha x (Herm uplo ap) = hermCreate $ do
    hp' <- Herm uplo `fmap` newCopy ap
    hermRank1UpdateM_ alpha x hp'
    return hp'

-- | @hermRank2Update alpha x y a@ returns
-- @alpha * x * y^H + conj(alpha) * y * x^H + a@.
hermRank2Update :: (BLAS2 e)
                => e -> Vector e -> Vector e -> Herm Packed e
                -> Herm Packed e
hermRank2Update alpha x y (Herm uplo ap) = hermCreate $ do
    hp' <- Herm uplo `fmap` newCopy ap
    hermRank2UpdateM_ alpha x y hp'
    return hp'

-- | @hermRank1UpdateM_ alpha x a@ sets
-- @a := alpha * x * x^H + a@.
hermRank1UpdateM_ :: (RVector v, BLAS2 e)
                  => Double -> v e -> Herm (STPacked s) e -> ST s ()
hermRank1UpdateM_ alpha x (Herm uplo a) = do
    nx <- V.getDim x
    na <- getDim a
    let n = nx

    when ((not . and) [ nx == n, na == n ]) $ error $
        printf ("hermRank1UpdateM_ _ <vector with dim %d>"
                 ++ " (Herm _ <packed matrix with dim %d>):"
                 ++ " invalid dimensions") nx na

    unsafeIOToST $
        V.unsafeWith x $ \px ->
        unsafeWith a $ \pa ->
            BLAS.hpr uplo n alpha px 1 pa

-- | @hermRank2UpdateM_ alpha x y a@ sets
-- @a := alpha * x * y^H + conj(alpha) * y * x^H + a@.
hermRank2UpdateM_ :: (RVector v1, RVector v2, BLAS2 e)
                  => e -> v1 e -> v2 e -> Herm (STPacked s) e -> ST s ()
hermRank2UpdateM_ alpha x y (Herm uplo a) = do
    nx <- V.getDim x
    ny <- V.getDim y
    na <- getDim a
    let n = nx

    when ((not . and) [ nx == n, ny == n, na == n ]) $ error $
        printf ("hermRank2UpdateM_ _ <vector with dim %d>"
                 ++ " <vector with dim %d>"
                 ++ " (Herm _ <packed matrix with dim %d>):"
                 ++ " invalid dimensions") nx ny na

    unsafeIOToST $
        V.unsafeWith x $ \px ->
        V.unsafeWith y $ \py ->
        unsafeWith a $ \pa ->
            BLAS.hpr2 uplo n alpha px 1 py 1 pa

-- | @hermMulVector a x@ returns @a * x@.
hermMulVector :: (BLAS2 e)
                    => Herm Packed e
                    -> Vector e
                    -> Vector e
hermMulVector a x =
    V.create $ do
        y <- V.new_ (V.dim x)
        hermMulVectorTo y a x
        return y

-- | @hermMulVectorWithScale alpha a x@ retunrs @alpha * a * x@.
hermMulVectorWithScale :: (BLAS2 e)
                       => e
                       -> Herm Packed e
                       -> Vector e
                       -> Vector e
hermMulVectorWithScale alpha a x =
    V.create $ do
        y <- V.new_ (V.dim x)
        hermMulVectorWithScaleTo y alpha a x
        return y

-- | @addHermMulVectorWithScales alpha a x y@
-- returns @alpha * a * x + beta * y@.
addHermMulVectorWithScales :: (BLAS2 e)
                           => e
                           -> Herm Packed e
                           -> Vector e
                           -> e
                           -> Vector e
                           -> Vector e
addHermMulVectorWithScales alpha a x beta y =
    V.create $ do
        y' <- V.newCopy y
        addHermMulVectorWithScalesM_ alpha a x beta y'
        return y'

-- | @hermMulVectorTo dst a x@ sets @dst := a * x@.
hermMulVectorTo :: (RPacked p, RVector v, BLAS2 e)
                => STVector s e
                -> Herm p e
                -> v e
                -> ST s ()
hermMulVectorTo dst = hermMulVectorWithScaleTo dst 1

-- | @hermMulVectorWithScaleTo dst alpha a x@
-- sets @dst := alpha * a * x@.
hermMulVectorWithScaleTo :: (RPacked p, RVector v, BLAS2 e)
                         => STVector s e
                         -> e
                         -> Herm p e
                         -> v e
                         -> ST s ()
hermMulVectorWithScaleTo dst alpha a x =
    addHermMulVectorWithScalesM_ alpha a x 0 dst

-- | @addHermMulVectorWithScalesM_ alpha a x beta y@
-- sets @y := alpha * a * x + beta * y@.
addHermMulVectorWithScalesM_ :: (RPacked p, RVector v, BLAS2 e)
                             => e
                             -> Herm p e
                             -> v e
                             -> e
                             -> STVector s e
                             -> ST s ()
addHermMulVectorWithScalesM_ alpha (Herm uplo a) x beta y = do
    na <- getDim a
    nx <- V.getDim x
    ny <- V.getDim y
    let n = ny

    when ((not . and) [ na == n
                      , nx == n
                      , ny == n
                      ]) $ error $
        printf ("addHermMulVectorWithScalesM_ _"
                ++ " (Herm %s <packed matrix with dim %d>)"
                ++ " %s <vector with dim %d>"
                ++ " _"
                ++ " <vector with dim %d>: dimension mismatch")
               (show uplo) na
               nx ny

    unsafeIOToST $
        unsafeWith a $ \pa ->
        V.unsafeWith x $ \px ->
        V.unsafeWith y $ \py ->
            BLAS.hpmv uplo n alpha pa px 1 beta py 1

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.AD.Newton
-- Copyright   :  (c) Edward Kmett 2010
-- License     :  BSD3
-- Maintainer  :  ekmett@gmail.com
-- Stability   :  experimental
-- Portability :  GHC only
--
-----------------------------------------------------------------------------

module Numeric.AD.Newton
    (
    -- * Newton's Method (Forward AD)
      findZero
    , inverse
    , fixedPoint
    , extremum
    -- * Gradient Ascent/Descent (Reverse AD)
    , gradientDescent
    , gradientAscent
    , conjugateGradientDescent
    , conjugateGradientAscent
    ) where

import Prelude hiding (all, mapM, sum)
import Data.Functor
import Data.Foldable (all, sum)
import Data.Traversable
import Numeric.AD.Types
import Numeric.AD.Mode.Forward (diff, diff')
import Numeric.AD.Mode.Reverse (grad, gradWith')
import Numeric.AD.Internal.Combinators
import Numeric.AD.Internal.Composition

-- | The 'findZero' function finds a zero of a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes constant
-- ("it converges"), no further elements are returned.
--
-- Examples:
--
-- >>> take 10 $ findZero (\x->x^2-4) 1
-- [1.0,2.5,2.05,2.000609756097561,2.0000000929222947,2.000000000000002,2.0]
--
-- >>> import Data.Complex
-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)
-- 0.0 :+ 1.0
findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
findZero f = go where
  go x = x : if x == xn then [] else go xn where
    (y,y') = diff' f x
    xn = x - y/y'
{-# INLINE findZero #-}

-- | The 'inverse' function inverts a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes
-- constant ("it converges"), no further elements are returned.
--
-- Example:
--
-- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10)
-- 10.0
inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
inverse f x0 y = findZero (\x -> f x - lift y) x0
{-# INLINE inverse  #-}

-- | The 'fixedPoint' function find a fixedpoint of a scalar
-- function using Newton's method; its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- If the stream becomes constant ("it converges"), no further
-- elements are returned.
--
-- >>> last $ take 10 $ fixedPoint cos 1
-- 0.7390851332151607
fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
fixedPoint f = findZero (\x -> f x - x)
{-# INLINE fixedPoint #-}

-- | The 'extremum' function finds an extremum of a scalar
-- function using Newton's method; produces a stream of increasingly
-- accurate results.  (Modulo the usual caveats.) If the stream
-- becomes constant ("it converges"), no further elements are returned.
--
-- >>> last $ take 10 $ extremum cos 1
-- 0.0
extremum :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
extremum f = findZero (diff (decomposeMode . f . composeMode))
{-# INLINE extremum #-}

-- | The 'gradientDescent' function performs a multivariate
-- optimization, based on the naive-gradient-descent in the file
-- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the
-- VLAD compiler Stalingrad sources.  Its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- It uses reverse mode automatic differentiation to compute the gradient.
gradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
gradientDescent f x0 = go x0 fx0 xgx0 0.1 (0 :: Int)
    where
        (fx0, xgx0) = gradWith' (,) f x0
        go x fx xgx !eta !i
            | eta == 0     = [] -- step size is 0
            | fx1 > fx     = go x fx xgx (eta/2) 0 -- we stepped too far
            | zeroGrad xgx = [] -- gradient is 0
            | otherwise    = x1 : if i == 10
                                  then go x1 fx1 xgx1 (eta*2) 0
                                  else go x1 fx1 xgx1 eta (i+1)
            where
                zeroGrad = all (\(_,g) -> g == 0)
                x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
                (fx1, xgx1) = gradWith' (,) f x1
{-# INLINE gradientDescent #-}

-- | Perform a gradient descent using reverse mode automatic differentiation to compute the gradient.
gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
gradientAscent f = gradientDescent (negate . f)
{-# INLINE gradientAscent #-}

-- | Perform a conjugate gradient descent using reverse mode automatic differentiation to compute the gradient.
conjugateGradientDescent :: (Traversable f, Fractional a, Ord a) =>
                (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
conjugateGradientDescent f x0 = go x0 d0 d0
  where
    dot x y = sum $ zipWithT (*) x y
    d0 = negate <$> grad f x0
    go xi ri di = xi : go xi1 ri1 di1
      where
        ai  = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> lift x + a * lift d) xi di) 0
        xi1 = zipWithT (\x d -> x + ai*d) xi di
        ri1 = negate <$> grad f xi1
        bi1 = max 0 $ dot ri1 (zipWithT (-) ri1 ri) / dot ri1 ri1
        -- bi1 = max 0 $ sum (zipWithT (\a b -> a * (a - b)) ri1 ri) / dot ri1 ri1
        di1 = zipWithT (\r d -> r * bi1*d) ri1 di
{-# INLINE conjugateGradientDescent #-}

-- | Perform a conjugate gradient ascent using reverse mode automatic differentiation to compute the gradient.
conjugateGradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
conjugateGradientAscent f = conjugateGradientDescent (negate . f)
{-# INLINE conjugateGradientAscent #-}

max_stars_repo_name: yairchu/ad

max_issues_repo_name: yairchu/ad

max_forks_repo_name: yairchu/ad


text:
-- | Reanimate Fuckery
module CalamityBot.Commands.Reanimate
    ( reanimateGroup,
    ) where

import Calamity.Commands as C
import Calamity.Commands.Context (FullContext)
import Calamity
import CalamityBot.Utils.Reanimate
import qualified Data.Text as T
import qualified Polysemy as P
import Reanimate
import Reanimate.Render (Format(RenderWebm, RenderGif), Raster(RasterAuto))

import           Data.Complex
import           Graphics.SvgTree
import qualified Graphics.SvgTree.Types
import           Linear.V2
import           Codec.Picture
import qualified Data.ByteString as B
import Network.Mime (defaultMimeLookup)
import Control.Lens
import Network.Wreq

protectText :: T.Text -> T.Text
protectText = T.concatMap (fromString . protectChar)

protectChar :: Char -> String
protectChar '#'  = "\\#"
protectChar '$'  = "\\$"
protectChar '%'  = "\\%"
protectChar '^'  = "\\^{}"
protectChar '&'  = "\\&"
protectChar '{'  = "\\{"
protectChar '}'  = "\\}"
protectChar '~'  = "\\~{}"
protectChar '\\' = "\\textbackslash{}"
protectChar '_'  = "\\_{}"
protectChar x = [x]

renderText :: T.Text -> Animation
renderText text =
  let rt = center . latex $ text
      rtRotated = mkAnimation 5 (\t -> rotateAroundCenter (t * 360) rt)
  in addStatic (mkBackground "white") rtRotated

findSVG :: [Attachment] -> Maybe Attachment
findSVG = listToMaybe . filter (\a -> "image/svg" `B.isPrefixOf` defaultMimeLookup (a ^. #filename))

usableTree :: Tree -> Bool
usableTree None = False
usableTree (UseTree _ _) = False
usableTree (DefinitionTree _) = False
usableTree _ = True

reanimateGroup :: BotC r => P.Sem (DSLState FullContext r) ()
reanimateGroup = void
  . help (const "Commands related to reanimate fuckery")
  . C.group "reanimate"
  $ do
    help (const "Render a message") $
      command @'[KleenePlusConcat T.Text] "render" \ctx msg -> do
        let anim = renderText msg
        r <- P.embed $ renderToMemory anim RasterAuto RenderWebm 400 400 30
        case r of
          Right s -> do
            void $ tell ctx (TFile "lol.webm" s)
          Left e ->
            putLBSLn $ "Failed with reason: " <> e

    help (const "Render a fourier thing of an svg") $
      command @'[] "renders" \ctx -> do
        case findSVG (ctx ^. #message . #attachments) of
          Just svg -> do
            r <- P.embed $ Network.Wreq.get (T.unpack $ svg ^. #url)
            let file = r ^. responseBody
            let Just doc = parseSvgFile "a.svg" (decodeUtf8 file)
            let tree = head . fromList . filter usableTree $ doc ^. Graphics.SvgTree.Types.documentElements
            let f = mkSVGLatex $ flipYAxis tree
            let anim = setDuration 20 $ scene $ do
                  _ <- newSpriteSVG $ mkBackgroundPixel (PixelRGBA8 252 252 252 0xFF)
                  play $ fourierA f (fromToS 0 5)      -- Rotate 15 times
                    & setDuration 30
                    & signalA (reverseS . powerS 2 . reverseS) -- Start fast, end slow
                    & pauseAtEnd 2
                  play $ fourierA f (constantS 0)       -- Don't rotate at all
                    & setDuration 10
                    & reverseA
                    & signalA (powerS 2)                       -- Start slow, end fast
                    & pauseAtEnd 2
            r <- P.embed $ renderToMemory anim RasterAuto RenderGif 480 360 15
            case r of
              Right s -> do
                void $ tell ctx (TFile "lol.gif" s)
              Left e ->
                putLBSLn $ "Failed with reason: " <> e
          Nothing ->
            void $ tell @T.Text ctx "Couldn't find an svg"

    help (const "Render a fourier thing") $
      command @'[KleenePlusConcat T.Text] "renderf" \ctx msg -> do
        let f = mkFourierLatex msg
        let anim = setDuration 20 $ scene $ do
              _ <- newSpriteSVG $ mkBackgroundPixel (PixelRGBA8 252 252 252 0xFF)
              play $ fourierA f (fromToS 0 5)      -- Rotate 15 times
                & setDuration 30
                & signalA (reverseS . powerS 2 . reverseS) -- Start fast, end slow
                & pauseAtEnd 2
              play $ fourierA f (constantS 0)       -- Don't rotate at all
                & setDuration 10
                & reverseA
                & signalA (powerS 2)                       -- Start slow, end fast
                & pauseAtEnd 2
        r <- P.embed $ renderToMemory anim RasterAuto RenderGif 480 360 15
        case r of
          Right s -> do
            void $ tell ctx (TFile "lol.gif" s)
          Left e ->
            putLBSLn $ "Failed with reason: " <> e

scaleToFit :: Double -> Double -> Tree -> Tree
scaleToFit w h t =
  let oh = svgHeight t
      ow = svgWidth t
  in scale (min (w / ow) (h / oh)) t

-- layer 2
fourierA :: Fourier -> (Double -> Double) -> Animation
fourierA f genPhi = animate $ \t ->
    let circles = setFourierLength (t* fourierLen f) f
        coeffs = fourierCoefficients $ rotateFourier (genPhi t) circles
    in mkGroup
    [ drawCircles coeffs
    , withStrokeColor "green" $
      withStrokeLineJoin JoinRound $
      withFillOpacity 0 $
      withStrokeWidth (defaultStrokeWidth*2) $
      mkLinePath $ mkFourierOutline circles
    , let x :+ y = sum coeffs in
      translate x y $ withFillColor "red" $ mkCircle (defaultStrokeWidth*3)
    ]

drawCircles :: [Complex Double] -> SVG
drawCircles [] = mkGroup []
drawCircles ( x :+ y : xs) =
  translate x y $ drawCircles' xs

drawCircles' :: [Complex Double] -> SVG
drawCircles' circles = mkGroup
    [ worker circles
    , withStrokeColor "black" $
      withStrokeLineJoin JoinRound $
      withFillOpacity 0 $
      mkLinePath [ (x, y) | x :+ y <- scanl (+) 0 circles ]]
  where
    worker [] = None
    worker (x :+ y : rest) =
      let radius = sqrt(x*x+y*y) in
      mkGroup
      [ withStrokeColor "dimgrey" $
        withFillOpacity 0 $
        mkCircle radius
      , translate x y $ worker rest ]

-- layer 1
newtype Fourier = Fourier {fourierCoefficients :: [Complex Double]}

mkSVGLatex :: Tree -> Fourier
mkSVGLatex t = mkFourier $ lineToPoints 500 $
  toLineCommands $ extractPath $ center $ scaleToFit 9.6 7.2 $ t

mkFourierLatex :: T.Text -> Fourier
mkFourierLatex t = mkFourier $ lineToPoints 500 $
  toLineCommands $ extractPath $ center $ scaleToFit 9.6 7.2 $ latex t

fourierLen :: Fourier -> Double
fourierLen f = sum $ map magnitude $ drop 1 $ take 500 $ fourierCoefficients f

pointAtFourier :: Fourier -> Complex Double
pointAtFourier = sum . fourierCoefficients

mkFourier :: [RPoint] -> Fourier
mkFourier points = Fourier $ findCoefficient 0 :
    concat [ [findCoefficient n, findCoefficient (-n)] | n <- [1..] ]
  where
    findCoefficient :: Int -> Complex Double
    findCoefficient n =
        sum [ toComplex point * exp (negate (fromIntegral n) * 2 *pi * i*t) * deltaT
            | (idx, point) <- zip [0::Int ..] points, let t = fromIntegral idx/nPoints ]
    i = 0 :+ 1
    toComplex (V2 x y) = x :+ y
    deltaT = recip nPoints
    nPoints = fromIntegral (length points)

setFourierLength :: Double -> Fourier -> Fourier
setFourierLength _ (Fourier []) = Fourier []
setFourierLength len0 (Fourier (first:lst)) = Fourier $ first : worker len0 lst
  where
    worker _len [] = []
    worker len (c:cs) =
      if magnitude c < len
        then c : worker (len - magnitude c) cs
        else [c * realToFrac (len / magnitude c)]

rotateFourier :: Double -> Fourier -> Fourier
rotateFourier phi (Fourier coeffs) =
    Fourier $ worker coeffs (0::Integer)
  where
    worker [] _ = []
    worker (x:rest) 0 = x : worker rest 1
    worker [left] n = worker [left,0] n
    worker (left:right:rest) n =
      let n' = fromIntegral n in
      left * exp (negate n' * 2 * pi * i * phi') :
      right * exp (n' * 2 * pi * i * phi') :
      worker rest (n+1)
    i = 0 :+ 1
    phi' = realToFrac phi

mkFourierOutline :: Fourier -> [(Double, Double)]
mkFourierOutline fourier =
    [ (x, y)
    | idx <- [0 .. granularity]
    , let x :+ y = pointAtFourier $ rotateFourier (idx/granularity) fourier
    ]
  where
    granularity = 500

max_stars_repo_name: simmsb/calamity-bot

max_issues_repo_name: simmsb/calamity-bot

max_forks_repo_name: simmsb/calamity-bot


text:
{-# LANGUAGE CPP #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Instances where

import Control.Monad(replicateM)

#if !MIN_VERSION_base(4,8,0)
import Control.Applicative((<$>), (<*>), pure)
#endif

import Test.Tasty.QuickCheck
import Text.Namelist.Types
import Data.Complex(Complex((:+)))
import Data.CaseInsensitive(mk)
import Data.Char(toUpper)
import Numeric(showFFloat)

arbitrarySafeDouble :: Gen Double
arbitrarySafeDouble = suchThat (arbitrary :: Gen Double) $ \d ->
    let (n, _:f) = break (== '.') $ showFFloat Nothing d ""
    in d == read (n ++ f) / 10 ** fromIntegral (length f)

arbitraryIndex :: Gen Index
arbitraryIndex = do
    i <- choose (1, maxBound)
    return $ Index i

arbitraryRange :: Gen Index
arbitraryRange = do
    a <- oneof [Just <$> choose (1, maxBound - 1),          pure Nothing]
    b <- oneof [Just <$> choose (maybe 1 succ a, maxBound), pure Nothing]
    s <- oneof [Just . unSmall <$> arbitrary,               pure Nothing]
    return $ Range a b s
  where
    unSmall (Small a) = a

instance Arbitrary Index where
    arbitrary = oneof
        [ arbitraryIndex
        , arbitraryRange
        ]

arbitraryName :: Gen String
arbitraryName = sized $ \i -> do
    a  <- elements alpha
    as <- replicateM (min i 30) (elements an_)
    return $ a : as
  where
    an_   = '_': ['0' .. '9'] ++ alpha
    alpha = map toUpper lower ++ lower
    lower = ['a'..'z']

arbitraryKeyName :: Gen Key
arbitraryKeyName = oneof
    [ Key . mk <$> arbitraryName
    , Indexed  <$> (mk <$> arbitraryName) <*> listOf1 arbitrary
    ]

arbitrarySub :: Gen Key
arbitrarySub = sized $ \depth ->
    if depth < 1
    then arbitraryKeyName
    else Sub <$> (mk <$> arbitraryName) <*> arbitrary

instance Arbitrary Key where
    arbitrary = oneof
        [ arbitraryKeyName
        , arbitrarySub
        ]

newtype Scalar = Scalar { unScalar :: Value }

instance Arbitrary Scalar where
    arbitrary = Scalar <$> oneof
        [ Integer <$> arbitrary
        , Real    <$> arbitrarySafeDouble
        , Complex <$> ((:+) <$> arbitrarySafeDouble <*> arbitrarySafeDouble)
        , Logical <$> arbitrary
        , String  <$> arbitrary
        ]

instance Arbitrary Value where
    arbitrary = oneof
        [ unScalar <$> arbitrary
        , pure Null
        , sized $ \i -> do
            as <- map unScalar <$> replicateM (max 1 i) arbitrary
            a  <- unScalar <$> arbitrary
            return $ Array (reverse $ a:as)
        , sized $ \i -> (:*) <$> pure i <*> (unScalar <$> arbitrary)
        ]

instance Arbitrary Pair where
    arbitrary = (:=) <$> arbitrary <*> arbitrary

scaleDouble :: (Double -> Double) -> Gen a -> Gen a
scaleDouble f = scale (round . f . fromIntegral)

instance Arbitrary Group where
    arbitrary = Group <$> (mk <$> arbitraryName) <*> arbitrary

newtype Namelist = Namelist { getNamelist :: [Group] }

instance Show Namelist where
    show (Namelist gs) = show gs

instance Arbitrary Namelist where
    arbitrary = Namelist <$> scaleDouble sqrt arbitrary

max_stars_repo_name: philopon/namelist-hs

max_issues_repo_name: philopon/namelist-hs

max_forks_repo_name: philopon/namelist-hs


text:
module Main where

import           Isumi.Math.FFT.Bruun

import           Control.DeepSeq            (NFData)
import           Control.Monad.Random.Class
import           Criterion.Main
import           Criterion.Types
import           Data.Complex
import qualified Data.Vector.Unboxed        as UV

main :: IO ()
main = do
  realInputs <- sampleInput
  complexInputs <- sampleInputC
  defaultMainWith config [
      bgroup "real input" [
        bgroup "fftDirect" $
          genBenchs (\xs -> fftDirect $ UV.map (:+0) xs) realInputs
        ,
        bgroup "fftBruun" $
          genBenchs fftBruun realInputs
      ],
      bgroup "complex input" [
        bgroup "fftDirect" $
          genBenchs fftDirect complexInputs
        ,
        bgroup "fftBruun" $
          genBenchs fftBruunC complexInputs
      ]
    ]

config = defaultConfig
  { reportFile = Just "dist/bench_report.html"
  }

genBenchs :: (UV.Unbox a, NFData b)
          => (UV.Vector a -> b) -> [UV.Vector a] -> [Benchmark]
genBenchs f = fmap (\xs -> bench (show . UV.length $ xs) $ nf f xs)

sampleInputC :: MonadRandom m => m [UV.Vector (Complex Double)]
sampleInputC = do
  reals <- sampleInput
  comps <- sampleInput
  pure $ zipWith (UV.zipWith (:+)) reals comps

sampleInput :: MonadRandom m => m [UV.Vector Double]
sampleInput = do
  flip traverse [0..11] $ \e -> do
    xs <- take (2^e) <$> getRandoms
    pure $ UV.fromList xs

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun