text:
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MagicHash                  #-}
{-# LANGUAGE MultiParamTypeClasses      #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE TypeInType                 #-}
{-# LANGUAGE UnboxedTuples              #-}
{-# LANGUAGE UndecidableInstances       #-}
{-# LANGUAGE ViewPatterns               #-}
{-# OPTIONS_GHC -fno-warn-orphans       #-}
module Numeric.Quaternion.Internal.QDouble
    ( QDouble, Quater (..)
    ) where

import qualified Control.Monad.ST                     as ST
import           Data.Coerce                          (coerce)
import           Numeric.Basics
import           Numeric.DataFrame.Internal.PrimArray
import qualified Numeric.DataFrame.ST                 as ST
import           Numeric.DataFrame.Type
import           Numeric.Quaternion.Internal
import           Numeric.Vector.Internal

type QDouble = Quater Double

deriving instance PrimBytes (Quater Double)
deriving instance PrimArray Double (Quater Double)

instance Quaternion Double where
    newtype Quater Double = QDouble (Vector Double 4)
    {-# INLINE packQ #-}
    packQ = coerce (vec4 :: Double -> Double -> Double -> Double -> Vector Double 4)
    {-# INLINE unpackQ# #-}
    unpackQ# = coerce (unpackV4# :: Vector Double 4 -> (# Double, Double, Double, Double #))
    {-# INLINE fromVecNum #-}
    fromVecNum (unpackV3# -> (# x, y, z #))  = packQ x y z
    {-# INLINE fromVec4 #-}
    fromVec4 = coerce
    {-# INLINE toVec4 #-}
    toVec4 = coerce
    {-# INLINE square #-}
    square (unpackQ# -> (# x, y, z, w #)) = (x * x) + (y * y) + (z * z) + (w * w)
    {-# INLINE im #-}
    im (unpackQ# -> (# x, y, z, _ #)) = packQ x y z 0.0
    {-# INLINE re #-}
    re (unpackQ# -> (# _, _, _, w #)) = packQ 0 0 0 w
    {-# INLINE imVec #-}
    imVec (unpackQ# -> (# x, y, z, _ #)) = vec3 x y z
    {-# INLINE taker #-}
    taker (unpackQ# -> (# _, _, _, w #)) = w
    {-# INLINE takei #-}
    takei (unpackQ# -> (# x, _, _, _ #)) = x
    {-# INLINE takej #-}
    takej (unpackQ# -> (# _, y, _, _ #)) = y
    {-# INLINE takek #-}
    takek (unpackQ# -> (# _, _, z, _ #)) = z
    {-# INLINE conjugate #-}
    conjugate (unpackQ# -> (# x, y, z, w #))
      = packQ (negate x) (negate y) (negate z) w
    {-# INLINE rotScale #-}
    rotScale (unpackQ# -> (# i, j, k, t #))
             (unpackV3# -> (# x, y, z #))
      = let l = t*t - i*i - j*j - k*k
            d = 2.0 * ( i*x + j*y + k*z)
            t2 = t * 2.0
        in vec3 (l*x + d*i + t2 * (z*j - y*k))
                (l*y + d*j + t2 * (x*k - z*i))
                (l*z + d*k + t2 * (y*i - x*j))
    {-# INLINE getRotScale #-}
    getRotScale a b = case (# unpackV3# a, unpackV3# b #) of
      (# _, (# 0, 0, 0 #) #) -> packQ 0 0 0 0
      (# (# 0, 0, 0 #), _ #) -> let x = (1 / 0 :: Double) in packQ x x x x
      (# (# a1, a2, a3 #), (# b1, b2, b3 #) #) ->
        let ma = sqrt (a1*a1 + a2*a2 + a3*a3)
            mb = sqrt (b1*b1 + b2*b2 + b3*b3)
            d  = a1*b1 + a2*b2 + a3*b3
            c  = sqrt (ma*mb + d)
            ma2 = ma * 1.4142135623730951 -- sqrt 2.0
            r  = recip (ma2 * c)
            c' = sqrt (mb / ma) -- ratio of a and b for corner cases
            r' = recip (sqrt ( negate (a1*b1 + a2*b2) ))
        in case unpackV3# (cross a b) of
          (# 0, 0, 0 #)
              -- if a and b face the same direction, q is fully real
            | d >= 0       -> packQ 0 0 0 c'
              -- if a and b face opposite directions, find an orthogonal vector
              -- prerequisites: w == 0  and  a·(x,y,z) == 0
              -- corner cases: only one vector component is non-zero
            | b1 == 0      -> packQ c' 0 0 0
              -- otherwise set the last component to zero,
              -- and get an orthogonal vector in 2D.
            | otherwise    -> packQ (-b2*r') (b1*r') 0 0
              -- NB: here we have some precision troubles
              --     when a and b are close to parallel and opposite.
          (# t1, t2, t3 #) -> packQ (t1 * r) (t2 * r) (t3 * r) (c / ma2)
    {-# INLINE axisRotation #-}
    axisRotation v a = case unpackV3# v of
      (# 0, 0, 0 #) -> packQ 0 0 0 (negateUnless (abs a < M_PI) 1)
      (# x, y, z #) ->
        let c = cos (a * 0.5)
            s = sin (a * 0.5)
                / sqrt (x*x + y*y + z*z)
        in packQ (x * s) (y * s) (z * s) c
    {-# INLINE qArg #-}
    qArg (unpackQ# -> (# x, y, z, w #)) = 2 * atan2 (sqrt (x*x + y*y + z*z)) w
    {-# INLINE fromMatrix33 #-}
    fromMatrix33 m = fromM 1
      (ix# 0# m) (ix# 1# m) (ix# 2# m)
      (ix# 3# m) (ix# 4# m) (ix# 5# m)
      (ix# 6# m) (ix# 7# m) (ix# 8# m)

    {-# INLINE fromMatrix44 #-}
    fromMatrix44 m = fromM (ix# 15# m)
      (ix# 0# m) (ix# 1# m) (ix# 2# m)
      (ix# 4# m) (ix# 5# m) (ix# 6# m)
      (ix# 8# m) (ix# 9# m) (ix# 10# m)

    {-# INLINE toMatrix33 #-}
    toMatrix33 (unpackQ# -> (# 0.0, 0.0, 0.0, w #))
      = let x = w * w
            f 0 = (# 3 :: Int , x #)
            f k = (# k-1, 0 #)
        in case gen# (CumulDims [9,3,1]) f 0 of
            (# _, m #) -> m -- diag (scalar (w * w))
    toMatrix33 (unpackQ# -> (# x', y', z', w' #)) =
      let x = scalar x'
          y = scalar y'
          z = scalar z'
          w = scalar w'
          x2 = x * x
          y2 = y * y
          z2 = z * z
          w2 = w * w
          l2 = x2 + y2 + z2 + w2
      in ST.runST $ do
        df <- ST.newDataFrame
        ST.writeDataFrameOff df 0 $ l2 - 2*(z2 + y2)
        ST.writeDataFrameOff df 1 $ 2*(x*y + z*w)
        ST.writeDataFrameOff df 2 $ 2*(x*z - y*w)
        ST.writeDataFrameOff df 3 $ 2*(x*y - z*w)
        ST.writeDataFrameOff df 4 $ l2 - 2*(z2 + x2)
        ST.writeDataFrameOff df 5 $ 2*(y*z + x*w)
        ST.writeDataFrameOff df 6 $ 2*(x*z + y*w)
        ST.writeDataFrameOff df 7 $ 2*(y*z - x*w)
        ST.writeDataFrameOff df 8 $ l2 - 2*(y2 + x2)
        ST.unsafeFreezeDataFrame df
    {-# INLINE toMatrix44 #-}
    toMatrix44 (unpackQ# -> (# 0.0, 0.0, 0.0, w #)) = ST.runST $ do
      df <- ST.newDataFrame
      mapM_ (flip (ST.writeDataFrameOff df) 0) [0..15]
      let w2 = scalar (w * w)
      ST.writeDataFrameOff df 0 w2
      ST.writeDataFrameOff df 5 w2
      ST.writeDataFrameOff df 10 w2
      ST.writeDataFrameOff df 15 1
      ST.unsafeFreezeDataFrame df
    toMatrix44 (unpackQ# -> (# x', y', z', w' #)) =
      let x = scalar x'
          y = scalar y'
          z = scalar z'
          w = scalar w'
          x2 = x * x
          y2 = y * y
          z2 = z * z
          w2 = w * w
          l2 = x2 + y2 + z2 + w2
      in ST.runST $ do
        df <- ST.newDataFrame
        ST.writeDataFrameOff df 0 $ l2 - 2*(z2 + y2)
        ST.writeDataFrameOff df 1 $ 2*(x*y + z*w)
        ST.writeDataFrameOff df 2 $ 2*(x*z - y*w)
        ST.writeDataFrameOff df 3 0
        ST.writeDataFrameOff df 4 $ 2*(x*y - z*w)
        ST.writeDataFrameOff df 5 $ l2 - 2*(z2 + x2)
        ST.writeDataFrameOff df 6 $ 2*(y*z + x*w)
        ST.writeDataFrameOff df 7 0
        ST.writeDataFrameOff df 8 $ 2*(x*z + y*w)
        ST.writeDataFrameOff df 9 $ 2*(y*z - x*w)
        ST.writeDataFrameOff df 10 $ l2 - 2*(y2 + x2)
        ST.writeDataFrameOff df 11 0
        ST.writeDataFrameOff df 12 0
        ST.writeDataFrameOff df 13 0
        ST.writeDataFrameOff df 14 0
        ST.writeDataFrameOff df 15 1
        ST.unsafeFreezeDataFrame df

{- Calculate quaternion from a 3x3 matrix.

   First argument is a constant; it is either 1 for a 3x3 matrix,
   or m44 for a 4x4 matrix. I just need to multiply all components by
   this number.

   Further NB for the formulae:

   d == square q == det m ** (1/3)
   t == trace m  == 4 w w - d
   m01 - m10 == 4 z w
   m20 - m02 == 4 y w
   m12 - m21 == 4 x w
   m01 + m10 == 4 x y
   m20 + m02 == 4 x z
   m12 + m21 == 4 y z
   m00 == + x x - y y - z z + w w
   m11 == - x x + y y - z z + w w
   m22 == - x x - y y + z z + w w
   4 x x == d + m00 - m11 - m22
   4 y y == d - m00 + m11 - m22
   4 z z == d - m00 - m11 + m22
   4 w w == d + m00 + m11 + m22
 -}
fromM :: Double
      -> Double -> Double -> Double
      -> Double -> Double -> Double
      -> Double -> Double -> Double
      -> QDouble
fromM c'
  m00 m01 m02
  m10 m11 m12
  m20 m21 m22
    | t > 0
      = let dd = sqrt ( d + t )
            is = c / dd
        in packQ ((m12 - m21)*is) ((m20 - m02)*is) ((m01 - m10)*is) (c*dd)
    | m00 > m11 && m00 > m22
      = let dd = sqrt ( d + m00 - m11 - m22 )
            is = c / dd
        in packQ (c*dd) ((m01 + m10)*is) ((m02 + m20)*is) ((m12 - m21)*is)
    | m11 > m22
      = let dd = sqrt ( d - m00 + m11 - m22 )
            is = c / dd
        in packQ ((m01 + m10)*is) (c*dd) ((m12 + m21)*is) ((m20 - m02)*is)
    | otherwise
      = let dd = sqrt ( d - m00 - m11 + m22 )
            is = c / dd
        in packQ ((m02 + m20)*is) ((m12 + m21)*is) (c*dd) ((m01 - m10)*is)

  where
    -- normalizing constant
    c = recip $ 2 * sqrt c'
    -- trace
    t = m00 + m11 + m22
    -- cubic root of determinant
    d = ( m00 * ( m11 * m22 - m12 * m21 )
        - m01 * ( m10 * m22 - m12 * m20 )
        + m02 * ( m10 * m21 - m11 * m20 )
        ) ** 0.33333333333333333333333333333333

instance Num QDouble where
    QDouble a + QDouble b
      = QDouble (a + b)
    {-# INLINE (+) #-}
    QDouble a - QDouble b
      = QDouble (a - b)
    {-# INLINE (-) #-}
    (unpackQ# -> (# a1, a2, a3, a4 #)) * (unpackQ# -> (# b1, b2, b3, b4 #))
      = packQ ((a4 * b1) + (a1 * b4) + (a2 * b3) - (a3 * b2))
              ((a4 * b2) - (a1 * b3) + (a2 * b4) + (a3 * b1))
              ((a4 * b3) + (a1 * b2) - (a2 * b1) + (a3 * b4))
              ((a4 * b4) - (a1 * b1) - (a2 * b2) - (a3 * b3))
    {-# INLINE (*) #-}
    negate (QDouble a) = QDouble (negate a)
    {-# INLINE negate #-}
    abs = packQ 0 0 0 . sqrt . square
    {-# INLINE abs #-}
    signum q@(unpackQ# -> (# x, y, z, w #))
      | qd == 0   = q
      | otherwise = case ix + iy + iz + iw + nn of
        0 -> packQ (x * l) (y * l) (z * l) (w * l)
        1 -> packQ (copysign 1 x) 0 0 0
        2 -> packQ 0 (copysign 1 y) 0 0
        4 -> packQ 0 0 (copysign 1 z) 0
        8 -> packQ 0 0 0 (copysign 1 w)
        _ -> packQ n n n n
      where
        n  = 0 / 0 :: Double
        qd = x*x + y*y + z*z + w*w
        ix = if isInfinite x then 1 else 0 :: Int
        iy = if isInfinite y then 2 else 0 :: Int
        iz = if isInfinite z then 4 else 0 :: Int
        iw = if isInfinite w then 8 else 0 :: Int
        nn = if isNaN x || isNaN y || isNaN z || isNaN w then 16 else 0 :: Int
        l  = recip (sqrt qd)
    {-# INLINE signum #-}
    fromInteger = packQ 0 0 0 . fromInteger
    {-# INLINE fromInteger #-}

instance Fractional QDouble where
    {-# INLINE recip #-}
    recip q@(unpackQ# -> (# x, y, z, w #)) = case negate (recip (square q)) of
      c -> packQ (x * c) (y * c) (z * c) (negate (w * c))
    {-# INLINE (/) #-}
    a / b = a * recip b
    {-# INLINE fromRational #-}
    fromRational = packQ 0 0 0 . fromRational

instance Floating QDouble where
    {-# INLINE pi #-}
    pi = packQ 0 0 0 M_PI
    {-# INLINE exp #-}
    exp (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0  = packQ x y z ew
      | otherwise = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv  = sqrt mv2
        ew  = exp w
        l   = ew * sin mv / mv
        arg = ew * cos mv
    {-# INLINE log #-}
    log = log' (Vec3 1 0 0)
    {-# INLINE sqrt #-}
    sqrt = sqrt' (Vec3 1 0 0)
    {-# INLINE sin #-}
    sin (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0  = packQ x y z (sin w)
      | otherwise = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv  = sqrt mv2
        l   = cos w * sinh mv / mv
        arg = sin w * cosh mv
    {-# INLINE cos #-}
    cos (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0  = packQ x y z (cos w)
      | otherwise = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv  = sqrt mv2
        l   = sin w * sinh mv / negate mv
        arg = cos w * cosh mv
    {-# INLINE tan #-}
    tan (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0       = packQ x y z (tan w)
      | isInfinite mv2 = signum (packQ x y z 0)
      | otherwise      = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv = sqrt mv2
        b = 2*mv
        a = 2*w
        sina = sin a
        eb = exp (-b)
        eb2 = eb*eb
        d = 1 + eb2 + 2 * eb * cos a
        rd = recip d
        pa = M_PI - abs a
        rd' = 2 / (b*b + pa*pa)
        (l, arg) =
          if d >= M_EPS
          then ((1 - eb2) * rd / mv, 2 * sina * eb * rd)
          else (2 * rd' , negate sina * rd')
    {-# INLINE sinh #-}
    sinh (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0  = packQ x y z (sinh w)
      | otherwise = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv  = sqrt mv2
        l   = cosh w * sin mv / mv
        arg = sinh w * cos mv
    {-# INLINE cosh #-}
    cosh (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0  = packQ x y z (cosh w)
      | otherwise = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv  = sqrt mv2
        l   = sinh w * sin mv / mv
        arg = cosh w * cos mv
    {-# INLINE tanh #-}
    tanh (unpackQ# -> (# x, y, z, w #))
      | mv2 == 0       = packQ x y z (tanh w)
      | isInfinite mv2 = packQ 0 0 0 (signum w)
      | otherwise      = packQ (x * l) (y * l) (z * l) arg
      where
        mv2 = (x * x) + (y * y) + (z * z)
        mv = sqrt mv2
        b = 2*w
        a = 2*mv
        eb = exp (- abs b)
        eb2 = eb*eb
        d = 1 + eb2 + 2 * eb * cos a
        rd = recip d
        pa = M_PI - a
        rd' = 2 / (b*b + pa*pa)
        (l, arg) =
          if d >= M_EPS
          then (2 * sin a * eb * rd / mv, copysign (1 - eb2) b * rd)
          else (2 * rd' , b * rd')
    {-# INLINE asin #-}
    -- The original formula:
    -- asin q = -i * log (i*q + sqrt (1 - q*q))
    -- below is a more numerically stable version.
    asin (unpackQ# -> (# x, y, z, w #))
      | v2 == 0   = if w2 <= 1
                    then packQ x y z (asin w)
                    else packQ l 0 0 arg
      | otherwise  = packQ (x*c) (y*c) (z*c) arg
      where
        v2 = (x * x) + (y * y) + (z * z)
        v = sqrt v2
        w2 = w*w
        w1qq = 0.5 *(1 - w2 + v2)       -- real part of (1 - q*q)/2
        l1qq = sqrt (w1qq*w1qq + w2*v2) -- length of (1 - q*q)/2
        sp2 = l1qq + w1qq
        sn2 = l1qq - w1qq
        sp = sqrt sp2
        sn = copysign (sqrt sn2) w
        -- choose a more stable (symbolically equiv) version
        dp = if 2 * v2 <= sp2 then sp - v else v2 / ((sp + v)*(sn2 + v2))
        dn = if 2 * w2 <= sn2 then w - sn else w2 / ((sn + w)*(sp2 + w2))
        (wD, vD) = case compare w1qq 0 of
            GT -> (dp, w * dp / sp)
            LT -> (v * dn / sn, dn)
            EQ -> (-v, w)
        l = -0.5 * log (wD*wD + vD*vD)
        c = l / v
        arg = atan2 vD wD
    {-# INLINE acos #-}
    acos q = M_PI_2 - asin q
    {-# INLINE atan #-}
    -- atan q = i / 2 * log ( (i + q) / (i - q) )
    atan (unpackQ# -> (# x, y, z, w #))
      | v2 == 0   = packQ x y z (atan w)
      | otherwise = packQ (x*c) (y*c) (z*c) arg
      where
        v2 = (x * x) + (y * y) + (z * z)
        v = sqrt v2
        w2 = w*w
        q2 = w2 + v2
        v' = v - 1
        mzero = w2 + v'*v'
        (c, arg) =
          if mzero == 0
          then ( sqrt maxFinite / v, 0)
          else ( 0.25 * (log (1 + q2 + 2*v) - log mzero) / v
               , 0.5 * atan2 (2*w) (1 - q2) )
    {-# INLINE asinh #-}
    -- The original formula:
    -- asinh q = log (q + sqrt (q*q + 1))
    -- below is a more numerically stable version.
    asinh (unpackQ# -> (# x, y, z, w #))
      | v2 == 0   = packQ x y z (asinh w)
      | otherwise = packQ (x*c) (y*c) (z*c) arg
      where
        v2 = (x * x) + (y * y) + (z * z)
        v = sqrt v2
        w2 = w*w
        w1qq = 0.5 *(1 + w2 - v2)       -- real part of (1 + q*q)/2
        l1qq = sqrt (w1qq*w1qq + w2*v2) -- length of (1 + q*q)/2
        sp2 = l1qq + w1qq
        sn2 = l1qq - w1qq
        sp = sqrt sp2
        sn = copysign (sqrt sn2) w
        -- choose a more stable (symbolically equiv) version
        dp = if 2 * w >= - sp then w + sp else w2 / ((sp - w)*(w2 + sn2))
        dn = if 2 * v <= - sn || sn >= 0
                              then v + sn else v2 / ((v - sn)*(v2 + sp2))
        (wD, vD) = case compare w1qq 0 of
            GT -> (dp, v * dp / sp)
            LT -> (w * dn / sn, dn)
            EQ -> (w, v)
        c = atan2 vD wD / v
        arg = 0.5 * log (wD*wD + vD*vD)
    {-# INLINE acosh #-}
    -- The original formula:
    -- asinh q = log (q + sqrt (q + 1) * sqrt (q - 1))
    -- below is a more numerically stable version.
    -- note, log (q + sqrt (q*q - 1)) would not work, because that would not
    -- be the principal value.
    acosh (unpackQ# -> (# x, y, z, w #))
      | v2 == 0   = packQ x y z (acosh w)
      | otherwise = packQ (x*c) (y*c) (z*c) arg
      where
        v2 = (x * x) + (y * y) + (z * z)
        v = sqrt v2
        w2 = w*w
        w1qq = 0.5 *(w2 - v2 - 1)       -- real part of (q*q - 1)/2
        l1qq = sqrt (w1qq*w1qq + w2*v2) -- length of (q*q - 1)/2
        sp2 = l1qq + w1qq
        sn2 = l1qq - w1qq
        sp = sqrt sp2
        sn = copysign (sqrt sn2) w
        -- choose a more stable (symbolically equiv) version
        dp = if 2 * w >= - sp then w + sp else w2 / ((w - sp)*(w2 + sn2))
        dn = if 2 * v <= - sn || sn >= 0
                              then v + sn else v2 / ((sn - v)*(v2 + sp2))
        (wD, vD) = case compare w1qq 0 of
            GT -> (dp, v * dp / sp)
            LT -> (w * dn / sn, dn)
            EQ -> (w, v)
        c = atan2 vD wD / v
        arg = 0.5 * log (wD*wD + vD*vD)
    {-# INLINE atanh #-}
    -- atanh q = 0.5 * log ( (1 + q) / (1 - q) )
    atanh (unpackQ# -> (# x, y, z, w #))
      | v2 ==  0  = packQ x y z (atanh w)
      | otherwise = packQ (x*c) (y*c) (z*c) (copysign arg w)
      where
        v2 = (x * x) + (y * y) + (z * z)
        v = sqrt v2
        w2 = w*w
        q2 = w2 + v2
        w' = abs w - 1
        c  = 0.5 * atan2 (2*v) (1 - q2) / v
        arg = if w' == 0
              then (1/0)
              else 0.25 * (log (1 + q2 + 2 * abs w) - log (v2 + w'*w'))

-- If q is negative real, provide a fallback axis to align log.
log' :: Vector Double 3 -> QDouble -> QDouble
log' r (unpackQ# -> (# x, y, z, w #))
  = case (x * x) + (y * y) + (z * z) of
    0.0 | w >= 0
           -> packQ 0 0 0 (log w)
        | Vec3 rx ry rz <- r
           ->  packQ (M_PI*rx) (M_PI*ry) (M_PI*rz) (log (negate w))
    mv2 -> case (# mv2 + w * w, sqrt mv2 #) of
      (# q2, mv #) -> case atan2 mv w / mv of
        l -> packQ (x * l) (y * l) (z * l) (0.5 * log q2)

-- If q is negative real, provide a fallback axis to align sqrt.
sqrt' :: Vector Double 3 -> QDouble -> QDouble
sqrt' r (unpackQ# -> (# x, y, z, w #))
  | v2 == 0 && w >= 0
    = packQ x y z (sqrt w)
  | v2 == 0
  , Vec3 rx ry rz <- r
  , sw <- sqrt (negate w)
    = packQ (sw*rx) (sw*ry) (sw*rz) 0
  | otherwise
    = packQ (x * c) (y * c) (z * c) arg
  where
    v2 = (x * x) + (y * y) + (z * z)
    mq = sqrt (v2 + w * w)
    arg = sqrt $ 0.5 * if w >= 0 then mq + w else v2 / (mq - w)
    c = 0.5 / arg

instance Eq QDouble where
    {-# INLINE (==) #-}
    QDouble a == QDouble b = a == b
    {-# INLINE (/=) #-}
    QDouble a /= QDouble b = a /= b

max_stars_repo_name: achirkin/fasttensor

max_issues_repo_name: achirkin/fasttensor

max_forks_repo_name: achirkin/fasttensor


text:
module Statistics.Covariance.Shrinkage.Strimmer (
    psmallSVD
  , nsmallSVD
  , posSVD
  , covShrink
  , corShrink
  , corCoef
  , varShrink
  , varCoef
  ) where

import Numeric.LinearAlgebra hiding ((<>))
import Numeric.Sum (kbn, sumVector)
import Statistics.Quantile (quantile, s)
import qualified Data.Vector.Storable as VS

covShrink :: Matrix Double -> Herm Double
covShrink x = trustSym $ (asColumn sc * unSym c) * asRow sc
  where
    sc = sqrt $ varShrink x
    c = corShrink x

psmallSVD :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
psmallSVD x = (u', d', v')
  where
    b = unSym $ mTm x
    (_, d, v) = svd b
    tol = 1e-6
    positive = find (> tol) d
    u' = (x <> v') * asRow (recip d')
    d' = sqrt . vector $ map (atIndex d) positive
    v' = v ¿ positive

nsmallSVD :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
nsmallSVD x = (u', d', v')
  where
    b = unSym $ mTm (tr x)
    (u, d, _) = svd b
    tol = 1e-6
    positive = find (> tol) d
    u' = u ¿ positive
    d' = sqrt . vector $ map (atIndex d) positive
    v' = (tr x <> u') * asRow (recip d')

posSVD :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
posSVD x = (u', d', v')
  where
    (u, d, v) = svd x
    tol = 1e-6
    positive = find (> tol) d
    u' = u ¿ positive
    d' = vector $ map (atIndex d) positive
    v' = v ¿ positive

fastSVD :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
fastSVD x
  | n > edgeRatio * p = psmallSVD x
  | edgeRatio * n < p = nsmallSVD x
  | otherwise = posSVD x
  where
    edgeRatio = 2
    n = rows x
    p = cols x

mpower :: Herm Double -> Double -> Herm Double
mpower x p = trustSym $ (evec * asRow (cmap exp $ scale p (cmap log eval))) <> tr evec
  where (eval, evec) = eigSH x

corShrink :: Matrix Double -> Herm Double
corShrink x
  | lambda == 1 || alpha == 0 = trustSym $ ident p
  | alpha == 1 = unitDiag $ scale (1 - lambda) r0
  | lambda == 0 = unitDiag . trustSym $ svdxsV <> unSym (mpower c2 alpha) <> tr svdxsV
  | otherwise = unitDiag . trustSym . scale (exp (alpha * log lambda)) $ ident p - (svdxsV <> f <> tr svdxsV)
  where
    unitDiag :: Herm Double -> Herm Double
    unitDiag powr = trustSym (diag (1 - takeDiag (unSym powr))) `add` powr
    n = rows x
    p = cols x
    alpha = 1 :: Double
    w2 = recip $ fromIntegral n -- w = 1 / n => sum(w^2) =  1 / n
    h1 = 1 / (1 - w2)
    xs = let (u, v) = meanCov x in (x - asRow u) / asRow (sqrt (takeDiag (unSym v)))
    lambda = corCoef x
    r0 = scale h1 $ mTm (xs / sqrt (fromIntegral n))
    (svdxsU, svdxsS, svdxsV) = fastSVD xs
    m = size svdxsS
    utwu = tr svdxsU <> svdxsU / fromIntegral n
    c0 = utwu * asRow svdxsS * asColumn svdxsS
    c1 = scale ((1 - lambda) * h1) c0
    c2 = sym c1
    f = ident m - unSym (mpower (trustSym $ unSym (scale (recip lambda) c2) + ident m) alpha)

corCoef :: Matrix Double -> Double
corCoef x
  | denominator == 0 = 1
  | otherwise = numerator / denominator * h1w2
  where
    n = rows x
    p = cols x
    xs = let (u, v) = meanCov x in (x - asRow u) / asRow (sqrt (takeDiag (unSym v)))
    w2 = recip $ fromIntegral n
    h1w2 = w2 / (1 - w2)
    sw = sqrt . recip $ fromIntegral n
    xsw = scale sw xs
    (xswsvdU, xswsvdS, xswsvdV) = fastSVD xsw
    colSums = vector . map (sumVector kbn) . toColumns
    -- colCumSums = fromColumns . map (VS.postscanl (+) 0) . toColumns
    rowCumSums = fromRows . map (VS.postscanl (+) 0) . toRows
    sE2R = sumVector kbn (flatten (xsw * ((xswsvdU * asRow (cb xswsvdS)) <> tr xswsvdV)))
      - sumVector kbn (sq (colSums (sq xsw)))
    xs2w = scale sw $ sq xs
    sER2 = 2 * (sumVector kbn . flatten
      $ fliprl (takeColumns (p - 1) xs2w) * (rowCumSums . fliprl $ dropColumns 1 xs2w))
    denominator = sE2R
    numerator = sER2 - sE2R

varShrink :: Matrix Double -> Vector Double
varShrink x = scalar (lambda * target) + scale (1 - lambda) v
  where
    lambda = varCoef x
    v = takeDiag . unSym . snd $ meanCov x
    target = quantile s 1 2 v

varCoef :: Matrix Double -> Double
varCoef x
  | denominator == 0 = 1
  | otherwise = max 0 . min 1 $ numerator / denominator * h1w2
  where
    n = rows x
    w2 = recip (fromIntegral n)
    h1 = 1 / (1 - w2)
    h1w2 = w2 / (1 - w2)
    xc = x - asRow (fst (meanCov x))
    v = scale h1 q1
    target = quantile s 1 2 v
    zz = xc * xc
    q1 = fst $ meanCov zz
    q2 = fst (meanCov (zz * zz) ) - q1 * q1
    numerator = sumVector kbn q2
    denominator = sumVector kbn . sq $ q1 - scalar (target / h1)

sq, cb :: Num a => a -> a
sq a = a * a
cb a = a * sq a

ext: hs

max_stars_repo_name: tsbattman/coviest

max_issues_repo_name: tsbattman/coviest

max_forks_repo_name: tsbattman/coviest


text:
-----------------------------------------------------------------------------
-- |
-- Module      : Berp.Base.SemanticTypes
-- Copyright   : (c) 2010 Bernie Pope
-- License     : BSD-style
-- Maintainer  : florbitous@gmail.com
-- Stability   : experimental
-- Portability : ghc
--
-- The core types used to represent the state of Python programs. We put
-- them all here in one file because they tend to be mutually recursive.
-- Using one file for such types tends to avoid problems with unbreakable
-- cycles in the Haskell module imports. Try not to put functions in here
-- (except perhaps type class instances).
--
-----------------------------------------------------------------------------

module Berp.Base.SemanticTypes
   ( Procedure, ControlStack (..), EvalState (..), Object (..), Eval, ObjectRef
   , HashTable, HashSet, ListArray, Arity, ModuleCache
   , initState )  where

import Data.Typeable (Typeable)
import Control.Monad.State.Strict (StateT)
import Control.Monad.Cont (ContT)
import Data.IntMap (IntMap)
import Data.Map as Map (Map, empty)
import Data.IORef (IORef)
import Data.Complex (Complex)
import Data.Array.IO (IOArray)
import Berp.Base.Identity (Identity)

data ControlStack
   = EmptyStack
   | ProcedureCall
     { procedure_return :: Object -> Eval Object
     , control_stack_tail :: ControlStack
     }
   | ExceptionHandler
     { exception_handler :: Maybe (Object -> Eval Object)
     , exception_finally :: Maybe (Eval Object)
     , control_stack_tail :: ControlStack
     }
   | WhileLoop
     { loop_start :: Eval Object
     , loop_end :: Eval Object
     , control_stack_tail :: ControlStack
     }
   | GeneratorCall
     { generator_yield :: Object -> Eval Object
     , generator_object :: Object
     , control_stack_tail :: ControlStack
     }

instance Show ControlStack where
   show EmptyStack = "EmptyStack"
   show (ProcedureCall {}) = "ProcedureCall"
   show (ExceptionHandler {}) = "ExceptionHandler"
   show (WhileLoop {}) = "WhileLoop"
   show (GeneratorCall {}) = "GeneratorCall"

data EvalState =
   EvalState
   { state_control_stack :: !ControlStack
   , state_builtins :: !HashTable
   , state_moduleCache :: !ModuleCache
   }

initState :: HashTable -> EvalState
initState builtins =
   EvalState
   { state_control_stack = EmptyStack
   , state_builtins = builtins
   , state_moduleCache = Map.empty
   }

type ModuleCache = Map String Object

type Eval a = StateT EvalState (ContT Object IO) a

type ObjectRef = IORef Object
type Procedure = [Object] -> Eval Object

-- XXX maybe this should be:
-- IORef (IntMap (IORef [(Object, Object)]))
-- or even:
-- IORef (IntMap (IORef [(Object, ObjectRef)]))
-- type HashTable = IORef (IntMap [(Object, Object)])

type HashTable = IORef (IntMap [(Object, ObjectRef)])
type HashSet = IORef (IntMap [Object])

type ListArray = IOArray Integer Object
type Arity = Int

{-
-- Here's another possible encoding of objects which is more abstract.
-- It would make it easier to add new object types, and they could be added
-- in separate modules. The problem is that it would make it slower for
-- detecting the kind of object we have, compared to the Alegabraic approach
-- which gives us a tag to match against.

class ObjectLike t where
   identity :: t -> Identity
   ...

data Object = forall t . (ObjectLike t, Typeable t) => O t

instance ObjectLike Object where
   identity (O x) = identity x
   ...
-}

-- XXX probably need Bound Methods.
data Object
   = Object
     { object_identity :: !Identity
     , object_type :: !Object -- type
     , object_dict :: !Object -- dictionary
     }
   | Type
     { object_identity :: !Identity
     , object_type :: Object  -- type, is this needed? Is the type of all types == type?
     , object_dict :: !Object  -- dictionary
     , object_bases :: !Object -- tuple
     , object_constructor :: !Procedure
     , object_type_name :: !Object -- string
     , object_mro :: !Object -- tuple. Method Resolution Order.
     }
   | Integer
     { object_identity :: !Identity
     , object_integer :: !Integer
     }
   | Float
     { object_identity :: !Identity
     , object_float :: !Double
     }
   | Complex
     { object_identity :: !Identity
     , object_complex :: Complex Double
     }
   | Bool
     { object_identity :: !Identity
     , object_bool :: !Bool
     }
   | Tuple
     { object_identity :: !Identity
     , object_tuple :: ![Object]
     , object_length :: !Int
     }
   -- We must allow the contents and the size of an array to change.
   | List
     { object_identity :: !Identity
     , object_list_elements :: !(IORef ListArray)
     , object_list_num_elements :: !(IORef Integer)
     }
   | Function
     { object_identity :: !Identity
     , object_procedure :: !Procedure
     , object_arity :: !Arity
     , object_dict :: !Object -- dictionary
     }
   | String
     { object_identity :: !Identity
     , object_string :: !String
     }
   | Dictionary
     { object_identity :: !Identity
     , object_hashTable :: !HashTable
     }
   | Set
     { object_identity :: !Identity
     , object_hashSet :: !HashSet
     }
   | Generator
     { object_identity :: !Identity
     , object_continuation :: !(IORef (Eval Object))
     , object_stack_context :: !(IORef (ControlStack -> ControlStack))
     }
     -- Modules probably need names and source information.
   | Module
     { object_identity :: !Identity
     , object_dict :: !Object -- dictionary
     }
   | None

-- For debugging only
instance Show Object where
   show obj@(Object {}) = "object of (" ++ show (object_type obj) ++ ")"
   show obj@(Type {}) = "type(" ++ show (object_type_name obj) ++ ")"
   show obj@(Integer {}) = "integer(" ++ show (object_integer obj) ++ ")"
   show obj@(Float {}) = "float(" ++ show (object_float obj) ++ ")"
   show obj@(Bool {}) = "bool(" ++ show (object_bool obj) ++ ")"
   show (Tuple {}) = "tuple"
   show (List {}) = "list"
   show (Function {}) = "function"
   show obj@(String {}) = "string(" ++ show (object_string obj) ++ ")"
   show (Dictionary {}) = "dictionary"
   show (Set {}) = "set"
   show (Generator {}) = "generator"
   show obj@(Complex {}) = "complex(" ++ show (object_complex obj) ++ ")"
   show (Module {}) = "module"
   show (None {}) = "None"

-- equality instance for objects
-- NOTE: use with care. This does not call the user defined equality
-- on the object. It only uses identity equality.

instance Eq Object where
   None {} == None {} = True
   obj1 == obj2 = object_identity obj1 == object_identity obj2

max_stars_repo_name: ppelleti/berp

max_issues_repo_name: ppelleti/berp

max_forks_repo_name: bjpop/berp


text:
{-# LANGUAGE TupleSections #-}

{-|

Copyright:
  This file is part of the package zxcvbn-hs. It is subject to the
  license terms in the LICENSE file found in the top-level directory
  of this distribution and at:

    https://code.devalot.com/sthenauth/zxcvbn-hs

  No part of this package, including this file, may be copied,
  modified, propagated, or distributed except according to the terms
  contained in the LICENSE file.

License: MIT

-}
module Text.Password.Strength.Internal.Keyboard (
  -- * Keyboard Pattern Matching
  KeyboardPattern,
  keyboardToken,
  keyboardPattern,
  keyboardEstimate
  ) where

--------------------------------------------------------------------------------
-- Library Imports:
import Control.Lens (Lens, (^.), _3)
import Data.Foldable (foldl')
import Numeric.SpecFunctions (choose)

--------------------------------------------------------------------------------
-- Project Imports:
import Text.Password.Strength.Internal.Adjacency
import Text.Password.Strength.Internal.Token
import Text.Password.Strength.Internal.Math (variations')

--------------------------------------------------------------------------------
-- | Information about a found pattern.
newtype KeyboardPattern =
  KeyboardPattern (Int, Int, Token, AdjacencyScore)
  deriving (Show)

--------------------------------------------------------------------------------
-- | Allow other code to access the token used in a pattern.
keyboardToken :: Lens KeyboardPattern KeyboardPattern Token Token
keyboardToken f (KeyboardPattern t) = KeyboardPattern <$> _3 f t

--------------------------------------------------------------------------------
-- | Helper function to check if a token forms a keyboard pattern.
keyboardPattern :: AdjacencyTable -> Token -> Maybe KeyboardPattern
keyboardPattern graph token = KeyboardPattern .
  (graph ^. totalChars, graph ^. averageNeighbors, token,) .
    foldl' scoreSequence mempty <$> findSequence (token ^. tokenChars) graph

--------------------------------------------------------------------------------
-- | Estimate the number of guesses needed for a keyboard pattern to
-- be cracked.
keyboardEstimate :: KeyboardPattern -> Integer
keyboardEstimate (KeyboardPattern (s, d, _, a)) =
    e3 * e2 (a ^. primaryLayer) (a ^. secondaryLayer)

  where
    -- Equation 3, section 4, page 163 (8/18)
    --
    -- Deviations from the paper or other implementations:
    --
    --   * There's a typo in the paper: min(T, i - i)
    --     but should be: min(T, i - 1)
    --
    --   * Another typo, i should start at 2 and not 1.
    --
    --   * The other implementations don't seem to divide the outer
    --     sum by two but the equation clearly does.
    e3 :: Integer
    e3 = max 1 . (`div` 2) . sum $ do
      i <- [2 .. (a ^. patternLength)]
      j <- [1 .. min (a ^. totalTurns) (i - 1)]
      pure $ floor (choose (i - 1) (j - 1)) * toInteger s * (toInteger d ^ j)

    -- Modified version of equation 2 for primary layer vs. secondary.
    e2 :: Int -> Int -> Integer
    e2 = variations'


text:
{-# LANGUAGE DisambiguateRecordFields #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE NamedFieldPuns #-}

module SimulatedAnnealing (
   State,
   energy,
   metric,
   perturb,
   Args(..),
   CoolingTemp(..),
   anneal
) where

import Metric (MetricSpace, Metric, metric)

import qualified Metric

import Prelude hiding (min)

import Control.Monad.Random (evalRand)
import Control.Monad.Random.Class (MonadRandom, getRandom)
import Data.Function ((&))
import Data.Functor ((<&>))
import Data.Int (Int32)
import Numeric.GSL.SimulatedAnnealing (SimulatedAnnealingParams(..), simanSolve)
import Numeric.LinearAlgebra.Data (Vector, (!))
import System.Random (mkStdGen, StdGen)

import Debug.Trace

class MetricSpace a => State a where
    energy :: a -> Double
    perturb :: MonadRandom m => Double -> a -> m a

data Args = Args {
    numTries :: Word,
    numItersPerTemp :: Word,
    maxStepSize :: Double,
    boltzmannK :: Double,
    coolingTemp :: CoolingTemp
}

data CoolingTemp = CoolingTemp {
    initial :: Double,
    rate :: Double,
    min :: Double
}

argsToParams :: Args -> SimulatedAnnealingParams
argsToParams Args {
    numTries,
    numItersPerTemp,
    maxStepSize,
    boltzmannK,
    coolingTemp = CoolingTemp {
        initial,
        rate,
        min
    }
} = SimulatedAnnealingParams {
    n_tries = fromIntegral numTries,
    iters_fixed_T = fromIntegral numItersPerTemp,
    step_size = maxStepSize,
    boltzmann_k = boltzmannK,
    cooling_t_initial = initial,
    cooling_mu_t = rate,
    cooling_t_min = min
}

data Annealer a = Annealer {
    seed :: Int,
    numRandoms :: Int,
    params :: SimulatedAnnealingParams,
    initialState :: a,
    energy' :: a -> Double,
    metric' :: a -> a -> Double,
    step :: Vector Double -> Double -> a -> a,
    maybeShow :: Maybe (a -> String)
}

rawAnneal :: Annealer a -> a
rawAnneal Annealer {
    seed,
    numRandoms,
    params,
    initialState,
    energy',
    metric',
    step,
    maybeShow
} = simanSolve
    seed
    numRandoms
    params
    initialState
    energy'
    metric'
    step
    maybeShow

anneal :: (State a, MonadRandom m) => a -> Args -> Maybe (a -> String) -> m a
anneal initialState args maybeShow = do
    seed <- getRandom
    Annealer {
        seed,
        numRandoms = 1,
        params = argsToParams args,
        initialState,
        energy' = energy,
        metric' = metric,
        step,
        maybeShow
    }
        & rawAnneal
        & return
      where
        -- It's too hard to thread the monad through the SA,
        -- b/c the FFI library isn't designed with monads.
        -- So I'm doing it manually here, seeding a new StdGen using mkStdGen
        -- with a seed from the MonadRandom.
        step vectors stepSize s = s'
          where
            seed = (vectors ! 0)
                * (fromIntegral (maxBound :: Int32))
                & truncate
            g = mkStdGen seed
            rand = perturb stepSize s
            s' = evalRand rand g

realAnneal :: (State a, MonadRandom m) => a -> Args -> m a
realAnneal initialState args = do
    undefined

max_stars_repo_name: kkysen/ParallelBoggle

max_issues_repo_name: kkysen/ParallelBoggle

max_forks_repo_name: kkysen/ParallelBoggle


text:
-- This code was automatically generated by lhs2tex --code, from the file
-- HSoM/SpectrumAnalysis.lhs.  (See HSoM/MakeCode.bat.)

module Euterpea.Music.Signal.SpectrumAnalysis where

import Euterpea

import Data.Complex
import Data.List
-- import Math.FFT
dft :: RealFloat a => [Complex a] -> [Complex a]
dft xs =
  let  lenI = length xs
       lenR = fromIntegral lenI
       lenC = lenR :+ 0
  in [  let i = -2 * pi * fromIntegral k / lenR
        in (1/lenC) * sum [  (xs!!n) * exp (0 :+ i * fromIntegral n)
                             | n <- [0,1..lenI-1] ]
        | k <- [0,1..lenI-1] ]
mkTerm :: Int -> Double -> [Complex Double]
mkTerm num n = let f = 2 * pi / fromIntegral num
               in map (:+ 0) [  sin (n * f * fromIntegral i) / n
                                | i <- [0,1..num-1] ]

mkxa, mkxb, mkxc :: Int-> [Complex Double]
mkxa num = mkTerm num 1
mkxb num = zipWith (+) (mkxa num) (mkTerm num 3)
mkxc num = zipWith (+) (mkxb num) (mkTerm num 5)
printComplexL :: [Complex Double] -> IO ()
printComplexL xs  =
  let  f (i,rl:+im) =
            do  putStr (spaces (3 - length (show i))  )
                putStr (show i       ++ ":  ("        )
                putStr (niceNum rl  ++ ", "           )
                putStr (niceNum im  ++ ")\n"          )
  in mapM_ f (zip [0..length xs - 1] xs)

niceNum :: Double -> String
niceNum d =
  let  d' = fromIntegral (round (1e10 * d)) / 1e10
       (dec, fra)  = break (== '.') (show d')
       (fra',exp)  = break (== 'e') fra
  in  spaces (3  - length dec) ++ dec ++ take 11 fra'
      ++ exp ++ spaces (12 - length fra' - length exp)

spaces :: Int -> String
spaces  n = take n (repeat ' ')
mkPulse :: Int -> [Complex Double]
mkPulse n = 100 : take (n-1) (repeat 0)
x1 num = let f = pi * 2 * pi / fromIntegral num
         in map (:+ 0) [  sin (f * fromIntegral i)
                          | i <- [0,1..num-1] ]
mkPolars :: [Complex Double] -> [Complex Double]
mkPolars = map ((\(m,p)-> m:+p) . polar)

max_stars_repo_name: kianwilcox/haskell-music

max_issues_repo_name: kianwilcox/haskell-music

max_forks_repo_name: kianwilcox/haskell-music


text:
module HVX.Internal.Solvers
  ( constStep
  , decNonSumStep
  , subgradMinimize
  , subgradMaximize
  , subgradLoop
  , ellipsoidMinimize
  , ellipsoidMaximize
  , ellipsoidLoop
  ) where

import Data.List
import Data.Ord
import Data.Maybe
import Numeric.LinearAlgebra

import HVX.Internal.Constraints
import HVX.Internal.DCP
import HVX.Internal.Matrix
import HVX.Internal.Primitives
import HVX.Internal.SymbolicSubgrad
import HVX.Internal.Util

class CanMinimize vex
instance CanMinimize Affine
instance CanMinimize Convex

class CanMaximize vex
instance CanMaximize Affine
instance CanMaximize Concave

-- | Constant step size.
constStep :: Double -> Int -> Double
constStep alpha _ = alpha

-- | 1/n step size.
decNonSumStep :: Double -> Int -> Double
decNonSumStep alpha itr = alpha / fromIntegral itr

-- | Use subgradient method to solve a minimization problem.
subgradMinimize :: CanMinimize vex =>
     Expr vex mon     -- ^ Objective to minimize.
  -> [Constraint]     -- ^ Constraints on minimization problem.
  -> (Int -> Double)  -- ^ Stepsize function.
  -> Int              -- ^ Number of iterations to run.
  -> Vars             -- ^ Variables and their starting values.
  -> (Vars, Double)   -- ^ Optimal variable values and optimal objective value.
subgradMinimize = subgradLoop 1

-- | Use subgradient method to solve a maximization problem.
subgradMaximize :: CanMaximize vex =>
     Expr vex mon     -- ^ Objective to maximize.
  -> [Constraint]     -- ^ Constraints on maximization problem.
  -> (Int -> Double)  -- ^ Stepsize function.
  -> Int              -- ^ Number of iterations to run.
  -> Vars             -- ^ Variables and their starting values.
  -> (Vars, Double)   -- ^ Optimal variable values and optimal objective value.
subgradMaximize objective constraints stepFun maxItr vars = (fst result, negate $ snd result)
  where result = subgradLoop 1 (EFun Neg objective) constraints stepFun maxItr vars

-- | Subgradient solver that minimizes an expression subject to constraints.
subgradLoop :: Int -> Expr vex mon -> [Constraint] -> (Int -> Double) -> Int -> Vars
  -> (Vars, Double)
subgradLoop itr objective constraints stepFun maxItr vars =
  if itr >= maxItr || vars == varsNext
    then (vars, evaluate objective vars @@> (0,0))
    else subgradLoop (itr + 1) objective constraints stepFun maxItr varsNext
      where varsNext = map (updateWithSubgrad objective constraints (stepFun itr) vars) vars

updateWithSubgrad :: Expr vex mon -> [Constraint] -> Double -> Vars -> (Var, Mat) -> (Var, Mat)
updateWithSubgrad objective constraints stepSize vars (varname, val) =
  if any isNaN (toList $ flatten g)
    then (varname, val)
    else (varname, val - scale stepSize g)
      where g = fromMaybe (trans $ jacobianWrtVar objective vars varname) $
                          getConstraintSubgrad constraints vars varname

type Soid = (Var, Mat, Mat)

-- | Use ellipsoid method to solve a minimization problem.
ellipsoidMinimize :: CanMinimize vex =>
     Expr vex mon            -- ^ Objective to minimize.
  -> [Constraint]            -- ^ Constraints on minimization problem.
  -> [(Var,Int)]             -- ^ Variable names and their sizes.
  -> Double                  -- ^ Tolerance.
  -> Double                  -- ^ Radius to seed initial sphere with.
  -> (Vars, Double, Double)  -- ^  Optimal variable values, lower bound, and upper bound.
ellipsoidMinimize objective constraints varsizes tol radius =
  ellipsoidLoop objective constraints tol soids (-infinity) infinity
    where soids = map (\(name,n) -> (name, zeroVec n, scale radius $ ident n)) varsizes

-- | Use ellipsoid method to solve a maximization problem.
ellipsoidMaximize :: CanMaximize vex =>
     Expr vex mon            -- ^ Objective to maximize.
  -> [Constraint]            -- ^ Constraints on maximization problem.
  -> [(Var,Int)]             -- ^ Variable names and their sizes.
  -> Double                  -- ^ Tolerance.
  -> Double                  -- ^ Radius to seed initial sphere with.
  -> (Vars, Double, Double)  -- ^  Optimal variable values, lower bound, and upper bound.
ellipsoidMaximize objective constraints varsizes tol radius = (optvars, negate ubound, negate lbound)
  where (optvars, lbound, ubound) = ellipsoidLoop (EFun Neg objective) constraints tol soids (-infinity) infinity
        soids = map (\(name,n) -> (name, zeroVec n, scale radius $ ident n)) varsizes

-- | Ellipsoid method solver that minimizes an expression subject to constraints.
ellipsoidLoop :: Expr vex mon -> [Constraint] -> Double -> [Soid] -> Double -> Double
  -> (Vars, Double, Double)
ellipsoidLoop objective constraints tol soids lbound ubound =
  if ubound - lbound < tol || sqrtgpgsum == 0.0
    then (centers, nlbound, nubound)
    else ellipsoidLoop objective constraints tol nsoids nlbound nubound
      where updates = map (updateEllipsoid objective constraints centers) soids
            nsoids = map fst updates
            sqrtgpgsum = sqrt . sum . map snd $ updates
            objval = evaluate objective centers @@> (0,0)
            nlbound = max lbound $ objval - sqrtgpgsum
            nubound = min ubound objval
            centers = map (\(name, val, _) -> (name, val)) soids

updateEllipsoid :: Expr vex mon -> [Constraint] -> Vars -> Soid -> (Soid, Double)
updateEllipsoid objective constraints vars soid =
  if gpg == 0.0
    then (soid, objgpg)
    else (nsoid, objgpg)
      where objg = trans $ jacobianWrtVar objective vars varname
            g = fromMaybe objg $ getConstraintSubgrad constraints vars varname
            n = fromIntegral $ rows val
            gpg = (@@> (0,0)) $ trans g <> p <> g
            objgpg = (@@> (0,0)) $ trans objg <> p <> objg
            gtilde = scale (1 / sqrt gpg) g
            nval = val - scale (1 / (n + 1)) (p <> gtilde)
            np = scale (n**2 / (n**2 - 1)) $
              p - scale (2/(n+1)) (p <> (gtilde <> trans gtilde) <> p)
            (varname, val, p) = soid
            nsoid = (varname, nval, np)

-- | Get the subgradient of the maximally violated constraint if there is a
-- violation. Returns nothing if all constraints are satisfied
getConstraintSubgrad :: [Expr vex mon] -> Vars -> String -> Maybe Mat
getConstraintSubgrad constraints vars var =
  if null constraints || maxVal < 0.0
    then Nothing
    else Just cg
      where cg = trans $ jacobianWrtVar maxConstraint vars var
            (maxConstraint, maxVal) = maximumBy (comparing snd) constraintViolations
            constraintViolations = map (\c -> (c, evaluate c vars @@> (0,0))) constraints

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx


text:
-- |
-- Module:      Math.NumberTheory.RootsOfUnity
-- Copyright:   (c) 2018 Bhavik Mehta
-- Licence:     MIT
-- Maintainer:  Bhavik Mehta <bhavikmehta8@gmail.com>
--
-- Implementation of roots of unity
--

module Math.NumberTheory.RootsOfUnity
(
-- * Roots of unity
   RootOfUnity (..)
-- ** Conversions
   , toRootOfUnity
   , toComplex )

where

import Data.Complex                                        (Complex(..), cis)
import Data.Semigroup                                      (Semigroup(..))
import Data.Ratio                                          ((%), numerator, denominator)

-- | A representation of <https://en.wikipedia.org/wiki/Root_of_unity roots of unity>: complex
-- numbers \(z\) for which there is \(n\) such that \(z^n=1\).
newtype RootOfUnity =
  RootOfUnity { -- | Every root of unity can be expressed as \(e^{2 \pi i q}\) for some
                -- rational \(q\) satisfying \(0 \leq q < 1\), this function extracts it.
                fromRootOfUnity :: Rational }
  deriving (Eq)

instance Show RootOfUnity where
  show (RootOfUnity q)
    | n == 0    = "1"
    | d == 1    = "-1"
    | n == 1    = "e^(πi/" ++ show d ++ ")"
    | otherwise = "e^(" ++ show n ++ "πi/" ++ show d ++ ")"
    where n = numerator (2*q)
          d = denominator (2*q)

-- | Given a rational \(q\), produce the root of unity \(e^{2 \pi i q}\).
toRootOfUnity :: Rational -> RootOfUnity
toRootOfUnity q = RootOfUnity ((n `rem` d) % d)
  where n = numerator q
        d = denominator q
        -- effectively q `mod` 1
  -- This smart constructor ensures that the rational is always in the range 0 <= q < 1.

-- | This Semigroup is in fact a group, so @'stimes'@ can be called with a negative first argument.
instance Semigroup RootOfUnity where
  RootOfUnity q1 <> RootOfUnity q2 = toRootOfUnity (q1 + q2)
  stimes k (RootOfUnity q) = toRootOfUnity (q * fromIntegral k)

instance Monoid RootOfUnity where
  mappend = (<>)
  mempty = RootOfUnity 0

-- | Convert a root of unity into an inexact complex number. Due to floating point inaccuracies,
-- it is recommended to avoid use of this until the end of a calculation. Alternatively, with
-- the [cyclotomic](http://hackage.haskell.org/package/cyclotomic-0.5.1) package, one can use
-- @[polarRat](https://hackage.haskell.org/package/cyclotomic-0.5.1/docs/Data-Complex-Cyclotomic.html#v:polarRat)
-- 1 . @'fromRootOfUnity' to convert to a cyclotomic number.
toComplex :: Floating a => RootOfUnity -> Complex a
toComplex (RootOfUnity t)
  | t == 1/2 = (-1) :+ 0
  | t == 1/4 = 0 :+ 1
  | t == 3/4 = 0 :+ (-1)
  | otherwise = cis . (2*pi*) . fromRational $ t

max_stars_repo_name: cjlarose/arithmoi

max_issues_repo_name: cjlarose/arithmoi

max_forks_repo_name: cjlarose/arithmoi


text:
module Vulkyrie.Vulkan.Default.RenderPass
  ( createPrivateAttachments
  , createRenderPass
  , createRenderPassBeginInfo
  ) where

import           Data.Bits
import           Graphics.Vulkan
import           Graphics.Vulkan.Core_1_0
import           Graphics.Vulkan.Marshal.Create
import           Graphics.Vulkan.Marshal.Create.DataFrame
import           Numeric.Vector

import           Vulkyrie.Program
import           Vulkyrie.Program.Foreign
import           Vulkyrie.Resource
import           Vulkyrie.Vulkan.Engine
import           Vulkyrie.Vulkan.Image

  -- The colorOut is manually passed as last attachment to the renderpass and has 1 sample.
  -- If samples = 1bit then colorOut is used as the color attachment, otherwise
  -- colorOut is used as the resolve attachment.
  --
  -- => order of attachments has to be: depth, color, maybe resolve

-- | Create attachments that are only internally needed by the renderpass, in the expected order.
createPrivateAttachments :: EngineCapability
                         -> VkExtent2D
                         -> VkFormat
                         -> VkSampleCountFlagBits
                         -> Resource ([(VkSemaphore, VkPipelineStageBitmask a)], [VkImageView])
createPrivateAttachments cap extent imgFormat samples = Resource $ do
  let msaaOn = samples /= VK_SAMPLE_COUNT_1_BIT
  fmap unzip . sequence $
    [auto $ createDepthAttImgView cap extent samples]
    <> [auto $ createColorAttImgView cap imgFormat extent samples | msaaOn]

-- | Attachment order: Private attachments first, followed by one color output attachment having 1 sample.
createRenderPass :: VkDevice
                 -> VkFormat
                 -> VkFormat
                 -> VkSampleCountFlagBits
                 -> VkImageLayout
                 -> MetaResource VkRenderPass
createRenderPass dev colorFormat depthFormat samples colorOutFinalLayout =
  let msaaOn = samples /= VK_SAMPLE_COUNT_1_BIT
      finalColorLayout =
        if msaaOn then
          VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL
        else
          colorOutFinalLayout
      depthAttachment = createVk @VkAttachmentDescription
        $  set @"flags" VK_ZERO_FLAGS
        &* set @"format" depthFormat
        &* set @"samples" samples
        &* set @"loadOp" VK_ATTACHMENT_LOAD_OP_CLEAR
        &* set @"storeOp" VK_ATTACHMENT_STORE_OP_DONT_CARE
        &* set @"stencilLoadOp" VK_ATTACHMENT_LOAD_OP_DONT_CARE
        &* set @"stencilStoreOp" VK_ATTACHMENT_STORE_OP_DONT_CARE
        &* set @"initialLayout" VK_IMAGE_LAYOUT_UNDEFINED
        &* set @"finalLayout" VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL

      colorAttachment = createVk @VkAttachmentDescription
        $  set @"flags" VK_ZERO_FLAGS
        &* set @"format" colorFormat
        &* set @"samples" samples
        &* set @"loadOp" VK_ATTACHMENT_LOAD_OP_CLEAR
        &* set @"storeOp" VK_ATTACHMENT_STORE_OP_STORE
        &* set @"stencilLoadOp" VK_ATTACHMENT_LOAD_OP_DONT_CARE
        &* set @"stencilStoreOp" VK_ATTACHMENT_STORE_OP_DONT_CARE
        &* set @"initialLayout" VK_IMAGE_LAYOUT_UNDEFINED
        &* set @"finalLayout" finalColorLayout

      -- this is needed when msaaOn
      resolveAttachment = createVk @VkAttachmentDescription
        $  set @"flags" VK_ZERO_FLAGS
        &* set @"format" colorFormat
        &* set @"samples" VK_SAMPLE_COUNT_1_BIT
        &* set @"loadOp" VK_ATTACHMENT_LOAD_OP_DONT_CARE
        &* set @"storeOp" VK_ATTACHMENT_STORE_OP_STORE
        &* set @"stencilLoadOp" VK_ATTACHMENT_LOAD_OP_DONT_CARE
        &* set @"stencilStoreOp" VK_ATTACHMENT_STORE_OP_DONT_CARE
        &* set @"initialLayout" VK_IMAGE_LAYOUT_UNDEFINED
        -- resolve attachment used => it is the color output
        &* set @"finalLayout" colorOutFinalLayout

      -- subpasses and attachment references
      depthAttachmentRef = createVk @VkAttachmentReference
        $  set @"attachment" 0
        &* set @"layout" VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL

      colorAttachmentRef = createVk @VkAttachmentReference
        $  set @"attachment" 1
        &* set @"layout" VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL

      resolveAttachmentRef = createVk @VkAttachmentReference
        $  set @"attachment" 2
        &* set @"layout" VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL

      subpass = createVk @VkSubpassDescription
        $  set @"pipelineBindPoint" VK_PIPELINE_BIND_POINT_GRAPHICS
        &* set @"colorAttachmentCount" 1
        &* setVkRef @"pColorAttachments" colorAttachmentRef
        &* setVkRef @"pDepthStencilAttachment" depthAttachmentRef
        &* ( if msaaOn then
               setVkRef @"pResolveAttachments" resolveAttachmentRef
             else
               set @"pResolveAttachments" VK_NULL
           )
        &* set @"pPreserveAttachments" VK_NULL
        &* set @"pInputAttachments" VK_NULL

      -- subpass dependencies
      dependency = createVk @VkSubpassDependency
        $  set @"srcSubpass" VK_SUBPASS_EXTERNAL
        &* set @"dstSubpass" 0
        &* set @"srcStageMask" VK_PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT
        &* set @"srcAccessMask" VK_ZERO_FLAGS
        &* set @"dstStageMask" VK_PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT
        &* set @"dstAccessMask"
            (   VK_ACCESS_COLOR_ATTACHMENT_READ_BIT
            .|. VK_ACCESS_COLOR_ATTACHMENT_WRITE_BIT )

      -- render pass
      rpCreateInfo = createVk @VkRenderPassCreateInfo
        $  set @"sType" VK_STRUCTURE_TYPE_RENDER_PASS_CREATE_INFO
        &* set @"pNext" VK_NULL
        &* setListCountAndRef @"attachmentCount" @"pAttachments"
            ([depthAttachment, colorAttachment] <> [resolveAttachment | msaaOn])
        &* set @"subpassCount" 1
        &* setVkRef @"pSubpasses" subpass
        &* set @"dependencyCount" 1
        &* setVkRef @"pDependencies" dependency

  in metaResource
       (\rp -> liftIO $ vkDestroyRenderPass dev rp VK_NULL) $
       withVkPtr rpCreateInfo $ \rpciPtr -> allocaPeek $
         runVk . vkCreateRenderPass dev rpciPtr VK_NULL

createRenderPassBeginInfo :: VkRenderPass -> VkFramebuffer -> VkExtent2D -> VkRenderPassBeginInfo
createRenderPassBeginInfo renderPass framebuffer extent =
  createVk @VkRenderPassBeginInfo
      $  set @"sType" VK_STRUCTURE_TYPE_RENDER_PASS_BEGIN_INFO
      &* set @"pNext" VK_NULL
      &* set @"renderPass" renderPass
      &* set @"framebuffer" framebuffer
      &* setVk @"renderArea"
          (  setVk @"offset"
              ( set @"x" 0 &* set @"y" 0 )
          &* set @"extent" extent
          )
      &* setListCountAndRef @"clearValueCount" @"pClearValues"
          -- This needs to fit the renderpass attachments. Clear values for
          -- attachments that don't use VK_ATTACHMENT_LOAD_OP_CLEAR in loadOp or
          -- stencilLoadOp are ignored.
          [ createVk @VkClearValue
             $ setVk @"depthStencil"
             -- this needs to fit pipeline settings regarding depth, including depthCompareOp
             $  set @"depth" 1.0
             &* set @"stencil" 0
          , createVk @VkClearValue
             $ setVk @"color"
             $ setVec @"float32" (vec4 0 0 0.2 1)
          ]

max_stars_repo_name: cjay/Vulkyrie

max_issues_repo_name: cjay/Vulkyrie

max_forks_repo_name: cjay/Vulkyrie


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE MagicHash #-}

#if __GLASGOW_HASKELL__ >= 701
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE Trustworthy #-}
#endif

#if __GLASGOW_HASKELL__ >= 705
{-# LANGUAGE PolyKinds #-}
#endif

#include "HsBaseConfig.h"

module Generics.Deriving.Eq (
  -- * Generic Eq class
    GEq(..)

  -- * Default definition
  , geqdefault

  -- * Internal Eq class
  , GEq'(..)

  ) where

import           Control.Applicative (Const, ZipList)

import           Data.Char (GeneralCategory)
import           Data.Int
import qualified Data.Monoid as Monoid (First, Last)
import           Data.Monoid (All, Any, Dual, Product, Sum)
import           Data.Version (Version)
import           Data.Word

import           Foreign.C.Error
import           Foreign.C.Types
import           Foreign.ForeignPtr (ForeignPtr)
import           Foreign.Ptr
import           Foreign.StablePtr (StablePtr)

import           Generics.Deriving.Base

import           GHC.Exts hiding (Any)

import           System.Exit (ExitCode)
import           System.IO (BufferMode, Handle, HandlePosn, IOMode, SeekMode)
import           System.IO.Error (IOErrorType)
import           System.Posix.Types

#if MIN_VERSION_base(4,4,0)
import           Data.Complex (Complex)
#endif

#if MIN_VERSION_base(4,7,0)
import           Data.Proxy (Proxy)
#endif

#if MIN_VERSION_base(4,8,0)
import           Data.Functor.Identity (Identity)
import           Data.Monoid (Alt)
import           Data.Void (Void)
import           Numeric.Natural (Natural)
#endif

#if MIN_VERSION_base(4,9,0)
import           Data.List.NonEmpty (NonEmpty)
import qualified Data.Semigroup as Semigroup (First, Last)
import           Data.Semigroup (Arg(..), Max, Min, WrappedMonoid)
#endif

--------------------------------------------------------------------------------
-- Generic show
--------------------------------------------------------------------------------

class GEq' f where
  geq' :: f a -> f a -> Bool

instance GEq' V1 where
  geq' _ _ = True

instance GEq' U1 where
  geq' _ _ = True

instance (GEq c) => GEq' (K1 i c) where
  geq' (K1 a) (K1 b) = geq a b

-- No instances for P or Rec because geq is only applicable to types of kind *

instance (GEq' a) => GEq' (M1 i c a) where
  geq' (M1 a) (M1 b) = geq' a b

instance (GEq' a, GEq' b) => GEq' (a :+: b) where
  geq' (L1 a) (L1 b) = geq' a b
  geq' (R1 a) (R1 b) = geq' a b
  geq' _      _      = False

instance (GEq' a, GEq' b) => GEq' (a :*: b) where
  geq' (a1 :*: b1) (a2 :*: b2) = geq' a1 a2 && geq' b1 b2

-- Unboxed types
instance GEq' UAddr where
  geq' (UAddr a1) (UAddr a2)     = isTrue# (eqAddr# a1 a2)
instance GEq' UChar where
  geq' (UChar c1) (UChar c2)     = isTrue# (eqChar# c1 c2)
instance GEq' UDouble where
  geq' (UDouble d1) (UDouble d2) = isTrue# (d1 ==## d2)
instance GEq' UFloat where
  geq' (UFloat f1) (UFloat f2)   = isTrue# (eqFloat# f1 f2)
instance GEq' UInt where
  geq' (UInt i1) (UInt i2)       = isTrue# (i1 ==# i2)
instance GEq' UWord where
  geq' (UWord w1) (UWord w2)     = isTrue# (eqWord# w1 w2)

#if !(MIN_VERSION_base(4,7,0))
isTrue# :: Bool -> Bool
isTrue# = id
#endif

class GEq a where
  geq :: a -> a -> Bool

#if __GLASGOW_HASKELL__ >= 701
  default geq :: (Generic a, GEq' (Rep a)) => a -> a -> Bool
  geq = geqdefault
#endif

geqdefault :: (Generic a, GEq' (Rep a)) => a -> a -> Bool
geqdefault x y = geq' (from x) (from y)

-- Base types instances
instance GEq () where
  geq = geqdefault

instance (GEq a, GEq b) => GEq (a, b) where
  geq = geqdefault

instance (GEq a, GEq b, GEq c) => GEq (a, b, c) where
  geq = geqdefault

instance (GEq a, GEq b, GEq c, GEq d) => GEq (a, b, c, d) where
  geq = geqdefault

instance (GEq a, GEq b, GEq c, GEq d, GEq e) => GEq (a, b, c, d, e) where
  geq = geqdefault

instance (GEq a, GEq b, GEq c, GEq d, GEq e, GEq f)
    => GEq (a, b, c, d, e, f) where
  geq = geqdefault

instance (GEq a, GEq b, GEq c, GEq d, GEq e, GEq f, GEq g)
    => GEq (a, b, c, d, e, f, g) where
  geq = geqdefault

instance GEq a => GEq [a] where
  geq = geqdefault

instance (GEq (f p), GEq (g p)) => GEq ((f :+: g) p) where
  geq = geqdefault

instance (GEq (f p), GEq (g p)) => GEq ((f :*: g) p) where
  geq = geqdefault

instance GEq (f (g p)) => GEq ((f :.: g) p) where
  geq = geqdefault

instance GEq All where
  geq = geqdefault

#if MIN_VERSION_base(4,8,0)
instance GEq (f a) => GEq (Alt f a) where
  geq = geqdefault
#endif

instance GEq Any where
  geq = geqdefault

#if !(MIN_VERSION_base(4,9,0))
instance GEq Arity where
  geq = geqdefault
#endif

#if MIN_VERSION_base(4,9,0)
instance GEq a => GEq (Arg a b) where
  geq (Arg a _) (Arg b _) = geq a b
#endif

instance GEq Associativity where
  geq = geqdefault

instance GEq Bool where
  geq = geqdefault

instance GEq BufferMode where
  geq = (==)

#if defined(HTYPE_CC_T)
instance GEq CCc where
  geq = (==)
#endif

instance GEq CChar where
  geq = (==)

instance GEq CClock where
  geq = (==)

#if defined(HTYPE_DEV_T)
instance GEq CDev where
  geq = (==)
#endif

instance GEq CDouble where
  geq = (==)

instance GEq CFloat where
  geq = (==)

#if defined(HTYPE_GID_T)
instance GEq CGid where
  geq = (==)
#endif

instance GEq Char where
  geq = (==)

#if defined(HTYPE_INO_T)
instance GEq CIno where
  geq = (==)
#endif

instance GEq CInt where
  geq = (==)

instance GEq CIntMax where
  geq = (==)

instance GEq CIntPtr where
  geq = (==)

instance GEq CLLong where
  geq = (==)

instance GEq CLong where
  geq = (==)

#if defined(HTYPE_MODE_T)
instance GEq CMode where
  geq = (==)
#endif

#if defined(HTYPE_NLINK_T)
instance GEq CNlink where
  geq = (==)
#endif

#if defined(HTYPE_OFF_T)
instance GEq COff where
  geq = (==)
#endif

#if MIN_VERSION_base(4,4,0)
instance GEq a => GEq (Complex a) where
  geq = geqdefault
#endif

instance GEq a => GEq (Const a b) where
  geq = geqdefault

#if defined(HTYPE_PID_T)
instance GEq CPid where
  geq = (==)
#endif

instance GEq CPtrdiff where
  geq = (==)

#if defined(HTYPE_RLIM_T)
instance GEq CRLim where
  geq = (==)
#endif

instance GEq CSChar where
  geq = (==)

#if defined(HTYPE_SPEED_T)
instance GEq CSpeed where
  geq = (==)
#endif

#if MIN_VERSION_base(4,4,0)
instance GEq CSUSeconds where
  geq = (==)
#endif

instance GEq CShort where
  geq = (==)

instance GEq CSigAtomic where
  geq = (==)

instance GEq CSize where
  geq = (==)

#if defined(HTYPE_SSIZE_T)
instance GEq CSsize where
  geq = (==)
#endif

#if defined(HTYPE_TCFLAG_T)
instance GEq CTcflag where
  geq = (==)
#endif

instance GEq CTime where
  geq = (==)

instance GEq CUChar where
  geq = (==)

#if defined(HTYPE_UID_T)
instance GEq CUid where
  geq = (==)
#endif

instance GEq CUInt where
  geq = (==)

instance GEq CUIntMax where
  geq = (==)

instance GEq CUIntPtr where
  geq = (==)

instance GEq CULLong where
  geq = (==)

instance GEq CULong where
  geq = (==)

#if MIN_VERSION_base(4,4,0)
instance GEq CUSeconds where
  geq = (==)
#endif

instance GEq CUShort where
  geq = (==)

instance GEq CWchar where
  geq = (==)

#if MIN_VERSION_base(4,9,0)
instance GEq DecidedStrictness where
  geq = geqdefault
#endif

instance GEq Double where
  geq = (==)

instance GEq a => GEq (Down a) where
  geq = geqdefault

instance GEq a => GEq (Dual a) where
  geq = geqdefault

instance (GEq a, GEq b) => GEq (Either a b) where
  geq = geqdefault

instance GEq Errno where
  geq = (==)

instance GEq ExitCode where
  geq = geqdefault

instance GEq Fd where
  geq = (==)

instance GEq a => GEq (Monoid.First a) where
  geq = geqdefault

#if MIN_VERSION_base(4,9,0)
instance GEq a => GEq (Semigroup.First a) where
  geq = geqdefault
#endif

instance GEq Fixity where
  geq = geqdefault

instance GEq Float where
  geq = (==)

instance GEq (ForeignPtr a) where
  geq = (==)

instance GEq (FunPtr a) where
  geq = (==)

instance GEq GeneralCategory where
  geq = (==)

instance GEq Handle where
  geq = (==)

instance GEq HandlePosn where
  geq = (==)

#if MIN_VERSION_base(4,8,0)
instance GEq a => GEq (Identity a) where
  geq = geqdefault
#endif

instance GEq Int where
  geq = (==)

instance GEq Int8 where
  geq = (==)

instance GEq Int16 where
  geq = (==)

instance GEq Int32 where
  geq = (==)

instance GEq Int64 where
  geq = (==)

instance GEq Integer where
  geq = (==)

instance GEq IntPtr where
  geq = (==)

instance GEq IOError where
  geq = (==)

instance GEq IOErrorType where
  geq = (==)

instance GEq IOMode where
  geq = (==)

instance GEq c => GEq (K1 i c p) where
  geq = geqdefault

instance GEq a => GEq (Monoid.Last a) where
  geq = geqdefault

#if MIN_VERSION_base(4,9,0)
instance GEq a => GEq (Semigroup.Last a) where
  geq = geqdefault
#endif

instance GEq (f p) => GEq (M1 i c f p) where
  geq = geqdefault

instance GEq a => GEq (Maybe a) where
  geq = geqdefault

#if MIN_VERSION_base(4,9,0)
instance GEq a => GEq (Max a) where
  geq = geqdefault

instance GEq a => GEq (Min a) where
  geq = geqdefault
#endif

#if MIN_VERSION_base(4,8,0)
instance GEq Natural where
  geq = (==)
#endif

#if MIN_VERSION_base(4,9,0)
instance GEq a => GEq (NonEmpty a) where
  geq = geqdefault
#endif

instance GEq Ordering where
  geq = geqdefault

instance GEq p => GEq (Par1 p) where
  geq = geqdefault

instance GEq a => GEq (Product a) where
  geq = geqdefault

#if MIN_VERSION_base(4,7,0)
instance GEq
# if MIN_VERSION_base(4,9,0)
             (Proxy s)
# else
             (Proxy (s :: *))
# endif
             where
  geq = geqdefault
#endif

instance GEq (Ptr a) where
  geq = (==)

instance GEq (f p) => GEq (Rec1 f p) where
  geq = geqdefault

instance GEq SeekMode where
  geq = (==)

instance GEq (StablePtr a) where
  geq = (==)

#if MIN_VERSION_base(4,9,0)
instance GEq SourceStrictness where
  geq = geqdefault

instance GEq SourceUnpackedness where
  geq = geqdefault
#endif

instance GEq a => GEq (Sum a) where
  geq = geqdefault

instance GEq (U1 p) where
  geq = geqdefault

instance GEq (UAddr p) where
  geq = geqdefault

instance GEq (UChar p) where
  geq = geqdefault

instance GEq (UDouble p) where
  geq = geqdefault

instance GEq (UFloat p) where
  geq = geqdefault

instance GEq (UInt p) where
  geq = geqdefault

instance GEq (UWord p) where
  geq = geqdefault

instance GEq Version where
  geq = (==)

#if MIN_VERSION_base(4,8,0)
instance GEq Void where
  geq = (==)
#endif

instance GEq Word where
  geq = (==)

instance GEq Word8 where
  geq = (==)

instance GEq Word16 where
  geq = (==)

instance GEq Word32 where
  geq = (==)

instance GEq Word64 where
  geq = (==)

instance GEq WordPtr where
  geq = (==)

#if MIN_VERSION_base(4,9,0)
instance GEq m => GEq (WrappedMonoid m) where
  geq = geqdefault
#endif

instance GEq a => GEq (ZipList a) where
  geq = geqdefault

#if MIN_VERSION_base(4,10,0)
instance GEq CBool where
  geq = (==)

# if defined(HTYPE_BLKSIZE_T)
instance GEq CBlkSize where
  geq = (==)
# endif

# if defined(HTYPE_BLKCNT_T)
instance GEq CBlkCnt where
  geq = (==)
# endif

# if defined(HTYPE_CLOCKID_T)
instance GEq CClockId where
  geq = (==)
# endif

# if defined(HTYPE_FSBLKCNT_T)
instance GEq CFsBlkCnt where
  geq = (==)
# endif

# if defined(HTYPE_FSFILCNT_T)
instance GEq CFsFilCnt where
  geq = (==)
# endif

# if defined(HTYPE_ID_T)
instance GEq CId where
  geq = (==)
# endif

# if defined(HTYPE_KEY_T)
instance GEq CKey where
  geq = (==)
# endif

# if defined(HTYPE_TIMER_T)
instance GEq CTimer where
  geq = (==)
# endif
#endif

max_stars_repo_name: treeowl/generic-deriving

max_issues_repo_name: treeowl/generic-deriving

max_forks_repo_name: treeowl/generic-deriving