{-# LANGUAGE NamedFieldPuns #-}

-- | Principal component analysis as optic residual ownership over harpie arrays.
--
-- Pipeline (see circuits @examples\/pca.md@):
--
-- @
-- center → dagger-compose (Gram / SVD of centered data) → spectral cut
--       → Optic (,) with major = focus, minor = owned residual
-- @
--
-- 'Dagger' from circuits promises reverse wires; here reverse is matrix
-- transpose and the spectral cut is hmatrix. The keep\/discard protocol is
-- 'Circuit.PCA.Optic.Lens' — hinge table cell A (own + finite residual).
--
-- 'circuits-ad' is a sibling dependency for reverse-mode work on surrounding
-- nets; this module does not yet differentiate through the spectral cut.
module Circuit.PCA
  ( -- * Model
    PCAModel (..),
    fit,
    fitSVD,

    -- * Project / reconstruct (optic face)
    scores,
    reconstruct,
    projectRows,

    -- * Single-row optic
    rowLens,
    viewMajor,
    setMajor,

    -- * Re-exports
    module Circuit.PCA.Optic,
    module Circuit.PCA.Lin,
  )
where

-- circuits-ad: reverse-mode sibling; kept linked for diff-through-PCA work.

import Circuit.Diff.Circuit ()
import Circuit.PCA.Lin
import Circuit.PCA.Optic
import Data.Vector.Storable qualified as VS
import Harpie.Array.Storable (Array)
import Harpie.Array.Storable qualified as A
import Numeric.LinearAlgebra qualified as LA
import Prelude

-- | Fitted PCA: column means, loadings (features × k), singular values.
data PCAModel = PCAModel
  { -- | Column means, shape @[p]@.
    PCAModel -> Array Double
pcaMean :: !(Array Double),
    -- | Principal axes as columns, shape @[p, k]@.
    PCAModel -> Array Double
pcaComponents :: !(Array Double),
    -- | Top-k singular values of the centered data matrix (length k).
    PCAModel -> Vector Double
pcaSingularValues :: !(VS.Vector Double),
    -- | Requested component count.
    PCAModel -> Int
pcaK :: !Int
  }
  deriving stock (Int -> PCAModel -> ShowS
[PCAModel] -> ShowS
PCAModel -> String
(Int -> PCAModel -> ShowS)
-> (PCAModel -> String) -> ([PCAModel] -> ShowS) -> Show PCAModel
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PCAModel -> ShowS
showsPrec :: Int -> PCAModel -> ShowS
$cshow :: PCAModel -> String
show :: PCAModel -> String
$cshowList :: [PCAModel] -> ShowS
showList :: [PCAModel] -> ShowS
Show)

-- | Fit PCA with @k@ components via thin SVD of centered data.
--
-- Data shape: samples × features. Uses right singular vectors as loadings
-- (same subspace as eigendecomposition of the feature Gram).
fit :: Int -> Array Double -> PCAModel
fit :: Int -> Array Double -> PCAModel
fit = Int -> Array Double -> PCAModel
fitSVD

-- | Explicit SVD path (BLAS via hmatrix).
fitSVD :: Int -> Array Double -> PCAModel
fitSVD :: Int -> Array Double -> PCAModel
fitSVD Int
k Array Double
x
  | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 = String -> PCAModel
forall a. HasCallStack => String -> a
error String
"Circuit.PCA.fitSVD: k must be >= 1"
  | Bool
otherwise =
      case Array Double -> Maybe (Matrix Double)
toMatrix Array Double
x of
        Maybe (Matrix Double)
Nothing -> String -> PCAModel
forall a. HasCallStack => String -> a
error String
"Circuit.PCA.fitSVD: expected rank-2 samples×features"
        Just Matrix Double
m ->
          let n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
LA.rows Matrix Double
m
              p :: Int
p = Matrix Double -> Int
forall t. Matrix t -> Int
LA.cols Matrix Double
m
              k' :: Int
k' = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
k (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
n Int
p)
              (Array Double
xc, Array Double
mu) = Array Double -> (Array Double, Array Double)
centerColumns Array Double
x
           in case Array Double -> Maybe (Matrix Double)
toMatrix Array Double
xc of
                Maybe (Matrix Double)
Nothing -> String -> PCAModel
forall a. HasCallStack => String -> a
error String
"Circuit.PCA.fitSVD: centered data not rank-2"
                Just Matrix Double
mc ->
                  -- thin SVD: mc == u <> diag s <> tr v; columns of v = right SVs
                  let (Matrix Double
_u, Vector Double
s, Matrix Double
v) = Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
forall t.
Field t =>
Matrix t -> (Matrix t, Vector Double, Matrix t)
LA.thinSVD Matrix Double
mc
                      vk :: Matrix Double
vk = Matrix Double
v Matrix Double -> (Extractor, Extractor) -> Matrix Double
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
LA.?? (Extractor
LA.All, Int -> Extractor
LA.Take Int
k')
                      sTop :: Vector Double
sTop = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
VS.fromList (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
k' (Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
LA.toList Vector Double
s))
                   in PCAModel
                        { pcaMean :: Array Double
pcaMean = Array Double
mu,
                          pcaComponents :: Array Double
pcaComponents = Matrix Double -> Array Double
fromMatrix Matrix Double
vk,
                          pcaSingularValues :: Vector Double
pcaSingularValues = Vector Double
sTop,
                          pcaK :: Int
pcaK = Int
k'
                        }

meanRow :: Array Double -> LA.Vector Double
meanRow :: Array Double -> Vector Double
meanRow Array Double
v = Array Double -> Vector Double
forall t a. FromVector t a => t -> Vector a
A.asVector Array Double
v

-- | Project centered rows to k-dimensional scores. Shape: samples × k.
scores :: PCAModel -> Array Double -> Array Double
scores :: PCAModel -> Array Double -> Array Double
scores PCAModel {Array Double
pcaMean :: PCAModel -> Array Double
pcaMean :: Array Double
pcaMean, Array Double
pcaComponents :: PCAModel -> Array Double
pcaComponents :: Array Double
pcaComponents} Array Double
x =
  case (Array Double -> Maybe (Matrix Double)
toMatrix Array Double
x, Array Double -> Maybe (Matrix Double)
toMatrix Array Double
pcaComponents) of
    (Just Matrix Double
m, Just Matrix Double
vk) ->
      let mu :: Vector Double
mu = Array Double -> Vector Double
meanRow Array Double
pcaMean
          xc :: Matrix Double
xc = Matrix Double
m Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
- Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
LA.asRow Vector Double
mu
       in Matrix Double -> Array Double
fromMatrix (Matrix Double
xc Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
LA.<> Matrix Double
vk)
    (Maybe (Matrix Double), Maybe (Matrix Double))
_ -> String -> Array Double
forall a. HasCallStack => String -> a
error String
"Circuit.PCA.scores: rank-2 expected"

-- | Reconstruct data from scores: @scores × components† + mean@.
reconstruct :: PCAModel -> Array Double -> Array Double
reconstruct :: PCAModel -> Array Double -> Array Double
reconstruct PCAModel {Array Double
pcaMean :: PCAModel -> Array Double
pcaMean :: Array Double
pcaMean, Array Double
pcaComponents :: PCAModel -> Array Double
pcaComponents :: Array Double
pcaComponents} Array Double
sc =
  case (Array Double -> Maybe (Matrix Double)
toMatrix Array Double
sc, Array Double -> Maybe (Matrix Double)
toMatrix Array Double
pcaComponents) of
    (Just Matrix Double
s, Just Matrix Double
vk) ->
      let mu :: Vector Double
mu = Array Double -> Vector Double
meanRow Array Double
pcaMean
          xh :: Matrix Double
xh = Matrix Double
s Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
LA.<> Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix Double
vk
       in Matrix Double -> Array Double
fromMatrix (Matrix Double
xh Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
LA.asRow Vector Double
mu)
    (Maybe (Matrix Double), Maybe (Matrix Double))
_ -> String -> Array Double
forall a. HasCallStack => String -> a
error String
"Circuit.PCA.reconstruct: rank-2 expected"

-- | Full project-through-model: @reconstruct . scores@.
projectRows :: PCAModel -> Array Double -> Array Double
projectRows :: PCAModel -> Array Double -> Array Double
projectRows PCAModel
model Array Double
x = PCAModel -> Array Double -> Array Double
reconstruct PCAModel
model (PCAModel -> Array Double -> Array Double
scores PCAModel
model Array Double
x)

-- | Product optic on one ambient row.
--
-- Focus = major scores (@k@); residual = reconstruction error in ambient
-- space (owned). @put@ replaces scores and re-adds the previous residual.
rowLens ::
  PCAModel ->
  Lens (VS.Vector Double) (VS.Vector Double) (VS.Vector Double) (VS.Vector Double)
rowLens :: PCAModel
-> Lens
     (Vector Double) (Vector Double) (Vector Double) (Vector Double)
rowLens PCAModel {Array Double
pcaComponents :: PCAModel -> Array Double
pcaComponents :: Array Double
pcaComponents, Array Double
pcaMean :: PCAModel -> Array Double
pcaMean :: Array Double
pcaMean} =
  (Vector Double -> (Vector Double, Vector Double -> Vector Double))
-> Lens
     (Vector Double) (Vector Double) (Vector Double) (Vector Double)
forall s a b t. (s -> (a, b -> t)) -> Lens s t a b
fromClassical ((Vector Double -> (Vector Double, Vector Double -> Vector Double))
 -> Lens
      (Vector Double) (Vector Double) (Vector Double) (Vector Double))
-> (Vector Double
    -> (Vector Double, Vector Double -> Vector Double))
-> Lens
     (Vector Double) (Vector Double) (Vector Double) (Vector Double)
forall a b. (a -> b) -> a -> b
$ \Vector Double
row ->
    let meanV :: Vector Double
meanV = Array Double -> Vector Double
forall t a. FromVector t a => t -> Vector a
A.asVector Array Double
pcaMean
        xc :: Vector Double
xc = (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VS.zipWith (-) Vector Double
row Vector Double
meanV
     in case Array Double -> Maybe (Matrix Double)
toMatrix Array Double
pcaComponents of
          Maybe (Matrix Double)
Nothing -> String -> (Vector Double, Vector Double -> Vector Double)
forall a. HasCallStack => String -> a
error String
"Circuit.PCA.rowLens: components not rank-2"
          Just Matrix Double
vk ->
            let sc :: Vector Double
sc = Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
LA.flatten (Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
LA.asRow Vector Double
xc Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
LA.<> Matrix Double
vk)
                recon :: Vector Double
recon = (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VS.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Vector Double
meanV (Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
LA.flatten (Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
LA.asRow Vector Double
sc Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
LA.<> Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix Double
vk))
                err :: Vector Double
err = (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VS.zipWith (-) Vector Double
row Vector Double
recon
             in ( Vector Double
sc,
                  \Vector Double
sc' ->
                    let recon' :: Vector Double
recon' = (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VS.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Vector Double
meanV (Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
LA.flatten (Vector Double -> Matrix Double
forall a. Storable a => Vector a -> Matrix a
LA.asRow Vector Double
sc' Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
LA.<> Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
LA.tr Matrix Double
vk))
                     in (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VS.zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Vector Double
recon' Vector Double
err
                )

-- | View major scores of one row.
viewMajor :: PCAModel -> VS.Vector Double -> VS.Vector Double
viewMajor :: PCAModel -> Vector Double -> Vector Double
viewMajor PCAModel
model Vector Double
row = Lens
  (Vector Double) (Vector Double) (Vector Double) (Vector Double)
-> Vector Double -> Vector Double
forall s a. Lens s s a a -> s -> a
view (PCAModel
-> Lens
     (Vector Double) (Vector Double) (Vector Double) (Vector Double)
rowLens PCAModel
model) Vector Double
row

-- | Set major scores of one row (residual error held from @row@).
setMajor :: PCAModel -> VS.Vector Double -> VS.Vector Double -> VS.Vector Double
setMajor :: PCAModel -> Vector Double -> Vector Double -> Vector Double
setMajor PCAModel
model Vector Double
sc Vector Double
row = Lens
  (Vector Double) (Vector Double) (Vector Double) (Vector Double)
-> Vector Double -> Vector Double -> Vector Double
forall s t a b. Lens s t a b -> b -> s -> t
set (PCAModel
-> Lens
     (Vector Double) (Vector Double) (Vector Double) (Vector Double)
rowLens PCAModel
model) Vector Double
sc Vector Double
row