text:
import System.Environment (getArgs)
import qualified Graphics.Image as G
import qualified Graphics.Image.Interface as GI
import qualified Graphics.Image.Interface.Vector as GIV
import qualified Data.HashMap as H
import qualified Data.Hashable as Ha
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU
import Statistics.LinearRegression (linearRegressionTLS)
-- use mutable vectors, unboxed preferably
import Data.Either (fromRight, partitionEithers)
import Data.List (groupBy, nub, (\\), sort, sortBy)

edgeSensitivity :: Double
edgeSensitivity = 0.3             -- [0,1]

noiseThreshold :: Int
noiseThreshold = 30               -- nT >= 0

vThreshold = 2
hThreshold = 0.5

data Orientation = Vertical | Horizontal | Suspicious deriving (Eq,Show)

data PixelCluster = PXC { pixelList :: V.Vector (Int,Int),
                          boundingBox :: ((Int,Int),(Int,Int)),
                          vhRatio :: Double,
                          slopeNintersect :: (Double,Double),
                          orientation :: Orientation } deriving (Eq,Show)

toPixelCluster :: [(Int,Int)] -> PixelCluster
toPixelCluster pixList = PXC pixelList boundingBox vhRatio slopeNintersect orientation
  where
    pixelList = V.fromList pixList
    boundingBox = ((\(x,y)-> ((minimum x, minimum y),(maximum x, maximum y))) $ unzip pixList)
    vhRatio = (\((minX,minY),(maxX,maxY))-> (fromIntegral $ (maxY-minY))/(fromIntegral $ (maxX-minX))) $ boundingBox
    slopeNintersect = uncurry linearRegressionTLS (V.unzip $ V.map (\(x,y)->
                                                                      ((fromIntegral x :: Double),(fromIntegral y :: Double))) pixelList)
    orientation = if vhRatio <= hThreshold then Horizontal else (if vhRatio >= vThreshold then Vertical else Suspicious)

main :: IO ()
main = do
  argList <- getArgs
  let imgName = head argList
  origImg <- G.readImageY G.VU imgName
  let origImgM = GI.toManifest origImg  -- Convert to manifest array
      edgeImg = toEdgeImage origImgM
      edgeHM = VU.foldl (\hm (pos,_)-> H.insert pos (1::Int) hm) (H.empty) (GIV.filter (\pixelY-> pixelY == (G.PixelY 1)) edgeImg)
      edgeSetGC = getAllConnected edgeHM
      properSet = filter (\x->(length x)>noiseThreshold) edgeSetGC
      pixelClusterList = map toPixelCluster properSet
      verticalPXC = filter (\pxc-> orientation pxc == Vertical) pixelClusterList
      horizontalPXC = filter (\pxc-> orientation pxc == Horizontal) pixelClusterList
      clusterH = map (\x-> concat $ map (V.toList.pixelList) x) $ getClusters horizontalPXC Horizontal
      clusterV = map (\x-> concat $ map (V.toList.pixelList) x) $ getClusters verticalPXC Vertical
      rectPxList = foldl (\l k-> l ++ (makeRect k)) [] $
                   map (\e-> (\(x,y)-> ((minimum x, minimum y),(maximum x, maximum y))) $
                         unzip e) (clusterV++clusterH++properSet)
      activePixelHM = H.fromList $ zip (concat properSet ++ rectPxList) (repeat 0) :: H.Map (Int,Int) Int
      activePixelImage = G.makeImageR G.VU (G.dims origImg) (\tup-> if (H.member tup activePixelHM)
                                                then (G.PixelY (1::Double))
                                                else (G.PixelY (0::Double)))
  G.writeImage ((takeWhile (/= '.') imgName) ++ "Output.jpg") activePixelImage

distanceThreshold = 60

distanceX pxc1 pxc2 = (\(((minX1,minY1),(maxX1,maxY1)),((minX2,minY2),(maxX2,maxY2)))->
                         minimum [abs (maxX1 - minX2), abs (minX1 - maxX2)]) $
                      (boundingBox pxc1, boundingBox pxc2)

distanceY pxc1 pxc2 = (\(((minX1,minY1),(maxX1,maxY1)),((minX2,minY2),(maxX2,maxY2)))->
                         minimum [abs (maxY1 - minY2), abs (minY1 - maxY2)]) $
                      (boundingBox pxc1, boundingBox pxc2)

getClusters :: [PixelCluster] -> Orientation -> [[PixelCluster]]
getClusters pxcList orientation = getClusters' pxcList orientation []

getClusters' pxcList orientation clusterList = if null pxcList then clusterList
                                               else getClusters' newpxcList orientation (clusterList ++ [newCluster])
  where
    distance pxc1 pxc2 = if orientation == Vertical then distanceX pxc1 pxc2 else distanceY pxc1 pxc2
    newCluster = filter (\pxc-> distance (head pxcList) pxc < distanceThreshold) pxcList
    newpxcList = filter (\pxc-> distance (head pxcList) pxc > distanceThreshold) pxcList

toEdgeImage :: G.Image G.VU G.Y Double -> G.Image G.VU G.Y Double
toEdgeImage origImgM = G.imap (\(r,c) pixelY-> if (any id $
                                  map (\pos -> (>)
                                        (pixelY - (GI.handleBorderIndex (GI.Fill 1) (GI.dims origImgM) (GI.index origImgM) pos))
                                        (G.PixelY edgeSensitivity))
                                   [(r-1,c),(r+1,c),(r,c-1),(r,c+1)])
                              then G.PixelY (1::Double)
                              else G.PixelY (0::Double)) origImgM              -- Convert to bilevel edge Image

