| Maths, Stats and Functional Progamming | ||
| Blogs by Idontgetoutmuch | ||
| Latest | Previous | About |
Suppose you have some array-like data stored in NetCDF such as surface sea temperature at given positions (latitude and longitude) over time. It’s possible to do APL-like programming in Haskell using Naperian functors. APL-like programming seems like a good fit for analysing such data but how does one get such data from a NetCDF file into an array in which each dimension is typed with the size of the data in that dimension?
Data.NetCDF.Store contains the NcStore store class. So all we need to is implement this class for the appropriate structure in the Naperian package.
Before we start having fun with pointers, let’s try implementing a simpler but similar class.
First we need some extensions and to import the necessary modules.
> {-# LANGUAGE DataKinds #-}
> {-# LANGUAGE FlexibleInstances #-}
> {-# LANGUAGE TypeFamilies #-}
> {-# LANGUAGE ExplicitForAll #-}
> {-# LANGUAGE TypeOperators #-}
> {-# LANGUAGE FlexibleContexts #-}
> {-# LANGUAGE InstanceSigs #-}
> {-# LANGUAGE ScopedTypeVariables #-}
> {-# LANGUAGE OverloadedLists #-}
> {-# LANGUAGE Rank2Types #-}
> module NaperianNetCDF where
> import Naperian
> import Data.Proxy
> import GHC.TypeLits
> import Data.List.Split
> import Data.Maybe
> import Data.Foldable
> import Data.NetCDF.Store
> import Data.NetCDF
> import qualified Data.Vector.Storable as SV
> import GHC.ForeignPtr
> import Foreign.Storable
> import Foreign.C
> import GHC.Ptr
> import Control.Monad
Here’s the class definition.
> class ListLike fs where
> toListLike :: fs e -> [e]
> unToListLike :: [e] -> fs e
The base instance for the hypercuboid of the empty list of types is straightforward.
> instance ListLike (Hyper '[]) where
> toListLike (Scalar c) = [c]
> unToListLike [c] = Scalar c
> unToListLike cs = error $ "ToListLike: wrong length = " ++ show (Prelude.length cs)
The instance for recursion also seems simple in retrospect but I took many wrong turns before getting here.
> instance (KnownNat n, ListLike (Hyper fs), Shapely fs) =>
> ListLike (Hyper ((Vector n) : fs)) where
> toListLike (Prism c) = concat $ toListLike $ fmap toList c
> unToListLike :: forall e . [e] -> Hyper ((Vector n) : fs) e
> unToListLike cs = Prism $ unToListLike us
> where
> us :: [Vector n e]
> us = map (fromJust . fromList) $
> chunksOf n cs
>
> n = fromIntegral $ natVal (Proxy :: Proxy n)
Let’s check things are working as intended.
> roundTrip :: Hyper '[Vector 5, Vector 2, Vector 3] Int ->
> Hyper '[Vector 5, Vector 2, Vector 3] Int
> roundTrip = unToListLike . toListLike
> v5 :: Hyper '[Vector 5, Vector 2, Vector 3] Int
> v5 = Prism (Prism (Prism (Scalar a)))
> a :: Vector 3 (Vector 2 (Vector 5 Int))
> a = fromJust $ fromList $ [b, b1, b2]
> b :: Vector 2 (Vector 5 Int)
> b = [ [ 1, 2, 3, 4, 5 ]
> , [ 6, 7, 8, 9, 10 ]
> ]
> b1 :: Vector 2 (Vector 5 Int)
> b1 = [ [ 11, 12, 13, 14, 15 ]
> , [ 16, 17, 18, 19, 20 ]
> ]
> b2 :: Vector 2 (Vector 5 Int)
> b2 = [ [ 21, 22, 23, 24, 25 ]
> , [ 26, 27, 28, 29, 30 ]
> ]
*NaperianNetCDF> roundTrip v5
<<<1,2,3,4,5>,<6,7,8,9,10>>,
<<11,12,13,14,15>,<16,17,18,19,20>>,
<<21,22,23,24,25>,<26,27,28,29,30>>>
*NaperianNetCDF> v5
<<<1,2,3,4,5>,<6,7,8,9,10>>,
<<11,12,13,14,15>,<16,17,18,19,20>>,
<<21,22,23,24,25>,<26,27,28,29,30>>>Now let’s implement an instance for the type class that we really want.
Again the base instance is straightforward.
> instance NcStore (Hyper '[]) where
> type instance NcStoreExtraCon (Hyper '[]) a = ()
> toForeignPtr = fst . SV.unsafeToForeignPtr0 . SV.fromList . elements
> fromForeignPtr p _ = Scalar . head . SV.toList $ SV.unsafeFromForeignPtr0 p 1
> smap = error "smap: NcStore (Hyper '[])"
The instance for recursion is not much harder.
> instance (KnownNat n, Shapely fs, NcStore (Hyper fs)) =>
> NcStore (Hyper ((Vector n) : fs)) where
> type NcStoreExtraCon (Hyper ((Vector n) : fs)) e = NcStoreExtraCon (Hyper fs) (Vector n e)
> toForeignPtr :: forall e . Storable e =>
> Hyper ((Vector n) : fs) e -> ForeignPtr e
> toForeignPtr = fst . SV.unsafeToForeignPtr0 . SV.fromList . elements
> fromForeignPtr :: forall e . (Storable e , NcStoreExtraCon (Hyper fs) (Vector n e)) =>
> ForeignPtr e -> [Int] -> Hyper ((Vector n) : fs) e
> fromForeignPtr p _ = Prism ws
> where
> q :: ForeignPtr (Vector n e)
> q = castForeignPtr p
> ws :: Hyper fs (Vector n e)
> ws = fromForeignPtr q undefined
> smap = error "smap: NcStore (Hyper ((Vector n) : fs))"
This won’t compile complaining that there is no instance for Storable for Vector n a. So let’s create such an instance.
> instance (KnownNat n, Storable e) => Storable (Vector n e) where
> poke p v = zipWithM_ poke qs (toList v)
> where
> qs = map (\i -> plusPtr q (l * i)) [0 .. n - 1]
> n = fromIntegral $ natVal (Proxy :: Proxy n)
> l = sizeOf (undefined :: e)
> q = castPtr p
> peek p = do vs <- mapM peek qs
> return $ fromJust $ fromList vs
> where
> q :: Ptr e
> q = castPtr p
> qs = map (\i -> plusPtr q (l * i)) [0 .. n - 1]
> l = sizeOf (undefined :: e)
> n = fromIntegral $ natVal (Proxy :: Proxy n)
> sizeOf _ = n * sizeOf (undefined :: e)
> where
> n = fromIntegral $ natVal (Proxy :: Proxy n)
> alignment _ = alignment (undefined :: e)
One thing I have glossed over is the fact that the NcStore typeclass has a constraint. We can’t ignore this in our instance. Effectively we are saying that the constraint NcStoreExtraCon (Hyper ((Vector n) : fs)) e is satisfied if NcStoreExtraCon (Hyper fs) (Vector n e) is satisifed. This constraint also has to be in the defintion of fromForeignPtr so that the constraint on the call of fromForeignPtr at the smaller type is satisfied. Ultimately these constraints all evaluate (at compile time) to ().
Now we can do some analysis of our data knowing that at compile time it has consistent dimensions.
> f :: (KnownNat m, KnownNat n) =>
> Hyper '[Vector m, Vector n] Double ->
> Hyper '[Vector m, Vector n] Double ->
> IO ()
> f x y = do
> print "Maximum difference in latitude bounds"
> print $ foldrH max (read "-Infinity" :: Double) $
> fmap abs $
> foldrH (-) 0 $ transposeH x
> print "Maximum difference in longitude bounds"
> print $ foldrH max (read "-Infinity" :: Double) $
> fmap abs $
> foldrH (-) 0 $ transposeH y
> return ()
Of course we need a way of tying together the dynamic nature of external data with our statically consistent program.
> type NaperianRet2 m n a = IO (Either NcError (Hyper '[Vector m, Vector n] a))
> withNc :: NcInfo NcRead ->
> (forall m n . (KnownNat m, KnownNat n) =>
> Hyper '[Vector m, Vector n] Double ->
> Hyper '[Vector m, Vector n] Double ->
> IO ()) ->
> IO ()
> withNc nc f = do
> case ncVar nc "lat_bnds" of
> Nothing -> error "Missing lat_bnds"
> Just latBnds -> case ncVar nc "lon_bnds" of
> Nothing -> error "Missing lon_bnds"
> Just lonBnds -> do
> let dims = ncVarDims latBnds
> let ls = map ncDimLength dims
> case someNatVal (fromIntegral (ls!!0)) of
> Nothing -> error "static / dynamic mismatch"
> Just (SomeNat (_ :: Proxy aa)) ->
> case someNatVal (fromIntegral (ls!!1)) of
> Nothing -> error "static / dynamic mismatch"
> Just (SomeNat (_bb :: Proxy bb)) -> do
> elat <- get nc latBnds :: NaperianRet2 aa bb CDouble
> case elat of
> Left ncErr1 -> error $ show ncErr1
> Right lat -> do
> elon <- get nc lonBnds :: NaperianRet2 aa bb CDouble
> case elon of
> Left ncErr2 -> error $ show ncErr2
> Right lon -> f (fmap realToFrac lat) (fmap realToFrac lon)
> main :: IO ()
> main = do
> enc <- openFile "tos_O1_2001-2002.nc"
> case enc of
> Left e -> error $ show e
> Right nc -> withNc nc f
And finally get our answer (for what it’s worth):
"Maximum difference in latitude bounds"
85.0
"Maximum difference in longitude bounds"
170.0