text:
{-# language TypeFamilies #-}
{-# language FlexibleInstances #-}
{-# language FlexibleContexts #-}
{-# language UndecidableInstances #-}
{-# language DerivingVia #-}
{-# language GeneralizedNewtypeDeriving #-}

module Geometry.Polygon
  (
    PiecewiseLinear (..)
  , isDegenerate, isNondegenerate
  , Polyline (..)
  , mkPolyline, trivialPolyline
  , Polygonal, Polygon (..)
  , innerAngles
  , mkPolygon, trivialPolygon, closePolyline
  , Triangle (..)
  , mkTriangle
  , Rectangle (..)
  , mkRectangle
  , boxRectangle
  , RightTriangle (..), isRightTriangle, hypotenuse, catets
  )
where

import Data.Complex
import Data.Foldable
import Data.List.Extra (minimumOn)
import Data.Monoid
import Data.Maybe
import Data.Fixed

import Geometry.Base
import Geometry.Point
import Geometry.Line

-- | A class for polylines and polygons.
class (Trans p, Manifold p) => PiecewiseLinear p where
  {-# MINIMAL vertices | asPolyline  #-}
  -- | A list of polyline vertices.
  vertices :: p -> [Cmp]
  vertices = vertices . asPolyline

  -- | A representation of the instance as a polyline.
  asPolyline :: p -> Polyline
  asPolyline = Polyline . vertices

  -- | The number of vertices of a polyline.
  verticesNumber :: p -> Int
  verticesNumber p = length (vertices p)

  -- | The list of segments of a polyline.
  segments :: p -> [Segment]
  segments p = case vertices p of
    [] -> []
    [x] -> [Segment (x,x)]
    vs -> Segment <$> zip vs (tail vs)

  -- | The indexed selector for vertices of a polyline.
  vertex :: Int -> p -> Cmp
  vertex i p = vs !! j
    where vs = vertices p
          n = verticesNumber p
          j = (0 `max` i) `min` (n - 1)

  -- | The indexed selector for sides of a polyline.
  side :: Int -> p -> Segment
  side i p = segments p !! j
    where j = (0 `max` i) `min` (n - 2)
          n = verticesNumber p

  -- | The list of angles of a polyline.
  vertexAngles :: p -> [Direction]
  vertexAngles p = case segments p of
    [] -> []
    [x] -> []
    s -> zipWith (\a b -> angle a - angle b) s (tail s)

innerAngles t = flip <$> vertexAngles t
  where flip a | deg a > 180 = 360 - a
               | otherwise = a

instance PiecewiseLinear p => PiecewiseLinear (Maybe p) where
  vertices = maybe mempty vertices
  asPolyline = maybe mempty asPolyline
  verticesNumber = maybe 0 verticesNumber
  segments = maybe mempty segments
  vertex i = maybe 0 (vertex i)
  side i = maybe (asCmp 0) (side i)
  vertexAngles = maybe mempty vertexAngles

-- | A predicate. Returns `True` if any of polyline's segment has zero length.
isDegenerate :: PiecewiseLinear p => p -> Bool
isDegenerate p = any isTrivial (segments  p) ||
                 any (0 ~=) (vertexAngles  p)

-- | Returns @True@ for non degenerate polyline or polygon.
isNondegenerate :: PiecewiseLinear p => p -> Bool
isNondegenerate = not . isDegenerate

interpolation :: PiecewiseLinear p => p -> Double -> Maybe Cmp
interpolation p x = param' <$> find interval tbl
  where
    interval ((a, b), _) = a <= x && x <= b
    param' ((a, b), s) = s @-> ((x - a)/(b-a))
    tbl = zip (zip ds (tail ds)) $ segments p
    ds = scanl (+) 0 $ unit <$> segments p

-- | Instantiates as a single segment polyline.
instance PiecewiseLinear Segment where
  vertices s = let (p1,p2) = refPoints s in [p1,p2]

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

-- | Representation of a polygonal chain as a list of vertices.
-- Vector representation of a polyline is given by the first segment of a polyline.
-- It is used for alignment or creation of polylines.
--
-- >>> cmp $ mkTriangle [0, 1:+0, 0:+1]
-- 1.0 :+ 0.0
--
-- >>> asCmp (1 :+ 2) :: Triangle
-- <Triangle (0.0,0.0) (1.0,2.0) (-1.2320508075688772,1.8660254037844393)>
--
-- >>> asCmp 1 :: Rectangle
-- <Rectangle (0.0,0.0) (1.0,0.0) (1.0,1.0) (0.0,1.0)>
--
newtype Polyline = Polyline [Cmp]

instance PiecewiseLinear Polyline where
  vertices (Polyline vs) = vs
  asPolyline = id

-- | Creates an empty polyline.
trivialPolyline :: Polyline
trivialPolyline = mempty

-- | The main polyline constructor.
mkPolyline :: Affine a => [a] -> Polyline
mkPolyline [] = Polyline [0, 0]
mkPolyline [x] = Polyline [cmp x, cmp x]
mkPolyline vs = Polyline $ cmp <$> vs

instance Show Polyline where
  show p = concat ["<Polyline ", n, ">"]
    where vs = vertices p
          n = if length vs < 5
              then unwords $ show . xy <$> vs
              else "-" <> show (length vs) <> "-"

instance Eq Polyline where
  p1 == p2 = vertices p1 ~= vertices p2

-- | Addition of polylines, excluding zero segments.
instance Semigroup Polyline where
  p1 <> p2 = Polyline vs
    where
      vs | null v1 = v2
         | null v2 = v1
         | otherwise = v1 <> m v2
        where
          v1 = vertices p1
          v2 = vertices p2
          m = if last v1 == head v2 then tail else id

instance Monoid Polyline where
  mempty = Polyline []

instance Affine Polyline where
  cmp p = case segments p of
            [] -> 0
            (x:_) -> cmp x
  asCmp x = mkPolyline [0, x]

instance Trans Polyline where
  transform t p = Polyline $ transform t <$> vertices p

instance Manifold Polyline where
  type Domain Polyline = Cmp
  param p t | t < 0 = param (asLine (side 0 p)) t
            | t > 1 = param (asLine (side (verticesNumber p - 1) p)) t
            | otherwise = fromJust $  interpolation p (t * unit p)

  project p pt = (x0 + (project s pt * unit s)) / unit p
    where
      ss = segments p
      ds = scanl (+) 0 $ unit <$> ss
      (x0, s) = minimumOn (\(_,s) -> distanceTo pt s) $ zip ds ss
      x = (x0 + (project s pt * unit s)) / unit p

  isContaining p x = any (`isContaining` x) (segments p)
  unit p = sum $ unit <$> segments p

instance Curve Polyline where
  tangent p t =  asCmp $ (p @-> (t + dt)) - (p @-> (t - dt))
    where dt = 1e-5
  normal p t = tangent p t # rotate (-90)

instance Figure Polyline where
  isTrivial p = length (vertices p) < 2
  refPoint p = if isNontrivial p
               then head $ vertices p
               else 0
  box p = foldMap (box . mkPoint) (vertices p)

------------------------------------------------------------
-- | A class for distinguishing `Polygon`s within `PiecewiseLinear` objects.
class PiecewiseLinear p => Polygonal p where

-- | Representation of a closed polygon as a list of vertices.
newtype Polygon = Polygon [Cmp]
  deriving ( Eq
           , Trans
           , Affine
           , Figure
           ) via Polyline

instance Polygonal Polygon where

-- | Creates an empty polyline.
trivialPolygon :: Polygon
trivialPolygon = Polygon []

-- | The main polyline constructor.
mkPolygon :: Affine a => [a] -> Polygon
mkPolygon = Polygon . fmap cmp

-- | Turns a polyline to a polygon with the same vertices.
closePolyline :: Polyline -> Polygon
closePolyline p = Polygon $
              if last vs == head vs then (init vs) else vs
  where vs = vertices p

instance PiecewiseLinear Polygon where
  asPolyline p = Polyline $ take (length vs + 1) (cycle vs)
    where vs = vertices p
  vertices (Polygon vs) = vs
  verticesNumber p = length (vertices p)
  segments = segments . asPolyline

  vertex i p = vs !! (i `mod` n)
    where vs = vertices (asPolyline p)
          n = verticesNumber p

  side i p = segments p !! (i `mod` n)
    where n = verticesNumber p

  vertexAngles p = case segments p of
    [] -> []
    [x] -> [0,0]
    s -> zipWith (\a b -> 180 - (angle b - angle a)) (last s : s) s
    where inner a | a > 180 = 360 - a
                  | a < 0 = 360 + a
                  | otherwise = a

instance Show Polygon where
  show p = concat ["<Polygon ", n, ">"]
    where vs = vertices p
          n = if length vs < 5
              then unwords $ show . xy <$> vs
              else "-" <> show (length vs) <> "-"

instance Manifold Polygon where
  type Domain Polygon = Cmp
  param p t = fromJust $ interpolation (asPolyline p) $ (t `mod'` 1) * unit p
  project = project . asPolyline
  isContaining = isContaining . asPolyline
  unit = unit . asPolyline

instance Curve Polygon where
  tangent p t =  asCmp $ (p @-> (t + dt)) - (p @-> (t - dt))
    where dt = 1e-5
  normal p t = tangent p t # rotate (-90)

instance ClosedCurve Polygon where
  location p pt = case foldMap go (segments p') of
                    (Any True, _) -> OnCurve
                    (_, Sum n) | odd n -> Inside
                    _ -> Outside
    where
      p' = p # translate' (negate (cmp pt))
      go s | y0 * y1 == 0       = (Any True, mempty)
           | y0 == y1           = mempty
           | x == 0             = (Any True, mempty)
           | x > 0 && y0*y1 < 0 = (mempty, Sum 1)
           | otherwise          = mempty
        where
          (x0:+y0, x1:+y1) = refPoints s
          x = (x0*y1-x1*y0)/(y1-y0)

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

-- | Representation of a triangle as a list of vertices.
newtype Triangle = Triangle [Cmp]
  deriving ( Figure
           , Manifold
           , Curve
           , ClosedCurve
           , Trans
           , Eq
           , PiecewiseLinear
           , Polygonal
           ) via Polygon

{- | The main triangle constructor.

-}
mkTriangle :: Affine a => a -> a -> a -> Triangle
mkTriangle a b c = Triangle [cmp a, cmp b, cmp c]

instance Show Triangle where
  show t = concat ["<Triangle ", ss, ">"]
    where ss = unwords $ show . xy <$> vertices t

instance Affine Triangle where
  cmp = cmp . asPolyline
  asCmp = Triangle . scanl (+) 0 . take 2 . iterate (rotate 120)

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

-- | Representation of a rectangle as a list of vertices.
newtype Rectangle = Rectangle [Cmp]
  deriving ( Figure
           , Manifold
           , Curve
           , ClosedCurve
           , Trans
           , Eq
           , PiecewiseLinear
           , Polygonal
           ) via Polygon

-- | The main rectangle constructor, uses two side lengths and guarantees
-- that polygon will have right angles.
mkRectangle :: Double -> Double -> Rectangle
mkRectangle a b = (asCmp 1 :: Rectangle) # scaleX a # scaleY b

instance Show Rectangle where
  show t = concat ["<Rectangle ", ss, ">"]
    where ss = unwords $ show . xy <$> vertices t

instance Affine Rectangle where
  cmp = cmp . asPolyline
  asCmp = Rectangle . scanl (+) 0 . take 3 . iterate (rotate 90)

{- | Returns a figures' box as a rectangle.

> let fs = scanl1 (<+>) [ G $ segment (-1,1) (1,2)
>                       , G $ point (1, 3)
>                       , G $ circle 1 (3, 4)
>                       , G $ aTriangle ]
>     withBox f = f <+> boxRectangle f #: thin <> dotted
> in rowSep (space 1) $ withBox <$> fs

<< figs/boxRectangle.svg >>
-}
boxRectangle :: Figure f => f -> Rectangle
boxRectangle f = Rectangle [ p4, p3, p2, p1 ]
  where ((p4,p3),(p1,p2)) = corner f

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

{- | A wrapper for the right triangles.

>>> aRightTriangle # scaleX 2
RightTriangle <Triangle (0.0,0.0) (2.0,0.0) (1.0,1.0)>
-}
newtype RightTriangle = RightTriangle Triangle
  deriving ( Eq
           , Show
           , Affine
           , Trans
           , Figure
           , Manifold
           , Curve
           , ClosedCurve
           , Polygonal
           , PiecewiseLinear)

{- | Returns `True` if polygon is a right triangle. -}
isRightTriangle :: PiecewiseLinear p => p -> Bool
isRightTriangle t = any (90 ~=) $ vertexAngles t

{- | Returns a hypotenuse for a rigth triangle.

>>> hypotenuse aRightTriangle
Segment (1.0 :+ 0.0, 1.0 :+ 0.0)

-}
hypotenuse :: RightTriangle -> Segment
hypotenuse (RightTriangle t) = side 1 t

{- | Returns a pair of catets for a rigth triangle.

>>> catets aRightTriangle
(Segment (0.0 :+ 0.0, 1.0 :+ 0.0),Segment (0.0 :+ 1.0, 0.0 :+ 0.0))

-}
catets :: RightTriangle -> (Segment, Segment)
catets (RightTriangle t) = (side 0 t, side 2 t)

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
{-# LANGUAGE BangPatterns, CPP #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleContexts #-}

module Mandelbrot(mandelBench) where
{- Code taken from: https://github.com/simonmar/monad-par/blob/master/examples/src/mandel/mandel.hs-}

import Control.Monad
import Control.DeepSeq
import Control.Exception
import Data.Complex
import System.Environment
import System.IO
#ifdef ALIST
import Control.Monad.Par.AList as A
#endif

#ifdef PARSCHED
import PARSCHED
#else
import Control.Monad.Par
#endif

import qualified Data.Vector.Unboxed as V

#ifdef WRITE_IMAGE
import Codec.Picture  -- JuicyPixels
import qualified Data.Vector.Storable as V
#endif

#ifdef NEW_GENERIC
import qualified Data.Par as C
import qualified Data.Par.Range as C
import Control.Par.Class
parMapReduceRangeThresh :: (NFData a, ParFuture p, FutContents p a)
                        => Int -> C.InclusiveRange -> (Int -> p a) -> (a -> a -> p a) -> a -> p a
parMapReduceRangeThresh = C.parMapReduceThresh
#else
import Control.Monad.Par.Combinator as C
#endif

import Control.Monad.SD as Ohua
import Criterion
import Criterion.Main
import Control.Monad.Identity (runIdentity)

mandel :: Int -> Complex Double -> Int
mandel max_depth c = loop 0 0
  where
   fn = magnitude
   loop i !z
    | i == max_depth = i
    | fn(z) >= 2.0   = i
    | otherwise      = loop (i+1) (z*z + c)

threshold = 1

data VecTree = Leaf (V.Vector Int)
	     | MkNode VecTree VecTree
  deriving (Eq,Show)

#if 0
instance V.Unbox a => NFData (V.Vector a) where
  rnf v = rnf (V.length v)
#endif

instance NFData VecTree where
  rnf (Leaf v) = rnf v
  rnf (MkNode v1 v2) = rnf v1 `seq` rnf v2

leftmost (Leaf v)     = v
leftmost (MkNode l _) = leftmost l

#ifdef ALIST
runMandel :: Double -> Double -> Double -> Double -> Int -> Int -> Int -> Par (AList [Int])
runMandel minX minY maxX maxY winX winY max_depth = do
  -- This version does a ROW at a time in parallel:
  A.parBuildThreshM threshold (C.InclusiveRange 0 (winY-1)) $ \y ->
       do let l = [ mandelStep y x | x <- [0.. winX-1] ]
          deepseq l (return l)

#else
-- RRN: Comment to handle LVish types too:
-- runMandel :: Double -> Double -> Double -> Double -> Int -> Int -> Int -> Par VecTree
runMandel minX minY maxX maxY winX winY max_depth runner = do
  -- Auto-partitioning version.  A bit worse:
  -- C.parMapReduceRange (C.InclusiveRange 0 (winY-1))
  runner threshold (C.InclusiveRange 0 (winY-1))
     (\y ->
       do
          let vec = V.generate winX (\x -> mandelStep y x)
          seq (vec V.! 0) $
           return (Leaf vec))
     (\ a b -> return$ MkNode a b)
     (Leaf V.empty)
#endif
  where
    mandelStep i j = mandel max_depth z
        where z = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
                  ((fromIntegral i * c_scale) / fromIntegral winX + minX)
    r_scale  =  maxY - minY  :: Double
    c_scale =   maxX - minX  :: Double

-- pure version of the above
runMandelP minX minY maxX maxY winX winY max_depth runner =
  -- Auto-partitioning version.  A bit worse:
  -- C.parMapReduceRange (C.InclusiveRange 0 (winY-1))
  runner threshold (C.InclusiveRange 0 (winY-1))
     (\y ->
        let vec = V.generate winX (\x -> mandelStep y x)
        in seq (vec V.! 0) $ (Leaf vec))
     (\ a b ->  MkNode a b)
     (Leaf V.empty)
  where
    mandelStep i j = mandel max_depth z
        where z = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
                  ((fromIntegral i * c_scale) / fromIntegral winX + minX)
    r_scale  =  maxY - minY  :: Double
    c_scale =   maxX - minX  :: Double

-- runMandelP_ minX minY maxX maxY winX winY max_depth runner =
--   -- Auto-partitioning version.  A bit worse:
--   -- C.parMapReduceRange (C.InclusiveRange 0 (winY-1))
--   runner threshold (C.InclusiveRange 0 (winY-1))
--      (\y ->
--         let vec = V.generate winX (\x -> mandelStep y x)
--         in seq (vec V.! 0) $ (Leaf vec))
--      (\ a b ->  MkNode a b)
--      (Leaf V.empty)
--   where
--     mandelStep i j = mandel max_depth z
--         where z = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
--                   ((fromIntegral i * c_scale) / fromIntegral winX + minX)
--     r_scale  =  maxY - minY  :: Double
--     c_scale =   maxX - minX  :: Double

#ifdef WRITE_IMAGE
-- makeImage :: Integer -> Integer -> Int -> A.AList [Int] -> Image PixelRGB8
makeImage :: Int -> Int -> Int -> A.AList [Int] -> Image PixelRGB8
makeImage x y depth ls =
  -- This is a quite silly series of intermediate structures:
  Image x y (V.fromList $ concat $
	     Prelude.map prettyRGB $
	     concat$ A.toList ls)
 where
   prettyRGB s =
      let t = (fromIntegral (depth - s)) in
--      PixelRGB8 (fromIntegral s) t t
      [fromIntegral s, t, t]
#endif

simple x y depth = runMandel (-2) (-2) 2 2 x y depth
simpleP :: Monad m => Int -> Int -> Int -> (
    Int -- ^ threshold
    -> InclusiveRange -- ^ range over which to calculate
    -> (Int -> VecTree) -- ^ compute one result
    -> (VecTree -> VecTree -> VecTree) -- ^ combine two results (associative)
    -> VecTree -- ^ initial result
    -> m VecTree) -> m VecTree
simpleP x y depth = runMandelP (-2) (-2) 2 2 x y depth

-- simpleP_ :: Int -> Int -> Int -> (Int -- ^ threshold
--     -> InclusiveRange -- ^ range over which to calculate
--     -> (Int -> VecTree) -- ^ compute one result
--     -> (VecTree -> VecTree -> VecTree) -- ^ combine two results (associative)
--     -> VecTree -- ^ initial result
--     -> IO VecTree) -> IO VecTree
-- simpleP_ x y depth = runMandelP_ (-2) (-2) 2 2 x y depth

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

#ifdef ALIST
-- A meaningless checksum.  This match the C++ CnC benchmark:
mandelCheck :: AList [Int] -> Int -> Int -> Int
mandelCheck als max_col max_depth = loop 0 als 0
 where
 -- Sum positions of the pixels that are equal to max_depth:
 loop i als !sum | A.null als = sum
 loop i als !sum = loop (i+1) (A.tail als)
		        (loop2 i 0 (A.head als) sum)
 loop2 i j []    !sum = sum
 loop2 i j (h:t) !sum | h == max_depth = loop2 i (j+1) t (sum + i*max_col + j)
		      | otherwise      = loop2 i (j+1) t  sum
#endif

 -- This kind of checksum is much simpler:
checkSum :: VecTree -> Int
checkSum (Leaf vec) = V.foldl (+) 0 vec
checkSum (MkNode v1 v2) = checkSum v1 + checkSum v2

-- main = do args <- getArgs
--
--           let (x,y,depth) =
-- 		case args of
-- 		 [] ->
--     --                 runPar $ simple 3 3 3
--                        -- Defaults recommended by Simon Marlow
-- 		       -- (1024,1024,256)
--                        -- Note, this should do something very *quick* when called with no args:
--                        (3,3,3)
--
-- 		 [x,y,depth] ->
--     --		   simple (read x) (read y) (read depth)
-- 		       (read x, read y, read depth)
--
-- 		 -- [minX,minY,maxX,maxY,winX,winY,depth] ->
-- 		 --    runPar $
-- 		 -- 	    runMandel (read minX) (read minY)
-- 		 -- 		      (read maxX) (read maxY)
-- 		 -- 		      (read winX) (read winY) (read depth)
--
--           let pixels = runPar$ simple x y depth parMapReduceRangeThresh
--
--           putStrLn$ "About to poke pixels to launch runPar:"
--           evaluate pixels
--           putStrLn$ "Should be done with runpar at this point."
--
-- #ifdef WRITE_IMAGE
-- 	  writePng "mandel_image.png" (makeImage (fromIntegral x) (fromIntegral y) depth pixels)
-- 	  putStrLn$ "File written."
--           putStrLn$ "Checksum " ++ show (checkSumpixels y depth)
-- #endif
--           putStrLn$ "Spot check, length of leftmost leaf: " ++ show (V.length (leftmost pixels))
--           return ()

parMapReduceRangeThreshSeq:: Int -- ^ threshold
  -> InclusiveRange -- ^ range over which to calculate
  -> (Int -> VecTree) -- ^ compute one result
  -> (VecTree -> VecTree -> VecTree) -- ^ combine two results (associative)
  -> VecTree -- ^ initial result
  -> IO VecTree
parMapReduceRangeThreshSeq t (InclusiveRange mi ma) fn binop init =
  return $ loop mi ma
  where
    loop mi ma
      | ma - mi <= t =
        let mapred a b = let x = fn b
                             result = a `binop` x
                          in result
         in foldl mapred init [mi .. ma]
      | otherwise =
        let mid = mi + ((ma - mi) `quot` 2)
            r = loop (mid + 1) ma
            l = loop mi mid
         in l `binop` r

mandelBench =
  bgroup
    "mandel-bench"
    [ bench "sequential" (nfIO (simpleP x y depth parMapReduceRangeThreshSeq))
    , bench "ohua" (nfIO $ simpleP x y depth Ohua.mapReduceRangeThresh)
    , bench "par" (nf (runPar . (simple x y depth)) C.parMapReduceRangeThresh)
    ]
  where
    {- Values taken from: https://github.com/simonmar/monad-par/blob/master/examples/src/run_benchmark.hs-}
    -- desktop values
    -- x = 1024
    -- y = 1024
    -- depth = 256
    -- server values
    x = 1024
    y = 1024
    depth = 256
{-
[2012.03.08] {Updated to use Vector}
Looks like this was out of date in that it was still using AList.
AList has never performed well enough, and sure enough, this test was
falling to ~50% productivity at >16 threads.
However, simply replacing AList with vector is no good (14.4s for 512
512 2048 as opposed to 15.5).  It is of course a quadratic copy at
that point.
-}

ext: hs

max_stars_repo_name: ohua-dev/statefulness

max_issues_repo_name: ohua-dev/statefulness

max_forks_repo_name: ohua-dev/statefulness


text:
-- | Module providing various functions popular in quantum calculations, but
-- not implemented in "Numeric.LinearAlgebra".
module Hoqus.MtxFun where

import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra

-- | Trace of a matrix. This function ignores the non-sqare part of the matrix.
trace :: Matrix C -> C
trace = sum.toList.takeDiag

-- | Logarithm of a matrix. Based on 'Numeric.LinearAlgebra.matFunc' function
-- from "Numeric.LinearAlgebra" package.
logm :: Matrix C -> Matrix C
logm = matFunc log

max_stars_repo_name: jmiszczak/hoqus

max_issues_repo_name: jmiszczak/hoqus

max_forks_repo_name: jmiszczak/hoqus


text:
module Numeric.Network where
import Data.Kind
import Dependent.Size
import Numeric.Layer
import Numeric.Vector.Sized
import Numeric.Loss

infixr 0 :~>

type family NetworkOutput (layers :: [Type]) where
  NetworkOutput '[x] = Outputs x
  NetworkOutput (x ': y ': rest) = NetworkOutput (y ': rest)

type family NetworkInput (layers :: [Type]) where
  NetworkInput (x ': _) = Inputs x

data Network (layers :: [Type]) where
  Last :: Layer x => x -> Network '[x]
  (:~>) :: (Layer x, Outputs x ~ Inputs y) => x -> Network (y:xs) -> Network (x : y : xs)

instance (Show x) => Show (Network '[x]) where
  showsPrec n (Last g) = showsPrec n g
instance (Show x, Show (Network (y:xs))) => Show (Network (x:y:xs)) where
  showsPrec n (layer :~> net) = showsPrec n layer . showString " :~> " . showsPrec n net

infixr 0 :|>
data Tapes (layers :: [Type]) where
  LastTape :: Layer x => Tape x -> Tapes '[x]
  (:|>) :: (Layer x, Outputs x ~ Inputs y) => Tape x -> Tapes (y:xs) -> Tapes (x:y:xs)

instance (Show (Tape x)) => Show (Tapes '[x]) where
  showsPrec n (LastTape g) = showsPrec n g
instance (Show (Tape x), Show (Tapes (y:xs))) => Show (Tapes (x:y:xs)) where
  showsPrec n (tape :|> tapes) = showsPrec n tape . showString " :|> " . showsPrec n tapes

infixr 0 :<|
data Gradients (layers :: [Type]) where
  LastGradient :: Layer x => Gradient x -> Gradients '[x]
  (:<|) :: (Layer x, Outputs x ~ Inputs y) => Gradient x -> Gradients (y:xs) -> Gradients (x : y : xs)

instance (Show (Gradient x)) => Show (Gradients '[x]) where
  showsPrec n (LastGradient g) = showsPrec n g
instance (Show (Gradient x), Show (Gradients (y:xs))) => Show (Gradients (x:y:xs)) where
  showsPrec n (grad :<| grads) = showsPrec n grad . showString " :<| " . showsPrec n grads

predict
  :: Network layers
  -> SizedArray (NetworkInput layers)
  -> SizedArray (NetworkOutput layers)
predict (Last layer ) x = fst (forward layer x)
predict (layer:~>net) x = predict net y where y = fst (forward layer x)

predict'
  :: Network layers
  -> SizedArray (NetworkInput layers)
  -> (SizedArray (NetworkOutput layers), Tapes layers)
predict' (Last layer) x = let (y, tape) = forward layer x in (y, LastTape tape)
predict' (layer:~>net) x =
  let (y, tape) = forward layer x
  in  let (z, tapes) = predict' net y in (z, tape :|> tapes)

backprop
  :: Network layers
  -> Tapes layers
  -> SizedArray (NetworkOutput layers)
  -> (Gradients layers, SizedArray (NetworkInput layers))
backprop (Last layer) (LastTape tape) dy =
  let (grad, dx) = backward layer tape dy in (LastGradient grad, dx)
backprop (layer:~>net) (tape:|>tapes) dz =
  let (grads, dy) = backprop net tapes dz
  in  let (grad, dx) = backward layer tape dy in (grad :<| grads, dx)

backpropagation
  :: (Loss (NetworkOutput layers) loss)
  => Network layers
  -> loss
  -> SizedArray (NetworkInput layers)
  -> SizedArray (NetworkOutput layers)
  -> (Gradients layers, SizedArray (NetworkInput layers))
backpropagation net loss x target = backprop net tapes dy
  where
    (y, tapes) = predict' net x
    dy         = lossDerivative loss y target

learn :: Network layers -> Gradients layers -> Network layers
learn (Last layer ) (LastGradient grad) = Last (applyGradient layer grad)
learn (layer:~>net) (grad:<|grads     ) = applyGradient layer grad :~> learn net grads

train
  :: (Loss (NetworkOutput layers) loss)
  => Network layers
  -> loss
  -> SizedArray (NetworkInput layers)
  -> SizedArray (NetworkOutput layers)
  -> Network layers
train net loss x t = f $ backpropagation net loss x t
  where f (gradients, _) = learn net gradients

-- A network is a layer
instance Layer (Network layers) where
  type Inputs (Network layers) = NetworkInput layers
  type Outputs (Network layers) = NetworkOutput layers
  type Gradient (Network layers) = Gradients layers
  type Tape (Network layers) = Tapes layers
  forward = predict'
  backward = backprop
  applyGradient = learn

max_stars_repo_name: mixed-signals/mixed-signals

max_issues_repo_name: mixed-signals/mixed-signals

max_forks_repo_name: mixed-signals/mixed-signals


text:
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}
module Simulation where

import AmpOp
    ( SignalState(..),
      State(TwoNodes),
      Output,
      Value,
      dc100,
      dc12,
      sen,
      ground,
      getSignalOutput,
      lm741,
      r1,
      r2,
      r3,
      r4,
      -- ampOpSpecial,
      ampOpBuffer,
      ampOpInverting,
      ampOpInvertingTest,
      ampOpInverting' )
import Circuit
import Signal
import Control.Monad.Fix
import Data.Complex

executeTimeSimulation :: SignalState -> (SignalState -> Either String (Circuit a SignalState)) -> Signal a -> [Metric] -> Either String [(Time, Value)]
executeTimeSimulation _ _ _ [] = Right []
executeTimeSimulation initialState buildCircuit input (t:ts) = buildCircuit initialState >>= \circuit -> do let result = circuit `simulate` input `at` t
                                                                                rest <- executeTimeSimulation result buildCircuit input ts
                                                                                Right $ (realPart t, realPart (getSignalOutput result `at` t)) : rest

executeTimeSimulation2 :: SignalState -> (SignalState -> Circuit a SignalState) -> Signal a -> [Metric] -> [(Metric, Output)]
executeTimeSimulation2 initialState buildCircuit input samples = fix calculate (initialState, [], samples)
    where calculate f (state, l, t:ts) = if null ts then next else f (newState, next, ts)
            where newState = buildCircuit state `simulate` input `at` t
                  next = l ++ [result]
                  result = (t, getSignalOutput newState `at` t)

iSignalState = SignalTwoNodes (dc100, ground)
iSignalState3 = SignalThreeNodes (dc100, dc100, dc100)
iSignalBuffer = SignalOneNode dc100
c0 = ampOpBuffer lm741
c1 = ampOpInverting lm741 r1 r2
c3 = ampOpInvertingTest lm741 r1 r2
-- cS = ampOpSpecial lm741 r1 r2 r3 r4

simulationSignal = executeTimeSimulation iSignalState c1 sen (map (:+ 0) [1..5])
simulationSignalTest = executeTimeSimulation2 iSignalState c3 sen (map (:+ 0) [1..5])
-- simulationSpecial = executeTimeSimulation iSignalState3 cS dc12 (map (:+ 0) [1])

iState = TwoNodes (100, 0)
c2 = ampOpInverting' lm741 r1 r2 12

simulation = c2 iState

max_stars_repo_name: FP-Modeling/fixingAnalog

max_issues_repo_name: FP-Modeling/fixingAnalog

max_forks_repo_name: FP-Modeling/fixingAnalog


text:
{-# LANGUAGE CPP                    #-}
{-# LANGUAGE ConstraintKinds        #-}
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE DefaultSignatures      #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE ParallelListComp       #-}
{-# LANGUAGE PartialTypeSignatures  #-}
{-# LANGUAGE PatternGuards          #-}
{-# LANGUAGE Rank2Types             #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TemplateHaskell        #-}
{-# LANGUAGE TypeOperators          #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeFamilies           #-}
{-# LANGUAGE TypeSynonymInstances   #-}
{-# LANGUAGE UndecidableInstances   #-} -- See below

-- experiment
-- {-# LANGUAGE PartialTypeSignatures #-}

{-# OPTIONS_GHC -Wall #-}
-- {-# OPTIONS_GHC -fno-warn-unused-imports #-} -- TEMP

-- #define TESTING

#ifdef TESTING
{-# OPTIONS_GHC -fno-warn-unused-binds   #-} -- TEMP
#endif

-- #define GenericPowFFT

----------------------------------------------------------------------
-- | Generic FFT
----------------------------------------------------------------------

module ConCat.FFT
  ( dft, FFT(..), DFTTy, genericFft, GFFT
  -- -- Temporary while debugging
  , twiddle, twiddles, omega, cis
  , o8sq
  ) where

import Prelude hiding (zipWith)

import Data.Complex
-- import Control.Applicative (liftA2)
import GHC.Generics hiding (C,S)

import Data.Key (Zip(..))
import Data.Pointed

#ifdef TESTING
import Data.Foldable (toList)
-- import Test.QuickCheck (quickCheck)
import Test.QuickCheck.All (quickCheckAll)
import ShapedTypes.ApproxEq
#endif

import ConCat.Misc (Unop,transpose,C,inGeneric1,inComp)
import ConCat.Sized
import ConCat.Scan (LScan,lproducts)
import ConCat.Pair
import ConCat.Free.LinearRow (($*))

{--------------------------------------------------------------------
    DFT
--------------------------------------------------------------------}

type AS  h = (Pointed h, Zip h, LScan h)
type ASZ h = (AS h, Sized h)

-- To resolve: Zip vs Applicative. traverse/transpose needs Applicative.

dft :: forall f a. (ASZ f, Foldable f, RealFloat a) => Unop (f (Complex a))
dft xs = omegas (size @f) $* xs
{-# INLINE dft #-}

omegas :: (AS f, AS g, RealFloat a) => Int -> g (f (Complex a))
omegas = fmap powers . powers . omega
-- omegas n = powers <$> powers (omega n)
-- omegas n = powers <$> powers (exp (- i * 2 * pi / fromIntegral n))
{-# INLINE omegas #-}

omegas' :: forall f g a. (ASZ f, ASZ g, RealFloat a) => g (f (Complex a))
omegas' = fmap powers (powers (cis (- 2 * pi / fromIntegral (size @(g :.: f)))))

-- i :: Num a => Complex a
-- i = 0 :+ 1

omega :: RealFloat a => Int -> Complex a
omega n = cis (- 2 * pi / fromIntegral n)
-- omega n = exp (0 :+ (- 2 * pi / fromIntegral n))
-- omega n = exp (- 2 * (0:+1) * pi / fromIntegral n)
{-# INLINE omega #-}

{--------------------------------------------------------------------
    FFT
--------------------------------------------------------------------}

type DFTTy f = forall a. RealFloat a => f (Complex a) -> FFO f (Complex a)

class FFT f where
  type FFO f :: * -> *
  fft :: DFTTy f
  -- default fft :: ( Generic1 f, Generic1 (FFO f)
  --                , FFT (Rep1 f), FFO (Rep1 f) ~ Rep1 (FFO f) ) => DFTTy f
  -- fft = genericFft
  -- Temporary hack to avoid newtype-like representation.
  fftDummy :: f a
  fftDummy = undefined

-- TODO: Eliminate FFO, in favor of fft :: Unop (f (Complex a)).
-- Use dft as spec.

twiddle :: forall g f a. (ASZ g, ASZ f, RealFloat a) => Unop (g (f (Complex a)))
-- twiddle = (zipWith.zipWith) (*) twiddles
twiddle = (zipWith.zipWith) (*) omegas'
-- twiddle = (zipWith.zipWith) (*) (omegas (size @(g :.: f)))
{-# INLINE twiddle #-}

twiddles :: forall g f a. (ASZ g, ASZ f, RealFloat a) => g (f (Complex a))
twiddles = omegas (size @(g :.: f))
{-# INLINE twiddles #-}

o8sq :: C
o8sq = omega (8 :: Int) ^ (2 :: Int)

-- Powers of x, starting x^0. Uses 'LScan' for log parallel time
powers :: (LScan f, Pointed f, Num a) => a -> f a
powers = fst . lproducts . point
{-# INLINE powers #-}

-- TODO: Consolidate with powers in TreeTest and rename sensibly. Maybe use
-- "In" and "Ex" suffixes to distinguish inclusive and exclusive cases.

{--------------------------------------------------------------------
    Generic support
--------------------------------------------------------------------}

instance FFT Par1 where
  type FFO Par1 = Par1
  fft = id

#if 0
inTranspose :: (Traversable f', Traversable g, Applicative g', Applicative f)
            => (f (g a) -> f' (g' a)) -> g (f a) -> g' (f' a)
inTranspose = transpose <-- transpose

ffts' :: ( FFT g, Traversable f, Traversable g
         , Applicative (FFO g), Applicative f, RealFloat a) =>
     g (f (Complex a)) -> FFO g (f (Complex a))
ffts' = transpose . fmap fft . transpose
#endif

#if 0

transpose :: g (f C)     -> f (g C)
fmap fft  :: f (g C)     -> f (FFO g C)
transpose :: f (FFO g C) -> FFO g (f C)

#endif

instance ( Zip f,  Traversable f , Traversable g
         , Applicative f, Applicative (FFO f), Applicative (FFO g), Zip (FFO g)
         , Pointed f, Traversable (FFO g), Pointed (FFO g)
         , FFT f, FFT g, LScan f, LScan (FFO g), Sized f, Sized (FFO g) )
      => FFT (g :.: f) where
  type FFO (g :.: f) = FFO f :.: FFO g
  fft = inComp (traverse fft . twiddle . traverse fft . transpose)
  -- fft = inComp (ffts' . transpose . twiddle . ffts')
  {-# INLINE fft #-}

#if 0
  fft = Comp1 . transpose . fmap fft . twiddle . transpose . fmap fft . transpose . unComp1
  fft = Comp1 . traverse fft . twiddle . traverse fft . transpose . unComp1

-- Types in fft for (g :. f):
  unComp1   :: (g :. f) a -> g  (f  a)
  transpose :: g  (f  a)  -> f  (g  a)
  fmap fft  :: f  (g  a)  -> f  (g' a)
  transpose :: f  (g' a)  -> g' (f  a)
  twiddle   :: g' (f  a)  -> g' (f  a)
  fmap fft  :: g' (f  a)  -> g' (f' a)
  transpose :: g' (f' a)  -> f' (g' a)
  Comp1     :: f' (g' a)  -> (f' :. g') a
#endif

#if 0

--   fft = inComp (ffts' . transpose . twiddle . ffts')

ffts'     :: g (f C)     -> FFO g (f C)
twiddle   :: FFO g (f C) -> FFO g (f C)
transpose :: FFO g (f C) -> f (FFO g C)
ffts'     :: f (FFO g C) -> FFO f (FFO g C)

#endif

-- -- Generalization of 'dft' to traversables. Note that liftA2 should
-- -- work zippily (unlike with lists).
-- dftT :: forall f a. (ASZ f, Traversable f, RealFloat a)
--      => Unop (f (Complex a))
-- dftT xs = out <$> indices
--  where
--    out k   = sum (liftA2 (\ n x -> x * ok^n) indices xs) where ok = om ^ k
--    indices = fst iota :: f Int
--    om      = omega (size @f)
-- {-# INLINE dftT #-}

-- | Generic FFT
genericFft :: ( Generic1 f, Generic1 (FFO f)
              , FFT (Rep1 f), FFO (Rep1 f) ~ Rep1 (FFO f) ) => DFTTy f
genericFft = inGeneric1 fft

type GFFT f = (Generic1 f, Generic1 (FFO f), FFT (Rep1 f), FFO (Rep1 f) ~ Rep1 (FFO f))

#define GenericFFT(f,g) instance GFFT (f) => FFT (f) where { type FFO (f) = (g); fft = genericFft }

-- #define GenericFFT(f,g) instance GFFT (f) => FFT (f) where { type FFO (f) = g; INLINE }

-- TODO: Replace Applicative with Zippable.
-- Can't, because Traversable needs Applicative.

-- Perhaps dftT isn't very useful. Its result and argument types match, unlike fft.

{--------------------------------------------------------------------
    Specialized FFT instances.
--------------------------------------------------------------------}

-- I put the specific instances here in order to avoid an import loop between
-- the LPow and RPow modules. I'd still like to find an elegant FFT that maps f
-- to f, and then move the instances to RPow and LPow.

-- Radix 2 butterfly
instance FFT Pair where
  type FFO Pair = Pair
  -- bogus "non-exhaustive" warning in ghc 8.0.2
  -- fft (a :# b) = (a + b) :# (a - b)
  fft = \ (a :# b) -> (a + b) :# (a - b)
  -- fft = dft
  {-# INLINE fft #-}

#ifdef TESTING

#if 0
twiddles :: (ASZ g, ASZ f, RealFloat a) => g (f (Complex a))

as :: f C
(<.> as) :: f C -> C
twiddles :: f (f C)
(<.> as) <$> twiddles :: f C
#endif

-- -- Binary dot product
-- infixl 7 <.>
-- (<.>) :: (Foldable f, Applicative f, Num a) => f a -> f a -> a
-- u <.> v = sum (liftA2 (*) u v)

{--------------------------------------------------------------------
    Simple, quadratic DFT (for specification & testing)
--------------------------------------------------------------------}

-- Adapted from Dave's definition
dftL :: RealFloat a => Unop [Complex a]
dftL xs = [ sum [ x * ok^n | x <- xs | n <- [0 :: Int ..] ]
          | k <- [0 .. length xs - 1], let ok = om ^ k ]
 where
   om = omega (length xs)

{--------------------------------------------------------------------
    Tests
--------------------------------------------------------------------}

-- > powers 2 :: LTree N2 Int
-- B (B (L ((1 :# 2) :# (4 :# 8))))
-- > powers 2 :: LTree N3 Int
-- B (B (B (L (((1 :# 2) :# (4 :# 8)) :# ((16 :# 32) :# (64 :# 128))))))

fftl :: (FFT f, Foldable (FFO f), RealFloat a) => f (Complex a) -> [Complex a]
fftl = toList . fft

type LTree = L.Pow Pair
type RTree = R.Pow Pair

type LC n = LTree n C
type RC n = RTree n C

p1 :: Pair C
p1 = 1 :# 0

tw1 :: LTree N1 (Pair C)
tw1 = twiddles

tw2 :: RTree N2 (Pair C)
tw2 = twiddles

tw3 :: RTree N2 (RTree N2 C)
tw3 = twiddles

tw3' :: [[C]]
tw3' = toList (toList <$> tw3)

-- Adapted from Dave's testing

-- test :: (FFT f, Foldable f, Foldable (FFO f)) => f C -> IO ()
-- test fx =
--   do ps "\nTesting input" xs
--      ps "Expected output" (dftL xs)
--      ps "Actual output  " (toList (fft fx))
--  where
--    ps label z = putStrLn (label ++ ": " ++ show z)
--    xs = toList fx

#if 0
t0 :: LC N0
t0 = L.fromList [1]

t1 :: LC N1
t1 = L.fromList [1, 0]

t2s :: [LC N2]
t2s = L.fromList <$>
        [ [1,  0,  0,  0]  -- Delta
        , [1,  1,  1,  1]  -- Constant
        , [1, -1,  1, -1]  -- Nyquist
        , [1,  0, -1,  0]  -- Fundamental
        , [0,  1,  0, -1]  -- Fundamental w/ 90-deg. phase lag
       ]

tests :: IO ()
tests = do test p1
           test t0
           test t1
           mapM_ test t2s
#endif

infix 4 ===
(===) :: Eq b => (a -> b) -> (a -> b) -> a -> Bool
(f === g) x = f x == g x

infix 4 =~=
(=~=) :: ApproxEq b => (a -> b) -> (a -> b) -> a -> Bool
(f =~= g) x = f x =~ g x

fftIsDftL :: (FFT f, Foldable f, Foldable (FFO f), RealFloat a, ApproxEq a) =>
             f (Complex a) -> Bool
fftIsDftL = toList . fft =~= dftL . toList

dftTIsDftL :: (ASZ f, Traversable f, RealFloat a, ApproxEq a) =>
              f (Complex a) -> Bool
dftTIsDftL = toList . dftT =~= dftL . toList

dftIsDftL :: (ASZ f, Foldable f, RealFloat a, ApproxEq a) =>
             f (Complex a) -> Bool
dftIsDftL = toList . dft =~= dftL . toList

-- -- TEMP:
-- dftDft :: (ASZ f, Traversable f, RealFloat a, ApproxEq a) =>
--           f (Complex a) -> ([Complex a], [Complex a])
-- dftDft xs = (toList . dft $ xs, dftL . toList $ xs)

{--------------------------------------------------------------------
    Properties to test
--------------------------------------------------------------------}

transposeTwiddleCommutes :: (ASZ g, Traversable g, ASZ f, (ApproxEq (f (g C))))
                         => g (f C) -> Bool
transposeTwiddleCommutes =
 twiddle . transpose =~= transpose . twiddle

prop_transposeTwiddle_L3P :: LTree N3 (Pair C) -> Bool
prop_transposeTwiddle_L3P = transposeTwiddleCommutes

prop_transposeTwiddle_R3P :: RTree N3 (Pair C) -> Bool
prop_transposeTwiddle_R3P = transposeTwiddleCommutes

-- dft tests fail. Hm!

-- prop_dft_R3 :: RTree N3 C -> Bool
-- prop_dft_R3 = dftIsDftL

-- prop_dft_L3 :: LTree N3 C -> Bool
-- prop_dft_L3 = dftIsDftL

prop_dftT_p :: Pair C -> Bool
prop_dftT_p = dftTIsDftL

prop_dftT_L3 :: LTree N3 C -> Bool
prop_dftT_L3 = dftTIsDftL

prop_dftT_R3 :: RTree N3 C -> Bool
prop_dftT_R3 = dftTIsDftL

prop_fft_p :: Pair C -> Bool
prop_fft_p = fftIsDftL

prop_fft_L1 :: LTree N1 C -> Bool
prop_fft_L1 = fftIsDftL

prop_fft_L2 :: LTree N2 C -> Bool
prop_fft_L2 = fftIsDftL

prop_fft_L3 :: LTree N3 C -> Bool
prop_fft_L3 = fftIsDftL

prop_fft_L4 :: LTree N4 C -> Bool
prop_fft_L4 = fftIsDftL

prop_fft_R1 :: RTree N1 C -> Bool
prop_fft_R1 = fftIsDftL

prop_fft_R2 :: RTree N2 C -> Bool
prop_fft_R2 = fftIsDftL

prop_fft_R3 :: RTree N3 C -> Bool
prop_fft_R3 = fftIsDftL

prop_fft_R4 :: RTree N4 C -> Bool
prop_fft_R4 = fftIsDftL

-- TH oddity
return []

runTests :: IO Bool
runTests = $quickCheckAll

-- end of tests
#endif

max_stars_repo_name: kenranunderscore/concat

max_issues_repo_name: con-kitty/concat

max_forks_repo_name: con-kitty/concat


text:
{-# LANGUAGE TypeFamilies #-}

module FPNLA.Operations.Parameters(
    -- * Elements
    Elt(..),
    -- * Strategies and contexts
    StratCtx(),

    -- * Result type
    -- | In BLAS it's common that operations in higher levels use operations in the lower levels, so, an operation in level three that by its signature manipulates matrices only, internally uses level two operations that manipulates vectors. In order to avoid the /show . read/ problem, the type of the vector (or any other internal data type) must appear in the signature of an operation.
    -- To solve the problem we use phantom types to pass the internally used types to the Haskell type system.
    ResM(),
    ResV(),
    ResS(),
    blasResultM,
    blasResultV,
    blasResultS,
    getResultDataM,
    getResultDataV,
    getResultDataS,

    -- * Miscellaneous
    TransType(..),
    UnitType(..),
    TriangType(..),
    unTransT,
    unUnitT,
    unTriangT,
    elemTrans_m,
    dimTrans_m,
    elemSymm,
    dimTriang,
    elemUnit_m,
    dimUnit_m,
    elemTransUnit_m,
    dimTransUnit_m,
    transTrans_m

) where

import FPNLA.Matrix(Matrix(..))
import Data.Complex (Complex, conjugate)
import Data.Tuple (swap)

-- | This class represents the elements that can be used in the BLAS operations.
-- The elements in BLAS are real or complex numbers, so we provide default instances for the Haskell 'Double', 'Float' and 'Complex' types.
class (Eq e, Floating e) => Elt e where
    -- | Returns the conjugate of a number. For real numbers it's the identity function and for complex numbers it's the common 'Complex.conjugate' function.
    getConjugate :: e -> e
    getConjugate = id

instance Elt Double
instance Elt Float
instance (RealFloat e) => Elt (Complex e) where
    getConjugate = conjugate

-- | This type family is used to represent the /context/ of an operation.
-- A particular implementation is a combination of an algorithm and a parallelism technique, and we call it a /strategy/. A particular strategy may need particular information to execute. For example, an operation that computes the matrix-matrix multiplication by splitting the matrices in blocks must require the size of the blocks.
-- With this context we allows to pass any additional information that the operation needs to execute as parameters, but maintaining a common signature.
-- The /s/ type parameter is the strategy so, there must exist a Haskell data type to represent a particular strategy.
type family StratCtx s :: *

-- | The 'ResM' data type is used as result of level three BLAS operations and returns a matrix /m/ of elements /e/ and contains the strategy /s/ and vector /v/ as phantom types.
data ResM s (v :: * -> *) m e = ResM { unResM :: m e } deriving (Show)
-- | The 'ResV' data type is used as result of level two BLAS operations and returns a vector /v/ of elements /e/ and contains the strategy /s/ as phantom types.
data ResV s v e = ResV { unResV :: v e } deriving (Show)
-- | The 'ResS' data type is used as result of level one BLAS operations and returns an scalar /e/ and contains the strategy /s/ as phantom types.
data ResS s e = ResS { unResS :: e } deriving (Show)

-- | Wrap a matrix into a 'ResM'.
blasResultM :: m e -> ResM s v m e
blasResultM = ResM
-- | Unwrap a matrix from a 'ResM'.
getResultDataM :: ResM s v m e -> m e
getResultDataM = unResM

-- | Wrap a vector into a 'ResV'.
blasResultV :: v e -> ResV s v e
blasResultV = ResV
-- | Unwrap a vector from a 'ResV'.
getResultDataV :: ResV s v e -> v e
getResultDataV = unResV

-- | Wrap a scalar into a 'ResS'.
blasResultS :: e -> ResS s e
blasResultS = ResS
-- | Unwrap a scalar from a 'ResS'.
getResultDataS :: ResS s e -> e
getResultDataS = unResS

-- | Indicates if a matrix must be considered as normal, transposed or transposed conjugated.
-- This is part of the common flags in the BLAS operation signatures and it's useful to work with a transposed matrix without really computing the transposed matrix.
data TransType m = Trans m | NoTrans m | ConjTrans m deriving (Eq, Show)
-- | Indicates if a matrix must be considered as unitary or not. An unitary matrix is a matrix that contains ones in the diagonal.
-- This is part of the common flags in the BLAS operation signatures.
data UnitType m = Unit m | NoUnit m deriving (Eq, Show)
-- | Indicates that a matrix is symmetric and with which triangular part of the matrix the operation is going to work ('Upper' or 'Lower').
-- The operation only will see the indicated part of the matrix and should not try to access the other part.
-- This is part of the common flags in the BLAS operation signatures.
data TriangType m = Lower m | Upper m deriving (Eq, Show)

-- | Given a data type flagged by a TransType, returns a pair containing the TransType constructor and the data type.
unTransT :: TransType a -> (b -> TransType b, a)
unTransT (Trans a) = (Trans, a)
unTransT (NoTrans a) = (NoTrans, a)
unTransT (ConjTrans a) = (ConjTrans, a)

-- | Given a data type flagged by a UnitType, returns a pair containing the UnitType constructor and the data type.
unUnitT :: UnitType a -> (b -> UnitType b, a)
unUnitT (Unit a) = (Unit, a)
unUnitT (NoUnit a) = (NoUnit, a)

-- | Given a data type flagged by a TriangType, returns a pair containing the TriangType constructor and the data type.
unTriangT :: TriangType a -> (b -> TriangType b, a)
unTriangT (Lower a) = (Lower, a)
unTriangT (Upper a) = (Upper, a)

-- | Given an /i,j/ position and a TransType flagged matrix, returns the element in that position without computing the transpose.
elemTrans_m :: (Elt e, Matrix m e) => Int -> Int -> TransType (m e) -> e
elemTrans_m i j (NoTrans m) = elem_m i j m
elemTrans_m i j (Trans m) = elem_m j i m
elemTrans_m i j (ConjTrans m) = getConjugate $ elem_m j i m

-- | Given a TransType flagged matrix, returns the dimension of the matrix without computing the transpose.
dimTrans_m :: (Matrix m e) => TransType (m e) -> (Int, Int)
dimTrans_m (NoTrans m) = dim_m m
dimTrans_m (ConjTrans m) = swap $ dim_m m
dimTrans_m (Trans m) = swap $ dim_m m

-- | Given an /i,j/ position and a TransType flagged matrix, returns the element in that position only accessing the part indicated by the TransType.
elemSymm :: (Matrix m e) => Int -> Int -> TriangType (m e) -> e
elemSymm i j (Upper m)
    | i > j = elem_m j i m
    | otherwise = elem_m i j m
elemSymm i j (Lower m)
    | i > j = elem_m i j m
    | otherwise = elem_m j i m

-- | Given a TransType flagged matrix, returns the dimension of the matrix.
dimTriang :: (Matrix m e) => TriangType (m e) -> (Int, Int)
dimTriang = dim_m . snd . unTriangT

-- | Given an /i,j/ position and a UnitType flagged matrix, returns the element in that position. If the matrix is flagged as Unit and /i == j/ (the element is in the diagonal) returns one.
elemUnit_m :: (Elt e, Matrix m e) => Int -> Int -> UnitType (m e) -> e
elemUnit_m i j (Unit m)
    | i == j = 1
    | otherwise = elem_m i j m
elemUnit_m i j (NoUnit m) = elem_m i j m

-- | Given a UnitType flagged matrix, returns the dimension of the matrix.
dimUnit_m :: (Matrix m e) => UnitType (m e) -> (Int, Int)
dimUnit_m (Unit m) = dim_m m
dimUnit_m (NoUnit m) = dim_m m

-- | Given an /i,j/ position and a TransType-UnitType flagged matrix, returns the element in that position without computing the transpose.
elemTransUnit_m :: (Elt e, Matrix m e) => Int -> Int -> TransType (UnitType (m e)) -> e
elemTransUnit_m i j (NoTrans pmA) = elemUnit_m i j pmA
elemTransUnit_m i j (Trans pmA) = elemUnit_m j i pmA
elemTransUnit_m i j (ConjTrans pmA) = getConjugate $ elemUnit_m j i pmA

-- | Given a TransType-UnitType flagged matrix, returns the dimension of the matrix.
dimTransUnit_m :: Matrix m e => TransType (UnitType (m e)) -> (Int, Int)
dimTransUnit_m = dimUnit_m . snd . unTransT

-- | Given a TransType flagged matrix, computes and returns its transpose.
transTrans_m :: (Elt e, Matrix m e) => TransType (m e) -> m e
transTrans_m (NoTrans m) = m
transTrans_m (ConjTrans m) = map_m getConjugate $ transpose_m m
transTrans_m (Trans m) = transpose_m m

max_stars_repo_name: mauroblanco/fpnla

max_issues_repo_name: mauroblanco/fpnla

max_forks_repo_name: mauroblanco/fpnla


text:
module Calibrate where

import Data.WAVE
import Data.Complex
import Data.List
import Numeric.FFT

import Util
import ExtProg

-- Calibrated or uncalibrated voice: calibrated means we know the fundamental frequency
-- for pitch 50, uncalibrated means we use the given voice name and text to get it.

type Calibration = Either (String, String, String) Double

-- Determine the fundamental frequency for a voice at the given pitch. Get the spectrum,
-- smoothen it and get the first peak.

voiceFundamental :: String -> String -> Int -> String -> IO (Maybe Double)

voiceFundamental exec text pitch voice = do
  p <- voiceSpectrumPitchText exec text pitch voice
  let smpp = filter ((> 0.05) . snd) $ smoothen 3 p
      lmpp = locMax smpp
      fund = case lmpp of
        [] -> Nothing
        _ -> Just $ fst $ head lmpp
  return fund

-- The most general case: set text, calibration pitch, and voice name.

voiceSpectrumPitchText :: String -> String -> Int -> String -> IO [(Double, Double)]

voiceSpectrumPitchText exec text pitch voice = do
  spd <- findSpeed exec 120 voice pitch 0.25 text (1, 300) >>= return . fst
  wav <- runEspeak exec 120 voice pitch spd text >>= getSample
  let fft = filter ((< 2000) . fst) $ waveFFT wav
  return fft

-- Run a FFT on a signal given WAVE, returning pairs of frequency-normalized magnitude
-- sorted by decreasing magnitude

waveFFT :: WAVE -> [(Double, Double)]

waveFFT wav =
  let sampldbl = map (sampleToDouble . head) $ waveSamples wav
      sampcmplx = map (:+ 0) sampldbl
      freqdom = map magnitude $ fft sampcmplx
      rate = waveFrameRate $ waveHeader wav
      freqs = map (\i -> fromIntegral (i * fromIntegral rate) / fromIntegral (length freqdom))
                  [0 .. length freqdom - 1]
      spectrum = zip freqs (normalize freqdom)
      sortDescSnd = sortBy $ \(x1, y1) (x2, y2) -> compare y2 y1
  in  sortDescSnd spectrum

-- Find a minimal range of speed that approximates the given time in seconds

findSpeed :: String -> Int -> String -> Int -> Rational -> String -> (Int, Int) -> IO (Int, Int)

findSpeed exec ampl voice pitch tsec utter (lo, hi) =
  rangeSearch 1 (lo, hi) $ \mid -> do
    wav <- runEspeak exec ampl voice pitch mid utter >>= getSample
    return $ compare (soundLength wav) tsec

max_stars_repo_name: dmgolubovsky/hsespeak

max_issues_repo_name: dmgolubovsky/hsespeak

max_forks_repo_name: dmgolubovsky/hsespeak


text:
-- This file is auto-generated.  Do not edit directly.
-- |
--   Stability: Experimental
--
--   Generic interface to Blas using unsafe foreign calls.  Refer to the GHC
--   documentation for more information regarding appropriate use of safe and unsafe
--   foreign calls.
--
--   The functions here are named in a similar fashion to the original Blas interface, with
--   the type-dependent letter(s) removed.  Some functions have been merged with
--   others to allow the interface to work on both real and complex numbers.  If you can't a
--   particular function, try looking for its corresponding complex equivalent (e.g.
--   @symv@ is a special case of 'hemv' applied to real numbers).
--
--   Note: although complex versions of @rot@ and @rotg@ exist in many implementations,
--   they are not part of the official Blas standard and therefore not included here.  If
--   you /really/ need them, submit a ticket so we can try to come up with a solution.
--
--   The documentation here is still incomplete.  Consult the
--   <http://netlib.org/blas/#_blas_routines official documentation> for more
--   information.
--
--   Notation:
--
--     * @⋅@ denotes dot product (without any conjugation).
--     * @*@ denotes complex conjugation.
--     * @⊤@ denotes transpose.
--     * @†@ denotes conjugate transpose (Hermitian conjugate).
--
--   Conventions:
--
--     * All scalars are denoted with lowercase Greek letters
--     * All vectors are denoted with lowercase Latin letters and are
--       assumed to be column vectors (unless transposed).
--     * All matrices are denoted with uppercase Latin letters.
--
--   /Since: 0.1.1/
module Blas.Specialized.ComplexDouble.Unsafe (
  -- * Level 1: vector-vector operations
  -- ** Basic operations
    swap
  , scal
  , copy
  , axpy
  , dotu
  , dotc
  -- ** Norm operations
  , nrm2
  , asum
  , iamax
  -- * Level 2: matrix-vector operations
  -- ** Multiplication
  , gemv
  , gbmv
  , hemv
  , hbmv
  , hpmv
  -- ** Triangular operations
  , trmv
  , tbmv
  , tpmv
  , trsv
  , tbsv
  , tpsv
  -- ** Rank updates
  , geru
  , gerc
  , her
  , hpr
  , her2
  , hpr2
  -- * Level 3: matrix-matrix operations
  -- ** Multiplication
  , gemm
  , symm
  , hemm
  -- ** Rank updates
  , syrk
  , herk
  , syr2k
  , her2k
  -- ** Triangular operations
  , trmm
  , trsm
  ) where
import Prelude (Double, Int, IO)
import Foreign (Ptr)
import Data.Complex (Complex((:+)))
import FFI (getReturnValue)
import Blas.Primitive.Types (Order, Transpose, Uplo, Diag, Side)
import qualified Blas.Primitive.Unsafe as C

-- | Swap two vectors:
--
--   > (x, y) ← (y, x)
swap :: Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
swap = C.zswap

-- | Multiply a vector by a scalar.
--
--   > x ← α x
scal :: Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
scal n (alpha :+ 0) = C.zdscal n alpha
scal n alpha = C.zscal n alpha

-- | Copy a vector into another vector:
--
--   > y ← x
copy :: Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
copy = C.zcopy

-- | Add a scalar-vector product to a vector.
--
--   > y ← α x + y
axpy :: Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
axpy = C.zaxpy

-- | Calculate the bilinear dot product of two vectors:
--
--   > x ⋅ y ≡ ∑[i] x[i] y[i]
dotu :: Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO (Complex Double)
dotu a b c d e = getReturnValue (C.zdotu_sub a b c d e)

-- | Calculate the sesquilinear dot product of two vectors.
--
--   > x* ⋅ y ≡ ∑[i] x[i]* y[i]
dotc :: Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO (Complex Double)
dotc a b c d e = getReturnValue (C.zdotc_sub a b c d e)

-- | Calculate the Euclidean (L²) norm of a vector:
--
--   > ‖x‖₂ ≡ √(∑[i] x[i]²)
nrm2 :: Int
     -> Ptr (Complex Double)
     -> Int
     -> IO Double
nrm2 = C.dznrm2

-- | Calculate the Manhattan (L¹) norm, equal to the sum of the magnitudes of the elements:
--
--   > ‖x‖₁ = ∑[i] |x[i]|
asum :: Int
     -> Ptr (Complex Double)
     -> Int
     -> IO Double
asum = C.dzasum

-- | Calculate the index of the element with the maximum magnitude (absolute value).
iamax :: Int
      -> Ptr (Complex Double)
      -> Int
      -> IO Int
iamax = C.izamax

-- | Perform a general matrix-vector update.
--
--   > y ← α T(A) x + β y
gemv :: Order
     -> Transpose
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
gemv = C.zgemv

-- | Perform a general banded matrix-vector update.
--
--   > y ← α T(A) x + β y
gbmv :: Order
     -> Transpose
     -> Int
     -> Int
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
gbmv = C.zgbmv

-- | Perform a hermitian matrix-vector update.
--
--   > y ← α A x + β y
hemv :: Order
     -> Uplo
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
hemv = C.zhemv

-- | Perform a hermitian banded matrix-vector update.
--
--   > y ← α A x + β y
hbmv :: Order
     -> Uplo
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
hbmv = C.zhbmv

-- | Perform a hermitian packed matrix-vector update.
--
--   > y ← α A x + β y
hpmv :: Order
     -> Uplo
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Ptr (Complex Double)
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
hpmv = C.zhpmv

-- | Multiply a triangular matrix by a vector.
--
--   > x ← T(A) x
trmv :: Order
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
trmv = C.ztrmv

-- | Multiply a triangular banded matrix by a vector.
--
--   > x ← T(A) x
tbmv :: Order
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
tbmv = C.ztbmv

-- | Multiply a triangular packed matrix by a vector.
--
--   > x ← T(A) x
tpmv :: Order
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Ptr (Complex Double)
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
tpmv = C.ztpmv

-- | Multiply an inverse triangular matrix by a vector.
--
--   > x ← T(A⁻¹) x
trsv :: Order
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
trsv = C.ztrsv

-- | Multiply an inverse triangular banded matrix by a vector.
--
--   > x ← T(A⁻¹) x
tbsv :: Order
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
tbsv = C.ztbsv

-- | Multiply an inverse triangular packed matrix by a vector.
--
--   > x ← T(A⁻¹) x
tpsv :: Order
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Ptr (Complex Double)
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
tpsv = C.ztpsv

-- | Perform an unconjugated rank-1 update of a general matrix.
--
--   > A ← α x y⊤ + A
geru :: Order
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
geru = C.zgeru

-- | Perform a conjugated rank-1 update of a general matrix.
--
--   > A ← α x y† + A
gerc :: Order
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
gerc = C.zgerc

-- | Perform a rank-1 update of a Hermitian matrix.
--
--   > A ← α x y† + A
her :: Order
    -> Uplo
    -> Int
    -> Double
    -> Ptr (Complex Double)
    -> Int
    -> Ptr (Complex Double)
    -> Int
    -> IO ()
her = C.zher

-- | Perform a rank-1 update of a Hermitian packed matrix.
--
--   > A ← α x y† + A
hpr :: Order
    -> Uplo
    -> Int
    -> Double
    -> Ptr (Complex Double)
    -> Int
    -> Ptr (Complex Double)
    -> IO ()
hpr = C.zhpr

-- | Perform a rank-2 update of a Hermitian matrix.
--
--   > A ← α x y† + y (α x)† + A
her2 :: Order
     -> Uplo
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
her2 = C.zher2

-- | Perform a rank-2 update of a Hermitian packed matrix.
--
--   > A ← α x y† + y (α x)† + A
hpr2 :: Order
     -> Uplo
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> IO ()
hpr2 = C.zhpr2

-- | Perform a general matrix-matrix update.
--
--   > C ← α T(A) U(B) + β C
gemm :: Order -- ^ Layout of all the matrices.
     -> Transpose -- ^ The operation @T@ to be applied to @A@.
     -> Transpose -- ^ The operation @U@ to be applied to @B@.
     -> Int -- ^ Number of rows of @T(A)@ and @C@.
     -> Int -- ^ Number of columns of @U(B)@ and @C@.
     -> Int -- ^ Number of columns of @T(A)@ and number of rows of @U(B)@.
     -> Complex Double -- ^ Scaling factor @α@ of the product.
     -> Ptr (Complex Double) -- ^ Pointer to a matrix @A@.
     -> Int -- ^ Stride of the major dimension of @A@.
     -> Ptr (Complex Double) -- ^ Pointer to a matrix @B@.
     -> Int -- ^ Stride of the major dimension of @B@.
     -> Complex Double -- ^ Scaling factor @β@ of the original @C@.
     -> Ptr (Complex Double) -- ^ Pointer to a mutable matrix @C@.
     -> Int -- ^ Stride of the major dimension of @C@.
     -> IO ()
gemm = C.zgemm

-- | Perform a symmetric matrix-matrix update.
--
--   > C ← α A B + β C    or    C ← α B A + β C
--
--   where @A@ is symmetric.  The matrix @A@ must be in an unpacked format, although the
--   routine will only access half of it as specified by the @'Uplo'@ argument.
symm :: Order -- ^ Layout of all the matrices.
     -> Side -- ^ Side that @A@ appears in the product.
     -> Uplo -- ^ The part of @A@ that is used.
     -> Int -- ^ Number of rows of @C@.
     -> Int -- ^ Number of columns of @C@.
     -> Complex Double -- ^ Scaling factor @α@ of the product.
     -> Ptr (Complex Double) -- ^ Pointer to a symmetric matrix @A@.
     -> Int -- ^ Stride of the major dimension of @A@.
     -> Ptr (Complex Double) -- ^ Pointer to a matrix @B@.
     -> Int -- ^ Stride of the major dimension of @B@.
     -> Complex Double -- ^ Scaling factor @α@ of the original @C@.
     -> Ptr (Complex Double) -- ^ Pointer to a mutable matrix @C@.
     -> Int -- ^ Stride of the major dimension of @C@.
     -> IO ()
symm = C.zsymm

-- | Perform a Hermitian matrix-matrix update.
--
--   > C ← α A B + β C    or    C ← α B A + β C
--
--   where @A@ is Hermitian.  The matrix @A@ must be in an unpacked format, although the
--   routine will only access half of it as specified by the @'Uplo'@ argument.
hemm :: Order -- ^ Layout of all the matrices.
     -> Side -- ^ Side that @A@ appears in the product.
     -> Uplo -- ^ The part of @A@ that is used.
     -> Int -- ^ Number of rows of @C@.
     -> Int -- ^ Number of columns of @C@.
     -> Complex Double -- ^ Scaling factor @α@ of the product.
     -> Ptr (Complex Double) -- ^ Pointer to a Hermitian matrix @A@.
     -> Int -- ^ Stride of the major dimension of @A@.
     -> Ptr (Complex Double) -- ^ Pointer to a matrix @B@.
     -> Int -- ^ Stride of the major dimension of @B@.
     -> Complex Double -- ^ Scaling factor @α@ of the original @C@.
     -> Ptr (Complex Double) -- ^ Pointer to a mutable matrix @C@.
     -> Int -- ^ Stride of the major dimension of @C@.
     -> IO ()
hemm = C.zhemm

-- | Perform a symmetric rank-k update.
--
--   > C ← α A A⊤ + β C    or    C ← α A⊤ A + β C
syrk :: Order
     -> Uplo
     -> Transpose
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
syrk = C.zsyrk

-- | Perform a Hermitian rank-k update.
--
--   > C ← α A A† + β C    or    C ← α A† A + β C
herk :: Order
     -> Uplo
     -> Transpose
     -> Int
     -> Int
     -> Double
     -> Ptr (Complex Double)
     -> Int
     -> Double
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
herk = C.zherk

-- | Perform a symmetric rank-2k update.
--
--   > C ← α A B⊤ + α* B A⊤ + β C    or    C ← α A⊤ B + α* B⊤ A + β C
syr2k :: Order
      -> Uplo
      -> Transpose
      -> Int
      -> Int
      -> Complex Double
      -> Ptr (Complex Double)
      -> Int
      -> Ptr (Complex Double)
      -> Int
      -> Complex Double
      -> Ptr (Complex Double)
      -> Int
      -> IO ()
syr2k = C.zsyr2k

-- | Perform a Hermitian rank-2k update.
--
--   > C ← α A B† + α* B A† + β C    or    C ← α A† B + α* B† A + β C
her2k :: Order
      -> Uplo
      -> Transpose
      -> Int
      -> Int
      -> Complex Double
      -> Ptr (Complex Double)
      -> Int
      -> Ptr (Complex Double)
      -> Int
      -> Double
      -> Ptr (Complex Double)
      -> Int
      -> IO ()
her2k = C.zher2k

-- | Perform a triangular matrix-matrix multiplication.
--
--   > B ← α T(A) B    or    B ← α B T(A)
--
--   where @A@ is triangular.
trmm :: Order
     -> Side
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
trmm = C.ztrmm

-- | Perform an inverse triangular matrix-matrix multiplication.
--
--   > B ← α T(A⁻¹) B    or    B ← α B T(A⁻¹)
--
--   where @A@ is triangular.
trsm :: Order
     -> Side
     -> Uplo
     -> Transpose
     -> Diag
     -> Int
     -> Int
     -> Complex Double
     -> Ptr (Complex Double)
     -> Int
     -> Ptr (Complex Double)
     -> Int
     -> IO ()
trsm = C.ztrsm

max_stars_repo_name: Rufflewind/blas-hs

max_issues_repo_name: Rufflewind/blas-hs

max_forks_repo_name: Rufflewind/blas-hs


text:
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Strict              #-}
module Pinwheel.FourierSeries2D where

import           Control.Monad                  as M
import           Control.Monad.IO.Class
import           Control.Monad.Trans.Resource
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX as A
import           Data.Array.IArray              as IA
import           Data.Array.Repa                as R
import           Data.Complex
import           Data.Conduit                   as C
import           Data.Conduit.List              as CL
import           Data.List                      as L
import           Data.Vector.Generic            as VG
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.FourierSeries2D
import           Math.Gamma
import           Pinwheel.Base
import           Pinwheel.List
import           Utils.BLAS
import           Utils.Distribution
import           Utils.List
import           Utils.Parallel                 hiding ((.|))
import           Utils.SimpsonRule
import           Utils.Time
import DFT.Plan
import Filter.Utils
import Control.Concurrent.Async
import Numeric.GSL.Special.Bessel

pinwheelFourierSeries ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     , VG.Vector vector (Complex e)
     , NFData (vector (Complex e))
     )
  => [Int]
  -> [PTX]
  -> Int
  -> Int
  -> e
  -> e
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> Int
  -> Int
  -> IO (IA.Array (Int, Int) (vector (Complex e)))
pinwheelFourierSeries deviceIDs ptxs numR2Freqs numPoints delta periodR2 phiFreq rhoFreq thetaFreq rFreq sigma numBatchR2 numBatchPinwheelFreqs = do
  let idxs =
        [ (radialFreq, angularFreq)
        | radialFreq <- pinwheelFreqs rhoFreq rFreq
        , angularFreq <- pinwheelFreqs phiFreq thetaFreq
        ]
      centerR2Freq = div numR2Freqs 2
  pinwheels <-
    computeUnboxedP .
    R.traverse
      (fromListUnboxed (Z :. (L.length idxs)) idxs)
      (\(Z :. freqs) -> (Z :. freqs :. numR2Freqs :. numR2Freqs)) $ \f (Z :. pFreq :. xFreq :. yFreq) ->
      fourierMellin
        sigma
        (snd $ f (Z :. pFreq))
        (fst $ f (Z :. pFreq))
        ( delta * fromIntegral (xFreq - centerR2Freq)
        , delta * fromIntegral (yFreq - centerR2Freq))
  -- pinwheels <-
  --   M.mapM
  --     (\idx ->
  --        fmap
  --          (CuMat (numR2Freqs ^ 2) (L.length idx) .
  --           CuVecHost . VU.convert . toUnboxed) .
  --        computeP .
  --        R.traverse
  --          (fromListUnboxed (Z :. (L.length idx)) idx)
  --          (\(Z :. freqs) -> (Z :. numR2Freqs :. numR2Freqs :. freqs)) $ \f (Z :. xFreq :. yFreq :. pFreq) ->
  --          fourierMellin
  --            sigma
  --            (snd $ f (Z :. pFreq))
  --            (fst $ f (Z :. pFreq))
  --            ( fromIntegral (xFreq - centerR2Freq)
  --            , fromIntegral (yFreq - centerR2Freq))) .
  --   divideListN numBatchPinwheelFreqs $
  --   idxs
  -- printCurrentTime "Compute Fourier Series... "
  -- arr <-
  --   computeFourierSeriesR2Stream
  --     deviceIDs
  --     ptxs
  --     numR2Freqs
  --     numPoints
  --     periodR2
  --     delta
  --     numBatchR2
  --     pinwheels
  -- printCurrentTime "Compute Fourier Series Done"
  let (radialLB, radialUB) = pinwheelFreqsBound rhoFreq rFreq
      (angularLB, angularUB) = pinwheelFreqsBound phiFreq thetaFreq
  return .
    listArray ((radialLB, angularLB), (radialUB, angularUB)) .
    parMap
      rdeepseq
      (\i ->
         VG.convert . toUnboxed . computeS . R.slice pinwheels $  -- arr $
         (Z :. i :. All :. All)) $
    [0 .. (L.length idxs - 1)]

-- Stream
{-# INLINE source #-}
source :: Int -> Int -> Int -> Int -> ConduitT () (Int, Int) (ResourceT IO) ()
source phiFreq rhoFreq thetaFreq rFreq =
  CL.sourceList
    [ (radialFreq, angularFreq)
    | radialFreq <- pinwheelFreqs rhoFreq rFreq
    , angularFreq <- pinwheelFreqs phiFreq thetaFreq
    ]

conduit ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     , VG.Vector vector (Complex e)
     )
  => [Int]
  -> Int
  -> Int
  -> e
  -> Int
  -> e
  -> [[CuMat (Complex e)]]
  -> ConduitT (Int, Int) [vector (Complex e)] (ResourceT IO) ()
conduit deviceIDs numR2Freqs numPoints periodR2 batchSize sigma inverseR2Harmonics = do
  liftIO $ printCurrentTime "conduit"
  idx <- CL.take batchSize
  unless
    (L.null idx)
    (do let centerR2Freq = div numR2Freqs 2
        pinwheels <-
          fmap
            (CuMat (L.length idx) (numR2Freqs ^ 2) .
             CuVecHost . VU.convert . toUnboxed) .
          liftIO .
          computeP .
          R.traverse
            (fromListUnboxed (Z :. (L.length idx)) idx)
            (\(Z :. freqs) -> (Z :. freqs :. numR2Freqs :. numR2Freqs)) $ \f (Z :. pFreq :. xFreq :. yFreq) ->
            fourierMellin
              sigma
              (snd $ f (Z :. pFreq))
              (fst $ f (Z :. pFreq))
              ( fromIntegral (xFreq - centerR2Freq)
              , fromIntegral (yFreq - centerR2Freq))
        arr <-
          liftIO $
          computeFourierSeriesR2
            deviceIDs
            numR2Freqs
            numPoints
            periodR2
            inverseR2Harmonics
            [pinwheels]
        yield .
          L.map
            (\i ->
               VG.convert . toUnboxed . computeS . R.slice arr $
               (Z :. i :. All :. All)) $
          [0 .. (L.length idx - 1)])

pinwheelFourierSeriesStream ::
     ( Storable e
     , CUBLAS (Complex e)
     , Unbox e
     , RealFloat e
     , A.Elt e
     , A.Elt (Complex e)
     , Floating (A.Exp e)
     , A.FromIntegral Int e
     , VG.Vector vector (Complex e)
     )
  => [Int]
  -> [PTX]
  -> Int
  -> Int
  -> e
  -> e
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> Int
  -> Int
  -> Int
  -> IO (IA.Array (Int, Int) (vector (Complex e)))
pinwheelFourierSeriesStream deviceIDs ptxs numR2Freqs numPoints delta periodR2 phiFreq rhoFreq thetaFreq rFreq sigma numBatchR2 batchSize numBatchPinwheelFreqs = do
  inverseR2Harmonics <-
    createInverseHarmonicMatriesGPU
      ptxs
      numBatchR2
      numPoints
      numR2Freqs
      periodR2
      delta
  let (radialLB, radialUB) = pinwheelFreqsBound rhoFreq rFreq
      (angularLB, angularUB) = pinwheelFreqsBound rhoFreq rFreq
  fmap (listArray ((radialLB, angularLB), (radialUB, radialUB)) . L.concat) .
    runConduitRes $
    source phiFreq rhoFreq thetaFreq rFreq .|
    conduit
      deviceIDs
      numR2Freqs
      numPoints
      periodR2
      batchSize
      sigma
      inverseR2Harmonics .|
    CL.consume

-- 2D spatial domain to frequency domain:
-- transform: cis (-freq * x)
-- inverse transform: cis (freq * x)

-- Analytical solution type 1
-- The pinwheel harmonics are cis (-freq * x)
{-# INLINE analyticalFourierSeriesFunc1 #-}
analyticalFourierSeriesFunc1 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e))
  => Int
  -> Int
  -> e
  -> e
  -> e
  -> e
  -> e
  -> Complex e
analyticalFourierSeriesFunc1 angularFreq radialFreq sigma periodR2 periodEnv phi rho =
  let radialConst = 2 * pi / log periodEnv
   in pi * ((0 :+ 1) ^ abs angularFreq) *
      cis (fromIntegral (-angularFreq) * phi) *
      ((periodR2 / (pi * rho) :+ 0) **
       ((2 + sigma) :+ (radialConst * fromIntegral (-radialFreq)))) *
      gamma
        (((2 + fromIntegral (abs angularFreq) + sigma) :+
          (radialConst * fromIntegral (-radialFreq))) /
         2) /
      gamma
        (((fromIntegral (abs angularFreq) - sigma) :+
          (radialConst * fromIntegral radialFreq)) /
         2)

-- The pinwheel is in the frequency domain, the function computes the Fourier series
{-# INLINE analyticalFourierSeries1 #-}
analyticalFourierSeries1 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierSeries1 numPoints delta angularFreq radialFreq sigma periodR2 periodEnv =
  let center = div numPoints 2
  in fromFunction (Z :. numPoints :. numPoints) $ \(Z :. i :. j) ->
       let x = fromIntegral $ i - center
           y = fromIntegral $ j - center
           rho = sqrt $ x ^ 2 + y ^ 2
           phi = atan2 y x
       in if rho == 0
            then 0
            else analyticalFourierSeriesFunc1
                   angularFreq
                   radialFreq
                   sigma
                   periodR2
                   periodEnv
                   phi
                   (rho * delta)

-- The pinwheel is in the spatial domain, the function computes the Fourier coefficients
{-# INLINE analyticalFourierCoefficients1 #-}
analyticalFourierCoefficients1 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierCoefficients1 numFreqs delta angularFreq radialFreq sigma periodR2 periodEnv =
  let c = ((-1) ^ (abs angularFreq)) :+ 0
  in R.map (* c) $
     analyticalFourierSeries1 numFreqs delta angularFreq radialFreq sigma periodR2 periodEnv

-- Analytical solution type 2
-- The pinwheel harmonics are cis (freq * x)
{-# INLINE analyticalFourierSeriesFunc2 #-}
analyticalFourierSeriesFunc2 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e))
  => Int
  -> Int
  -> e
  -> e
  -> e
  -> e
  -> e
  -> Complex e
analyticalFourierSeriesFunc2 angularFreq radialFreq sigma periodR2 periodEnv phi rho =
  let radialConst = 2 * pi / (log periodEnv)
  in pi * ((0 :+ 1) ^ (abs angularFreq)) *
     ((cis (fromIntegral angularFreq * phi))) *
     ((periodR2 / (pi * rho) :+ 0) **
      ((2 + sigma) :+ (radialConst * fromIntegral radialFreq))) *
     (gamma $
      ((2 + fromIntegral (abs angularFreq) + sigma) :+
       (radialConst * fromIntegral radialFreq)) /
      2) /
     (gamma $
      ((fromIntegral (abs angularFreq) - sigma) :+
       (radialConst * fromIntegral (-radialFreq))) /
      2)

-- The pinwheel is in the frequency domain, the function computes the Fourier series
{-# INLINE analyticalFourierSeries2 #-}
analyticalFourierSeries2 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierSeries2 numPoints delta angularFreq radialFreq sigma periodR2 periodEnv =
  let center = div numPoints 2
      arr =
        fromFunction (Z :. numPoints :. numPoints) $ \(Z :. i :. j) ->
          let x = delta * (fromIntegral $ i - center)
              y = delta * (fromIntegral $ j - center)
              rho = sqrt $ x ^ 2 + y ^ 2
              phi = atan2 y x
          in if rho == 0 -- || (rho > period / 2) -- || phi /= 0
               then 0
               else analyticalFourierSeriesFunc2
                      angularFreq
                      radialFreq
                      sigma
                      periodR2
                      periodEnv
                      phi
                      rho
     -- if angularFreq == 0
     --   then fromFunction (Z :. numPoints :. numPoints) $ \_ -> 0
     --   else
  in arr

-- The pinwheel is in the spatial domain, the function computes the Fourier coefficients
{-# INLINE analyticalFourierCoefficients2 #-}
analyticalFourierCoefficients2 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierCoefficients2 numFreqs delta angularFreq radialFreq sigma periodR2 periodEnv =
  let c = ((-1) ^ (abs angularFreq)) / periodR2 :+ 0
  in R.map (* c) $
     analyticalFourierSeries2
       numFreqs
       delta
       angularFreq
       radialFreq
       sigma
       periodR2
       periodEnv

-- Analytical solution type 3
-- The pinwheel harmonics are cis (-freq * x)
-- The envelope is r^2 * r^\alpha, where \alpha \in (-2,-0.5)
{-# INLINE analyticalFourierSeriesFunc3 #-}
analyticalFourierSeriesFunc3 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e))
  => Int
  -> Int
  -> e
  -> e
  -> e
  -> e
  -> e
  -> Complex e
analyticalFourierSeriesFunc3 angularFreq radialFreq sigma periodR2 periodEnv phi rho =
  let radialConst = 2 * pi / log periodEnv
   in ((fromIntegral angularFreq ^ 2 :+ 0) -
       ((sigma + 2) :+ (radialConst * fromIntegral (-radialFreq))) ^ 2) /
      (8 * pi) *
      ((0 :+ 1) ^ abs angularFreq) *
      cis (fromIntegral (-angularFreq) * phi) /
      (rho ^ 2 :+ 0) *
      ((periodR2 / (pi * rho) :+ 0) **
       ((2 + sigma) :+ (radialConst * fromIntegral (-radialFreq)))) *
      gamma
        (((2 + fromIntegral (abs angularFreq) + sigma) :+
          (radialConst * fromIntegral (-radialFreq))) /
         2) /
      gamma
        (((fromIntegral (abs angularFreq) - sigma) :+
          (radialConst * fromIntegral radialFreq)) /
         2)

{-# INLINE analyticalFourierSeries3 #-}
analyticalFourierSeries3 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierSeries3 numPoints delta angularFreq radialFreq sigma periodR2 periodEnv =
  let center = div numPoints 2
  in fromFunction (Z :. numPoints :. numPoints) $ \(Z :. i :. j) ->
       let x = fromIntegral $ i - center
           y = fromIntegral $ j - center
           rho = sqrt $ x ^ 2 + y ^ 2
           phi = atan2 y x
       in if rho == 0
            then 0
            else analyticalFourierSeriesFunc3
                   angularFreq
                   radialFreq
                   sigma
                   periodR2
                   periodEnv
                   phi
                   (rho * delta)

{-# INLINE analyticalFourierCoefficients3 #-}
analyticalFourierCoefficients3 ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierCoefficients3 numFreqs delta angularFreq radialFreq sigma periodR2 periodEnv =
  let c = ((-1) ^ (abs angularFreq)) :+ 0
  in R.map (* c) $
     analyticalFourierSeries3 numFreqs delta angularFreq radialFreq sigma periodR2 periodEnv

pinwheelFourierCoefficientsAnatical ::
     ( Unbox e
     , RealFloat e
     , Gamma (Complex e)
     , VG.Vector vector (Complex e)
     , NFData (vector (Complex e))
     )
  => Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> IA.Array (Int, Int) (vector (Complex e))
pinwheelFourierCoefficientsAnatical numR2Freqs phiFreq rhoFreq thetaFreq rFreq sigma periodR2 periodEnv =
  let idxs =
        [ (radialFreq, angularFreq)
        | radialFreq <- pinwheelFreqs rhoFreq rFreq
        , angularFreq <- pinwheelFreqs phiFreq thetaFreq
        ]
      centerR2Freq = div numR2Freqs 2
      (radialLB, radialUB) = pinwheelFreqsBound rhoFreq rFreq
      (angularLB, angularUB) = pinwheelFreqsBound phiFreq thetaFreq
      pinwheels =
        parMap
          rdeepseq
          (\(radialFreq, angularFreq) ->
             VG.convert .
             toUnboxed .
             computeS -- .
             -- R.zipWith
             --   (*)
             --   (fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i :. j) ->
             --      gaussian2D
             --        (fromIntegral $ i - centerR2Freq)
             --        (fromIntegral $ j - centerR2Freq)
             --        (fromIntegral $ div centerR2Freq 1))
            $
             analyticalFourierCoefficients2
               numR2Freqs
               1
               (angularFreq)
               (radialFreq)
               sigma
               periodR2
               periodEnv)
          idxs
  in listArray ((radialLB, angularLB), (radialUB, angularUB)) pinwheels

pinwheelFourierCoefficientsAnaticalList ::
     ( Unbox e
     , RealFloat e
     , Gamma (Complex e)
     , VG.Vector vector (Complex e)
     , NFData (vector (Complex e))
     )
  => Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> [vector (Complex e)]
pinwheelFourierCoefficientsAnaticalList numR2Freqs phiFreq rhoFreq thetaFreq rFreq sigma periodR2 periodEnv =
  let idxs =
        [ (radialFreq, angularFreq)
        | radialFreq <- [-rhoFreq .. rhoFreq]
        , angularFreq <- [-phiFreq .. phiFreq]
        ]
      centerR2Freq = div numR2Freqs 2
      pinwheels =
        parMap
          rdeepseq
          (\(radialFreq, angularFreq) ->
             VG.convert . toUnboxed . computeS $
             analyticalFourierCoefficients2
               numR2Freqs
               1
               (angularFreq)
               (radialFreq)
               sigma
               periodR2
               periodEnv)
          idxs
  in pinwheels

pinwheelFourierCoefficientsAnaticalList1 ::
     ( Unbox e
     , RealFloat e
     , Gamma (Complex e)
     , VG.Vector vector (Complex e)
     , NFData (vector (Complex e))
     )
  => Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> [vector (Complex e)]
pinwheelFourierCoefficientsAnaticalList1 numR2Freqs phiFreq rhoFreq thetaFreq rFreq sigma periodR2 periodEnv =
  let idxs =
        [ (radialFreq, angularFreq)
        | radialFreq <- [-rhoFreq .. rhoFreq]
        , angularFreq <- [-phiFreq .. phiFreq]
        ]
      pinwheels =
        parMap
          rdeepseq
          (\(radialFreq, angularFreq) ->
             VG.convert . toUnboxed . computeS $
             analyticalFourierCoefficients2
               numR2Freqs
               1
               (angularFreq)
               (radialFreq)
               sigma
               periodR2
               periodEnv)
          idxs
  in pinwheels

pinwheelFourierCoefficientsAnaticalList2 ::
     ( Unbox e
     , RealFloat e
     , Gamma (Complex e)
     , VG.Vector vector (Complex e)
     , NFData (vector (Complex e))
     )
  => Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> [vector (Complex e)]
pinwheelFourierCoefficientsAnaticalList2 numR2Freqs phiFreq rhoFreq thetaFreq rFreq sigma periodR2 periodEnv =
  let idxs =
        [ (radialFreq, angularFreq)
        | radialFreq <- [rFreq,(rFreq - 1) .. -rFreq]
        , angularFreq <- [thetaFreq,(thetaFreq - 1) .. -thetaFreq]
        ]
      centerR2Freq = div numR2Freqs 2
      radialConstant = 2 * pi / (log periodEnv)
      pinwheels =
        parMap
          rdeepseq
          (\(radialFreq, angularFreq) ->
             VG.convert . toUnboxed . computeS $
             fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. x :. y) ->
               let xFreq = fromIntegral $ x - centerR2Freq
                   yFreq = fromIntegral $ y - centerR2Freq
                   rho = sqrt $ xFreq ^ 2 + yFreq ^ 2
                   phi = atan2 yFreq xFreq
               in if rho == 0
                    then 0
                    else ((0 :+ (-1)) ^ 1) *
                         ((periodR2 / (pi * rho) :+ 0) **
                          (0 :+ fromIntegral radialFreq * radialConstant)) *
                         (cis $ phi * fromIntegral angularFreq))
          idxs
  in pinwheels

{-# INLINE analyticalFourierSeries2' #-}
analyticalFourierSeries2' ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierSeries2' numPoints _ angularFreq radialFreq sigma periodR2 periodEnv =
  let center = div numPoints 2
      arr =
        fromFunction (Z :. numPoints :. numPoints) $ \(Z :. i :. j) ->
          let x' = i - center
              y' = j - center
          in if x' == 0 && y' == 0
               then 0
               else let x =
                          if x' == 0
                            then 0
                            else if x' > 0
                                   then 1 / (fromIntegral $ center + 1 - x')
                                   else (-1) / (fromIntegral $ center + 1 + x')
                        y =
                          if y' == 0
                            then 0
                            else if y' > 0
                                   then 1 / (fromIntegral $ center + 1 - y')
                                   else (-1) / (fromIntegral $ center + 1 + y')
                        rho = sqrt $ x ^ 2 + y ^ 2
                        phi = atan2 y x
                    in analyticalFourierSeriesFunc2
                         angularFreq
                         radialFreq
                         sigma
                         periodR2
                         periodEnv
                         phi
                         rho
  in arr

{-# INLINE analyticalFourierCoefficients2' #-}
analyticalFourierCoefficients2' ::
     (Eq e, Fractional e, RealFloat e, Gamma (Complex e), Unbox e)
  => Int
  -> e
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> R.Array D DIM2 (Complex e)
analyticalFourierCoefficients2' numFreqs delta angularFreq radialFreq sigma periodR2 periodEnv =
  let c = (-1) ^ (abs angularFreq) :+ 0
  in R.map (* c) $
     analyticalFourierSeries2'
       numFreqs
       delta
       angularFreq
       radialFreq
       sigma
       periodR2
       periodEnv

pinwheelFourierCoefficientsAnatical' ::
     ( Unbox e
     , RealFloat e
     , Gamma (Complex e)
     , VG.Vector vector (Complex e)
     , NFData (vector (Complex e))
     )
  => Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> e
  -> e
  -> IA.Array (Int, Int) (vector (Complex e))
pinwheelFourierCoefficientsAnatical' numR2Freqs phiFreq rhoFreq thetaFreq rFreq sigma periodR2 periodEnv =
  let idxs =
        [ (radialFreq, angularFreq)
        | radialFreq <- pinwheelFreqs rhoFreq rFreq
        , angularFreq <- pinwheelFreqs phiFreq thetaFreq
        ]
      centerR2Freq = div numR2Freqs 2
      (radialLB, radialUB) = pinwheelFreqsBound rhoFreq rFreq
      (angularLB, angularUB) = pinwheelFreqsBound phiFreq thetaFreq
      pinwheels =
        parMap
          rdeepseq
          (\(radialFreq, angularFreq) ->
             VG.convert .
             toUnboxed .
             computeS -- .
             -- R.zipWith
             --   (*)
             --   (fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i :. j) ->
             --      gaussian2D
             --        (fromIntegral $ i - centerR2Freq)
             --        (fromIntegral $ j - centerR2Freq)
             --        (fromIntegral $ div centerR2Freq 2))
            $
             analyticalFourierCoefficients2'
               numR2Freqs
               1
               (angularFreq)
               (radialFreq)
               sigma
               periodR2
               periodEnv)
          idxs
  in listArray ((radialLB, angularLB), (radialUB, angularUB)) pinwheels

{-# INLINE idealHighPassFilter #-}
idealHighPassFilter ::
     (VG.Vector vector (Complex Double))
  => Double
  -> Int
  -> vector (Complex Double)
idealHighPassFilter radius numR2Freq =
  let a = 1 / (radius * 2)
      r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
      sinc x =
        if x == 0
          then 1
          else (sin (pi * x)) / (pi * x)
  in VG.fromList
       [ if xFreq == 0 && yFreq == 0
         then (1 - 1 / a ^ 2) :+ 0
         else ((-1) * (sinc (xFreq / a)) * (sinc (yFreq / a)) / (a ^ 2)) :+ 0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE idealHighPassFilter1 #-}
idealHighPassFilter1 :: Int -> VS.Vector (Complex Double)
idealHighPassFilter1 numR2Freq =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
  in VS.fromList
       [ if xFreq == 0 && yFreq == 0
         then 0
         else (-1) :+ 0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE gaussianHighPassFilter #-}
gaussianHighPassFilter :: Double -> Int -> VS.Vector (Complex Double)
gaussianHighPassFilter alpha numR2Freq =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
  in VS.fromList
       [ if xFreq == 0 && yFreq == 0
         then (1 - (pi / alpha)) :+ 0
         else ((-1) * (pi / alpha) *
               exp ((pi ^ 2 * (xFreq ^ 2 + yFreq ^ 2)) / (-alpha))) :+
              0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE gaussianHighPassFilter1 #-}
gaussianHighPassFilter1 :: Double -> Double -> Int -> VS.Vector (Complex Double)
gaussianHighPassFilter1 radius alpha numR2Freq =
  let a = 1 / (radius * 2)
      r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
      sinc x =
        if x == 0
          then 1
          else (sin (pi * x)) / (pi * x)
  in VS.fromList
       [ (sinc (xFreq / a) * sinc (yFreq / a) / (a ^ 2) -
          (pi / alpha) * exp ((pi ^ 2 * (xFreq ^ 2 + yFreq ^ 2)) / (-alpha))) :+
       0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE idealLowPassFilter #-}
idealLowPassFilter ::
     (VG.Vector vector (Complex Double))
  => Double
  -> Double
  -> Int
  -> vector (Complex Double)
idealLowPassFilter radius periodR2 numR2Freq =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
   in VG.fromList
        [ let rho = 2 * pi * sqrt (xFreq ^ 2 + yFreq ^ 2) / periodR2
           in if rho == 0
                then 0
                else radius / rho  * bessel_J1 (radius * rho) / periodR2 / (2 * pi)^2 :+ 0
        | xFreq <- r2Freqs
        , yFreq <- r2Freqs
        ]

{-# INLINE gaussianLowPassFilter #-}
gaussianLowPassFilter :: Double -> Int -> VS.Vector (Complex Double)
gaussianLowPassFilter std numR2Freq =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
  in VS.fromList
       [ ((std ^ 2) * exp ((-pi) * (xFreq ^ 2 + yFreq ^ 2) * (std ^ 2))) :+
       0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE laplacianLowPassFilter #-}
laplacianLowPassFilter :: Double -> Int -> VS.Vector (Complex Double)
laplacianLowPassFilter a numR2Freq =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
  in VS.fromList
       [ let s2 = xFreq ^ 2 + yFreq ^ 2
         in if xFreq == 0 && yFreq == 0
              then (1 - a / ((4 * pi * pi * s2 + a ^ 2) ** 1.5)) :+ 0
              else (-a) / ((4 * pi * pi * s2 + a ^ 2) ** 1.5) :+ 0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE laplacianHighPassFilter #-}
laplacianHighPassFilter :: Int -> VS.Vector (Complex Double)
laplacianHighPassFilter numR2Freq =
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
  in VS.fromList
       [ if xFreq == 0 && yFreq == 0
         then 8 -- -3.33
         else if xFreq == 0 && abs yFreq == 1
                then -1 -- 0.67
                else if abs xFreq == 1 && yFreq == 0
                       then -1 -- 0.67
                       else if abs xFreq == 1 && abs yFreq == 1
                              then -1 -- 0.17
                              else 0
       | xFreq <- r2Freqs
       , yFreq <- r2Freqs
       ]

{-# INLINE convolveFrequency #-}
convolveFrequency ::
     (VG.Vector vector (Complex Double))
  => DFTPlan
  -> Int
  -> vector (Complex Double)
  -> IA.Array (Int, Int) (vector (Complex Double))
  -> IO (IA.Array (Int, Int) (vector (Complex Double)))
convolveFrequency plan numR2Freq filter' arr = do
  let dftPlanID = DFTPlanID DFT1DG [numR2Freq, numR2Freq] [0, 1]
      idftPlanID = DFTPlanID IDFT1DG [numR2Freq, numR2Freq] [0, 1]
      vecs = L.map VG.convert . IA.elems $ arr
      filter =
        VU.convert .
        toUnboxed .
        computeUnboxedS .
        makeFilter2D . fromUnboxed (Z :. numR2Freq :. numR2Freq) . VG.convert $
        filter'
  dftF <- dftExecute plan dftPlanID filter
  vecsF <- mapConcurrently (dftExecute plan dftPlanID) vecs
  outputs <-
    mapConcurrently (dftExecute plan idftPlanID . VS.zipWith (*) dftF) vecsF
  return . listArray (bounds arr) . L.map VG.convert $ outputs

{-# INLINE convolveFrequency1 #-}
convolveFrequency1 ::
     (VG.Vector vector (Complex Double))
  => DFTPlan
  -> Int
  -> vector (Complex Double)
  -> IA.Array (Int, Int) (vector (Complex Double))
  -> IO (IA.Array (Int, Int) (vector (Complex Double)))
convolveFrequency1 plan numR2Freq filter' arr = do
  let dftPlanID = DFTPlanID DFT1DG [numR2Freq, numR2Freq] [0, 1]
      idftPlanID = DFTPlanID IDFT1DG [numR2Freq, numR2Freq] [0, 1]
      vecs = L.map VG.convert . IA.elems $ arr
      filter =
        VU.convert .
        toUnboxed .
        computeUnboxedS .
        makeFilter2D . fromUnboxed (Z :. numR2Freq :. numR2Freq) . VG.convert $
        filter'
  dftF <- dftExecute plan dftPlanID filter
  vecsF <- mapConcurrently (dftExecute plan dftPlanID) vecs
  outputs <-
    mapConcurrently (dftExecute plan idftPlanID . VS.zipWith (*) dftF) vecsF
  return .
    listArray (bounds arr) . L.map VG.convert . L.zipWith (VG.zipWith (-)) vecs $
    outputs

{-# INLINE centerHollow #-}
centerHollow ::
     (VG.Vector vector (Complex Double))
  => Int
  -> IA.Array (Int, Int) (vector (Complex Double))
  -> IA.Array (Int, Int) (vector (Complex Double))
centerHollow numR2Freq arr = centerHollowVector numR2Freq <$> arr

{-# INLINE centerHollowVector #-}
centerHollowVector ::
     (VG.Vector vector (Complex Double))
  => Int
  -> vector (Complex Double)
  -> vector (Complex Double)
centerHollowVector numR2Freq vec =
  let s = VG.sum vec / (fromIntegral (numR2Freq ^ 2) :+ 0)
   in VG.map (\x -> x - s) vec

{-# INLINE centerHollowArray #-}
centerHollowArray ::
     (R.Source s (Complex e), Unbox e, RealFloat e)
  => Int
  -> R.Array s DIM2 (Complex e)
  -> R.Array D DIM2 (Complex e)
centerHollowArray numR2Freq arr =
  let s = sumAllS arr / (fromIntegral (numR2Freq ^ 2) :+ 0)
   in R.map (\x -> x - s) $ arr

{-# INLINE centerHollowArray' #-}
centerHollowArray' ::
     (R.Source s (Complex e), Unbox e, RealFloat e)
  => Int
  -> R.Array s DIM2 (Complex e)
  -> R.Array D DIM2 (Complex e)
centerHollowArray' numR2Freq arr =
  let s = sumAllS arr
      c = div numR2Freq 2
   in R.traverse arr id $ \f idx@(Z :. i :. j) ->
        if i == c && j == c
          then (-s)
          else f idx

envelopIntegral ::
     Double -> Double -> Double -> Double -> Double -> Int -> Complex Double
envelopIntegral a b delta s period radialFreq =
  let m = round $ (b - a) / delta
      n =
        if odd m
          then m
          else m - 1
      -- weights = VU.fromList $ weightsSimpsonRule n
      vec =
        VU.generate n $ \i ->
          let x = (a + fromIntegral i * delta)
          -- in   ((x :+ 0) ** (s :+ freq)) * ((pi :+ 0) ** (0 :+ freq)) *
          --      (gamma ((0.5 :+ 0) * (1 :+ (-freq)))) /
          --      (gamma ((0.5 :+ 0) * (1 :+ freq)))
          in analyticalFourierSeriesFunc2 0 radialFreq s period period 0 x
  -- in (delta / 3 :+ 0) * (VU.sum . VU.zipWith (*) vec $ weights)
  in VU.sum vec

envelopIntegral2D ::
     Int -> Double -> Double -> Double -> Int -> Complex Double
envelopIntegral2D numPoints delta s period radialFreq =
  let center = div numPoints 2
      -- weights = VU.fromList $ weightsSimpsonRule n
      arr = analyticalFourierSeries2 numPoints delta 0 radialFreq s period period
  in sumAllS arr

printEnvelopIntegral ::
     Double -> Double -> Double -> Double -> Double -> Int -> IO ()
printEnvelopIntegral a b delta s period radialFreq =
  let m = round $ (b - a) / delta
      n =
        if odd m
          then m
          else m - 1
      weights = VU.fromList $ weightsSimpsonRule n
      vec =
        VU.generate n $ \i ->
          let x = a + fromIntegral i * delta
          in analyticalFourierSeriesFunc2 0 radialFreq s period period 0 x
  in print $ VU.zipWith (*) weights vec

{-# INLINE createFrequencyArray #-}
createFrequencyArray ::
     (RealFloat e) => Int -> (e -> e -> Complex e) -> R.Array D DIM2 (Complex e)
createFrequencyArray numR2Freqs f =
  fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i :. j) ->
    let x = fromIntegral $ i - div numR2Freqs 2
        y = fromIntegral $ j - div numR2Freqs 2
        rho = sqrt $ x ^ 2 + y ^ 2
        phi = atan2 y x
     in f phi rho

{-# INLINE createFrequencyArrayZeroCenter #-}
createFrequencyArrayZeroCenter ::
     (RealFloat e) => Int -> (e -> e -> Complex e) -> R.Array D DIM2 (Complex e)
createFrequencyArrayZeroCenter numR2Freqs f =
  fromFunction (Z :. numR2Freqs :. numR2Freqs) $ \(Z :. i :. j) ->
    let x = fromIntegral $ i - div numR2Freqs 2
        y = fromIntegral $ j - div numR2Freqs 2
        rho = sqrt $ x ^ 2 + y ^ 2
        phi = atan2 y x
     in if rho == 0
          then 0
          else f phi rho

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field