text:
{-# LANGUAGE TypeOperators, DeriveDataTypeable, FlexibleInstances, ScopedTypeVariables, ImplicitParams, ViewPatterns #-}
module Target.Regression 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, unormal, unit)
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

tmax :: Double
tmax = 1.000

dt :: Double
dt = 1.000e-2

(.==.) :: (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 = \((opts::[(T.Text,T.Text)])) -> \(Plot (pos::[(T.Text,T.Text)]) (plr::Prob ([Radian]))) -> Plot (append opts pos) plr

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}

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)))

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 = \(((_arg0)::Double)) -> \(((_arg1)::Bool)) -> case (_arg0,_arg1) of {((p::Double),True ) -> log p; ((p::Double),False ) -> log (1-p)}

bernoulli01LogPdf :: (Double -> ((Int -> Double)))
bernoulli01LogPdf = \(((_arg0)::Double)) -> \(((_arg1)::Int)) -> case (_arg0,_arg1) of {((p::Double),1) -> log p; ((p::Double),0) -> log (1-p)}

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))))

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

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))))

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

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

regress :: Prob ([((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))])
regress = (normal 0 1)>>=(\((offset::Double)) -> (gamma 1 1)>>=(\((sigma::Double)) -> (normal 0 1)>>=(\((slope::Double)) -> repeat 50 ((normal 0 1)>>=(\((w::Double)) -> (normal (offset+(slope*w)) sigma)>>=(\((y::Double)) -> return (R.X R.:& Y R.:= (y) R.:& W R.:= (w))))))))

regress1 :: Prob ([((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))])
regress1 = repeat 50 ((uniform 0.000 1.000)>>=(\((w::Double)) -> (normal ((1*w)-0.500) 0.100)>>=(\((y::Double)) -> return (R.X R.:& Y R.:= (y) R.:& W R.:= (w)))))

data W = W deriving Show
instance R.Name W where
   name = W
data Y = Y deriving Show
instance R.Name Y where
   name = Y
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
data Offset = Offset deriving Show
instance R.Name Offset where
   name = Offset
data Sigma = Sigma deriving Show
instance R.Name Sigma where
   name = Sigma
data Slope = Slope deriving Show
instance R.Name Slope where
   name = Slope