makeRect ((xMin,yMin),(xMax,yMax)) = (map (\x-> (x,yMin)) [xMin..xMax])++(map (\x->(x,yMax)) [xMin..xMax])++(map (\y->(xMax,y)) [yMin..yMax])++(map (\y->(xMin,y)) [yMin..yMax])

getAllConnected :: H.Map (Int,Int) Int -> [[(Int,Int)]]
getAllConnected edgeHM = getAllConnected' edgeHM []

getAllConnected' :: H.Map (Int,Int) Int -> [[(Int,Int)]] -> [[(Int,Int)]]
getAllConnected' edgeHM lineList = if (getNewPos edgeHM) == invalidPos
                                   then lineList
                                   else getAllConnected' newEdgeHM (lineList ++ [newLineList])
  where
    (newLineList, newEdgeHM) = getConnected (getNewPos edgeHM) edgeHM []

getConnected :: (Int,Int) -> H.Map (Int,Int) Int -> [(Int,Int)] -> ([(Int,Int)],H.Map (Int,Int) Int)
getConnected currentPos visitedHM connectedPositionsList = if null foundNeighbours
                                                           then (connectedPositionsList, newVisitedHM)
                                                           else (foundNeighbours ++ newConnectedPositionsList, newVisitedHM2)
  where
    rw = fst currentPos
    cl = snd currentPos
    possibleNeighbours = [(rw-1,cl-1),(rw-1,cl),(rw-1,cl+1),(rw,cl+1),(rw+1,cl+1),(rw+1,cl),(rw+1,cl-1),(rw,cl-1)]
    foundNeighbours = filter (\possiblePos->(H.member possiblePos visitedHM)&&(H.lookup possiblePos visitedHM==(Just 1)))
                      $ possibleNeighbours ++ [currentPos]
    newVisitedHM = foldl (\hm pos'-> H.adjust (\a-> 0) pos' hm) visitedHM foundNeighbours
    (newConnectedPositionsList, newVisitedHM2) = foldl (\(posList, hm) pos'-> (\(pL,h)-> (nub $ (pL ++ posList),h)) $ getConnected pos' hm posList) ([], newVisitedHM) foundNeighbours

getNewPos :: H.Map (Int,Int) Int -> (Int,Int)
getNewPos edgeHM = validPosCheck $ filter (\key->(H.member key edgeHM)&&(H.lookup key edgeHM==Just 1)) $ H.keys edgeHM

invalidPos = (-100,-100)

validPosCheck [] = invalidPos
validPosCheck xs = head xs

getNeigh :: (Int,Int) -> (Int,Int) -> H.Map (Int,Int) Int -> [(Int,Int)]
getNeigh pos prevPos hm = take 1 $ filter (\x->(H.member x hm) && (H.lookup x hm == Just 1)) directionList
  where
    directionalPoints = [(rw-1,cl-1),(rw-1,cl),(rw-1,cl+1),(rw,cl+1),(rw+1,cl+1),(rw+1,cl),(rw+1,cl-1),(rw,cl-1)]
    rw = fst pos
    cl = snd pos
    relativeDirection = getRelativeDirection pos prevPos
    d1 = (relativeDirection + 3)
    p1Direction = (\[x,y,z]-> [y,x,z]) $ take 3 $ drop d1 $ [0..7] ++ [0..7]
    p2Direction = [directionCircle!!(8+(p1Direction!!1) -1), directionCircle!!(p1Direction!!2 +1)]
    p3Direction = [directionCircle!!(8+(p2Direction!!0) -1), directionCircle!!(p2Direction!!1 +1)]
    directionList = map (\x-> directionalPoints!!x) $ p1Direction ++ p2Direction ++ p3Direction
    directionCircle = concat $ replicate 3 [0..7]

getRelativeDirection (rw,cl) prevPos | (rw-1,cl-1) == prevPos = 0
                                     | (rw-1,cl) == prevPos = 1
                                     | (rw-1,cl+1) == prevPos = 2
                                     | (rw,cl+1) == prevPos = 3
                                     | (rw+1,cl+1) == prevPos = 4
                                     | (rw+1,cl) == prevPos = 5
                                     | (rw+1,cl-1) == prevPos = 6
                                     | (rw,cl-1) == prevPos = 7
                                     | otherwise = 0

dfsEdge xs = recDFSEdge hashMap [(head xs)] []
  where
    hashMap = H.fromList $ zip xs (repeat 1)

recDFSEdge hm [] xs = (xs,hm)
recDFSEdge hm tup_list xs = recDFSEdge hm1 new_tup_list (xs ++ tup_list)
  where
    hm1 = H.adjust (\a-> 0) (head tup_list) hm
    new_tup_list = getNeigh (head tup_list) (if null xs then (0,0) else last xs) hm1

mySlice :: Int -> Int -> V.Vector a -> V.Vector a
mySlice s l xs = if (s < 0)
                 then (mySlice (V.length xs + s) l xs)
                 else (if ((s+l) > (V.length xs))
                       then (V.++) (mySlice s (V.length xs - s) xs) (mySlice 0 (l - (V.length xs - s)) xs)
                       else (V.slice s l xs))

fitLine :: V.Vector (Double, Double) -> (Int, Int) -> (Double, Double)
fitLine vec (s,l) = linearRegressionTLS xVec yVec
  where
    newVec = mySlice s l vec
    xVec = fst $ V.unzip newVec
    yVec = snd $ V.unzip newVec

