haskell bindings

5 files changed, 161 insertions, 1 deletions

diff --git a/Numeric/CODA/PCA.hs b/Numeric/CODA/PCA.hs
new file mode 100644
index 0000000..01fa3f1
--- /dev/null
+++ b/Numeric/CODA/PCA.hs
@@ -0,0 +1,92 @@
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell #-}
+
+-- |
+-- Module      : Numeric.CODA.PCA
+-- Description : Principal composition on untransformed compositional data (CoDA)
+-- Copyright   : (c) Justin Bedo 2023
+-- License     : MIT
+-- Maintainer  : cu@cua0.org
+-- Stability   : experimental
+--
+-- This module decomposes a (non-negative) matrix into @k@ components and computes
+-- an embedding in log space for futher analysis.
+module Numeric.CODA.PCA (PCA (..), pca) where
+
+import Data.Maybe
+import qualified Data.Vector.Storable as V
+import qualified Data.Vector.Storable.Mutable as VM
+import Foreign.C.Types
+import qualified Language.C.Inline as C
+import Numeric.LinearAlgebra.Data
+import System.IO.Unsafe
+
+-- | Datatype to hold the results of a decomposition
+data PCA = PCA
+  { -- | The decomposed row factors (n × k)
+    rowFactors :: Matrix Double,
+    -- | The decomposed column factors (m × k)
+    colFactors :: Matrix Double,
+    -- | Log-space embedding (n × m)
+    embedding :: Matrix Double,
+    -- | Final Bregman loss after optimisation
+    loss :: Double
+  }
+  deriving (Show)
+
+C.include "<../../../pca.h>"
+
+-- | Decomposes a non-negative matrix into the specified number components on the CoDA simplex
+-- via Bregman divergences
+pca ::
+  -- | An optional quantile to use as the reference; defaults to 0.5
+  Maybe Double ->
+  -- | number of components to decompose into
+  Int ->
+  -- | matrix to decompose
+  Matrix Double ->
+  PCA
+pca q k x = unsafePerformIO $ do
+  let (n, m) = size x
+      k' = fromIntegral k
+      n' = fromIntegral n
+      m' = fromIntegral m
+      q' = realToFrac $ fromMaybe 0.5 q
+  b <- VM.new (n * k)
+  c <- VM.new (m * k)
+  y <- VM.new (n * m)
+  l <- VM.new 1
+  V.unsafeWith (V.map realToFrac $ flatten x) $ \xptr ->
+    VM.unsafeWith l $ \lptr ->
+      VM.unsafeWith b $ \bptr ->
+        VM.unsafeWith c $ \cptr ->
+          VM.unsafeWith y $ \yptr ->
+            [C.block|void {
+           struct futhark_context_config *cfg = futhark_context_config_new();
+           struct futhark_context *ctx = futhark_context_new(cfg);
+           struct futhark_f64_2d *Xf = futhark_new_f64_2d(ctx, $(double *xptr), $(int n'), $(int m'));
+           struct futhark_f64_2d *Yf;
+           struct futhark_f64_2d *Bf;
+           struct futhark_f64_2d *Uf;
+
+           futhark_entry_pcaWithQuantile(ctx, &Bf, &Uf, &Yf, $(double *lptr), $(double q'), $(int k'), Xf);
+           futhark_context_sync(ctx);
+           futhark_values_f64_2d(ctx, Bf, $(double *bptr));
+           futhark_values_f64_2d(ctx, Uf, $(double *cptr));
+           futhark_values_f64_2d(ctx, Yf, $(double *yptr));
+
+           futhark_free_f64_2d(ctx, Xf);
+           futhark_free_f64_2d(ctx, Bf);
+           futhark_free_f64_2d(ctx, Uf);
+           futhark_free_f64_2d(ctx, Yf);
+           futhark_context_free(ctx);
+           futhark_context_config_free(cfg);
+        }|]
+  bF <- V.freeze b
+  cF <- V.freeze c
+  yF <- V.freeze y
+  lF <- VM.read l 0
+  pure $ PCA (toMatrix k bF) (toMatrix k cF) (toMatrix m yF) (realToFrac lF)
+  where
+    toMatrix d = reshape d . V.map realToFrac
diff --git a/coda-pca.cabal b/coda-pca.cabal
new file mode 100644
index 0000000..4fedff7
--- /dev/null
+++ b/coda-pca.cabal
@@ -0,0 +1,27 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.34.7.
+--
+-- see: https://github.com/sol/hpack
+
+name:           coda-pca
+version:        0.1
+synopsis:       PCA on CODA
+maintainer:     Justin Bedo <cu@cua0.org>
+license:        MIT
+license-file:   LICENSE
+build-type:     Simple
+
+library
+  exposed-modules:
+      Numeric.CODA.PCA
+  other-modules:
+      Paths_coda_pca
+  c-sources:
+      pca.c
+  build-depends:
+      base
+    , hmatrix
+    , inline-c
+    , vector
+  default-language: Haskell2010
diff --git a/flake.nix b/flake.nix
index bd1f1ca..6db06c1 100644
--- a/flake.nix
+++ b/flake.nix
@@ -12,7 +12,15 @@
     in rec {
       packages = rec {
         lib = pkgs.callPackage ./. {};
-        R = pkgs.callPackage ./R { inherit (pkgs) R; inherit (pkgs.rPackages) buildRPackage; coda-pca = lib; };
+        R = pkgs.callPackage ./R {
+          inherit (pkgs) R;
+          inherit (pkgs.rPackages) buildRPackage;
+          coda-pca = lib;
+        };
+        haskell = pkgs.callPackage ./haskell.nix {
+          inherit (pkgs.haskellPackages) mkDerivation base inline-c vector hmatrix;
+          inherit (pkgs) lib;
+        };
       };
       defaultPackage = packages.lib;
     });
diff --git a/haskell.nix b/haskell.nix
new file mode 100644
index 0000000..56f10d4
--- /dev/null
+++ b/haskell.nix
@@ -0,0 +1,16 @@
+{
+  mkDerivation,
+  base,
+  inline-c,
+  vector,
+  hmatrix,
+  lib,
+}:
+mkDerivation {
+  pname = "coda-pca";
+  version = "0.1";
+  src = ./.;
+  libraryHaskellDepends = [base inline-c vector hmatrix];
+  description = "PCA on CODA";
+  license = lib.licenses.mit;
+}
diff --git a/package.yaml b/package.yaml
new file mode 100644
index 0000000..dd4fac9
--- /dev/null
+++ b/package.yaml
@@ -0,0 +1,17 @@
+name: coda-pca
+version: 0.1
+synopsis: PCA on CODA
+maintainer: Justin Bedo <cu@cua0.org>
+license: MIT
+
+dependencies:
+  - base
+  - inline-c
+  - vector
+  - hmatrix
+
+library:
+  exposed-modules:
+    - Numeric.CODA.PCA
+  c-sources:
+    - pca.c