target = do
  ((fakedata::[((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))]))::[((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))] <- sample (regress1::Prob ([((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))]))
  let prims :: ((R.X R.:& Inisam R.::: Prob ([Double]) R.:& VToRec R.::: ((Vector Double) -> ((R.X R.:& Slope R.::: Double R.:& Sigma R.::: Double R.:& Offset R.::: Double) (Id KindStar))) R.:& Postgrad R.::: ((Vector Double) -> ((Double,(Vector Double)))) R.:& Posterior R.::: (([Double]) -> Double)) (Id KindStar))
      prims = let {(final0::[((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))]) = fakedata;
 posterior = runP1 (bindP headP (\(offset) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(sigma) -> bindP headP (\(slope) -> bindP (forP final0 (\((R.X R.:& Y R.:= (y::Double) R.:& W R.:= (w::Double))) -> returnP ())) (\(final0_pars) -> returnP (((normalLogPdf 0 1) offset)+(((gammaLogPdf 1 1) sigma)+(((normalLogPdf 0 1) slope)+((sum (map (\(((R.X R.:& Y R.:= (y::Double) R.:& W R.:= (w::Double)),(()))) -> ((normalLogPdf 0 1) w)+(((normalLogPdf (offset+(slope*w)) sigma) y)+0)) (zip final0 final0_pars)))+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)>>=(\((_offset_pos)) -> let {offset = _v@>_offset_pos;
 } in (((post_incr ((\(x) -> ((log 1)-(0.500*(log ((2.000*pi)*1))))-(((x-0)**2)/(2*1))) offset))>>(grad_incr _offset_pos (0-((0+((2*1)*((2*(offset-0))*(1-0))))/((2*1)*(2*1))))))>>(return ()))>>((count_incr 1)>>=(\((_sigma_pos)) -> let {sigma = exp (_v@>_sigma_pos);
 } in (((post_incr ((\(x) -> ((((1-1)*(log x))+((0.000-x)/1))-(1*(log 1)))-(gammaln 1)) sigma))>>(grad_incr _sigma_pos (((\((_x)) -> (exp _x)*1) (_v@>_sigma_pos))*((((0+((1-1)*((1/sigma)*1)))+((0+(1*(0-1)))/(1*1)))-0)-0))))>>(return ()))>>((count_incr 1)>>=(\((_slope_pos)) -> let {slope = _v@>_slope_pos;
 } in (((post_incr ((\(x) -> ((log 1)-(0.500*(log ((2.000*pi)*1))))-(((x-0)**2)/(2*1))) slope))>>(grad_incr _slope_pos (0-((0+((2*1)*((2*(slope-0))*(1-0))))/((2*1)*(2*1))))))>>(return ()))>>((forM final0 (\((R.X R.:& Y R.:= (y::Double) R.:& W R.:= (w::Double))) -> ((post_incr ((\(x) -> ((log 1)-(0.500*(log ((2.000*pi)*1))))-(((x-0)**2)/(2*1))) w))>>(return ()))>>(((post_incr ((\(x) -> ((log 1)-(0.500*(log ((2.000*pi)*sigma))))-(((x-(offset+(slope*w)))**2)/(2*sigma))) y))>>((((return ())>>(grad_incr _slope_pos (1*(0-((0+((2*sigma)*((2*(y-(offset+(slope*w))))*(0-(0+(0+(w*1)))))))/((2*sigma)*(2*sigma)))))))>>(grad_incr _sigma_pos (((\((_x)) -> (exp _x)*1) (_v@>_sigma_pos))*((0-(0+(0.500*((1/((2.000*pi)*sigma))*(0+((2.000*pi)*1))))))-((0+(0-(((y-(offset+(slope*w)))**2)*(0+(2*1)))))/((2*sigma)*(2*sigma)))))))>>(grad_incr _offset_pos (1*(0-((0+((2*sigma)*((2*(y-(offset+(slope*w))))*(0-(1+0)))))/((2*sigma)*(2*sigma))))))))>>(return (R.X R.:& Y R.:= (y) R.:& W R.:= (w))))))>>=(\((final0::[((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))])) -> (readSTRef postref)>>=(\(postval) -> (VS.unsafeFreeze gradref)>>=(\(gradval) -> return ((postval,gradval)))))))))))))));
 vToRec = \(theta) -> runP (toList theta) (bindP headP (\(offset) -> bindP (fmapP (\((_x)) -> exp _x) headP) (\(sigma) -> bindP headP (\(slope) -> bindP (forP final0 (\((R.X R.:& Y R.:= (y::Double) R.:& W R.:= (w::Double))) -> returnP ())) (\(final0_pars) -> returnP (R.X R.:& Slope R.:= (slope) R.:& Sigma R.:= (sigma) R.:& Offset R.:= (offset)))))));
 inisam = (normal 0 1)>>=(\(offset) -> (gamma 1 1)>>=(\(sigma) -> (normal 0 1)>>=(\(slope) -> (repeat (length final0) ((normal 0 1)>>=(\(w) -> (normal (offset+(slope*w)) sigma)>>=(\(y) -> return (((R.X R.:& Y R.:= (y) R.:& W R.:= (w)),(concat [])))))))>>=(\(final0_ret) -> let {(final0::[((R.X R.:& Y R.::: Double R.:& W R.::: Double) (Id KindStar))]) = map fst final0_ret;
 final0_iniret = map snd final0_ret;
 } in return (concat ((offset:[]):(((log sigma):[]):((slope:[]):((concat final0_iniret):[])))))))));
 } 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.:& Slope R.::: Double R.:& Sigma R.::: Double R.:& Offset 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:
-- |
-- Tests for data serialization instances
module Tests.Serialization where

import Data.Binary (Binary,decode,encode)
import Data.Aeson  (FromJSON,ToJSON,Result(..),toJSON,fromJSON)
import Data.Typeable

import Statistics.Distribution.Beta           (BetaDistribution)
import Statistics.Distribution.Binomial       (BinomialDistribution)
import Statistics.Distribution.CauchyLorentz
import Statistics.Distribution.ChiSquared     (ChiSquared)
import Statistics.Distribution.Exponential    (ExponentialDistribution)
import Statistics.Distribution.FDistribution  (FDistribution)
import Statistics.Distribution.Gamma          (GammaDistribution)
import Statistics.Distribution.Geometric
import Statistics.Distribution.Hypergeometric
import Statistics.Distribution.Laplace        (LaplaceDistribution)
import Statistics.Distribution.Normal         (NormalDistribution)
import Statistics.Distribution.Poisson        (PoissonDistribution)
import Statistics.Distribution.StudentT
import Statistics.Distribution.Transform      (LinearTransform)
import Statistics.Distribution.Uniform        (UniformDistribution)
import Statistics.Types

import Test.Tasty            (TestTree, testGroup)
import Test.Tasty.QuickCheck (testProperty)
import Test.QuickCheck         as QC

import Tests.Helpers
import Tests.Orphanage ()

tests :: TestTree
tests = testGroup "Test for data serialization"
  [ serializationTests (T :: T (CL Float))
  , serializationTests (T :: T (CL Double))
  , serializationTests (T :: T (PValue Float))
  , serializationTests (T :: T (PValue Double))
  , serializationTests (T :: T (NormalErr Double))
  , serializationTests (T :: T (ConfInt   Double))
  , serializationTests' "T (Estimate NormalErr Double)" (T :: T (Estimate NormalErr Double))
  , serializationTests' "T (Estimate ConfInt Double)" (T :: T (Estimate ConfInt   Double))
  , serializationTests (T :: T (LowerLimit Double))
  , serializationTests (T :: T (UpperLimit Double))
    -- Distributions
  , serializationTests (T :: T BetaDistribution        )
  , serializationTests (T :: T CauchyDistribution      )
  , serializationTests (T :: T ChiSquared              )
  , serializationTests (T :: T ExponentialDistribution )
  , serializationTests (T :: T GammaDistribution       )
  , serializationTests (T :: T LaplaceDistribution     )
  , serializationTests (T :: T NormalDistribution      )
  , serializationTests (T :: T UniformDistribution     )
  , serializationTests (T :: T StudentT                )
  , serializationTests (T :: T (LinearTransform NormalDistribution))
  , serializationTests (T :: T FDistribution           )
  , serializationTests (T :: T BinomialDistribution       )
  , serializationTests (T :: T GeometricDistribution      )
  , serializationTests (T :: T GeometricDistribution0     )
  , serializationTests (T :: T HypergeometricDistribution )
  , serializationTests (T :: T PoissonDistribution        )
  ]

serializationTests
  :: (Eq a, Typeable a, Binary a, Show a, Read a, ToJSON a, FromJSON a, Arbitrary a)
  => T a -> TestTree
serializationTests t = serializationTests' (typeName t) t

-- Not all types are Typeable, unfortunately
serializationTests'
  :: (Eq a, Binary a, Show a, Read a, ToJSON a, FromJSON a, Arbitrary a)
  => String -> T a -> TestTree
serializationTests' name t = testGroup ("Tests for: " ++ name)
  [ testProperty "show/read" (p_showRead t)
  , testProperty "binary"    (p_binary   t)
  , testProperty "aeson"     (p_aeson    t)
  ]

p_binary :: (Eq a, Binary a) => T a -> a -> Bool
p_binary _ a = a == (decode . encode) a

p_showRead :: (Eq a, Read a, Show a) => T a -> a -> Bool
p_showRead _ a = a == (read . show) a

p_aeson :: (Eq a, ToJSON a, FromJSON a) => T a -> a -> Bool
p_aeson _ a = Data.Aeson.Success a == (fromJSON . toJSON) a

ext: hs

max_stars_repo_name: haskell/statistics

max_issues_repo_name: haskell/statistics

max_forks_repo_name: haskell/statistics


text:
module Vote
  ( processVote,
    reprocessVotes,
    encryptBallot,
    decryptBallot,
  )
where

import Conduit (foldlC, mapC)
import Control.Lens
import Control.Monad.Trans.Reader (ReaderT)
import Data.Binary.Get (getWord32be, runGet)
import Data.Binary.Put (putWord32be, runPut)
import qualified Data.ByteString as B
import qualified Data.ByteString.Char8 as BC
import qualified Data.Map.Strict as M
import qualified Data.Text as T
import Data.Time (NominalDiffTime, UTCTime, addUTCTime, getCurrentTime)
import Data.Time.Clock.POSIX (posixSecondsToUTCTime, utcTimeToPOSIXSeconds)
import Database.Persist.Sql (SqlBackend)
import Import hiding (toList)
import Model.IsaacVersion
import Numeric.LinearAlgebra
import Ranks
import qualified Web.ClientSession as WS

processVote :: (MonadIO m) => IsaacVersion -> Int -> Int -> UTCTime -> Text -> Text -> ReaderT SqlBackend m ()
processVote ver w l timestamp voter rawBallot = do
  Entity w' _ <- getBy404 $ UniqueItem ver w
  Entity l' _ <- getBy404 $ UniqueItem ver l
  insert400_ $ Ballot rawBallot timestamp
  insert400_ $ Vote ver w' l' timestamp voter
  return ()

reprocessVotes :: (MonadIO m, MonadResource m) => ReaderT SqlBackend m ()
reprocessVotes = do
  items <- M.fromList . map (\(Entity a b) -> (a, b)) <$> selectList [] []
  bms <- runConduit $ selectSource [] [] .| mapC entityVal .| beatMatrix
  recalculatePairs items bms
  now <- liftIO getCurrentTime
  deleteWhere [BallotTimestamp <=. (- validTime) `addUTCTime` now]
  return ()

type Matchup = (Key Item, Key Item)

type BeatMatrix = M.Map Matchup Int

beatMatrix :: Monad m => ConduitT Vote o m (M.Map IsaacVersion BeatMatrix)
beatMatrix = foldlC upd M.empty
  where
    upd bm v = bm & at ver . non M.empty . at (w, l) . non 0 +~ 1
      where
        w = v ^. voteWinner
        l = v ^. voteLoser
        ver = v ^. voteVersion

ranks :: BeatMatrix -> [Key Item] -> [(Key Item, R)]
ranks bm items =
  zip items . toList $ ilsrPairwise bm' 0.01
  where
    n = length items
    itemToIdx = M.fromList (zip items [0 ..])
    bm' = assoc (n, n) 0 (convR <$> M.toList bm)
    convR ((w, l), s) = ((itemToIdx M.! w, itemToIdx M.! l), fromIntegral s)

recalculatePairs :: (MonadIO m, MonadResource m) => M.Map (Key Item) Item -> M.Map IsaacVersion BeatMatrix -> ReaderT SqlBackend m ()
recalculatePairs items bms =
  ifor_ bms $ \ver bm ->
    for_ (ranks bm (itemsFor ver items)) $ \(itemId, rating) ->
      update
        itemId
        [ ItemRating =. rating,
          ItemVotes =. sumOf (ifolded . ifiltered (\(i1, i2) _ -> i1 == itemId || i2 == itemId)) bm
        ]
  where
    itemsFor ver = findIndicesOf ifolded (\i -> i ^. itemVersion == ver)

encodeBallot :: UTCTime -> Int -> Int -> B.ByteString
encodeBallot expiry winner loser = toStrict . runPut $ do
  putWord32be (truncate . utcTimeToPOSIXSeconds $ expiry)
  putWord32be (fromIntegral winner)
  putWord32be (fromIntegral loser)

decodeBallot :: B.ByteString -> (UTCTime, Int, Int)
decodeBallot = go . fromStrict
  where
    go = runGet $ do
      expiry <- getWord32be
      winner <- getWord32be
      loser <- getWord32be
      return
        ( posixSecondsToUTCTime (fromIntegral expiry),
          fromIntegral winner,
          fromIntegral loser
        )

-- 1 hour
validTime :: NominalDiffTime
validTime = 3600

encryptBallot :: UTCTime -> Int -> Int -> Handler Text
encryptBallot now winner loser = do
  key <- getsYesod appBallotKey
  out <- liftIO $ WS.encryptIO key (encodeBallot (validTime `addUTCTime` now) winner loser)
  return . T.pack . BC.unpack $ out

decryptBallot :: UTCTime -> Text -> Handler (Int, Int)
decryptBallot now b = do
  key <- getsYesod appBallotKey
  (expiry, winner, loser) <-
    case decodeBallot <$> (WS.decrypt key . BC.pack . T.unpack) b of
      Just r -> pure r
      Nothing -> error "invalid ballot!"
  if expiry < now
    then error "expired ballot!"
    else return (winner, loser)

max_stars_repo_name: jsza/isaacranks

max_issues_repo_name: jsza/isaacranks

max_forks_repo_name: jsza/isaacranks


text:
{-#LANGUAGE ForeignFunctionInterface#-}
-- |This  module provides slow but functional means for exporting images from and to
--  CArrays, which can easily be passed into foreign functions.
module CV.Conversions (
     -- Arrays of Double
     copyCArrayToImage
    ,copy8UCArrayToImage
    ,copyImageToCArray
     -- Arrays of Float
    ,copyFCArrayToImage
    ,copyImageToFCArray
    -- * Complex arrays
    ,copyComplexCArrayToImage
    ,copyImageToComplexCArray
    -- * Copying
    ,copyImageToExistingCArray
    -- * Acquiring images from pointers
    ,unsafe8UC_RGBFromPtr
    ,unsafe8UC_BGRFromPtr
    ,unsafe8UC_MONOFromPtr
    ,acquireImageSlowF'
    ,acquireImageSlow'
    ,acquireImageSlow8URGB'
    ) where

import Data.Complex as C

import CV.Image
import Data.Word

import Data.Array.CArray
import Data.Array.IArray

import Foreign.C.Types
import Foreign.Ptr
import System.IO.Unsafe

unsafe8UC_RGBFromPtr :: (Int,Int) -> Ptr Word8 -> IO (Image RGB D8)
unsafe8UC_RGBFromPtr (w,h) ptr = S `fmap`  creatingBareImage (acquireImageSlow8URGB' w h ptr)

unsafe8UC_BGRFromPtr :: (Int,Int) -> Ptr Word8 -> IO (Image RGB D8)
unsafe8UC_BGRFromPtr (w,h) ptr = S `fmap`  creatingBareImage (acquireImageSlow8UBGR' w h ptr)

unsafe8UC_MONOFromPtr :: (Int,Int) -> Ptr Word8 -> IO (Image GrayScale D8)
unsafe8UC_MONOFromPtr (w,h) ptr = S `fmap`  creatingBareImage (acquireImageSlow8U' w h ptr)

-- |Copy the contents of a CArray into CV.Image type.
copy8UCArrayToImage :: CArray (Int,Int) Word8 -> Image GrayScale D8
copy8UCArrayToImage carr = S $ unsafePerformIO $
                          creatingBareImage (withCArray carr (acquireImageSlow8U' w h))
    where
     ((sx,sy),(ex,ey)) = bounds carr
     (w,h) = (fromIntegral $ ex-sx+1, fromIntegral $ ey-sy+1)
--
-- |Copy the contents of a CArray into CV.Image type.
copyCArrayToImage :: CArray (Int,Int) Double -> Image GrayScale D32
copyCArrayToImage carr = S $ unsafePerformIO $
                          creatingBareImage (withCArray carr (acquireImageSlow' w h))
    where
     ((sx,sy),(ex,ey)) = bounds carr
     (w,h) = (fromIntegral $ ex-sx+1, fromIntegral $ ey-sy+1)

-- |Copy CArray of floats to image
copyFCArrayToImage :: CArray (Int,Int) Float -> Image GrayScale D32
copyFCArrayToImage carr = S $ unsafePerformIO $
                          creatingBareImage (withCArray carr (acquireImageSlowF' w h))
    where
     ((sx,sy),(ex,ey)) = bounds carr
     (w,h) = (fromIntegral $ ex-sx+1, fromIntegral $ ey-sy+1)

-- |Copy D32 grayscale image to CArray
copyImageToFCArray :: Image GrayScale D32 -> CArray (Int,Int) Float
copyImageToFCArray (S img) = unsafePerformIO $
         withBareImage img $ \cimg ->
          createCArray ((0,0),(w-1,h-1)) (exportImageSlowF' cimg) --({#call exportImageSlow#} cimg)
    where
     (w,h) = getSize img

-- |Copy the real part of an array to image
copyComplexCArrayToImage :: CArray (Int,Int) (C.Complex Double) -> Image GrayScale D32
copyComplexCArrayToImage carr = S $ unsafePerformIO $
                          creatingBareImage (withCArray carr (acquireImageSlowComplex' w h))
    where
     ((sx,sy),(ex,ey)) = bounds carr
     (w,h) = (fromIntegral $ ex-sx+1, fromIntegral $ ey-sy+1)

-- |Copy the contents of a CV.Image into a CArray.
copyImageToCArray :: Image GrayScale D32 -> CArray (Int,Int) Double
copyImageToCArray (S img) = unsafePerformIO $
         withBareImage img $ \cimg ->
          createCArray ((0,0),(w-1,h-1)) (exportImageSlow' cimg) --({#call exportImageSlow#} cimg)
    where
     (w,h) = getSize img
-- |Copy the contents of CV.Image into a pre-existing CArray.
--
copyImageToExistingCArray (S img) arr =
         withBareImage img $ \cimg ->
          withCArray arr $ \carr -> (exportImageSlow' cimg carr) --({#call exportImageSlow#} cimg)
    where
     (w,h) = getSize img

-- |Copy image as a real part of a complex CArray
copyImageToComplexCArray :: Image GrayScale D32 -> CArray (Int,Int) (C.Complex Double)
copyImageToComplexCArray (S img) = unsafePerformIO $
         withBareImage img $ \cimg ->
          createCArray ((0,0),(w-1,h-1)) (exportImageSlowComplex' cimg) --({#call exportImageSlow#} cimg)
    where
     (w,h) = getSize img

foreign import ccall safe "CV/cvWrapLeo.h exportImageSlow"
  exportImageSlow' :: ((Ptr (BareImage)) -> ((Ptr Double) -> (IO ())))

foreign import ccall safe "CV/cvWrapLeo.h exportImageSlowF"
  exportImageSlowF' :: ((Ptr (BareImage)) -> ((Ptr Float) -> (IO ())))

foreign import ccall safe "CV/cvWrapLeo.h exportImageSlowComplex"
  exportImageSlowComplex' :: ((Ptr (BareImage)) -> ((Ptr (C.Complex Double)) -> (IO ())))

foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlow"
  acquireImageSlow' :: (Int -> (Int -> ((Ptr Double) -> (IO (Ptr (BareImage))))))

foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlowF"
  acquireImageSlowF' :: (Int -> (Int -> ((Ptr Float) -> (IO (Ptr (BareImage))))))

foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlow8URGB"
  acquireImageSlow8URGB' :: (Int -> (Int -> ((Ptr Word8) -> (IO (Ptr (BareImage))))))

foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlow8UBGR"
  acquireImageSlow8UBGR' :: (Int -> (Int -> ((Ptr Word8) -> (IO (Ptr (BareImage))))))

foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlow8U"
  acquireImageSlow8U' :: (Int -> (Int -> ((Ptr Word8) -> (IO (Ptr (BareImage))))))

foreign import ccall safe "CV/cvWrapLeo.h acquireImageSlowComplex"
  acquireImageSlowComplex' :: (Int -> (Int -> ((Ptr (C.Complex Double)) -> (IO (Ptr (BareImage))))))

max_stars_repo_name: maaleske/CV

max_issues_repo_name: maaleske/CV

max_forks_repo_name: maaleske/CV


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
module Data.Network.Activation
  (ActivationFunction, unAct, act, name, activate, derivate, logistic, softmax, softmaxMat, relu) where

import Numeric.AD.Mode.Kahn
import Numeric.LinearAlgebra
import Data.Network.Utils(Differentiable, squeeze, derivate, (.*))

data ActivationFunction a = Act {unAct :: (a -> (a,a)), name :: String}

act :: Floating a => (forall s. AD s (Kahn a) -> (AD s (Kahn a))) -> String -> ActivationFunction a
act f = Act (diff' f)

activate :: ActivationFunction a -> a -> a
activate af = fst . (unAct af)

instance Differentiable a a (ActivationFunction a) where
  derivate af = snd . (unAct af)

-- derivate :: ActivationFunction a -> a -> a
-- derivate af = snd . (unAct af)

logistic :: Floating a => ActivationFunction a
logistic = let f x =  1 / (1 + exp (-x))
           in  act f "logistic"

deltaf :: Floating a => Int -> Int -> a
deltaf i j = if (i == j) then 1 else 0

softmax :: (Fractional e, Linear e Vector, Floating (Vector e), Container Vector e) => ActivationFunction (Vector e)
softmax = let f x = let n = sumElements (exp x)
                    in (recip n) .* (exp x)
              f' x = 1 -- Softmax is used only for last layer, so derivative is not needed.
          in  Act (\x -> (f x, f' x)) "softmax"

-- Matrix version for batch computing
softmaxMat :: (Floating e, Container Matrix e, Num (Matrix e), Numeric e, Floating (Matrix e), Floating (Vector e)) =>
  ActivationFunction (Matrix e)
softmaxMat = let f x = let n = fromColumns (replicate (cols x) $ squeeze (exp x))
                       in (recip n) * (exp x)
                 f' x = 1 -- Softmax is used only for last layer, so derivative is not needed.
          in  Act (\x -> (f x, f' x)) "softmax"

relu = let f x = max x 0
       in  act f "relu"

-- We can add activation functions and their derivative using the "diff" function from Numeric.AD

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE UndecidableInstances #-}

module Numeric.PRIMME.Types
  ( PrimmeException (..),
    BlasDatatype (..),
    BlasRealPart,
    BlasDatatypeTag (..),
    MBlock (..),
    Block (..),
    PrimmeOperator,
    PrimmeTarget (..),
    PrimmeMethod (..),
    PrimmeInfo (..),
    PrimmeStats (..),
    PrimmeEventInfo (..),
    PrimmeMonitor (..),
    PrimmeOptions (..),
    PrimmeInt,
    Cprimme_params,
    CmatrixMatvec,
    CmonitorFun,
  )
where

import Control.Exception.Safe
import Control.Monad.ST (RealWorld)
import Data.Complex
import Data.Int (Int32, Int64)
import Data.Kind
import Data.Text
import Data.Vector.Storable (MVector, Vector)
import Foreign.C.Types (CChar, CDouble, CInt, CLong)
import Foreign.Ptr (Ptr)
import Foreign.Storable

#if !defined(PRIMME_INT_SIZE) || PRIMME_INT_SIZE == 64
type PrimmeInt = CLong -- Int64
#elif PRIMME_INT_SIZE == 0
type PrimmeInt = CInt
#elif PRIMME_INT_SIZE == 32
type PrimmeInt = Int32
#else
#error "invalid value of PRIMME_INT_SIZE"
#endif

data Cprimme_params

-- | Type of blocked matrix-vector product accepted by PRIMME C library.
-- @
-- void (*matrixMatvec)
--       ( void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize,
--         struct primme_params *primme, int *ierr);
-- @
type CmatrixMatvec = Ptr () -> Ptr PrimmeInt -> Ptr () -> Ptr PrimmeInt -> Ptr CInt -> Ptr Cprimme_params -> Ptr CInt -> IO ()

-- | Type of logging function.
-- @
-- void (*monitorFun)(void *basisEvals, int *basisSize, int *basisFlags, int *iblock,
--                    int *blockSize, void *basisNorms, int *numConverged, void *lockedEvals,
--                    int *numLocked, int *lockedFlags, void *lockedNorms, int *inner_its,
--                    void *LSRes, const char *msg, double *time, primme_event *event,
--                    struct primme_params *primme, int *ierr)
-- @
type CmonitorFun =
  Ptr () ->
  Ptr CInt ->
  Ptr CInt ->
  Ptr CInt ->
  Ptr CInt ->
  Ptr () ->
  Ptr CInt ->
  Ptr () ->
  Ptr CInt ->
  Ptr CInt ->
  Ptr () ->
  Ptr CInt ->
  Ptr () ->
  Ptr CChar ->
  Ptr CDouble ->
  Ptr CInt ->
  Ptr Cprimme_params ->
  Ptr CInt ->
  IO ()

-- | Exceptions thrown by this package.
data PrimmeException
  = PrimmeMaximumIterationsReached
  | PrimmeUserFailure
  | PrimmeLapackFailure
  | PrimmeOtherFailure {-# UNPACK #-} !CInt
  deriving stock (Show)

instance Exception PrimmeException

-- | A GADT which allows to dispatch between BLAS types at runtime.
data BlasDatatypeTag :: (Type -> Type) where
  FloatTag :: BlasDatatypeTag Float
  DoubleTag :: BlasDatatypeTag Double
  ComplexFloatTag :: BlasDatatypeTag (Complex Float)
  ComplexDoubleTag :: BlasDatatypeTag (Complex Double)

type family BlasRealPart a where
  BlasRealPart Float = Float
  BlasRealPart Double = Double
  BlasRealPart (Complex Float) = Float
  BlasRealPart (Complex Double) = Double

-- | BLAS datatype.
class
  ( Storable a,
    Storable (BlasRealPart a),
    Floating a,
    Floating (BlasRealPart a),
    Show a,
    Show (BlasRealPart a)
  ) =>
  BlasDatatype a
  where
  blasTag :: proxy a -> BlasDatatypeTag a

instance BlasDatatype Float where blasTag _ = FloatTag

instance BlasDatatype Double where blasTag _ = DoubleTag

instance BlasDatatype (Complex Float) where blasTag _ = ComplexFloatTag

instance BlasDatatype (Complex Double) where blasTag _ = ComplexDoubleTag

-- | Mutable dense matrix in column-major order.
data MBlock s a = MBlock
  { -- | Shape of the matrix as @(number of rows, number of columns)@
    mutableBlockShape :: {-# UNPACK #-} !(Int, Int),
    -- | Stride along the second dimension (typically called @ldim@ or "leading
    -- dimension" in BLAS-like libraries). Since the matrix is in column-major
    -- order, stride along the first dimension is always 1.
    mutableBlockStride :: {-# UNPACK #-} !Int,
    -- | The actual data buffer of length @number of columns * mutableBlockStride@
    mutableBlockData :: {-# UNPACK #-} !(MVector s a)
  }

-- | Dense matrix in column-major order.
data Block a = Block
  { -- | Shape of the matrix as @(number of rows, number of columns)@
    blockShape :: {-# UNPACK #-} !(Int, Int),
    -- | Stride along the second dimension (typically called @ldim@ or "leading
    -- dimension" in BLAS-like libraries). Since the matrix is in column-major
    -- order, stride along the first dimension is always 1.
    blockStride :: {-# UNPACK #-} !Int,
    -- | The actual data buffer of length @#columns * blockStride@
    blockData :: {-# UNPACK #-} !(Vector a)
  }
  deriving stock (Show)

-- | Matrix-matrix product which defines an operator @A@.
type PrimmeOperator a =
  -- | Matrix @x@ to which we apply @A@
  Block a ->
  -- | Where to store @Ax@
  MBlock RealWorld a ->
  -- | Procedure is run in 'IO' since it will typically involve FFI calls.
  IO ()

-- | Which eigenpairs to target.
--
-- /Note:/ PRIMME C library also supports computing non-extremal eigenvalues.
-- Wrappers for it may be added in the future.
data PrimmeTarget
  = -- | Target smallest algebraic eigenvalues
    PrimmeSmallest
  | -- | Target largest algebraic eigenvalues
    PrimmeLargest
  deriving stock (Read, Show, Eq)

-- | Which method to use.
data PrimmeMethod
  = PrimmeDynamic
  | -- | GD+k
    PrimmeGDOlsenPlusK
  | -- | JDQMR
    PrimmeJDQMRETol
  deriving stock (Read, Show, Eq)

newtype PrimmeMonitor = PrimmeMonitor
  { unPrimmeMonitor :: forall a. BlasDatatype a => PrimmeInfo a -> IO Bool
  }

data PrimmeInfo a = PrimmeInfo (PrimmeEventInfo a) PrimmeStats

deriving stock instance BlasDatatype a => Show (PrimmeInfo a)

data PrimmeStats = PrimmeStats
  { pStatsNumOuterIterations :: !Int,
    pStatsNumRestarts :: !Int,
    pStatsNumMatvecs :: !Int,
    pStatsTimeMatvec :: !Double,
    pStatsTimeOrtho :: !Double,
    pStatsTimeDense :: !Double
  }
  deriving stock (Show)

data PrimmeEventInfo a
  = PrimmeOuterInfo !(Vector (BlasRealPart a)) !(Vector (BlasRealPart a))
  | PrimmeInnerInfo
  | PrimmeLockedInfo
  | PrimmeConvergedInfo !Int !(BlasRealPart a) !(BlasRealPart a)
  | PrimmeMessageInfo !Text

deriving stock instance BlasDatatype a => Show (PrimmeEventInfo a)

-- | Specify some more information about the eigenvalue problem at hand.
--
-- /Warning:/ it is of utmost importance to use the correct dimension! If your
-- 'PrimmeOperator' has a different dimension than specified here, stuff will
-- break with segmentation faults...
data PrimmeOptions = PrimmeOptions
  { -- | Dimension of the operator. Please, set it carefully!
    pDim :: !Int,
    -- | Number of eigenpairs to compute
    pNumEvals :: !Int,
    -- | Which eigenpairs to target
    pTarget :: !PrimmeTarget,
    -- | Which algorithm to use
    pMethod :: !PrimmeMethod,
    -- | Tolerance
    pEps :: !Double,
    -- | Maximal basis size. [See also](http://www.cs.wm.edu/~andreas/software/doc/appendix.html#c.primme_params.maxBasisSize)
    pMaxBasisSize :: !Int,
    -- | Maximum number of Ritz vectors kept after restarting the basis.
    -- [See also](http://www.cs.wm.edu/~andreas/software/doc/appendix.html#c.primme_params.minRestartSize)
    pMinRestartSize :: !Int,
    -- | Maximal block size. [See also](http://www.cs.wm.edu/~andreas/software/doc/appendix.html#c.primme_params.maxBlockSize)
    pMaxBlockSize :: !Int,
    -- | Logging
    pLogAction :: !(Either Int PrimmeMonitor)
  }

max_stars_repo_name: twesterhout/primme-hs

max_issues_repo_name: twesterhout/primme-hs

max_forks_repo_name: twesterhout/primme-hs


text:
module Language.Scheme.Parser where

import Language.Scheme.Types
import Language.Scheme.Error
import Language.Scheme.Error.Types

import Text.ParserCombinators.Parsec hiding (spaces)
import Numeric (readHex, readInt, readOct, readFloat)
import Control.Monad (liftM)
import Control.Applicative ((<$>), (<*>))
import Data.Complex (Complex((:+)))
import Data.Ratio ((%))
import Control.Monad.Except (throwError)

symbol :: Parser Char
symbol = oneOf "!#$%&|*+-/:<=>?@^_~"

spaces :: Parser ()
spaces = skipMany1 space

escapedChars :: Parser Char
escapedChars = do
  char '\\'
  x <- oneOf "\\\"nrt"
  return x

parseString :: Parser LispVal
parseString = do
  char '"'
  x <- many $ escapedChars <|> noneOf "\"\\"
  char '"'
  return $ String x

parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest <- many $ letter <|> digit <|> symbol
  let atom = first:rest
  return $ Atom atom

parseHash :: Parser LispVal
parseHash = do
  c <- hashchar
  case c of
    -- Will always be one of these values since
    -- hashchar parses only these characters
    't' -> return $ Bool True
    'f' -> return $ Bool False
    'd' -> parseNumber
    'x' -> parseNumberHex
    'o' -> parseNumberOct
    'b' -> parseNumberBin
    '\\' -> parseChar

parseChar :: Parser LispVal
parseChar = do
  c <- anyChar
  return $ Character c

getDouble :: String -> Double
getDouble = getValue . readFloat

toDouble :: LispVal -> Double
toDouble (Number (Real f))  = f
toDouble (Number (Integer n)) = fromIntegral n

parseReal :: Parser LispVal
parseReal = liftM (Number . Real . getDouble) float
    where float = (++) <$> digits <*> decimal
          decimal = (:) <$> char '.' <*> digits

parseRational :: Parser LispVal
parseRational = do
  numer <- digits
  char '/'
  denom <- digits
  return $ Number $ Rational ((read numer) % (read denom))

parseComplex :: Parser LispVal
parseComplex = do
  x <- (try parseReal <|> parseInteger)
  char '+'
  y <- (try parseReal <|> parseInteger)
  char 'i'
  return $ Number $ Complex (toDouble x :+ toDouble y)

digits :: Parser String
digits = many1 digit

hashchar :: Parser Char
hashchar = char '#' >> oneOf "boxdtf\\"

-- Gets the value from a read function such as readHex
-- Since the parser handles parsing we know it will always
-- produce a value if the parser successes
getValue :: [(a, String)] -> a
getValue [(x,_)] = x
getValue _       = error "Should not happen"

hexToNum :: (Num a, Eq a) => String -> a
hexToNum = getValue . readHex

octToNum :: (Num a, Eq a) => String -> a
octToNum = getValue . readOct

binToNum :: (Num a, Eq a) => String -> a
binToNum = getValue . readBin

isBinChar :: Char -> Bool
isBinChar '0' = True
isBinChar '1' = True
isBinChar _   = False

binCharToInt :: Char -> Int
binCharToInt '0' = 0
binCharToInt '1' = 1
binCharToInt _   = error "Not a binary character"

readBin :: (Num a, Eq a) => ReadS a
readBin = readInt 2 isBinChar binCharToInt

binDigits :: Parser String
binDigits = many1 $ oneOf "01"

hexDigits :: Parser String
hexDigits = many1 hexDigit

octDigits :: Parser String
octDigits = many1 octDigit

parseNumberBin :: Parser LispVal
parseNumberBin = liftM (Number . Integer . binToNum) binDigits

parseNumberHex :: Parser LispVal
parseNumberHex = liftM (Number . Integer . hexToNum) hexDigits

parseNumberOct :: Parser LispVal
parseNumberOct = liftM (Number . Integer . octToNum) octDigits

parseInteger :: Parser LispVal
parseInteger = liftM (Number . Integer . read) digits

parseNumber :: Parser LispVal
parseNumber = try parseReal <|> try parseRational <|> try parseComplex <|> parseInteger

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
  head <- endBy parseExpr spaces
  tail <- char '.' >> spaces >> parseExpr
  return $ DottedList head tail

parseQuoted :: Parser LispVal
parseQuoted = do
  char '\''
  x <- parseExpr
  return $ List [Atom "quote", x]

parseExpr :: Parser LispVal
parseExpr = parseHash
            <|> parseAtom
            <|> parseString
            <|> parseNumber
            <|> parseQuoted
            <|> do char '('
                   x <- try parseList <|> parseDottedList
                   char ')'
                   return x

readExpr :: String -> ThrowsError LispVal
readExpr input = case parse parseExpr "lisp" input of
                   Left err -> throwError $ Parser err
                   Right val  -> return val


text:
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RecordWildCards   #-}

{-# OPTIONS_GHC -fno-warn-orphans #-}

module Internal.CG(
    cgSolve, cgSolve',
    CGState(..), R, V
) where

import           Control.Arrow       ((***))
import           Internal.Algorithms (NormType (..), linearSolve, linearSolveLS, pnorm,
                                      relativeError)
import           Internal.Container
import           Internal.Element
import           Internal.IO
import           Internal.Matrix
import           Internal.Numeric
import           Internal.Sparse
import           Internal.Vector
import           Numeric.Vector      ()

{-
import Util.Misc(debug, debugMat)

(//) :: Show a => a -> String -> a
infix 0 // -- , ///
a // b = debug b id a

(///) :: V -> String -> V
infix 0 ///
v /// b = debugMat b 2 asRow v
-}

type V = Vector R

data CGState = CGState
    { cgp  :: Vector R  -- ^ conjugate gradient
    , cgr  :: Vector R  -- ^ residual
    , cgr2 :: R         -- ^ squared norm of residual
    , cgx  :: Vector R  -- ^ current solution
    , cgdx :: R         -- ^ normalized size of correction
    }

cg :: Bool -> (V -> V) -> (V -> V) -> CGState -> CGState
cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
  where
    ap1 = a p
    ap  | sym       = ap1
        | otherwise = at ap1
    pap | sym       = p <.> ap1
        | otherwise = norm2 ap1 ** 2
    alpha = r2 / pap
    dx = scale alpha p
    x' = x + dx
    r' = r - scale alpha ap
    r'2 = r' <.> r'
    beta = r'2 / r2
    p' = r' + scale beta p

    rdx = norm2 dx / max 1 (norm2 x)

conjugrad
  :: Bool -> GMatrix -> V -> V -> R -> R -> [CGState]
conjugrad sym a b = solveG sym (tr a !#>) (a !#>) (cg sym) b

solveG
    :: Bool
    -> (V -> V) -> (V -> V)
    -> ((V -> V) -> (V -> V) -> CGState -> CGState)
    -> V
    -> V
    -> R -> R
    -> [CGState]
solveG sym mat ma meth rawb x0' ϵb ϵx
    = takeUntil ok . iterate (meth mat ma) $ CGState p0 r0 r20 x0 1
  where
    a = if sym then ma else mat . ma
    b = if sym then rawb else mat rawb
    x0  = if x0' == 0 then konst 0 (dim b) else x0'
    r0  = b - a x0
    r20 = r0 <.> r0
    p0  = r0
    nb2 = b <.> b
    ok CGState {..}
        =  cgr2 <nb2*ϵb**2
        || cgdx < ϵx

takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil q xs = a++ take 1 b
  where
    (a,b) = break q xs

-- | Solve a sparse linear system using the conjugate gradient method with default parameters.
cgSolve
  :: Bool          -- ^ is symmetric
  -> GMatrix       -- ^ coefficient matrix
  -> Vector R      -- ^ right-hand side
  -> Vector R      -- ^ solution
cgSolve sym a b  = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0
  where
    n = max 10 (round $ sqrt (fromIntegral (dim b) :: Float))

-- | Solve a sparse linear system using the conjugate gradient method with default parameters.
cgSolve'
  :: Bool      -- ^ symmetric
  -> R         -- ^ relative tolerance for the residual (e.g. 1E-4)
  -> R         -- ^ relative tolerance for δx (e.g. 1E-3)
  -> Int       -- ^ maximum number of iterations
  -> GMatrix   -- ^ coefficient matrix
  -> Vector R  -- ^ initial solution
  -> Vector R  -- ^ right-hand side
  -> [CGState] -- ^ solution
cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es

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

instance Testable GMatrix
  where
    checkT _ = (ok,info)
      where
        sma = convo2 20 3
        x1 = vect [1..20]
        x2 = vect [1..40]
        sm = mkSparse sma
        dm = toDense sma

        s1 = sm !#> x1
        d1 = dm #> x1

        s2 = tr sm !#> x2
        d2 = tr dm #> x2

        sdia = mkDiagR 40 20 (vect [1..10])
        s3 =    sdia !#> x1
        s4 = tr sdia !#> x2
        ddia = diagRect 0 (vect [1..10])  40 20
        d3 = ddia #> x1
        d4 = tr ddia #> x2

        v = testb 40
        s5 = cgSolve False sm v
        d5 = denseSolve dm v

        symassoc = [((0,0),1.0),((1,1),2.0),((0,1),0.5),((1,0),0.5)]
        b = vect [3,4]
        d6 = flatten $ linearSolve (toDense symassoc) (asColumn b)
        s6 = cgSolve True (mkSparse symassoc) b

        info = do
            print sm
            disp (toDense sma)
            print s1; print d1
            print s2; print d2
            print s3; print d3
            print s4; print d4
            print s5; print d5
            print $ relativeError (pnorm Infinity) s5 d5
            print s6; print d6
            print $ relativeError (pnorm Infinity) s6 d6

        ok = s1==d1
          && s2==d2
          && s3==d3
          && s4==d4
          && relativeError (pnorm Infinity) s5 d5 < 1E-10
          && relativeError (pnorm Infinity) s6 d6 < 1E-10

        disp = putStr . dispf 2

        vect = fromList :: [Float] -> Vector Float

        convomat :: Int -> Int -> AssocMatrix
        convomat n k = [ ((i,j `mod` n),1) | i<-[0..n-1], j <- [i..i+k-1]]

        convo2 :: Int -> Int -> AssocMatrix
        convo2 n k = m1 ++ m2
          where
            m1 = convomat n k
            m2 = map (((+n) *** id) *** id) m1

        testb n = vect $ take n $ cycle ([0..10]++[9,8..1])

        denseSolve a = flatten . linearSolveLS a . asColumn

        -- mkDiag v = mkDiagR (dim v) (dim v) v

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
module DeltaSigma.PowerSpectrumDensity where

import Numeric.FFT
import Data.Complex

psd :: [Double] -> [Double]
psd xs = take n_half . fmap (/ (fromIntegral n_half))
        . fmap magnitude . fft . fmap (:+ 0) $ take n xs
  where
    n = length xs
    n_half = n `div` 2


text:
{-
  Programa em haskell para gerar Mandelbrot
  Trabalho 1 de Paradigmas de Programação
  Autoria: Natã Schmitt
-}

import Codec.Picture
import Text.Printf ( printf )
import Codec.FFmpeg.Juicy ( imageWriter )
import Codec.FFmpeg.Encode
import System.IO
import System.Exit ( exitFailure, exitSuccess )
import Data.WAVE
import Math.FFT ( dftRC )
import Data.Array.CArray
import qualified Data.Array.IArray
import Data.Complex
import Data.Maybe ( fromMaybe )
import Data.Char ( toLower )
import System.Process ( shell, callCommand )

type CComplex = (Double, Double)

-----------------------------------------------------------
-- Variáveis pré-definidas
-----------------------------------------------------------
musicPath, videoPath, imagesPath :: FilePath
width, height, maxIter, hueOffset, sensitivity, dropped :: Int
coordY, coordX, maxZoom, framerate, fixedZoom :: Double

musicPath = "./mandeloso.wav"
videoPath = "./mandelbrot.mp4"
imagesPath = "./anim/%d.png"

width = 1920          -- largura
height = 1080        -- altura
maxIter = 90        -- numero maximo de iterações
hueOffset = -150       -- caso queira que a cor inicial seja diferente

maxZoom = 1.6519276269095182e16     --Zoom maximo suportado

coordX =  0.36024044343761436323612524444954530848260780795858575048837581474019534605  -- Coordenadas X
coordY = -0.64131306106480317486037501517930206657949495228230525955617754306444857417  -- Coordenadas Y

framerate = 60        -- Framerate do vídeo

sensitivity = 20      -- Sensibilidade ao grave maior = menos sensível

estático :: Bool
estático = False       -- Quando True ele não fará zoom e se manterá estático em fixedZoom

fixedZoom = 17.5        -- Zoom caso estático for True

dropped = 0       -- Se quiser retomar uma animação pelo numero da imagem, use aqui

-----------------------------------------------------------
-- Funções diversas
-----------------------------------------------------------
fst3, snd3, thr :: (a, a, a) -> a
fst3 (x, _, _) = x
snd3 (_, x, _) = x
thr (_, _, x)  = x

genPath :: Int -> FilePath
genPath = printf imagesPath

exitFalha :: String -> IO ()
exitFalha s = do
  a <- putStrLn s
  exitFailure

toFloat :: Int -> Double
toFloat i = fromIntegral i :: Double

status :: Int -> Int -> IO ()
status total n = do
  putStr (printf "\r%.f%%\t- %d/%d Frames" (toFloat n/toFloat total*100) n total)

-----------------------------------------------------------
-- Criação de paletas de cor
-----------------------------------------------------------
filterPallete :: [(Double, Double, Double)] -> [(Pixel8,Pixel8,Pixel8)]
filterPallete = map (\(r,g,b) -> (round r, round g, round b))

laranja, azul, roxo, verde :: [(Pixel8, Pixel8, Pixel8)]
laranja = filterPallete $ zip3 [0..254.0] [0,0.72..184] (replicate 254 0.0)

verde = filterPallete $ zip3 [0,0.46..120] [0..255.0] (replicate 254 0.0)

azul = filterPallete $ zip3 (replicate 254 0.0) (replicate 254 0.0) [0..255.0]

roxo = filterPallete $ zip3 [0..254] (replicate 254 0.0) [0,0.59..150.0]

pallete :: [(Pixel8,Pixel8,Pixel8)]
pallete = azulInOut ++ laranjaInOut ++ verdeInOut ++ roxoInOut
  where laranjaInOut = laranja ++ reverse laranja
        verdeInOut = verde ++ reverse verde
        azulInOut  = azul ++ reverse azul
        roxoInOut  = roxo ++ reverse roxo

palleteR, palleteG, palleteB :: CArray Int Pixel8
(palleteR, palleteG, palleteB) = ( palVec [fst3 x | x <- pallete],
                                   palVec [snd3 x | x <- pallete],
                                   palVec [thr x | x <- pallete]
                                 )

maxi :: Int
maxi = length laranja*2 + length verde*2 + length azul*2 + length roxo*2

-----------------------------------------------------------
-- Cálculo do Mandelbrot
-----------------------------------------------------------
calcPoint :: CComplex -> CComplex -> Int -> (Int, Double)
calcPoint (cx,cy) (zx,zy) iter
  | iter < maxIter && sqrt(zx^2 + zy^2) < 4 = calcPoint (cx,cy) newZ (iter+1)
  | otherwise = (iter, magnitude (zx:+zy) )
  where newZ = (zx^2 - zy^2 + cx , 2*zx*zy + cy)

colorFromIter :: Int -> (Int, Double) -> PixelRGB8
colorFromIter hue iterPoint
  | iter == maxIter = PixelRGB8 0 0 0
  | otherwise = PixelRGB8 (palleteR!i) (palleteG!i) (palleteB!i)
  where (iter,mag) = iterPoint
        i = mod (truncate (color * toFloat maxi)) maxi
        color = toFloat ix / toFloat points -- Entre 0 e 1, representa um ciclo de cores
        ix = truncate (sqrt (toFloat iter + 1 - logBase 2 (logBase 2 mag))*200 + toFloat hue*4 + toFloat hueOffset ) `mod` points
        points = 2048

genIter :: Int -> Int -> Int -> Int -> (Int, Double)
genIter  x y frame db = calcPoint (xPos, yPos) (0,0) 0
  where xPos = (x'-w/2)/size+coordX
        yPos = (y'-h/2)/size+coordY
        (w, h) = (fromIntegral width, fromIntegral height)
        (x', y') = (fromIntegral x, fromIntegral y)
        size = 2** if frame == -1 then fixedZoom else exponent
        exponent = if exp > maxZoom then maxZoom else exp
        exp = fromIntegral (frame + db `div` 160)/(framerate*10)+7

-----------------------------------------------------------
-- Funções de matrizes/vetores
-----------------------------------------------------------
carray :: [Double] -> CArray (Int, Int) Double
carray = listArray ((0, 0), (1024,0))

cMatrix :: [Int] -> CArray (Int, Int) Int
cMatrix = listArray ((0, 0), (width,height))

cMatrixDouble :: [Double] -> CArray (Int, Int) Double
cMatrixDouble = listArray ((0, 0), (width,height))

matrixIter :: CArray (Int, Int) Int
matrixIter = cMatrix [fst $ genIter x y (-1) 0 | x <- [0..width], y <- [0..height]]

matrixMag :: CArray (Int, Int) Double
matrixMag = cMatrixDouble [snd $ genIter x y (-1) 0 | x <- [0..width], y <- [0..height]]

readIter :: Int -> Int -> Int
readIter x y = matrixIter ! (x,y)

readMag :: Int -> Int -> Double
readMag x y = matrixMag ! (x,y)

palVec :: [Pixel8] -> CArray Int Pixel8
palVec = listArray (0,maxi)

-----------------------------------------------------------
-- Funções FFT
-----------------------------------------------------------
cleanComplex :: CArray (Int,Int) (Complex Double) -> Int
cleanComplex c = sum stripped `div` 10
  where stripped = take 7 cleaned
        cleaned = map (\comp -> round (magnitude comp) `div` 100)  $ elems c

getAnimationTimings :: [Double] -> [Int] -> Int -> Int -> [((Int, Int), Int)]
getAnimationTimings samples rangeFrames spf dur = map (\(sz, db) -> ((sz, db), foldl1 (\n x -> n + x `div` 15) $ take sz [snd t | t <- newSt])) newSt
  where st = map (\x -> cleanComplex $ dftRC $ carray $ drop (spf*x+1) samples) rangeFrames
        newSt = map (\sz -> (sz, foldl1 (\n x -> n - (n-x) `div` 4) $ take sz st)) [2..dur * round framerate]

-----------------------------------------------------------
-- Criação da imagem
-----------------------------------------------------------
doAnim :: ((Int, Int),Int) -> Bool -> Maybe (Image PixelRGB8)
doAnim info static = Just $ generateImage genPixel width height
  where genPixel x y = colorFromIter hue $ if static then (readIter x y, readMag x y) else genIter x y frameN hueDB
        h = hueDB `div` (10*sensitivity)
        hue = mod (db `div` (20*sensitivity) + h + frameN `div` 60) maxi
        (path, frameN, db) = (genPath $ fst s, fromIntegral $ fst s, snd s)
        (s, hueDB) = info

doAnimSave :: ((Int, Int),Int) -> FilePath -> Bool -> IO ()
doAnimSave info path static = do

  let image = fromMaybe (generateImage (\x y -> PixelRGB8 0 0 0) 2 2 ) $ doAnim info static

  if imageWidth image < 3 then putStrLn "" else writePng path image

writeVideo :: [((Int, Int), Int)] -> Int -> Int -> Int ->  IO ()
writeVideo dbList spf sr total = do
  let listOfImage = map (`doAnim` estático) dbList ++ [Nothing]

  save <- imageWriter (EncodingParams (fromIntegral width) (fromIntegral height) (round framerate) Nothing Nothing "medium" Nothing) videoPath

  mapM_ (\(n,z) -> status total n >> save z) $ zip [1..] listOfImage

  callCommand $ printf "ffmpeg -i %s -i %s -map 0:v -map 1:a -c:v copy -shortest %s -y" videoPath musicPath ( reverse (drop 4 $ reverse videoPath) ++ "_audio.mp4")

-----------------------------------------------------------
-- Main
-----------------------------------------------------------
main :: IO ()
main = do
  _ <- hSetBuffering stdout NoBuffering

  p <- getWAVEFile musicPath --Lendo nosso arquivo

  let header = waveHeader p

  -- Esse número 32768 é quando estamos lendo word8 signed como unsigned
  -- então o que deveria ser 0 signed se torna 32768 unsigned (10000000 = -32768 signed but 10000000 unsiged = 32768)
  -- o FFT não se importa mt, math is beautiful
  let waveSmpAll = [fromIntegral (head l) /32768 | l <- waveSamples p] --todas nossas samples

  let sampleRate = waveFrameRate header  --Numero de samples e.g 44100

  let waveSmp = drop (sampleRate `div` 9) waveSmpAll -- Isso irá colocar nosso vídeo 100ms adiantado, melhor sincronismo

  let samples = fromMaybe 0 $ waveFrames header     -- Total de samples da nossa música

  _ <- if samples == 0 then exitFalha "Falha na leitura de samples do arquivo." else putStrLn "\nFile OK."

  let duração = samples `div` waveFrameRate header

  let samplesPerFrame = sampleRate `div` round framerate

  let total = duração * round framerate

  let rangeFrames = [0..total]

  _ <- putStrLn $ printf "%dx%d - %.f fps" width height framerate
  _ <- putStrLn $ printf "Sample Count: %d" samples
  _ <- putStrLn $ printf "SampleRate: %d\nduracao em segundos: %d\nNumero de frames: %d" sampleRate duração total

  {- Chamamos essa função para retornar uma lista de informações úteis para
     desenharmos o mandelbrot reagindo a música.
     [((a, b), c)] onde a é o numero do frame
     b é quantidade de grave que o fft disse que tem no frame atual
     c é a versão suavizada de b, objetivo dela é não ter picos extremos -}
  let dbList = drop dropped $ getAnimationTimings waveSmp rangeFrames samplesPerFrame duração

  _ <- putStr "Deseja salvar em vídeo? [Y/N] "
  opt <- getLine

  _ <- if null opt || 'y' /= toLower (head opt) && 'n' /= toLower (head opt)
    then exitFalha "Não entendi. Abortando."
    else putStrLn (if 'y' == toLower (head opt) then "Saída de vídeo" else "Saída em imagens")

  if 'n' == toLower (head opt)
    then mapM_ (\(n, x) -> status total n >> doAnimSave x (genPath n) estático) $ zip [dropped..] dbList

    else writeVideo dbList samplesPerFrame sampleRate total

  putStrLn "\nDone"


text:
{-# LANGUAGE FlexibleContexts #-}
module Statistics.Test.Internal (
    rank
  , rankUnsorted
  , splitByTags
  ) where

import Data.Ord
import           Data.Vector.Generic           ((!))
import qualified Data.Vector.Generic         as G
import qualified Data.Vector.Generic.Mutable as M
import Statistics.Function

-- Private data type for unfolding
data Rank v a = Rank {
      rankCnt :: {-# UNPACK #-} !Int        -- Number of ranks to return
    , rankVal :: {-# UNPACK #-} !Double     -- Rank to return
    , rankNum :: {-# UNPACK #-} !Double     -- Current rank
    , rankVec :: v a                        -- Remaining vector
    }

-- | Calculate rank of every element of sample. In case of ties ranks
--   are averaged. Sample should be already sorted in ascending order.
--
--   Rank is index of element in the sample, numeration starts from 1.
--   In case of ties average of ranks of equal elements is assigned
--   to each
--
-- >>> rank (==) (fromList [10,20,30::Int])
-- > fromList [1.0,2.0,3.0]
--
-- >>> rank (==) (fromList [10,10,10,30::Int])
-- > fromList [2.0,2.0,2.0,4.0]
rank :: (G.Vector v a, G.Vector v Double)
     => (a -> a -> Bool)        -- ^ Equivalence relation
     -> v a                     -- ^ Vector to rank
     -> v Double
rank eq vec = G.unfoldr go (Rank 0 (-1) 1 vec)
  where
    go (Rank 0 _ r v)
      | G.null v  = Nothing
      | otherwise =
          case G.length h of
            1 -> Just (r, Rank 0 0 (r+1) rest)
            n -> go Rank { rankCnt = n
                         , rankVal = 0.5 * (r*2 + fromIntegral (n-1))
                         , rankNum = r + fromIntegral n
                         , rankVec = rest
                         }
          where
            (h,rest) = G.span (eq $ G.head v) v
    go (Rank n val r v) = Just (val, Rank (n-1) val r v)
{-# INLINE rank #-}

-- | Compute rank of every element of vector. Unlike rank it doesn't
--   require sample to be sorted.
rankUnsorted :: ( Ord a
                , G.Vector v a
                , G.Vector v Int
                , G.Vector v Double
                , G.Vector v (Int, a)
                )
             => v a
             -> v Double
rankUnsorted xs = G.create $ do
    -- Put ranks into their original positions
    -- NOTE: backpermute will do wrong thing
    vec <- M.new n
    for 0 n $ \i ->
      M.unsafeWrite vec (index ! i) (ranks ! i)
    return vec
  where
    n = G.length xs
    -- Calculate ranks for sorted array
    ranks = rank (==) sorted
    -- Sort vector and retain original indices of elements
    (index, sorted)
      = G.unzip
      $ sortBy (comparing snd)
      $ indexed xs
{-# INLINE rankUnsorted #-}

-- | Split tagged vector
splitByTags :: (G.Vector v a, G.Vector v (Bool,a)) => v (Bool,a) -> (v a, v a)
splitByTags vs = (G.map snd a, G.map snd b)
  where
    (a,b) = G.unstablePartition fst vs
{-# INLINE splitByTags #-}