neighbourhoodToConsiderForSlope :: Double
neighbourhoodToConsiderForSlope = 0.3

getSlopePoints :: V.Vector (Int,Int) -> V.Vector (Double,Double)
getSlopePoints xs = V.map (\(s,l)-> fitLine dxs (s,l)) (V.generate (V.length xs) getRange)
  where
    dxs = V.map (\(x,y)-> ((fromIntegral x)::Double, (fromIntegral y)::Double)) xs
    l = V.length xs
    offset = floor $ (fromIntegral l) * (neighbourhoodToConsiderForSlope / 2)
    getRange i = (i-offset, offset * 2)

myClamp i l = if i >= l then i - l else i

getMid s e l = if mid >= l then mid - l else mid
  where
    mid = (if (s > e) then (ceiling $ (fromIntegral (s + e + l - 1)) / 2) else ceiling $ (fromIntegral (s + e)) / 2) :: Int

isAlmostPerpendicular m1 m2 = if (m <= (1.6708)) && (m >= (1.4708)) then True else False
  where
    m = atan $ (m1 - m2) / (1 - m1 * m2)

detectShifts :: V.Vector (Double,Double) -> V.Vector Int
detectShifts slopeIntersectVec = V.imap (\i (c,m)-> if isAlmostPerpendicular m (snd $ (V.!) slopeIntersectVec (myClamp (i+offset) l))
                                                    then (getMid i (myClamp (i+offset) l) l) else -1) slopeIntersectVec
  where
    l = V.length slopeIntersectVec
    offset = floor $ (fromIntegral l) * neighbourhoodToConsiderForSlope

-- SCRATCHED CODE --

--edgeSet = filter (\x->(length x)>10) $ snd $ foldl (\(hm1,lds) _ -> (\(nlds,hm2)->(hm2,(nlds:lds))) $ if (null $ activeEdge hm1) then ([],H.empty) else (recDFSEdge hm1 (activeEdge hm1) [])) (hm,[]) [1..(H.size hm)]

-- newFxn g t = do
--   GI.write g t (G.PixelY (1::Double))
--   return g

-- g2 = G.imap (\(r,c) p-> if (any id $ map (\i -> (>) (p - (GI.handleBorderIndex (GI.Fill 1) (GI.dims g1) (GI.index g1) i)) (G.PixelY 0.5)) [(r-1,c),(r+1,c),(r,c-1),(r,c+1)]) then G.PixelY (1::Double) else G.PixelY (0::Double)) g1

-- e=[0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]

-- rangeAdd hm x = H.adjust (\y-> y+1) ((fromIntegral $ floor (x*10))/10 :: Double, (fromIntegral $ ceiling (x*10))/10 :: Double) hm

-- rangeAdd hm x = H.adjust (\y-> y+1) (c-0.1,c) hm
--   where
--     c = (fromIntegral $ ceiling (x*10))/10 :: Double

-- clamp mn mx val = max mn (min val mx)

-- matrixUpdate :: (a -> a) -> (Int,Int) -> V.Vector (V.Vector a) -> V.Vector (V.Vector a)
-- matrixUpdate f pos matrix = V.accum (\x y-> V.accum (\a b-> b a) x [(snd pos,y)]) matrix [(fst pos, f)]

-- threshold (r,g,b) (r1,g1,b1) = (<) 220 $ abs $ (r+g+b)-(r1+g1+b1)
-- isDarker (r,g,b) (r1,g1,b1) = (<) (div (r+g+b) 3) (div (r1+g1+b1) 3)

-- getEdgeIndices :: [[(Int,Int,(Int,Int,Int))]] -> [(Int,Int)]
-- getEdgeIndices s = concat.map (concat.map isDisc) $ s
--   where
--     (mr,mc) = (length s,length $ head s)
--     isDisc (rw,cl,p1@(r,g,b)) = map (\(x,y,p)-> if isDarker p p1 then (x,y) else (rw,cl)) $ filter (\(_,_,p2)-> threshold p1 p2 ) $ neighIndicesWP rw cl
--     neighIndicesWP rw cl = map (\(x,y) -> (s!!x)!!y) $ neighIndices rw cl
--     neighIndices rw cl = filter (\(x,y)-> x<mr && x>=0 && y<mc && y>=0) [(rw-1,cl),(rw+1,cl),(rw,cl-1),(rw,cl+1),(rw-1,cl-1),(rw+1,cl+1),(rw-1,cl+1),(rw+1,cl-1)]

-- getPts (x,y) = concat $ map (\i-> map (\j-> (x*5+i,y*5+j)) [0..4]) [0..4]
-- fun (x,y) = (G.index g ((x*5)+2,(y*5)+2))
-- g2 = G.makeImage ((\(x,y)-> (div x 5,div y 5)) (G.dims g)) fun :: GI.Image G.VU G.Y Double

-- g <- readFile "properSet"
-- g1 = read  g :: [[(Double,Double)]]
-- g2 = head $ g1
-- g31 = V.fromList $ fst (unzip g2)
-- g32 = V.fromList $ snd (unzip g2)
-- linearRegressionTLS g31 g32

max_stars_repo_name: arkouda/GridDetectionCV

max_issues_repo_name: arkouda/GridDetectionCV

max_forks_repo_name: arkouda/GridDetectionCV


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-|
Module      : Grenade.Layers.LRN
Description : LRN layer
Maintainer  : Aprova Varshney
License     : BSD2
Stability   : experimental

LRN is a non trainable layer which normalises the output of a neuron based on the values in a
neighborhood around it, especially useful when non-normalised activation functions such as ReLU are utilised.
-}

