diff options
| -rw-r--r-- | Numeric/CODA/PCA.hs | 92 | ||||
| -rw-r--r-- | coda-pca.cabal | 27 | ||||
| -rw-r--r-- | flake.nix | 10 | ||||
| -rw-r--r-- | haskell.nix | 16 | ||||
| -rw-r--r-- | package.yaml | 17 |
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 @@ -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 |