module Grenade.Layers.LRN
  (
  -- * Layer Definition
    LRN(..)
  ) where

import           Control.DeepSeq
import           GHC.TypeLits
import           Data.Proxy
import           Data.Serialize

import qualified Numeric.LinearAlgebra          as LA
import           Numeric.LinearAlgebra.Static

import           Grenade.Core
import           Grenade.Types
import           Grenade.Onnx
import           Data.Kind (Type)
import           Numeric.LinearAlgebra (Matrix)
import           Grenade.Utils.Symbols
import           Data.ProtoLens.Labels ()

data LRN :: Symbol -> Symbol -> Symbol -> Nat -> Type where
  LRN :: ( KnownSymbol a, ValidRealNum a
         , KnownSymbol b, ValidRealNum b
         , KnownSymbol k, ValidRealNum k
         , KnownNat n
         ) => LRN a b k n

instance NFData (LRN a b k n) where
  rnf LRN = ()

instance Show (LRN a b k n) where
  show LRN
    = "LRN "
    ++ (symbolVal (Proxy :: Proxy a)) ++ " "
    ++ (symbolVal (Proxy :: Proxy b)) ++ " "
    ++ (symbolVal (Proxy :: Proxy k)) ++ " "
    ++ show (natVal (Proxy :: Proxy n))

instance UpdateLayer (LRN a b k n) where
  type Gradient (LRN a b k n) = ()
  runUpdate _ x _ = x
  reduceGradient _ = ()

instance ( KnownSymbol a, ValidRealNum a
         , KnownSymbol b, ValidRealNum b
         , KnownSymbol k, ValidRealNum k
         , KnownNat n
         ) => RandomLayer (LRN a b k n) where
  createRandomWith _ _ = return LRN

instance ( KnownSymbol a, ValidRealNum a
         , KnownSymbol b, ValidRealNum b
         , KnownSymbol k, ValidRealNum k
         , KnownNat n
         ) => Serialize (LRN a b k n) where
  put _ = return ()
  get   = return LRN

instance ( KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (LRN a b k n) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols channels) where

  type Tape (LRN a b k n) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols channels)
    = (L (inputRows * channels) inputCols, L (inputRows * channels) inputCols)

  runForwards LRN (S3D input) = ((input, sums), S3D normalized)
    where
      normalized :: L (inputRows * channels) inputCols
      normalized = build buildNormalizedMatrix

      ex = extract input :: Matrix RealNum
      channels = fromIntegral $ natVal (Proxy :: Proxy channels)  :: Int
      rows     = fromIntegral $ natVal (Proxy :: Proxy inputRows) :: Int

      n  = natVal (Proxy :: Proxy n)
      a  = (read $ symbolVal (Proxy :: Proxy a)) :: RealNum
      b  = (read $ symbolVal (Proxy :: Proxy b)) :: RealNum
      k  = (read $ symbolVal (Proxy :: Proxy k)) :: RealNum

      sums  = build buildSums
      sums' = extract sums

      buildSums :: RealNum -> RealNum -> RealNum
      buildSums i j = k + a * summation
        where
          i' = floor i :: Int
          j' = floor j :: Int
          ch = i' `div` rows :: Int
          ro = i' `mod` rows :: Int
          co = j'

          sub = floor ((fromIntegral n) / 2     :: RealNum)
          add = floor ((fromIntegral n - 1) / 2 :: RealNum)
          lower = maximum [0, ch - sub]
          upper = minimum [channels - 1, ch + add]
          summation = sum [ (ex `LA.atIndex` (q * rows + ro, co)) ** 2 | q <- [lower..upper]]

      -- Calculates the normalised values for each cell of the input based
      -- on the cross-channel sums of the surrounding cells
      buildNormalizedMatrix :: RealNum -> RealNum -> RealNum
      buildNormalizedMatrix i j = val
        where
          i' = floor i :: Int
          j' = floor j :: Int

          f = ex `LA.atIndex` (i', j')
          val = f / den

          den  = den' ** b
          den' = sums' `LA.atIndex` (i', j')

  runBackwards LRN (input, sums) (S3D err) = ((), S3D err')
    where
      ex    = extract input :: Matrix RealNum
      exr   = extract err   :: Matrix RealNum
      sums' = extract sums  :: Matrix RealNum

      channels = fromIntegral $ natVal (Proxy :: Proxy channels)  :: Int
      rows     = fromIntegral $ natVal (Proxy :: Proxy inputRows) :: Int

      n  = natVal (Proxy :: Proxy n)
      a  = (read $ symbolVal (Proxy :: Proxy a)) :: RealNum
      b  = (read $ symbolVal (Proxy :: Proxy b)) :: RealNum

      err' = build buildErrorMatrix

      buildErrorMatrix :: RealNum -> RealNum -> RealNum
      buildErrorMatrix i j = res
        where
          i' = floor i :: Int
          j' = floor j :: Int
          ch = i' `div` rows
          ro = i' `mod` rows
          co = j'

          sub = floor ((fromIntegral n) / 2     :: RealNum)
          add = floor ((fromIntegral n - 1) / 2 :: RealNum)
          lower = maximum [0, ch - sub]
          upper = minimum [channels - 1, ch + add]

          s  = (ex `LA.atIndex` (i', j')) * s'
          s' = sum [ (ex `LA.atIndex` (q * rows + ro, co)) * (exr `LA.atIndex` (q * rows + ro, co)) | q <- [lower..upper] ]
          c  = sums' `LA.atIndex` (i', j')

          res = (exr `LA.atIndex` (i', j')) * (c ** (-b)) - 2 * b * a * (c ** (-b - 1)) * s

instance OnnxOperator (LRN a b k n) where
  onnxOpTypeNames _ = ["LRN"]

instance ( KnownSymbol a, ValidRealNum a
         , KnownSymbol b, ValidRealNum b
         , KnownSymbol k, ValidRealNum k
         , KnownNat n
         )  => OnnxLoadable (LRN a b k n) where
  loadOnnxNode _ node = do
    hasMatchingRealNum node (Proxy :: Proxy a) "alpha"
    hasMatchingRealNum node (Proxy :: Proxy b) "beta"
    hasMatchingRealNum node (Proxy :: Proxy k) "bias"
    hasMatchingInt     node (Proxy :: Proxy n) "size"
    return LRN

ext: hs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-- | This module defines the 'write-xx' style primitives.

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
module Primitives.WriteLib
  ( writeSH, writeSimpleSH, writeSharedSH, displaySH
  , showErrIO

  , ushowString
  ) where

import Val
import Primitives.Vector (vectorElemsPH)
import Primitives.String (stringSH, unwrapStringPH)

import System.Mem.StableName (StableName, makeStableName)

import Control.Monad.Reader (ReaderT(runReaderT), asks)
import Control.Monad.State
import qualified Data.HashSet as S
import qualified Data.HashMap.Strict as M
import GHC.Show (showLitChar)
import Data.Char (isPrint)
import Data.Complex (Complex(..))
import Data.Ratio (numerator, denominator)

writeSH :: Val -> IO String
writeSH v = do
  labeling <- labelDatum CyclicOnly v
  runShowS <$> withWriteComponents labeling (writeShowS v)

writeSimpleSH :: Val -> IO String
writeSimpleSH v =
  runShowS <$> withWriteComponents M.empty (writeShowS v)

writeSharedSH :: Val -> IO String
writeSharedSH v = do
  labeling <- labelDatum FullSharing v
  runShowS <$> withWriteComponents labeling (writeShowS v)

displaySH :: Val -> IO String
displaySH v = do
  labeling <- labelDatum CyclicOnly v
  runShowS <$> withDisplayComponents labeling (writeShowS v)

-- | Run a 'Write' action with the given labeling, using
-- component writers for the 'writeSH' family.
withWriteComponents :: Labeling -> Write a -> IO a
withWriteComponents = runWrite charHH stringHH symbolHH

-- | Run a 'Write' action with the given labeling, using
-- component writers for the 'displaySH' function.
withDisplayComponents :: Labeling -> Write a -> IO a
withDisplayComponents = runWrite writeCharHH writeStringHH writeSymbolHH

runShowS :: ShowS -> String
runShowS f = f ""

-- | Show a 'LispErr' in IO. Forms inside the error are shown with
-- 'write-shared'.
showErrIO :: LispErr -> IO String
showErrIO = fmap (prefix ++) . mkMsgFor
  where
    prefix = "Error: "
    mkMsgFor (UnboundVar msg name)
      = pure $ msg ++ " unbound symbol" `colonAnd` symbolAsString name
    mkMsgFor (EvaluateDuringInit name) = pure $
      symbolAsString name ++ " referred to itself during initialization"
    mkMsgFor (SetImmutable tyname) = pure $ "can't set immutable " ++ tyname
    mkMsgFor (BadSpecialForm form) =
      ("bad form" `colonAnd`) <$> writeSharedSH form
    mkMsgFor (NotFunction msg form) =
      (msg `colonAnd`) <$> writeSharedSH form
    mkMsgFor (NumArgs arity act) =
      ((expectedMsg arity ++ ", found values") ++) . concat
      <$> mapM (fmap (' ':) . writeSharedSH) act
      where expectedMsg a = "expected " ++ showArity a ++ " arg" ++ p
              where p | usePluralForArity a = "s"
                      | otherwise = ""
    mkMsgFor (TypeMismatch exp act) = (expectedMsg exp ++) <$> writeSharedSH act
      where expectedMsg s =
              "invalid type: expected " ++ s ++ ", found "
    mkMsgFor CircularList = pure "circular list"
    mkMsgFor EmptyBody = pure "attempt to define function with no body"
    mkMsgFor (Parser parseErr) = pure $ "parser error at " ++ show parseErr
    mkMsgFor (Condition _ obj) = writeSharedSH obj
    mkMsgFor (Default msg) = pure msg
    mkMsgFor Quit = pure "quit invoked"

    colonAnd s r = s ++ ": " ++ r

-------------------------------------------------------------------------------
-- HH versions of Scheme primitive functions
-- These are NOT the conversion functions used by writeSH!
-------------------------------------------------------------------------------

writeCharHH :: Char -> ShowS
writeCharHH c = (c:)

writeStringHH :: Text -> ShowS
writeStringHH t = showString $ unpack t

writeSymbolHH :: Symbol -> ShowS
writeSymbolHH t = showString $ unpack t

-------------------------------------------------------------------------------
-- HH conversion functions used by writeSH and friends
-------------------------------------------------------------------------------

charHH :: Char -> ShowS
charHH = shows . Char

stringHH :: Text -> ShowS
-- ushowString behaves like shows, not showString.
stringHH t = ushowString $ unpack t

-- TODO [r7rs]
-- these should be escaped according to r7rs, but we don't want
-- them to be escaped at the REPL prompt. 'repl-write'?
symbolHH :: Symbol -> ShowS
symbolHH t = showString $ unpack t

-------------------------------------------------------------------------------
-- Core implementation
-------------------------------------------------------------------------------

-- | stablename -> (label, used?)
type FirstPassMap = M.HashMap (StableName Val) Bool
type Labeling = M.HashMap (StableName Val) Int
type FirstPass = StateT FirstPassMap IO ()
data Sharing = FullSharing | CyclicOnly

labelDatum :: Sharing -> Val -> IO Labeling
labelDatum s v = fpmToLabeling <$>
  execStateT (labelDatumWorker s v) M.empty

labelDatumWorker :: Sharing -> Val -> FirstPass
labelDatumWorker s !p
  | pairSH p = labelPair s p
  | vectorSH p = labelVector s p
labelDatumWorker _s _notNested = pure ()

-- | Calling this on a non-pair will panic!
labelPair :: Sharing -> Val -> FirstPass
labelPair sharing pair = do
  name <- lift $ makeStableName pair
  keepGoing <- state $ \m -> case M.lookup name m of
        Just False -> (False, M.adjust (const True) name m)
        Just True  -> (False, m)
        Nothing -> (True, M.insert name False m)
  when keepGoing $ do
    car <- carPS pair
    labelDatumWorker sharing car
    cdr <- cdrPS pair
    labelDatumWorker sharing cdr
    case sharing of
      FullSharing -> pure ()
      CyclicOnly -> modify $ deleteUnused name

-- | See 'labelPair'.
labelVector :: Sharing -> Val -> FirstPass
labelVector sharing vec = do
  name <- lift $ makeStableName vec
  keepGoing <- state $ \m -> case M.lookup name m of
    Just False -> (False, M.adjust (const True) name m)
    Just True  -> (False, m)
    Nothing    -> (True, M.insert name False m)
  when keepGoing $ do
    vecElems <- vectorElemsPH vec
    mapM_ (labelDatumWorker sharing) vecElems
    case sharing of
      FullSharing -> pure ()
      CyclicOnly  -> modify $ deleteUnused name

deleteUnused :: StableName Val -> FirstPassMap -> FirstPassMap
deleteUnused = M.alter f
  where
    f Nothing = Nothing
    f (Just False) = Nothing
    f (Just True)  = Just True

fpmToLabeling :: FirstPassMap -> Labeling
fpmToLabeling fpm = evalState (M.traverseWithKey go used) 0
  where
    go _k _v = state $ \i -> (i, i+1)
    used = M.filter id fpm

type WriteState = S.HashSet (StableName Val)
type ComponentWriter a = a -> ShowS
data WriteEnv = WriteEnv
  { labeling :: Labeling
  , writeChar :: ComponentWriter Char
  , writeString :: ComponentWriter Text
  , writeSymbol :: ComponentWriter Symbol
  }
type Write a = ReaderT WriteEnv (StateT WriteState IO) a

runWrite :: ComponentWriter Char   -- ^ how to print characters
         -> ComponentWriter Text   -- ^ how to print (Scheme) strings
         -> ComponentWriter Symbol -- ^ how to print symbols
         -> Labeling               -- ^ mapping that witnesses sharing in v
         -> Write a                -- ^ Write action that computes v -> a
         -> IO a
runWrite wc wt ws lbls m =
  flip evalStateT S.empty $ runReaderT m $ WriteEnv lbls wc wt ws

writeShowS :: Val -> Write ShowS
writeShowS v
  | pairSH v = writePair v
  | vectorSH v = writeVector v
  | stringSH v = do
    ws <- asks writeString
    ws <$> unwrapStringPH v

-- we have to write symbols in all their glory because
-- literal (quoted) symbols need to be printed correctly.
-- If they aren't printed correctly at a prompt, a user
-- would merely be confused - but if they aren't written
-- correctly to a file, then the expander can't expand itself.
writeShowS (Symbol s) = do
  ws <- asks writeSymbol
  return $ ws s

writeShowS (Char c) = do
  wc <- asks writeChar
  return $ wc c

writeShowS (Number n) = pure $ writeNumber n
  where writeNumber (Real (Bignum i)) = shows i
        writeNumber (Real (Ratnum r))
          | denominator r == 1 = shows $ numerator r
          | otherwise = shows (numerator r)
                      . showChar '/'
                      . shows (denominator r)
        writeNumber (Real (Flonum f))
          | isInfinite f = showString $ sign : "inf.0"
          | isNaN f      = showString "+nan.0"
          | otherwise    = shows f
          where sign = if f < 0 then '-' else '+'
        writeNumber (Complex (r :+ i))
          | isExactZero $ Real i = writeNumber (Real r)
          | isExactZero $ Real r = writeNumber (Real i) . showChar 'i'

          -- minus sign will be given since i is negative
          | i < 0 || isNegativeZero i = writeNumber (Real r)
                                      . writeNumber (Real i)
                                      . showChar 'i'
          | otherwise = writeNumber (Real r)
                      . (if isInfinite i || isNaN i then id else showChar '+')
                      . writeNumber (Real i)
                      . showChar 'i'

writeShowS (MultipleValues vs) = fmt <$> writeList (makeImmutableList vs)
  where
    fmt vsStr = showString "#<values: " . vsStr . showChar '>'

writeShowS (Exception e) = showString <$> liftIO (showErrIO e)
writeShowS (Error msg irritants) = (prefix (unpack msg) .) <$> writeIrritants
  where
    prefix "" = showString ""
    prefix s  = showString s . showChar ' '
    writeIrritants = case irritants of
      [] -> pure $ showString ""
      _  -> writeList $ makeImmutableList irritants

writeShowS v = pure $ shows v

writePair :: Val -> Write ShowS
writePair p = do
  name <- liftIO $ makeStableName p
  mlbl <- askLabel name
  case mlbl of
    Nothing -> showParen True <$> writeList p
    Just lbl -> writeLabeledPair lbl name p

writeVector :: Val -> Write ShowS
writeVector v = do
  name <- liftIO $ makeStableName v
  elems <- vectorElemsPH v
  mlbl <- askLabel name
  case mlbl of
    Nothing -> (showChar '#' .) . showParen True <$> writeVectorElems elems
    Just lbl -> writeLabeledVector lbl name elems

writeLabeledPair :: Int -> StableName Val -> Val -> Write ShowS
writeLabeledPair lbl name p = do
  defined <- labelDefined name
  if defined
    then return $ showLabelUse lbl
    else do
      defineLabel name -- mark defined
      let pre = showLabelDef lbl
      (pre .) . showParen True <$> writeList p

writeLabeledVector :: Int -> StableName Val -> [Val] -> Write ShowS
writeLabeledVector lbl name p = do
  defined <- labelDefined name
  if defined
    then return $ showLabelUse lbl
    else do
      defineLabel name
      let pre = showLabelDef lbl . showChar '#'
      (pre .) . showParen True <$> writeVectorElems p

showLabelUse, showLabelDef :: Int -> ShowS
showLabelUse lbl = showChar '#' . shows lbl . showChar '#'
showLabelDef lbl = showChar '#' . shows lbl . showChar '='

-- we have to be careful here to correctly handle the case where the
-- cdr needs to be labeled. If the cdr needs to be labeled, we have
-- to use dotted notation immediately. This list ends, and we print
-- the cdr (which may still need its label defined, in which case it
-- will appear as a labeled list).
writeList :: Val -> Write ShowS
writeList pair = do
  car <- liftIO $ carPS pair
  showsCar <- writeShowS car
  (showsCar .) <$> (cdrPS pair >>= writeList1)

writeList1 :: Val -> Write ShowS
writeList1 p | pairSH p = do
  name <- liftIO $ makeStableName p
  mlbl <- askLabel name
  case mlbl of
    Nothing -> (showChar ' ' .) <$> writeList p
    Just lbl -> (showString " . " .) <$>
      writeLabeledPair lbl name p
writeList1 Nil = pure id
writeList1 dot = (showString " . " .) <$> writeShowS dot

writeVectorElems :: [Val] -> Write ShowS
writeVectorElems [] = pure id
writeVectorElems [val] = writeShowS val
writeVectorElems (v:vs) = spaceSep <$> writeShowS v <*> writeVectorElems vs
  where
    spaceSep s1 s2 = s1 . showChar ' ' . s2

labelDefined :: StableName Val -> Write Bool
labelDefined n = gets $ S.member n

defineLabel :: StableName Val -> Write ()
defineLabel n = modify $ S.insert n

askLabel :: StableName Val -> Write (Maybe Int)
askLabel name = asks $ M.lookup name . labeling

-- | So this is funny.
--
-- Haskell's 'Show' instance for 'Char' defines 'showList', this is what
-- is used to 'show' Strings. If you chase the source a bit, you'll find that
-- it calls 'showLitChar' - a function that does not propogate Unicode
-- characters! Which is pretty silly. This leads to the awkward situation where
-- we print #\λ correctly, but not (list->string '(#\λ)).
--
-- This function corrects the case of 'showLitChar'
-- that escapes Unicode code points.
ushowString :: String -> ShowS
ushowString s = showChar '\"' . showLitString s . showChar '\"'
  where
    showLitString [] s = s
    showLitString (c:cs) s = ushowLitChar c (showLitString cs s)

    ushowLitChar c s
      | c > '\DEL' && isPrint c = showChar c s
      | otherwise = showLitChar c s

max_stars_repo_name: JKTKops/Thcheme

max_issues_repo_name: JKTKops/Thcheme

max_forks_repo_name: JKTKops/Thcheme


text:
{-# LANGUAGE DataKinds,EmptyDataDecls #-}

module HLearn.Metrics.Mahalanobis.Lego
    where

import Control.DeepSeq
import Data.List
import qualified Data.Semigroup as SG
import qualified Data.Foldable as F
import qualified Data.Vector.Generic as VG

import Foreign.Storable
import Numeric.LinearAlgebra hiding ((<>))
import qualified Numeric.LinearAlgebra as LA

import HLearn.Algebra
import HLearn.Metrics.Mahalanobis
import HLearn.Models.Distributions.Multivariate.MultiNormalFast

-------------------------------------------------------------------------------
-- data types

type Lego reg (eta::Frac) dp = AddUnit1 (Lego' reg eta) dp

data Lego' reg (eta::Frac) dp = Lego'
    { b :: !(Matrix (Ring dp))
    , c :: !(Matrix (Ring dp))
    , x :: Matrix (Ring dp)
    , multinorm :: MultiNormal dp
    }

deriving instance (Element (Ring dp), Show (Ring dp)) => Show (Lego' reg eta dp)

instance NFData dp => NFData (Lego' reg eta dp) where
    rnf lego = deepseq b
             $ deepseq c
             $ rnf x

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

data RegMatrix_ = Identity_ | InvCovar_

class RegMatrix m where
    regmatrix :: m -> RegMatrix_

data Identity
instance RegMatrix Identity where
    regmatrix _ = Identity_

data InvCovar
instance RegMatrix InvCovar where
    regmatrix _ = InvCovar_

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

mkLego :: forall reg eta dp.
    ( SingI eta
    , MatrixField dp
    , RegMatrix reg
    ) => Matrix (Ring dp) -> Matrix (Ring dp) -> MultiNormal dp -> Lego' reg eta dp
mkLego b c multinorm = Lego'
    { b = b
    , c = c
    , multinorm = multinorm
    , x = inv $ completeTheSquare (scale (-1) b') (scale (-eta) c)
    }
    where
        m = case (regmatrix (undefined::reg)) of
            Identity_ -> ident (rows b)
            InvCovar_ -> inv $ covar multinorm
        b' = scale (1/2) (m `add` scale eta b)
        eta = fromSing (sing :: Sing eta)

-- | given an equation Y^2 + BY + YB + C, solve for Y
completeTheSquare ::
    ( Storable r
    , Element r
    , Container Vector r
    , LA.Product r
    , Field r
    , r ~ Double
    ) => Matrix r -> Matrix r -> Matrix r
-- completeTheSquare b c = (sqrtm' ((b LA.<> b) `sub` c)) `sub` b
completeTheSquare b c = (sqrtm' ((scale (1/4) $ b LA.<> b) `sub` c)) `sub` b
    where
        sqrtm' m = cmap realPart $ matFunc sqrt $ cmap (:+0) m

-------------------------------------------------------------------------------
-- algebra

instance (SingI eta, RegMatrix reg, MatrixField dp) => SG.Semigroup (Lego' reg eta dp) where
    lego1 <> lego2 = mkLego  b' c' mn'
        where
            b' = b lego1 `add` b lego2
            c' = c lego1 `add` c lego2
            mn' = multinorm lego1 <> multinorm lego2

instance HasRing dp => HasRing (Lego' reg eta dp) where
    type Ring (Lego' reg eta dp) = Ring dp

instance Num r => HasRing (Vector r) where
    type Ring (Vector r) = r

-------------------------------------------------------------------------------
-- training

instance
    ( SingI eta
    , RegMatrix reg
    , MatrixField (Vector r)
    , r ~ Ring (dp r), F.Foldable dp
    ) => HomTrainer (Lego reg eta (dp r))
        where
    type Datapoint (Lego reg eta (dp r)) = (r,dp r)

    train1dp (y,dp) = UnitLift1 $ mkLego (scale y zzT) (zzT LA.<> zzT) (train1dp dp)
        where
            zzT = asColumn dp' LA.<> asRow dp'
            dp' = fromList $ F.toList dp

-------------------------------------------------------------------------------
-- metric

instance
    ( VG.Vector dp r
    , Ring (dp r) ~ r
    , RegMatrix reg
    , LA.Product r
    ) =>  MkMahalanobis (Lego reg eta (dp r))
        where
    type MetricDatapoint (Lego reg eta (dp r)) = dp r
    mkMahalanobis (UnitLift1 lego) dp = Mahalanobis
        { rawdp = dp
        , moddp = VG.fromList $ toList $ flatten $ (x lego) LA.<> asColumn v
        }
        where
            v = fromList $ VG.toList dp

-------------------------------------------------------------------------------
-- test

xs =
    [ [1,2]
    , [2,3]
    , [3,4]
    ]
    :: [[Double]]

mb = (2><2)[-6.5,-10,-10,-14]      :: Matrix Double
mc = (2><2)[282,388,388,537] :: Matrix Double

max_stars_repo_name: silky/HLearn

max_issues_repo_name: silky/HLearn

max_forks_repo_name: silky/HLearn


text:
module Main where

import Lib
import Mandelbrot
import ParseInput

import System.IO
import Data.Complex
import Data.Array.Repa hiding ((++))
import Data.Array.Repa.IO.Matrix
import Data.Array.Repa.IO.BMP
import Data.Word
import System.Console.CmdArgs
import System.Environment (getArgs)

mandelbrotZoom :: String -> Complex Double -> Double -> Int -> Int -> Double -> Int -> Int -> IO ()
mandelbrotZoom _ _ _ _ _ _ 0 _ = return ()
mandelbrotZoom filebase centre size res n scale t maxt = do
    cs <- computeP $ genMandelbrot res n size centre
    writeImageToBMP (filebase++"_"++(show (maxt-t))++".bmp") cs
    putStrLn $ "Finished image "++(show (maxt-t))
    mandelbrotZoom filebase centre (size*scale) res n scale (t-1) maxt

mainSingle :: IO ()
mainSingle = do
  [i',j',size',res',n',filename] <- getArgs
  let res = read res'      :: Int
      n = read n'          :: Int
      size = read size'    :: Double
      i = read i'          :: Double
      j = read j'          :: Double
      centre = i :+ j
      cs' = genMandelbrot res n size centre
  cs <- computeP cs' :: IO (Array U DIM2 (Word8,Word8,Word8))
  writeImageToBMP filename cs

mainZoom :: InOpts -> IO ()
mainZoom inopts = do
  --[i',j',size',res',n',filebase,t'] <- getArgs
  let res = res_ inopts
      n = iters_ inopts
      size = width_ inopts
      i = i_ inopts
      j = j_ inopts
      t = zooms_ inopts
      s = scale_ inopts
      filebase = filepath_ inopts
      centre = i :+ j
  mandelbrotZoom filebase centre size res n s t t

main :: IO ()
main = mainZoom =<< (cmdArgs $ defopts)

max_stars_repo_name: ahgibbons/HaskMandelbrot

max_issues_repo_name: ahgibbons/HaskMandelbrot

max_forks_repo_name: ahgibbons/HaskMandelbrot