-- |
-- Module      : Math.LinearMap.Category
-- Copyright   : (c) Justus Sagemüller 2016
-- License     : GPL v3
-- 
-- Maintainer  : (@) jsag $ hvl.no
-- Stability   : experimental
-- Portability : portable
-- 


{-# LANGUAGE CPP                  #-}
{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE StandaloneDeriving   #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE UnicodeSyntax        #-}
{-# LANGUAGE TupleSections        #-}
{-# LANGUAGE ConstraintKinds      #-}
{-# LANGUAGE ExplicitNamespaces   #-}

module Math.LinearMap.Category (
            -- * Linear maps
            -- $linmapIntro

            -- ** Function implementation
              LinearFunction (..), type (-+>)(), Bilinear
            , lfun, (-+$>)
            -- ** Tensor implementation
            , LinearMap (..), type (+>)()
            , (⊕), (>+<)
            , adjoint
            -- ** Dual vectors
            -- $dualVectorIntro
            , (<.>^), (-+|>)
            -- * Tensor spaces
            , Tensor (..), type (⊗)(), (⊗)
            -- ** Symmetric
            , SymmetricTensor(..), squareV, squareVs
            , type (⊗〃+>)(), currySymBilin
            -- * Norms
            -- $metricIntro
            , Norm(..), Seminorm
            , spanNorm
            , euclideanNorm
            , (|$|)
            , normSq
            , (<$|)
            , scaleNorm
            , normSpanningSystem
            , normSpanningSystem'
            -- ** Variances
            , Variance, spanVariance, (|&>), varianceSpanningSystem
            , dualNorm, dualNorm', dependence
            -- ** Utility
            , densifyNorm, wellDefinedNorm
            -- * Solving linear equations
            , (\$), pseudoInverse, roughDet
            , linearRegressionW, linearRegression
            , LinearRegressionResult
            , linearFit_χν², linearFit_bestModel, linearFit_modelUncertainty 
            -- * Eigenvalue problems
            , eigen
            , constructEigenSystem
            , roughEigenSystem
            , finishEigenSystem
            , Eigenvector(..)
            -- * The classes of suitable vector spaces
            -- $classesExplanation
            -- ** Deriving from basis
            , module Math.LinearMap.Category.Instances.Deriving
            -- ** The classes
            , module Data.VectorSpace
            , LSpace
            , TensorSpace (..)
            , LinearSpace (..)
            -- ** Orthonormal systems
            , SemiInner (..), cartesianDualBasisCandidates, embedFreeSubspace
            -- ** Finite baseis
            , FiniteDimensional (..)
            -- * Utility
            -- ** Linear primitives
            , addV, scale, inner, flipBilin, bilinearFunction
            -- ** Tensors with basis decomposition
            , (.⊗)
            -- ** Hilbert space operations
            , (·), DualSpace, riesz, coRiesz, showsPrecAsRiesz, (.<)
            -- ** Standard decompositions
            , TensorDecomposable(..), RieszDecomposable(..)
            , tensorDecomposeShowsPrec, rieszDecomposeShowsPrec
            -- ** Constraint synonyms
            , HilbertSpace, SimpleSpace, RealSpace
            , Num'(..)
            , Fractional'
            , RealFrac', RealFloat', LinearShowable
            -- ** Double-dual, scalar-scalar etc. identity
            , ClosedScalarWitness(..), TrivialTensorWitness(..)
            , ScalarSpaceWitness(..), DualSpaceWitness(..), LinearManifoldWitness(..)
            , DualFinitenessWitness(..)
            -- ** Misc
            , relaxNorm, transformNorm, transformVariance
            , findNormalLength, normalLength
            , summandSpaceNorms, sumSubspaceNorms
            , sharedNormSpanningSystem, sharedSeminormSpanningSystem
            , sharedSeminormSpanningSystem'
            , convexPolytopeHull
            , symmetricPolytopeOuterVertices
            ) where

import Math.LinearMap.Category.Class
import Math.LinearMap.Category.Instances
import Math.LinearMap.Category.Instances.Deriving
import Math.LinearMap.Asserted
import Math.VectorSpace.Docile
import Math.LinearMap.Category.TensorQuot

import Data.Tree (Tree(..), Forest)
import Data.List (sortBy, foldl')
import qualified Data.Set as Set
import Data.Set (Set)
import Data.Ord (comparing)
import Data.List (maximumBy)
import Data.Maybe (catMaybes)
import Data.Foldable (toList)
import Data.Semigroup

import Data.VectorSpace
import Data.Basis

import Prelude ()
import qualified Prelude as Hask

import Control.Category.Constrained.Prelude hiding ((^))
import Control.Arrow.Constrained

import Linear ( V0(V0), V1(V1), V2(V2), V3(V3), V4(V4)
              , _x, _y, _z, _w )
import Data.VectorSpace.Free
import Math.VectorSpace.ZeroDimensional
import qualified Linear.Matrix as Mat
import qualified Linear.Vector as Mat
import Control.Lens ((^.))

import qualified Data.Vector.Unboxed as UArr

import Numeric.IEEE

import qualified GHC.Exts as GHC
import qualified Data.Type.Coercion as GHC

-- $linmapIntro
-- This library deals with linear functions, i.e. functions @f :: v -> w@
-- that fulfill
-- 
-- @
-- f $ μ 'Data.VectorSpace.^*' u 'Data.AdditiveGroup.^+^' v ≡ μ ^* f u ^+^ f v    ∀ u,v :: v;  μ :: 'Scalar' v
-- @
-- 
-- Such functions form a cartesian monoidal category (in maths called 
-- <https://en.wikipedia.org/wiki/Category_of_modules#Example:_the_category_of_vector_spaces VectK>).
-- This is implemented by 'Control.Arrow.Constrained.PreArrow', which is the
-- preferred interface for dealing with these mappings. The basic
-- “matrix operations” are then:
-- 
-- * Identity matrix: 'Control.Category.Constrained.id'
-- * Matrix addition: 'Data.AdditiveGroup.^+^' (linear maps form an ordinary vector space)
-- * Matrix-matrix multiplication: 'Control.Category.Constrained.<<<'
--     (or '>>>' or 'Control.Category.Constrained..')
-- * Matrix-vector multiplication: 'Control.Arrow.Constrained.$'
-- * Vertical matrix concatenation: 'Control.Arrow.Constrained.&&&'
-- * Horizontal matrix concatenation: '⊕' (aka '>+<')
-- 
-- But linear mappings need not necessarily be implemented as matrices:


-- $dualVectorIntro
-- A @'DualVector' v@ is a linear functional or
-- <https://en.wikipedia.org/wiki/Linear_form linear form> on the vector space @v@,
-- i.e. it is a linear function from the vector space into its scalar field.
-- However, these functions form themselves a vector space, known as the dual space.
-- In particular, the dual space of any 'InnerSpace' is isomorphic to the
-- space itself.
-- 
-- (More precisely: the continuous dual space of a
-- <https://en.wikipedia.org/wiki/Hilbert_space Hilbert space> is isomorphic to
-- that Hilbert space itself; see the 'riesz' isomorphism.)
-- 
-- As a matter of fact, in many applications, no distinction is made between a
-- space and its dual. Indeed, we have for the basic 'LinearSpace' instances
-- @'DualVector' v ~ v@, and '<.>^' is simply defined as a scalar product.
-- In this case, a general 'LinearMap' is just a tensor product / matrix.
-- 
-- However, scalar products are often not as natural as they are made to look:
-- 
-- * A scalar product is only preserved under orthogonal transformations.
--   It is not preserved under scalings, and certainly not under general linear
--   transformations. This is very important in applications such as relativity
--   theory (here, people talk about /covariant/ vs /contravariant/ tensors),
--   but also relevant for more mundane
--   <http://hackage.haskell.org/package/manifolds manifolds> like /sphere surfaces/:
--   on such a surface, the natural symmetry transformations do generally
--   not preserve a scalar product you might define.
-- 
-- * There may be more than one meaningful scalar product. For instance,
--   the <https://en.wikipedia.org/wiki/Sobolev_space Sobolev space> of weakly
--   differentiable functions also permits the
--   <https://en.wikipedia.org/wiki/Square-integrable_function 𝐿²> scalar product
--   – each has different and useful properties.
-- 
-- Neither of this is a problem if we keep the dual space a separate type.
-- Effectively, this enables the type system to prevent you from writing code that
-- does not behave natural (i.e. that depends on a concrete choice of basis / scalar
-- product).
-- 
-- For cases when you do have some given notion of orientation/scale in a vector space
-- and need it for an algorithm, you can always provide a 'Norm', which is essentially
-- a reified scalar product.
-- 
-- Note that @DualVector (DualVector v) ~ v@ in any 'LSpace': the /double-dual/
-- space is /naturally/ isomorphic to the original space, by way of
-- 
-- @
-- v '<.>^' dv  ≡  dv '<.>^' v
-- @



-- | A linear map that simply projects from a dual vector in @u@ to a vector in @v@.
-- 
-- @
-- (du '-+|>' v) u  ≡  v '^*' (du '<.>^' u)
-- @
infixr 7 -+|>
(-+|>) :: ( EnhancedCat f (LinearFunction s)
          , LSpace u, LSpace v, Scalar u ~ s, Scalar v ~ s
          , Object f u, Object f v )
             => DualVector u -> v -> f u v
DualVector u
du-+|> :: DualVector u -> v -> f u v
-+|>v
v = LinearFunction s u v -> f u v
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr (LinearFunction s u v -> f u v)
-> ((u -> v) -> LinearFunction s u v) -> (u -> v) -> f u v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (u -> v) -> LinearFunction s u v
forall s v w. (v -> w) -> LinearFunction s v w
LinearFunction ((u -> v) -> f u v) -> (u -> v) -> f u v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ (v
vv -> s -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*) (s -> v) -> (u -> s) -> u -> v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (DualVector u
duDualVector u -> u -> Scalar u
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^)




-- $metricIntro
-- A norm is a way to quantify the magnitude/length of different vectors,
-- even if they point in different directions.
-- 
-- In an 'InnerSpace', a norm is always given by the scalar product,
-- but there are spaces without a canonical scalar product (or situations
-- in which this scalar product does not give the metric you want). Hence,
-- we let the functions like 'constructEigenSystem', which depend on a norm
-- for orthonormalisation, accept a 'Norm' as an extra argument instead of
-- requiring 'InnerSpace'.

-- | A seminorm defined by
-- 
-- @
-- ‖v‖ = √(∑ᵢ ⟨dᵢ|v⟩²)
-- @
-- 
-- for some dual vectors @dᵢ@. If given a complete basis of the dual space,
-- this generates a proper 'Norm'.
-- 
-- If the @dᵢ@ are a complete orthonormal system, you get the 'euclideanNorm'
-- (in an inefficient form).
spanNorm :: ∀ v . LSpace v => [DualVector v] -> Seminorm v
spanNorm :: [DualVector v] -> Seminorm v
spanNorm = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \[DualVector v]
dvs -> LinearFunction (Scalar (DualVector v)) v (DualVector v)
-> Seminorm v
forall v. (v -+> DualVector v) -> Norm v
Norm (LinearFunction (Scalar (DualVector v)) v (DualVector v)
 -> Seminorm v)
-> ((v -> DualVector v)
    -> LinearFunction (Scalar (DualVector v)) v (DualVector v))
-> (v -> DualVector v)
-> Seminorm v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (v -> DualVector v)
-> LinearFunction (Scalar (DualVector v)) v (DualVector v)
forall s v w. (v -> w) -> LinearFunction s v w
LinearFunction ((v -> DualVector v) -> Seminorm v)
-> (v -> DualVector v) -> Seminorm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ \v
v -> [DualVector v] -> DualVector v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [DualVector v
dv DualVector v -> Scalar v -> DualVector v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^* (DualVector v
dvDualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v) | DualVector v
dv <- [DualVector v]
dvs]

spanVariance :: ∀ v . LSpace v => [v] -> Variance v
spanVariance :: [v] -> Variance v
spanVariance = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness -> [v] -> Variance v
forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm

-- | Modify a norm in such a way that the given vectors lie within its unit ball.
--   (Not /optimally/ – the unit ball may be bigger than necessary.)
relaxNorm :: ∀ v . SimpleSpace v => Norm v -> [v] -> Norm v
relaxNorm :: Norm v -> [v] -> Norm v
relaxNorm = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \Norm v
me [v]
vs -> let vs' :: [v]
vs' = Norm v -> [v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem' Norm v
me
                     in Norm (DualVector v) -> Norm v
forall v. SimpleSpace v => Norm v -> Variance v
dualNorm (Norm (DualVector v) -> Norm v)
-> ([v] -> Norm (DualVector v)) -> [v] -> Norm v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [v] -> Norm (DualVector v)
forall v. LSpace v => [v] -> Variance v
spanVariance ([v] -> Norm v) -> [v] -> Norm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [v]
vs' [v] -> [v] -> [v]
forall a. [a] -> [a] -> [a]
++ [v]
vs

-- | Scale the result of a norm with the absolute of the given number.
-- 
-- @
-- scaleNorm μ n |$| v = abs μ * (n|$|v)
-- @
-- 
-- Equivalently, this scales the norm's unit ball by the reciprocal of that factor.
scaleNorm :: ∀ v . LSpace v => Scalar v -> Norm v -> Norm v
scaleNorm :: Scalar v -> Norm v -> Norm v
scaleNorm = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness -> \Scalar v
μ (Norm v -+> DualVector v
n) -> (v -+> DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm ((v -+> DualVector v) -> Norm v) -> (v -+> DualVector v) -> Norm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Scalar v
μScalar v -> Int -> Scalar v
forall a. Num a => a -> Int -> a
^Int
2 Scalar (v -+> DualVector v)
-> (v -+> DualVector v) -> v -+> DualVector v
forall v. VectorSpace v => Scalar v -> v -> v
*^ v -+> DualVector v
n

-- | A positive (semi)definite symmetric bilinear form. This gives rise
--   to a <https://en.wikipedia.org/wiki/Norm_(mathematics) norm> thus:
-- 
--   @
--   'Norm' n '|$|' v = √(n v '<.>^' v)
--   @
--   
--   Strictly speaking, this type is neither strong enough nor general enough to
--   deserve the name 'Norm': it includes proper 'Seminorm's (i.e. @m|$|v ≡ 0@ does
--   not guarantee @v == zeroV@), but not actual norms such as the ℓ₁-norm on ℝⁿ
--   (Taxcab norm) or the supremum norm.
--   However, 𝐿₂-like norms are the only ones that can really be formulated without
--   any basis reference; and guaranteeing positive definiteness through the type
--   system is scarcely practical.
newtype Norm v = Norm {
    Norm v -> v -+> DualVector v
applyNorm :: v -+> DualVector v
  }

-- | A “norm” that may explicitly be degenerate, with @m|$|v ⩵ 0@ for some @v ≠ zeroV@.
type Seminorm v = Norm v

-- | @(m<>n|$|v)^2 ⩵ (m|$|v)^2 + (n|$|v)^2@
instance LSpace v => Semigroup (Norm v) where
  Norm v -+> DualVector v
m <> :: Norm v -> Norm v -> Norm v
<> Norm v -+> DualVector v
n = (v -+> DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm ((v -+> DualVector v) -> Norm v) -> (v -+> DualVector v) -> Norm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v -+> DualVector v
m(v -+> DualVector v) -> (v -+> DualVector v) -> v -+> DualVector v
forall v. AdditiveGroup v => v -> v -> v
^+^v -+> DualVector v
n
-- | @mempty|$|v ≡ 0@
instance LSpace v => Monoid (Seminorm v) where
  mempty :: Seminorm v
mempty = (v -+> DualVector v) -> Seminorm v
forall v. (v -+> DualVector v) -> Norm v
Norm v -+> DualVector v
forall v. AdditiveGroup v => v
zeroV
  mappend :: Seminorm v -> Seminorm v -> Seminorm v
mappend = Seminorm v -> Seminorm v -> Seminorm v
forall a. Semigroup a => a -> a -> a
(<>)

-- | A multidimensional variance of points @v@ with some distribution can be
--   considered a norm on the dual space, quantifying for a dual vector @dv@ the
--   expectation value of @(dv<.>^v)^2@.
type Variance v = Norm (DualVector v)

-- | The canonical standard norm (2-norm) on inner-product / Hilbert spaces.
euclideanNorm :: HilbertSpace v => Norm v
euclideanNorm :: Norm v
euclideanNorm = (v -+> DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm v -+> DualVector v
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id

-- | The norm induced from the (arbitrary) choice of basis in a finite space.
--   Only use this in contexts where you merely need /some/ norm, but don't
--   care if it might be biased in some unnatural way.
adhocNorm :: FiniteDimensional v => Norm v
adhocNorm :: Norm v
adhocNorm = (v -+> DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm v -+> DualVector v
forall v. FiniteDimensional v => v -+> DualVector v
uncanonicallyToDual

-- | A proper norm induces a norm on the dual space – the “reciprocal norm”.
--   (The orthonormal systems of the norm and its dual are mutually conjugate.)
--   The dual norm of a seminorm is undefined.
dualNorm :: SimpleSpace v => Norm v -> Variance v
dualNorm :: Norm v -> Variance v
dualNorm = [v] -> Variance v
forall v. LSpace v => [v] -> Variance v
spanVariance ([v] -> Variance v) -> (Norm v -> [v]) -> Norm v -> Variance v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Norm v -> [v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem'

-- | 'dualNorm' in the opposite direction. This is actually self-inverse;
--    with 'dualSpaceWitness' you can replace each with the other direction.
dualNorm' :: ∀ v . SimpleSpace v => Variance v -> Norm v
dualNorm' :: Variance v -> Norm v
dualNorm' = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
     DualSpaceWitness v
DualSpaceWitness -> [DualVector v] -> Norm v
forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm ([DualVector v] -> Norm v)
-> (Variance v -> [DualVector v]) -> Variance v -> Norm v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Variance v -> [DualVector v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem'

transformNorm :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w)
                             => (v+>w) -> Norm w -> Norm v
transformNorm :: (v +> w) -> Norm w -> Norm v
transformNorm = case ( DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v
                     , DualSpaceWitness w
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness w ) of
    (DualSpaceWitness v
DualSpaceWitness, DualSpaceWitness w
DualSpaceWitness)
            -> \v +> w
f (Norm w -+> DualVector w
m) -> LinearFunction (Scalar (DualVector v)) v (DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm (LinearFunction (Scalar (DualVector v)) v (DualVector v) -> Norm v)
-> (LinearMap (Scalar (DualVector w)) v (DualVector v)
    -> LinearFunction (Scalar (DualVector v)) v (DualVector v))
-> LinearMap (Scalar (DualVector w)) v (DualVector v)
-> Norm v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearMap (Scalar (DualVector w)) v (DualVector v)
-> LinearFunction (Scalar (DualVector v)) v (DualVector v)
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr (LinearMap (Scalar (DualVector w)) v (DualVector v) -> Norm v)
-> LinearMap (Scalar (DualVector w)) v (DualVector v) -> Norm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ (LinearFunction
  (Scalar (DualVector w))
  (LinearMap (Scalar v) v (DualVector (DualVector w)))
  (LinearMap (Scalar (DualVector w)) (DualVector w) (DualVector v))
forall v w.
(LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) =>
(v +> DualVector w) -+> (w +> DualVector v)
adjoint LinearFunction
  (Scalar (DualVector w))
  (LinearMap (Scalar v) v (DualVector (DualVector w)))
  (LinearMap (Scalar (DualVector w)) (DualVector w) (DualVector v))
-> LinearMap (Scalar v) v (DualVector (DualVector w))
-> LinearMap (Scalar (DualVector w)) (DualVector w) (DualVector v)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
LinearMap (Scalar v) v (DualVector (DualVector w))
f) LinearMap (Scalar (DualVector w)) (DualVector w) (DualVector v)
-> LinearMap (Scalar (DualVector w)) v (DualVector w)
-> LinearMap (Scalar (DualVector w)) v (DualVector v)
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (LinearFunction (Scalar w) w (DualVector w)
-> LinearFunction
     (Scalar w)
     (LinearMap (Scalar (DualVector w)) v w)
     (LinearMap (Scalar (DualVector w)) v (DualVector w))
forall (f :: * -> *) (r :: * -> * -> *) (t :: * -> * -> *) a b.
(Functor f r t, Object r a, Object t (f a), Object r b,
 Object t (f b)) =>
r a b -> t (f a) (f b)
fmap LinearFunction (Scalar w) w (DualVector w)
w -+> DualVector w
m LinearFunction
  (Scalar w)
  (LinearMap (Scalar (DualVector w)) v w)
  (LinearMap (Scalar (DualVector w)) v (DualVector w))
-> LinearMap (Scalar (DualVector w)) v w
-> LinearMap (Scalar (DualVector w)) v (DualVector w)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
LinearMap (Scalar (DualVector w)) v w
f)

transformVariance :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w)
                        => (v+>w) -> Variance v -> Variance w
transformVariance :: (v +> w) -> Variance v -> Variance w
transformVariance = case ( DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v
                     , DualSpaceWitness w
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness w ) of
    (DualSpaceWitness v
DualSpaceWitness, DualSpaceWitness w
DualSpaceWitness)
            -> \v +> w
f (Norm DualVector v -+> DualVector (DualVector v)
m) -> LinearFunction
  (Scalar w) (DualVector w) (DualVector (DualVector w))
-> Variance w
forall v. (v -+> DualVector v) -> Norm v
Norm (LinearFunction
   (Scalar w) (DualVector w) (DualVector (DualVector w))
 -> Variance w)
-> (LinearMap (Scalar w) (DualVector w) w
    -> LinearFunction
         (Scalar w) (DualVector w) (DualVector (DualVector w)))
-> LinearMap (Scalar w) (DualVector w) w
-> Variance w
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearMap (Scalar w) (DualVector w) w
-> LinearFunction
     (Scalar w) (DualVector w) (DualVector (DualVector w))
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr (LinearMap (Scalar w) (DualVector w) w -> Variance w)
-> LinearMap (Scalar w) (DualVector w) w -> Variance w
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
LinearMap (Scalar w) v w
f LinearMap (Scalar w) v w
-> LinearMap (Scalar w) (DualVector w) v
-> LinearMap (Scalar w) (DualVector w) w
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (LinearFunction (Scalar w) (DualVector v) v
-> LinearFunction
     (Scalar w)
     (LinearMap (Scalar w) (DualVector w) (DualVector v))
     (LinearMap (Scalar w) (DualVector w) v)
forall (f :: * -> *) (r :: * -> * -> *) (t :: * -> * -> *) a b.
(Functor f r t, Object r a, Object t (f a), Object r b,
 Object t (f b)) =>
r a b -> t (f a) (f b)
fmap LinearFunction (Scalar w) (DualVector v) v
DualVector v -+> DualVector (DualVector v)
m LinearFunction
  (Scalar w)
  (LinearMap (Scalar w) (DualVector w) (DualVector v))
  (LinearMap (Scalar w) (DualVector w) v)
-> LinearMap (Scalar w) (DualVector w) (DualVector v)
-> LinearMap (Scalar w) (DualVector w) v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearFunction
  (Scalar (DualVector w))
  (LinearMap (Scalar v) v (DualVector (DualVector w)))
  (LinearMap (Scalar w) (DualVector w) (DualVector v))
forall v w.
(LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) =>
(v +> DualVector w) -+> (w +> DualVector v)
adjoint LinearFunction
  (Scalar (DualVector w))
  (LinearMap (Scalar v) v (DualVector (DualVector w)))
  (LinearMap (Scalar w) (DualVector w) (DualVector v))
-> LinearMap (Scalar v) v (DualVector (DualVector w))
-> LinearMap (Scalar w) (DualVector w) (DualVector v)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
LinearMap (Scalar v) v (DualVector (DualVector w))
f)

infixl 6 ^%
(^%) :: (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
v
v ^% :: v -> Norm v -> v
^% Norm v -+> DualVector v
m = v
v v -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ Scalar v -> Scalar v
forall a. Floating a => a -> a
sqrt ((v -+> DualVector v
m(v -+> DualVector v) -> v -> DualVector v
forall s v w. LinearFunction s v w -> v -> w
-+$>v
v)DualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v)

-- | The unique positive number whose norm is 1 (if the norm is not constant zero).
findNormalLength :: ∀ s . RealFrac' s => Norm s -> Maybe s
findNormalLength :: Norm s -> Maybe s
findNormalLength (Norm s -+> DualVector s
m) = case ( ClosedScalarWitness s
forall s. Num' s => ClosedScalarWitness s
closedScalarWitness :: ClosedScalarWitness s
                                 , s -+> DualVector s
m(s -+> DualVector s) -> s -> DualVector s
forall s v w. LinearFunction s v w -> v -> w
-+$>s
1 ) of
   (ClosedScalarWitness s
ClosedScalarWitness, DualVector s
o) | DualVector s
o DualVector s -> DualVector s -> Bool
forall a. Ord a => a -> a -> Bool
> DualVector s
0      -> s -> Maybe s
forall a. a -> Maybe a
Just (s -> Maybe s) -> (DualVector s -> s) -> DualVector s -> Maybe s
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. DualVector s -> s
forall a. Floating a => a -> a
sqrt (DualVector s -> Maybe s) -> DualVector s -> Maybe s
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ DualVector s -> DualVector s
forall a. Fractional a => a -> a
recip DualVector s
o
                            | Bool
otherwise  -> Maybe s
forall a. Maybe a
Nothing

-- | Unsafe version of 'findNormalLength', only works reliable if the norm
--   is actually positive definite.
normalLength :: ∀ s . RealFrac' s => Norm s -> s
normalLength :: Norm s -> s
normalLength (Norm s -+> DualVector s
m) = case ( ClosedScalarWitness s
forall s. Num' s => ClosedScalarWitness s
closedScalarWitness :: ClosedScalarWitness s
                             , s -+> DualVector s
m(s -+> DualVector s) -> s -> DualVector s
forall s v w. LinearFunction s v w -> v -> w
-+$>s
1 ) of
   (ClosedScalarWitness s
ClosedScalarWitness, DualVector s
o) | DualVector s
o DualVector s -> DualVector s -> Bool
forall a. Ord a => a -> a -> Bool
>= DualVector s
0     -> DualVector s -> s
forall a. Floating a => a -> a
sqrt (DualVector s -> s) -> DualVector s -> s
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ DualVector s -> DualVector s
forall a. Fractional a => a -> a
recip DualVector s
o
                            | DualVector s
o DualVector s -> DualVector s -> Bool
forall a. Ord a => a -> a -> Bool
< DualVector s
0      -> [Char] -> s
forall a. HasCallStack => [Char] -> a
error [Char]
"Norm fails to be positive semidefinite."
                            | Bool
otherwise  -> [Char] -> s
forall a. HasCallStack => [Char] -> a
error [Char]
"Norm yields NaN."

infixr 0 <$|, |$|
-- | “Partially apply” a norm, yielding a dual vector
--   (i.e. a linear form that accepts the second argument of the scalar product).
-- 
-- @
-- ('euclideanNorm' '<$|' v) '<.>^' w  ≡  v '<.>' w
-- @
-- 
--   See also '|&>'.
(<$|) :: LSpace v => Norm v -> v -> DualVector v
Norm v -+> DualVector v
m <$| :: Norm v -> v -> DualVector v
<$| v
v = v -+> DualVector v
m(v -+> DualVector v) -> v -> DualVector v
forall s v w. LinearFunction s v w -> v -> w
-+$>v
v

-- | The squared norm. More efficient than '|$|' because that needs to take
--   the square root.
normSq :: LSpace v => Seminorm v -> v -> Scalar v
normSq :: Seminorm v -> v -> Scalar v
normSq (Norm v -+> DualVector v
m) v
v = (v -+> DualVector v
m(v -+> DualVector v) -> v -> DualVector v
forall s v w. LinearFunction s v w -> v -> w
-+$>v
v)DualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v

-- | Use a 'Norm' to measure the length / norm of a vector.
-- 
-- @
-- 'euclideanNorm' |$| v  ≡  √(v '<.>' v)
-- @
(|$|) :: (LSpace v, Floating (Scalar v)) => Seminorm v -> v -> Scalar v
|$| :: Seminorm v -> v -> Scalar v
(|$|) Seminorm v
m = Scalar v -> Scalar v
forall a. Floating a => a -> a
sqrt (Scalar v -> Scalar v) -> (v -> Scalar v) -> v -> Scalar v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Seminorm v -> v -> Scalar v
forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Seminorm v
m

infixl 1 |&>
-- | Flipped, “ket” version of '<$|'.
-- 
-- @
-- v '<.>^' (w |&> 'euclideanNorm')  ≡  v '<.>' w
-- @
(|&>) :: LSpace v => DualVector v -> Variance v -> v
DualVector v
dv |&> :: DualVector v -> Variance v -> v
|&> Norm DualVector v -+> DualVector (DualVector v)
m = Coercion v (DualVector (DualVector v))
-> Coercion (DualVector (DualVector v)) v
forall k (a :: k) (b :: k). Coercion a b -> Coercion b a
GHC.sym Coercion v (DualVector (DualVector v))
forall v. LinearSpace v => Coercion v (DualVector (DualVector v))
coerceDoubleDual Coercion (DualVector (DualVector v)) v
-> DualVector (DualVector v) -> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ DualVector v -+> DualVector (DualVector v)
m(DualVector v -+> DualVector (DualVector v))
-> DualVector v -> DualVector (DualVector v)
forall s v w. LinearFunction s v w -> v -> w
-+$>DualVector v
dv


-- | 'spanNorm' / 'spanVariance' are inefficient if the number of vectors
--   is similar to the dimension of the space, or even larger than it.
--   Use this function to optimise the underlying operator to a dense
--   matrix representation.
densifyNorm :: ∀ v . LSpace v => Norm v -> Norm v
densifyNorm :: Norm v -> Norm v
densifyNorm = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \(Norm v -+> DualVector v
m) -> (v -+> DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm ((v -+> DualVector v) -> Norm v)
-> (LinearMap (Scalar v) v (DualVector v) -> v -+> DualVector v)
-> LinearMap (Scalar v) v (DualVector v)
-> Norm v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearMap (Scalar v) v (DualVector v) -> v -+> DualVector v
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr (LinearMap (Scalar v) v (DualVector v) -> Norm v)
-> LinearMap (Scalar v) v (DualVector v) -> Norm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearFunction
  (Scalar v)
  (v -+> DualVector v)
  (LinearMap (Scalar v) v (DualVector v))
forall v w.
(LinearSpace v, TensorSpace w, Scalar v ~ Scalar w) =>
(v -+> w) -+> (v +> w)
sampleLinearFunction LinearFunction
  (Scalar v)
  (v -+> DualVector v)
  (LinearMap (Scalar v) v (DualVector v))
-> (v -+> DualVector v) -> LinearMap (Scalar v) v (DualVector v)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v -+> DualVector v
m

-- | Like 'densifyNorm', but also perform a “sanity check” to eliminate NaN etc. problems.
wellDefinedNorm :: ∀ v . LinearSpace v => Norm v -> Maybe (Norm v)
wellDefinedNorm :: Norm v -> Maybe (Norm v)
wellDefinedNorm = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \(Norm v -+> DualVector v
m) -> (v -+> DualVector v) -> Norm v
forall v. (v -+> DualVector v) -> Norm v
Norm ((v -+> DualVector v) -> Norm v)
-> Maybe (v -+> DualVector v) -> Maybe (Norm v)
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> (v -+> DualVector v) -> Maybe (v -+> DualVector v)
forall v. TensorSpace v => v -> Maybe v
wellDefinedVector v -+> DualVector v
m

data OrthonormalSystem v = OrthonormalSystem {
      OrthonormalSystem v -> Norm v
orthonormalityNorm :: Norm v
    , OrthonormalSystem v -> [v]
orthonormalVectors :: [v]
    }

orthonormaliseFussily :: ∀ v . (LSpace v, RealFloat (Scalar v))
                           => Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily :: Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily = DualSpaceWitness v -> Scalar v -> Norm v -> [v] -> [v]
onf DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness
 where onf :: DualSpaceWitness v -> Scalar v -> Norm v -> [v] -> [v]
       onf :: DualSpaceWitness v -> Scalar v -> Norm v -> [v] -> [v]
onf DualSpaceWitness v
DualSpaceWitness Scalar v
fuss Norm v
me = [(v, DualVector v)] -> [v] -> [v]
go []
        where go :: [(v, DualVector v)] -> [v] -> [v]
go [(v, DualVector v)]
_ [] = []
              go [(v, DualVector v)]
ws (v
v₀:[v]
vs)
                | Scalar v
mvd Scalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
> Scalar v
fuss  = let μ :: Scalar v
μ = Scalar v
1Scalar v -> Scalar v -> Scalar v
forall a. Fractional a => a -> a -> a
/Scalar v -> Scalar v
forall a. Floating a => a -> a
sqrt Scalar v
mvd
                                    v :: v
v = v
vdv -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ
                                in v
v v -> [v] -> [v]
forall a. a -> [a] -> [a]
: [(v, DualVector v)] -> [v] -> [v]
go ((v
v,DualVector v
dvdDualVector v -> Scalar v -> DualVector v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ)(v, DualVector v) -> [(v, DualVector v)] -> [(v, DualVector v)]
forall a. a -> [a] -> [a]
:[(v, DualVector v)]
ws) [v]
vs
                | Bool
otherwise   = [(v, DualVector v)] -> [v] -> [v]
go [(v, DualVector v)]
ws [v]
vs
               where vd :: v
vd = [(v, DualVector v)] -> LinearFunction (Scalar v) v v
forall v. LSpace v => [(v, DualVector v)] -> v -+> v
orthogonalComplementProj' [(v, DualVector v)]
ws LinearFunction (Scalar v) v v -> v -> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v
v₀
                     dvd :: DualVector v
dvd = Norm v -> v -+> DualVector v
forall v. Norm v -> v -+> DualVector v
applyNorm Norm v
me LinearFunction (Scalar v) v (DualVector v) -> v -> DualVector v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v
vd
                     mvd :: Scalar v
mvd = DualVector v
dvdDualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
vd

orthogonalComplementProj' :: LSpace v => [(v, DualVector v)] -> (v-+>v)
orthogonalComplementProj' :: [(v, DualVector v)] -> v -+> v
orthogonalComplementProj' [(v, DualVector v)]
ws = (v -> v) -> v -+> v
forall s v w. (v -> w) -> LinearFunction s v w
LinearFunction ((v -> v) -> v -+> v) -> (v -> v) -> v -+> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ \v
v₀
             -> (v -> (v, DualVector v) -> v) -> v -> [(v, DualVector v)] -> v
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (\v
v (v
w,DualVector v
dw) -> v
v v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
wv -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*(DualVector v
dwDualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v)) v
v₀ [(v, DualVector v)]
ws

orthogonalComplementProj :: LSpace v => Norm v -> [v] -> (v-+>v)
orthogonalComplementProj :: Norm v -> [v] -> v -+> v
orthogonalComplementProj (Norm v -+> DualVector v
m)
      = [(v, DualVector v)] -> v -+> v
forall v. LSpace v => [(v, DualVector v)] -> v -+> v
orthogonalComplementProj' ([(v, DualVector v)] -> v -+> v)
-> ([v] -> [(v, DualVector v)]) -> [v] -> v -+> v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (v -> (v, DualVector v)) -> [v] -> [(v, DualVector v)]
forall a b. (a -> b) -> [a] -> [b]
map (v -> v
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id (v -> v) -> (v -> DualVector v) -> v -> (v, DualVector v)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&& (v -+> DualVector v
m(v -+> DualVector v) -> v -> DualVector v
forall s v w. LinearFunction s v w -> v -> w
-+$>))


-- | A space in which you can use '·' both for scaling with a real number,
--   and as dot-product for obtaining such a number.
type RealSpace v = ( LinearSpace v, Scalar v ~ ℝ
                   , TensorQuot v ℝ, (v⨸ℝ) ~ DualVector v
                   , TensorQuot v v, (v⨸v) ~ ℝ )


data Eigenvector v = Eigenvector {
      Eigenvector v -> Scalar v
ev_Eigenvalue :: Scalar v -- ^ The estimated eigenvalue @λ@.
    , Eigenvector v -> v
ev_Eigenvector :: v       -- ^ Normalised vector @v@ that gets mapped to a multiple, namely:
    , Eigenvector v -> v
ev_FunctionApplied :: v   -- ^ @f $ v ≡ λ *^ v @.
    , Eigenvector v -> v
ev_Deviation :: v         -- ^ Deviation of @v@ to @(f$v)^/λ@. Ideally, this would of course be equal.
    , Eigenvector v -> Scalar v
ev_Badness :: Scalar v    -- ^ Squared norm of the deviation.
    }
deriving instance (Show v, Show (Scalar v)) => Show (Eigenvector v)

-- | Lazily compute the eigenbasis of a linear map. The algorithm is essentially
--   a hybrid of Lanczos/Arnoldi style Krylov-spanning and QR-diagonalisation,
--   which we don't do separately but /interleave/ at each step.
-- 
--   The size of the eigen-subbasis increases with each step until the space's
--   dimension is reached. (But the algorithm can also be used for
--   infinite-dimensional spaces.)
constructEigenSystem :: (LSpace v, RealFloat (Scalar v))
      => Norm v           -- ^ The notion of orthonormality.
      -> Scalar v           -- ^ Error bound for deviations from eigen-ness.
      -> (v-+>v)            -- ^ Operator to calculate the eigensystem of.
                            --   Must be Hermitian WRT the scalar product
                            --   defined by the given metric.
      -> [v]                -- ^ Starting vector(s) for the power method.
      -> [[Eigenvector v]]  -- ^ Infinite sequence of ever more accurate approximations
                            --   to the eigensystem of the operator.
constructEigenSystem :: Norm v -> Scalar v -> (v -+> v) -> [v] -> [[Eigenvector v]]
constructEigenSystem Norm v
me Scalar v
ε₀ v -+> v
f = ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v] -> [[Eigenvector v]]
forall a. (a -> a) -> a -> [a]
iterate (
                                             (Eigenvector v -> Eigenvector v -> Ordering)
-> [Eigenvector v] -> [Eigenvector v]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((Eigenvector v -> Scalar v)
-> Eigenvector v -> Eigenvector v -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Eigenvector v -> Scalar v)
 -> Eigenvector v -> Eigenvector v -> Ordering)
-> (Eigenvector v -> Scalar v)
-> Eigenvector v
-> Eigenvector v
-> Ordering
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$
                                               Scalar v -> Scalar v
forall a. Num a => a -> a
negate (Scalar v -> Scalar v)
-> (Eigenvector v -> Scalar v) -> Eigenvector v -> Scalar v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Scalar v -> Scalar v
forall a. Num a => a -> a
abs (Scalar v -> Scalar v)
-> (Eigenvector v -> Scalar v) -> Eigenvector v -> Scalar v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Eigenvector v -> Scalar v
forall v. Eigenvector v -> Scalar v
ev_Eigenvalue)
                                           ([Eigenvector v] -> [Eigenvector v])
-> ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v]
-> [Eigenvector v]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (v -> Eigenvector v) -> [v] -> [Eigenvector v]
forall a b. (a -> b) -> [a] -> [b]
map v -> Eigenvector v
asEV
                                           ([v] -> [Eigenvector v])
-> ([Eigenvector v] -> [v]) -> [Eigenvector v] -> [Eigenvector v]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Scalar v -> Norm v -> [v] -> [v]
forall v.
(LSpace v, RealFloat (Scalar v)) =>
Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily (Scalar v
1Scalar v -> Scalar v -> Scalar v
forall a. Fractional a => a -> a -> a
/Scalar v
4) Norm v
me
                                           ([v] -> [v]) -> ([Eigenvector v] -> [v]) -> [Eigenvector v] -> [v]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [Eigenvector v] -> [v]
newSys)
                                         ([Eigenvector v] -> [[Eigenvector v]])
-> ([v] -> [Eigenvector v]) -> [v] -> [[Eigenvector v]]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (v -> Eigenvector v) -> [v] -> [Eigenvector v]
forall a b. (a -> b) -> [a] -> [b]
map (v -> Eigenvector v
asEV (v -> Eigenvector v) -> (v -> v) -> v -> Eigenvector v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (v -> Norm v -> v
forall v. (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
^%Norm v
me))
 where newSys :: [Eigenvector v] -> [v]
newSys [] = []
       newSys (Eigenvector Scalar v
λ v
v v
fv v
dv Scalar v
ε : [Eigenvector v]
evs)
         | Scalar v
εScalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
>Scalar v
ε₀       = case [Eigenvector v] -> [v]
newSys [Eigenvector v]
evs of
                         []     -> [v
fvv -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v
λ, v
dvv -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v -> Scalar v
forall a. Floating a => a -> a
sqrt(Scalar v
εScalar v -> Scalar v -> Scalar v
forall a. Num a => a -> a -> a
+Scalar v
ε₀)]
                         v
vn:[v]
vns -> v
fvv -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v
λ v -> [v] -> [v]
forall a. a -> [a] -> [a]
: v
vn v -> [v] -> [v]
forall a. a -> [a] -> [a]
: v
dvv -> Scalar v -> v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v -> Scalar v
forall a. Floating a => a -> a
sqrt(Scalar v
εScalar v -> Scalar v -> Scalar v
forall a. Num a => a -> a -> a
+Scalar v
ε₀) v -> [v] -> [v]
forall a. a -> [a] -> [a]
: [v]
vns
         | Scalar v
εScalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
>=Scalar v
0       = v
v v -> [v] -> [v]
forall a. a -> [a] -> [a]
: [Eigenvector v] -> [v]
newSys [Eigenvector v]
evs
         | Bool
otherwise  = [Eigenvector v] -> [v]
newSys [Eigenvector v]
evs
       asEV :: v -> Eigenvector v
asEV v
v = Scalar v -> v -> v -> v -> Scalar v -> Eigenvector v
forall v. Scalar v -> v -> v -> v -> Scalar v -> Eigenvector v
Eigenvector Scalar v
λ v
v v
fv v
dv Scalar v
ε
        where λ² :: Scalar v
λ² = DualVector v
fv'DualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
fv
              λ :: Scalar v
λ = DualVector v
fv'DualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v
              ε :: Scalar v
ε = Norm v -> v -> Scalar v
forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm v
me v
dv
              fv :: v
fv = v -+> v
f (v -+> v) -> v -> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v
v
              fv' :: DualVector v
fv' = Norm v
meNorm v -> v -> DualVector v
forall v. LSpace v => Norm v -> v -> DualVector v
<$|v
fv
              dv :: v
dv | Scalar v
λ²Scalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
>Scalar v
0       = v
v v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
fvv -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*(Scalar v
λScalar v -> Scalar v -> Scalar v
forall a. Fractional a => a -> a -> a
/Scalar v
λ²) -- for stability reasons
                 | Bool
otherwise  = v
forall v. AdditiveGroup v => v
zeroV


finishEigenSystem :: ∀ v . (LSpace v, RealFloat (Scalar v))
                      => Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem :: Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem Norm v
me = [Eigenvector v] -> [Eigenvector v]
go ([Eigenvector v] -> [Eigenvector v])
-> ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v]
-> [Eigenvector v]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (Eigenvector v -> Eigenvector v -> Ordering)
-> [Eigenvector v] -> [Eigenvector v]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy ((Eigenvector v -> Scalar v)
-> Eigenvector v -> Eigenvector v -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing ((Eigenvector v -> Scalar v)
 -> Eigenvector v -> Eigenvector v -> Ordering)
-> (Eigenvector v -> Scalar v)
-> Eigenvector v
-> Eigenvector v
-> Ordering
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Scalar v -> Scalar v
forall a. Num a => a -> a
negate (Scalar v -> Scalar v)
-> (Eigenvector v -> Scalar v) -> Eigenvector v -> Scalar v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Eigenvector v -> Scalar v
forall v. Eigenvector v -> Scalar v
ev_Eigenvalue)
 where go :: [Eigenvector v] -> [Eigenvector v]
go [] = []
       go [Eigenvector v
v] = [Eigenvector v
v]
       go vs :: [Eigenvector v]
vs@[Eigenvector Scalar v
λ₀ v
v₀ v
fv₀ v
_dv₀ Scalar v
_ε₀, Eigenvector Scalar v
λ₁ v
v₁ v
fv₁ v
_dv₁ Scalar v
_ε₁]
          | Scalar v
λ₀Scalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
>Scalar v
λ₁      = [ v -> v -> Eigenvector v
asEV v
v₀' v
fv₀', v -> v -> Eigenvector v
asEV v
v₁' v
fv₁' ]
          | Bool
otherwise  = [Eigenvector v]
vs
        where
              v₀' :: v
v₀' = v
v₀v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₀ v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
v₁v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₁
              fv₀' :: v
fv₀' = v
fv₀v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₀ v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
fv₁v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₁
              
              v₁' :: v
v₁' = v
v₀v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₀ v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
v₁v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₁
              fv₁' :: v
fv₁' = v
fv₀v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₀ v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^+^ v
fv₁v -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₁
              
              fShift₁v₀ :: v
fShift₁v₀ = v
fv₀ v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ Scalar v
λ₁Scalar v -> v -> v
forall v. VectorSpace v => Scalar v -> v -> v
*^v
v₀
              
              (Scalar v
μ₀₀,Scalar v
μ₀₁) = (Scalar v, Scalar v) -> (Scalar v, Scalar v)
forall b. Floating b => (b, b) -> (b, b)
normalised ( Scalar v
λ₀ Scalar v -> Scalar v -> Scalar v
forall a. Num a => a -> a -> a
- Scalar v
λ₁
                                     , (Norm v
me Norm v -> v -> DualVector v
forall v. LSpace v => Norm v -> v -> DualVector v
<$| v
fShift₁v₀)DualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v₁ )
              (Scalar v
μ₁₀,Scalar v
μ₁₁) = (-Scalar v
μ₀₁, Scalar v
μ₀₀)
        
       go [Eigenvector v]
evs = [Eigenvector v]
lo'' [Eigenvector v] -> [Eigenvector v] -> [Eigenvector v]
forall a. [a] -> [a] -> [a]
++ [Eigenvector v]
upper'
        where l :: Int
l = [Eigenvector v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs
              lChunk :: Int
lChunk = Int
lInt -> Int -> Int
forall a. Integral a => a -> a -> a
`quot`Int
3
              ([Eigenvector v]
loEvs, ([Eigenvector v]
midEvs, [Eigenvector v]
hiEvs)) = ([Eigenvector v] -> ([Eigenvector v], [Eigenvector v]))
-> ([Eigenvector v], [Eigenvector v])
-> ([Eigenvector v], ([Eigenvector v], [Eigenvector v]))
forall (a :: * -> * -> *) d b c.
(Morphism a, ObjectPair a d b, ObjectPair a d c) =>
a b c -> a (d, b) (d, c)
second (Int -> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall a. Int -> [a] -> ([a], [a])
splitAt (Int -> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v]))
-> Int -> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
lChunk)
                                                    (([Eigenvector v], [Eigenvector v])
 -> ([Eigenvector v], ([Eigenvector v], [Eigenvector v])))
-> ([Eigenvector v], [Eigenvector v])
-> ([Eigenvector v], ([Eigenvector v], [Eigenvector v]))
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Int -> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
lChunk [Eigenvector v]
evs
              ([Eigenvector v]
lo',[Eigenvector v]
hi') = Int -> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
lChunk ([Eigenvector v] -> ([Eigenvector v], [Eigenvector v]))
-> ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v]
-> ([Eigenvector v], [Eigenvector v])
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [Eigenvector v] -> [Eigenvector v]
go ([Eigenvector v] -> ([Eigenvector v], [Eigenvector v]))
-> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [Eigenvector v]
loEvs[Eigenvector v] -> [Eigenvector v] -> [Eigenvector v]
forall a. [a] -> [a] -> [a]
++[Eigenvector v]
hiEvs
              ([Eigenvector v]
lo'',[Eigenvector v]
mid') = Int -> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
lChunk ([Eigenvector v] -> ([Eigenvector v], [Eigenvector v]))
-> ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v]
-> ([Eigenvector v], [Eigenvector v])
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [Eigenvector v] -> [Eigenvector v]
go ([Eigenvector v] -> ([Eigenvector v], [Eigenvector v]))
-> [Eigenvector v] -> ([Eigenvector v], [Eigenvector v])
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [Eigenvector v]
lo'[Eigenvector v] -> [Eigenvector v] -> [Eigenvector v]
forall a. [a] -> [a] -> [a]
++[Eigenvector v]
midEvs
              upper' :: [Eigenvector v]
upper'  = [Eigenvector v] -> [Eigenvector v]
go ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v] -> [Eigenvector v]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [Eigenvector v]
mid'[Eigenvector v] -> [Eigenvector v] -> [Eigenvector v]
forall a. [a] -> [a] -> [a]
++[Eigenvector v]
hi'
       
       asEV :: v -> v -> Eigenvector v
asEV v
v v
fv = Scalar v -> v -> v -> v -> Scalar v -> Eigenvector v
forall v. Scalar v -> v -> v -> v -> Scalar v -> Eigenvector v
Eigenvector Scalar v
λ v
v v
fv v
dv Scalar v
ε
        where λ :: Scalar v
λ = (Norm v
meNorm v -> v -> DualVector v
forall v. LSpace v => Norm v -> v -> DualVector v
<$|v
v)DualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
fv
              dv :: v
dv = v
vv -> Scalar v -> v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
λ v -> v -> v
forall v. AdditiveGroup v => v -> v -> v
^-^ v
fv
              ε :: Scalar v
ε = Norm v -> v -> Scalar v
forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm v
me v
dv Scalar v -> Scalar v -> Scalar v
forall a. Fractional a => a -> a -> a
/ Scalar v
λScalar v -> Int -> Scalar v
forall a. Num a => a -> Int -> a
^Int
2
       
       normalised :: (b, b) -> (b, b)
normalised (b
x,b
y) = (b
xb -> b -> b
forall a. Fractional a => a -> a -> a
/b
r, b
yb -> b -> b
forall a. Fractional a => a -> a -> a
/b
r)
        where r :: b
r = b -> b
forall a. Floating a => a -> a
sqrt (b -> b) -> b -> b
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ b
xb -> Int -> b
forall a. Num a => a -> Int -> a
^Int
2 b -> b -> b
forall a. Num a => a -> a -> a
+ b
yb -> Int -> b
forall a. Num a => a -> Int -> a
^Int
2

-- | Find a system of vectors that approximate the eigensytem, in the sense that:
--   each true eigenvalue is represented by an approximate one, and that is closer
--   to the true value than all the other approximate EVs.
-- 
--   This function does not make any guarantees as to how well a single eigenvalue
--   is approximated, though.
roughEigenSystem :: (FiniteDimensional v, IEEE (Scalar v))
        => Norm v
        -> (v+>v)
        -> [Eigenvector v]
roughEigenSystem :: Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem Norm v
me v +> v
f = [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go [v]
fBas Int
0 [[]]
 where go :: [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go [] Int
_ ([Eigenvector v]
_:[Eigenvector v]
evs:[[Eigenvector v]]
_) = [Eigenvector v]
evs
       go [] Int
_ ([Eigenvector v]
evs:[[Eigenvector v]]
_) = [Eigenvector v]
evs
       go (v
v:[v]
vs) Int
oldDim ([Eigenvector v]
evs:[[Eigenvector v]]
evss)
         | Norm v -> v -> Scalar v
forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm v
me v
vPerp Scalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
> Scalar v
fpε  = case [[Eigenvector v]]
evss of
             [Eigenvector v]
evs':[[Eigenvector v]]
_ | [Eigenvector v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
oldDim
               -> [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go (v
vv -> [v] -> [v]
forall a. a -> [a] -> [a]
:[v]
vs) ([Eigenvector v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs) [[Eigenvector v]]
evss
             [[Eigenvector v]]
_ -> let evss' :: [[Eigenvector v]]
evss' = [[Eigenvector v]] -> [[Eigenvector v]]
forall a. [a] -> [a]
tail ([[Eigenvector v]] -> [[Eigenvector v]])
-> ([v] -> [[Eigenvector v]]) -> [v] -> [[Eigenvector v]]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Norm v -> Scalar v -> (v -+> v) -> [v] -> [[Eigenvector v]]
forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> Scalar v -> (v -+> v) -> [v] -> [[Eigenvector v]]
constructEigenSystem Norm v
me Scalar v
fpε ((v +> v) -> v -+> v
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr v +> v
f)
                                ([v] -> [[Eigenvector v]]) -> [v] -> [[Eigenvector v]]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ (Eigenvector v -> v) -> [Eigenvector v] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map Eigenvector v -> v
forall v. Eigenvector v -> v
ev_Eigenvector ([[Eigenvector v]] -> [Eigenvector v]
forall a. [a] -> a
head ([[Eigenvector v]] -> [Eigenvector v])
-> [[Eigenvector v]] -> [Eigenvector v]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [[Eigenvector v]]
evss[[Eigenvector v]] -> [[Eigenvector v]] -> [[Eigenvector v]]
forall a. [a] -> [a] -> [a]
++[[Eigenvector v]
evs]) [v] -> [v] -> [v]
forall a. [a] -> [a] -> [a]
++ [v
vPerp]
                  in [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go [v]
vs ([Eigenvector v] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs) [[Eigenvector v]]
evss'
         | Bool
otherwise              = [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go [v]
vs Int
oldDim ([Eigenvector v]
evs[Eigenvector v] -> [[Eigenvector v]] -> [[Eigenvector v]]
forall a. a -> [a] -> [a]
:[[Eigenvector v]]
evss)
        where vPerp :: v
vPerp = Norm v -> [v] -> v -+> v
forall v. LSpace v => Norm v -> [v] -> v -+> v
orthogonalComplementProj Norm v
me (Eigenvector v -> v
forall v. Eigenvector v -> v
ev_Eigenvector(Eigenvector v -> v) -> [Eigenvector v] -> [v]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[Eigenvector v]
evs) (v -+> v) -> v -> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v
v
       fBas :: [v]
fBas = (v -> Norm v -> v
forall v. (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
^%Norm v
me) (v -> v) -> [v] -> [v]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> (SubBasis v, [v] -> [v]) -> [v] -> [v]
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd ((v +> v) -> (SubBasis v, [v] -> [v])
forall v w.
(FiniteDimensional v, LSpace w, Scalar w ~ Scalar v) =>
(v +> w) -> (SubBasis v, DList w)
decomposeLinMap v +> v
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id) []
       fpε :: Scalar v
fpε = Scalar v
forall a. IEEE a => a
epsilon Scalar v -> Scalar v -> Scalar v
forall a. Num a => a -> a -> a
* Scalar v
8

-- | Simple automatic finding of the eigenvalues and -vectors
--   of a Hermitian operator, in reasonable approximation.
-- 
--   This works by spanning a QR-stabilised Krylov basis with 'constructEigenSystem'
--   until it is complete ('roughEigenSystem'), and then properly decoupling the
--   system with 'finishEigenSystem' (based on two iterations of shifted Givens rotations).
--   
--   This function is a tradeoff in performance vs. accuracy. Use 'constructEigenSystem'
--   and 'finishEigenSystem' directly for more quickly computing a (perhaps incomplete)
--   approximation, or for more precise results.
eigen :: (FiniteDimensional v, HilbertSpace v, IEEE (Scalar v))
               => (v+>v) -> [(Scalar v, v)]
eigen :: (v +> v) -> [(Scalar v, v)]
eigen v +> v
f = (Eigenvector v -> (Scalar v, v))
-> [Eigenvector v] -> [(Scalar v, v)]
forall a b. (a -> b) -> [a] -> [b]
map (Eigenvector v -> Scalar v
forall v. Eigenvector v -> Scalar v
ev_Eigenvalue (Eigenvector v -> Scalar v)
-> (Eigenvector v -> v) -> Eigenvector v -> (Scalar v, v)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&& Eigenvector v -> v
forall v. Eigenvector v -> v
ev_Eigenvector)
   ([Eigenvector v] -> [(Scalar v, v)])
-> [Eigenvector v] -> [(Scalar v, v)]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v] -> [[Eigenvector v]]
forall a. (a -> a) -> a -> [a]
iterate (Norm v -> [Eigenvector v] -> [Eigenvector v]
forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem Norm v
forall v. HilbertSpace v => Norm v
euclideanNorm) (Norm v -> (v +> v) -> [Eigenvector v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem Norm v
forall v. HilbertSpace v => Norm v
euclideanNorm v +> v
f) [[Eigenvector v]] -> Int -> [Eigenvector v]
forall a. [a] -> Int -> a
!! Int
2


-- | Approximation of the determinant.
roughDet :: (FiniteDimensional v, IEEE (Scalar v)) => (v+>v) -> Scalar v
roughDet :: (v +> v) -> Scalar v
roughDet = Norm v -> (v +> v) -> [Eigenvector v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem Norm v
forall v. FiniteDimensional v => Norm v
adhocNorm ((v +> v) -> [Eigenvector v])
-> ([Eigenvector v] -> Scalar v) -> (v +> v) -> Scalar v
forall (k :: * -> * -> *) a b c.
(Category k, Object k a, Object k b, Object k c) =>
k a b -> k b c -> k a c
>>> (Eigenvector v -> Scalar v) -> [Eigenvector v] -> [Scalar v]
forall a b. (a -> b) -> [a] -> [b]
map Eigenvector v -> Scalar v
forall v. Eigenvector v -> Scalar v
ev_Eigenvalue ([Eigenvector v] -> [Scalar v])
-> ([Scalar v] -> Scalar v) -> [Eigenvector v] -> Scalar v
forall (k :: * -> * -> *) a b c.
(Category k, Object k a, Object k b, Object k c) =>
k a b -> k b c -> k a c
>>> [Scalar v] -> Scalar v
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product


orthonormalityError :: LSpace v => Norm v -> [v] -> Scalar v
orthonormalityError :: Norm v -> [v] -> Scalar v
orthonormalityError Norm v
me [v]
vs = Norm v -> v -> Scalar v
forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm v
me (v -> Scalar v) -> v -> Scalar v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Norm v -> [v] -> LinearFunction (Scalar v) v v
forall v. LSpace v => Norm v -> [v] -> v -+> v
orthogonalComplementProj Norm v
me [v]
vs LinearFunction (Scalar v) v v -> v -> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [v] -> v
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [v]
vs


normSpanningSystem :: SimpleSpace v
               => Seminorm v -> [DualVector v]
normSpanningSystem :: Seminorm v -> [DualVector v]
normSpanningSystem = ((v, DualVector v) -> DualVector v)
-> [(v, DualVector v)] -> [DualVector v]
forall a b. (a -> b) -> [a] -> [b]
map (v, DualVector v) -> DualVector v
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd ([(v, DualVector v)] -> [DualVector v])
-> (Seminorm v -> [(v, DualVector v)])
-> Seminorm v
-> [DualVector v]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Seminorm v -> [(v, DualVector v)]
forall v. SimpleSpace v => Seminorm v -> [(v, DualVector v)]
normSpanningSystems

normSpanningSystems :: SimpleSpace v
               => Seminorm v -> [(v, DualVector v)]
normSpanningSystems :: Seminorm v -> [(v, DualVector v)]
normSpanningSystems me :: Seminorm v
me@(Norm v -+> DualVector v
m)
     = [Maybe (v, DualVector v)] -> [(v, DualVector v)]
forall a. [Maybe a] -> [a]
catMaybes ([Maybe (v, DualVector v)] -> [(v, DualVector v)])
-> ([v] -> [Maybe (v, DualVector v)]) -> [v] -> [(v, DualVector v)]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. ((v, Maybe (DualVector v)) -> Maybe (v, DualVector v))
-> [(v, Maybe (DualVector v))] -> [Maybe (v, DualVector v)]
forall a b. (a -> b) -> [a] -> [b]
map (\(v
v,Maybe (DualVector v)
d)->(v
v,)(DualVector v -> (v, DualVector v))
-> Maybe (DualVector v) -> Maybe (v, DualVector v)
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>Maybe (DualVector v)
d) ([(v, Maybe (DualVector v))] -> [Maybe (v, DualVector v)])
-> ([v] -> [(v, Maybe (DualVector v))])
-> [v]
-> [Maybe (v, DualVector v)]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. Scalar v -> [(v, DualVector v)] -> [(v, Maybe (DualVector v))]
forall v.
(SemiInner v, RealFrac' (Scalar v)) =>
Scalar v -> [(v, DualVector v)] -> [(v, Maybe (DualVector v))]
orthonormaliseDuals Scalar v
0
         ([(v, DualVector v)] -> [(v, Maybe (DualVector v))])
-> ([v] -> [(v, DualVector v)])
-> [v]
-> [(v, Maybe (DualVector v))]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (v -> (v, DualVector v)) -> [v] -> [(v, DualVector v)]
forall a b. (a -> b) -> [a] -> [b]
map (v -> v
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id(v -> v) -> (v -> DualVector v) -> v -> (v, DualVector v)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&(v -+> DualVector v
m(v -+> DualVector v) -> v -> DualVector v
forall s v w. LinearFunction s v w -> v -> w
-+$>)) ([v] -> [(v, DualVector v)]) -> [v] -> [(v, DualVector v)]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Seminorm v -> [v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem' Seminorm v
me

normSpanningSystem' :: (FiniteDimensional v, IEEE (Scalar v))
               => Seminorm v -> [v]
normSpanningSystem' :: Seminorm v -> [v]
normSpanningSystem' Seminorm v
me = Scalar v -> Seminorm v -> [v] -> [v]
forall v.
(LSpace v, RealFloat (Scalar v)) =>
Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily Scalar v
0 Seminorm v
me ([v] -> [v]) -> [v] -> [v]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ SubBasis v -> [v]
forall v. FiniteDimensional v => SubBasis v -> [v]
enumerateSubBasis SubBasis v
forall v. FiniteDimensional v => SubBasis v
entireBasis

-- | Inverse of 'spanVariance'. Equivalent to 'normSpanningSystem' on the dual space.
varianceSpanningSystem :: ∀ v . SimpleSpace v => Variance v -> [v]
varianceSpanningSystem :: Variance v -> [v]
varianceSpanningSystem = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
                           DualSpaceWitness v
DualSpaceWitness -> Variance v -> [v]
forall v. SimpleSpace v => Seminorm v -> [DualVector v]
normSpanningSystem

-- | For any two norms, one can find a system of co-vectors that, with suitable
--   coefficients, spans /either/ of them: if @shSys = sharedNormSpanningSystem n₀ n₁@,
--   then
-- 
-- @
-- n₀ = 'spanNorm' $ fst<$>shSys
-- @
-- 
-- and
-- 
-- @
-- n₁ = 'spanNorm' [dv^*η | (dv,η)<-shSys]
-- @
-- 
-- A rather crude approximation ('roughEigenSystem') is used in this function, so do
-- not expect the above equations to hold with great accuracy.
sharedNormSpanningSystem :: SimpleSpace v
               => Norm v -> Seminorm v -> [(DualVector v, Scalar v)]
sharedNormSpanningSystem :: Norm v -> Norm v -> [(DualVector v, Scalar v)]
sharedNormSpanningSystem nn :: Norm v
nn@(Norm v -+> DualVector v
n) Norm v
nm
      = (v -> DualVector v) -> (v, Scalar v) -> (DualVector v, Scalar v)
forall (a :: * -> * -> *) b d c.
(Morphism a, ObjectPair a b d, ObjectPair a c d) =>
a b c -> a (b, d) (c, d)
first (v -+> DualVector v
n(v -+> DualVector v) -> v -> DualVector v
forall s v w. LinearFunction s v w -> v -> w
-+$>) ((v, Scalar v) -> (DualVector v, Scalar v))
-> [(v, Scalar v)] -> [(DualVector v, Scalar v)]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' Int
0 (Norm v
nn, Norm v -> Variance v
forall v. SimpleSpace v => Norm v -> Variance v
dualNorm Norm v
nn) Norm v
nm

sharedNormSpanningSystem' :: ∀ v . SimpleSpace v
               => Int -> (Norm v, Variance v) -> Seminorm v -> [(v, Scalar v)]
sharedNormSpanningSystem' :: Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' = DualSpaceWitness v
-> Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
snss DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness
 where snss :: DualSpaceWitness v -> Int -> (Norm v, Variance v)
                     -> Seminorm v -> [(v, Scalar v)]
       snss :: DualSpaceWitness v
-> Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
snss DualSpaceWitness v
DualSpaceWitness Int
nRefine (nn :: Norm v
nn@(Norm v -+> DualVector v
n), Norm DualVector v -+> DualVector (DualVector v)
n') (Norm v -+> DualVector v
m)
           = Eigenvector v -> [(v, Scalar v)]
forall a.
(Ord (Scalar a), Floating (Scalar a)) =>
Eigenvector a -> [(a, Scalar a)]
sep (Eigenvector v -> [(v, Scalar v)])
-> [Eigenvector v] -> [(v, Scalar v)]
forall (m :: * -> *) (k :: * -> * -> *) a b.
(Monad m k, Object k a, Object k b, Object k (m a), Object k (m b),
 Object k (m (m b))) =>
k a (m b) -> k (m a) (m b)
=<< ([Eigenvector v] -> [Eigenvector v])
-> [Eigenvector v] -> [[Eigenvector v]]
forall a. (a -> a) -> a -> [a]
iterate (Norm v -> [Eigenvector v] -> [Eigenvector v]
forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem Norm v
nn)
                        (Norm v -> LinearMap (Scalar v) v v -> [Eigenvector v]
forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem Norm v
nn (LinearMap (Scalar v) v v -> [Eigenvector v])
-> LinearMap (Scalar v) v v -> [Eigenvector v]
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearFunction (Scalar v) (DualVector v) v
-> LinearMap (Scalar v) (DualVector v) v
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr LinearFunction (Scalar v) (DualVector v) v
DualVector v -+> DualVector (DualVector v)
n' LinearMap (Scalar v) (DualVector v) v
-> LinearMap (Scalar v) v (DualVector v)
-> LinearMap (Scalar v) v v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearFunction (Scalar v) v (DualVector v)
-> LinearMap (Scalar v) v (DualVector v)
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr LinearFunction (Scalar v) v (DualVector v)
v -+> DualVector v
m) [[Eigenvector v]] -> Int -> [Eigenvector v]
forall a. [a] -> Int -> a
!! Int
nRefine
       sep :: Eigenvector a -> [(a, Scalar a)]
sep (Eigenvector Scalar a
λ a
v a
λv a
_ Scalar a
_)
         | Scalar a
λScalar a -> Scalar a -> Bool
forall a. Ord a => a -> a -> Bool
>=Scalar a
0       = [(a
v, Scalar a -> Scalar a
forall a. Floating a => a -> a
sqrt Scalar a
λ)]
         | Bool
otherwise  = []

-- | Like 'sharedNormSpanningSystem n₀ n₁', but allows /either/ of the norms
--   to be singular.
-- 
-- @
-- n₀ = 'spanNorm' [dv | (dv, Just _)<-shSys]
-- @
-- 
-- and
-- 
-- @
-- n₁ = 'spanNorm' $ [dv^*η | (dv, Just η)<-shSys]
--                 ++ [ dv | (dv, Nothing)<-shSys]
-- @
-- 
-- You may also interpret a @Nothing@ here as an “infinite eigenvalue”, i.e.
-- it is so small as an spanning vector of @n₀@ that you would need to scale it
-- by ∞ to use it for spanning @n₁@.
sharedSeminormSpanningSystem :: ∀ v . SimpleSpace v
               => Seminorm v -> Seminorm v -> [(DualVector v, Maybe (Scalar v))]
sharedSeminormSpanningSystem :: Seminorm v -> Seminorm v -> [(DualVector v, Maybe (Scalar v))]
sharedSeminormSpanningSystem Seminorm v
nn Seminorm v
nm
         = DualSpaceWitness v
-> (v, Scalar v) -> (DualVector v, Maybe (Scalar v))
finalise DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness
               ((v, Scalar v) -> (DualVector v, Maybe (Scalar v)))
-> [(v, Scalar v)] -> [(DualVector v, Maybe (Scalar v))]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Int -> (Seminorm v, Variance v) -> Seminorm v -> [(v, Scalar v)]
forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' Int
1 (Seminorm v
combined, Seminorm v -> Variance v
forall v. SimpleSpace v => Norm v -> Variance v
dualNorm Seminorm v
combined) Seminorm v
nn
 where combined :: Seminorm v
combined = Seminorm v -> Seminorm v
forall v. LSpace v => Norm v -> Norm v
densifyNorm (Seminorm v -> Seminorm v) -> Seminorm v -> Seminorm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Seminorm v
nnSeminorm v -> Seminorm v -> Seminorm v
forall a. Semigroup a => a -> a -> a
<>Seminorm v
nm
       finalise :: DualSpaceWitness v -> (v, Scalar v) -> (DualVector v, Maybe (Scalar v))
       finalise :: DualSpaceWitness v
-> (v, Scalar v) -> (DualVector v, Maybe (Scalar v))
finalise DualSpaceWitness v
DualSpaceWitness (v
v, Scalar v
μn)
           | Scalar v
μnScalar v -> Int -> Scalar v
forall a. Num a => a -> Int -> a
^Int
2 Scalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
> Scalar v
forall a. IEEE a => a
epsilon  = (DualVector v
v'DualVector v -> Scalar v -> DualVector v
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μn, Scalar v -> Maybe (Scalar v)
forall a. a -> Maybe a
Just (Scalar v -> Maybe (Scalar v)) -> Scalar v -> Maybe (Scalar v)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Scalar v -> Scalar v
forall a. Floating a => a -> a
sqrt (Scalar v -> Scalar v -> Scalar v
forall a. Ord a => a -> a -> a
max Scalar v
0 (Scalar v -> Scalar v) -> Scalar v -> Scalar v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Scalar v
1 Scalar v -> Scalar v -> Scalar v
forall a. Num a => a -> a -> a
- Scalar v
μnScalar v -> Int -> Scalar v
forall a. Num a => a -> Int -> a
^Int
2)Scalar v -> Scalar v -> Scalar v
forall a. Fractional a => a -> a -> a
/Scalar v
μn)
           | Bool
otherwise       = (DualVector v
v', Maybe (Scalar v)
forall a. Maybe a
Nothing)
        where v' :: DualVector v
v' = Seminorm v
combinedSeminorm v -> v -> DualVector v
forall v. LSpace v => Norm v -> v -> DualVector v
<$|v
v

-- | A system of vectors which are orthogonal with respect to both of the given
--   seminorms. (In general they are not /orthonormal/ to either of them.)
sharedSeminormSpanningSystem' :: ∀ v .  SimpleSpace v
               => Seminorm v -> Seminorm v -> [v]
sharedSeminormSpanningSystem' :: Seminorm v -> Seminorm v -> [v]
sharedSeminormSpanningSystem' Seminorm v
nn Seminorm v
nm
         = (v, Scalar v) -> v
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst ((v, Scalar v) -> v) -> [(v, Scalar v)] -> [v]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> Int -> (Seminorm v, Variance v) -> Seminorm v -> [(v, Scalar v)]
forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' Int
1 (Seminorm v
combined, Seminorm v -> Variance v
forall v. SimpleSpace v => Norm v -> Variance v
dualNorm Seminorm v
combined) Seminorm v
nn
 where combined :: Seminorm v
combined = Seminorm v -> Seminorm v
forall v. LSpace v => Norm v -> Norm v
densifyNorm (Seminorm v -> Seminorm v) -> Seminorm v -> Seminorm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Seminorm v
nnSeminorm v -> Seminorm v -> Seminorm v
forall a. Semigroup a => a -> a -> a
<>Seminorm v
nm


-- | Interpret a variance as a covariance between two subspaces, and
--   normalise it by the variance on @u@. The result is effectively
--   the linear regression coefficient of a simple regression of the vectors
--   spanning the variance.
dependence :: ∀ u v . (SimpleSpace u, SimpleSpace v, Scalar u~Scalar v)
                => Variance (u,v) -> (u+>v)
dependence :: Variance (u, v) -> u +> v
dependence = case ( DualSpaceWitness u
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness u
                  , DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v ) of
  (DualSpaceWitness u
DualSpaceWitness,DualSpaceWitness v
DualSpaceWitness)
        -> \(Norm DualVector (u, v) -+> DualVector (DualVector (u, v))
m) -> LinearFunction (Scalar v) (DualVector u) v
-> LinearFunction
     (Scalar v) (LinearMap (Scalar u) u (DualVector u)) (u +> v)
forall (f :: * -> *) (r :: * -> * -> *) (t :: * -> * -> *) a b.
(Functor f r t, Object r a, Object t (f a), Object r b,
 Object t (f b)) =>
r a b -> t (f a) (f b)
fmap ( LinearFunction (Scalar v) (u, v) v
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd LinearFunction (Scalar v) (u, v) v
-> LinearFunction (Scalar v) (DualVector u) (u, v)
-> LinearFunction (Scalar v) (DualVector u) v
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearFunction (Scalar v) (DualVector u, DualVector v) (u, v)
DualVector (u, v) -+> DualVector (DualVector (u, v))
m LinearFunction (Scalar v) (DualVector u, DualVector v) (u, v)
-> LinearFunction
     (Scalar v) (DualVector u) (DualVector u, DualVector v)
-> LinearFunction (Scalar v) (DualVector u) (u, v)
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (LinearFunction (Scalar v) (DualVector u) (DualVector u)
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
idLinearFunction (Scalar v) (DualVector u) (DualVector u)
-> LinearFunction (Scalar v) (DualVector u) (DualVector v)
-> LinearFunction
     (Scalar v) (DualVector u) (DualVector u, DualVector v)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&LinearFunction (Scalar v) (DualVector u) (DualVector v)
forall v. AdditiveGroup v => v
zeroV) )
              LinearFunction
  (Scalar v) (LinearMap (Scalar u) u (DualVector u)) (u +> v)
-> LinearMap (Scalar u) u (DualVector u) -> u +> v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ (DualVector u +> u) -> LinearMap (Scalar u) u (DualVector u)
forall u v.
(SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) =>
(u +> v) -> v +> u
pseudoInverse (LinearFunction (Scalar v) (DualVector u) u -> DualVector u +> u
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr (LinearFunction (Scalar v) (DualVector u) u -> DualVector u +> u)
-> LinearFunction (Scalar v) (DualVector u) u -> DualVector u +> u
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearFunction (Scalar v) (u, v) u
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst LinearFunction (Scalar v) (u, v) u
-> LinearFunction (Scalar v) (DualVector u) (u, v)
-> LinearFunction (Scalar v) (DualVector u) u
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearFunction (Scalar v) (DualVector u, DualVector v) (u, v)
DualVector (u, v) -+> DualVector (DualVector (u, v))
m LinearFunction (Scalar v) (DualVector u, DualVector v) (u, v)
-> LinearFunction
     (Scalar v) (DualVector u) (DualVector u, DualVector v)
-> LinearFunction (Scalar v) (DualVector u) (u, v)
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (LinearFunction (Scalar v) (DualVector u) (DualVector u)
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
idLinearFunction (Scalar v) (DualVector u) (DualVector u)
-> LinearFunction (Scalar v) (DualVector u) (DualVector v)
-> LinearFunction
     (Scalar v) (DualVector u) (DualVector u, DualVector v)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&LinearFunction (Scalar v) (DualVector u) (DualVector v)
forall v. AdditiveGroup v => v
zeroV))


summandSpaceNorms :: ∀ u v . (SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v)
                       => Norm (u,v) -> (Norm u, Norm v)
summandSpaceNorms :: Norm (u, v) -> (Norm u, Norm v)
summandSpaceNorms = case ( DualSpaceWitness u
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness u
                         , DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v ) of
  (DualSpaceWitness u
DualSpaceWitness,DualSpaceWitness v
DualSpaceWitness)
        -> \Norm (u, v)
nuv -> let spanSys :: [DualVector (u, v)]
spanSys = Norm (u, v) -> [DualVector (u, v)]
forall v. SimpleSpace v => Seminorm v -> [DualVector v]
normSpanningSystem Norm (u, v)
nuv
                   in ( Norm u -> Norm u
forall v. LSpace v => Norm v -> Norm v
densifyNorm (Norm u -> Norm u) -> Norm u -> Norm u
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [DualVector u] -> Norm u
forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm ((DualVector u, DualVector v) -> DualVector u
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst((DualVector u, DualVector v) -> DualVector u)
-> [(DualVector u, DualVector v)] -> [DualVector u]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[(DualVector u, DualVector v)]
[DualVector (u, v)]
spanSys)
                      , Norm v -> Norm v
forall v. LSpace v => Norm v -> Norm v
densifyNorm (Norm v -> Norm v) -> Norm v -> Norm v
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [DualVector v] -> Norm v
forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm ((DualVector u, DualVector v) -> DualVector v
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd((DualVector u, DualVector v) -> DualVector v)
-> [(DualVector u, DualVector v)] -> [DualVector v]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[(DualVector u, DualVector v)]
[DualVector (u, v)]
spanSys) )

sumSubspaceNorms :: ∀ u v . (LSpace u, LSpace v, Scalar u~Scalar v)
                         => Norm u -> Norm v -> Norm (u,v)
sumSubspaceNorms :: Norm u -> Norm v -> Norm (u, v)
sumSubspaceNorms = case ( DualSpaceWitness u
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness u
                         , DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v ) of
  (DualSpaceWitness u
DualSpaceWitness,DualSpaceWitness v
DualSpaceWitness)
        -> \(Norm u -+> DualVector u
nu) (Norm v -+> DualVector v
nv) -> LinearFunction
  (Scalar (DualVector u)) (u, v) (DualVector u, DualVector v)
-> Norm (u, v)
forall v. (v -+> DualVector v) -> Norm v
Norm (LinearFunction
   (Scalar (DualVector u)) (u, v) (DualVector u, DualVector v)
 -> Norm (u, v))
-> LinearFunction
     (Scalar (DualVector u)) (u, v) (DualVector u, DualVector v)
-> Norm (u, v)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ u -+> DualVector u
nu (u -+> DualVector u)
-> LinearFunction (Scalar (DualVector u)) v (DualVector v)
-> LinearFunction
     (Scalar (DualVector u)) (u, v) (DualVector u, DualVector v)
forall (a :: * -> * -> *) b b' c c'.
(Morphism a, ObjectPair a b b', ObjectPair a c c') =>
a b c -> a b' c' -> a (b, b') (c, c')
*** LinearFunction (Scalar (DualVector u)) v (DualVector v)
v -+> DualVector v
nv





instance (SimpleSpace v, Show (DualVector v)) => Show (Norm v) where
  showsPrec :: Int -> Norm v -> ShowS
showsPrec Int
p Norm v
n = Bool -> ShowS -> ShowS
showParen (Int
pInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
9) (ShowS -> ShowS) -> ShowS -> ShowS
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ ([Char]
"spanNorm "[Char] -> ShowS
forall a. [a] -> [a] -> [a]
++) ShowS -> ShowS -> ShowS
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [DualVector v] -> ShowS
forall a. Show a => a -> ShowS
shows (Norm v -> [DualVector v]
forall v. SimpleSpace v => Seminorm v -> [DualVector v]
normSpanningSystem Norm v
n)

type LinearShowable v = (Show v, RieszDecomposable v)



convexPolytopeHull :: ∀ v . SimpleSpace v => [v] -> [DualVector v]
convexPolytopeHull :: [v] -> [DualVector v]
convexPolytopeHull [v]
vs = case DualSpaceWitness v
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
         DualSpaceWitness v
DualSpaceWitness
             -> [DualVector v
dvDualVector v -> Scalar v -> DualVector v
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v
η | (DualVector v
dv,Scalar v
η) <- [(DualVector v, Scalar v)]
candidates, (v -> Bool) -> [v] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Scalar v -> Scalar v -> Bool
forall a. Ord a => a -> a -> Bool
<=Scalar v
η) (Scalar v -> Bool) -> (v -> Scalar v) -> v -> Bool
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (DualVector v
dvDualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^)) [v]
vs]
 where vrv :: Variance v
vrv = [v] -> Variance v
forall v. LSpace v => [v] -> Variance v
spanVariance [v]
vs
       nmv :: Norm v
nmv = Variance v -> Norm v
forall v. SimpleSpace v => Variance v -> Norm v
dualNorm' Variance v
vrv
       candidates :: [(DualVector v, Scalar v)]
       candidates :: [(DualVector v, Scalar v)]
candidates = [ (DualVector v
dv, DualVector v
dvDualVector v -> v -> Scalar v
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v) | v
v <- [v]
vs
                                   , let dv :: DualVector v
dv = Norm v
nmvNorm v -> v -> DualVector v
forall v. LSpace v => Norm v -> v -> DualVector v
<$|v
v ]

symmetricPolytopeOuterVertices :: ∀ v . SimpleSpace v => [DualVector v] -> [v]
symmetricPolytopeOuterVertices :: [DualVector v] -> [v]
symmetricPolytopeOuterVertices [DualVector v]
dvs
         = [ v -> [(DualVector v, v)] -> v
forall b.
(LinearSpace b, Fractional (Scalar b)) =>
b -> [(DualVector b, b)] -> b
seekExtreme v
forall v. AdditiveGroup v => v
zeroV [(DualVector v, v)]
group | [(DualVector v, v)]
group <- [[(DualVector v, v)]]
candidates ]
 where nmv :: Norm v
       nmv :: Norm v
nmv = [DualVector v] -> Norm v
forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm [DualVector v]
dvs
       vrv :: Variance v
vrv = Norm v -> Variance v
forall v. SimpleSpace v => Norm v -> Variance v
dualNorm Norm v
nmv
       withSomeVect :: [(DualVector v, v)]
       withSomeVect :: [(DualVector v, v)]
withSomeVect = [ (DualVector v
dv, v
v) | DualVector v
dv <- [DualVector v]
dvs
                                , let v :: v
v = DualVector v
dvDualVector v -> Variance v -> v
forall v. LSpace v => DualVector v -> Variance v -> v
|&>Variance v
vrv ]
       ([[(DualVector v, v)]]
candidates, [(DualVector v, v)]
_) = Int
-> Int
-> [(DualVector v, v)]
-> ([[(DualVector v, v)]], [(DualVector v, v)])
forall a. Int -> Int -> [a] -> ([[a]], [a])
multiSplit Int
d (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d) ([(DualVector v, v)]
 -> ([[(DualVector v, v)]], [(DualVector v, v)]))
-> ([(DualVector v, v)] -> [(DualVector v, v)])
-> [(DualVector v, v)]
-> ([[(DualVector v, v)]], [(DualVector v, v)])
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [[(DualVector v, v)]] -> [(DualVector v, v)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[(DualVector v, v)]] -> [(DualVector v, v)])
-> ([(DualVector v, v)] -> [[(DualVector v, v)]])
-> [(DualVector v, v)]
-> [(DualVector v, v)]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. [(DualVector v, v)] -> [[(DualVector v, v)]]
forall a.
SimpleSpace a =>
[(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions ([(DualVector v, v)]
 -> ([[(DualVector v, v)]], [(DualVector v, v)]))
-> [(DualVector v, v)]
-> ([[(DualVector v, v)]], [(DualVector v, v)])
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [(DualVector v, v)]
withSomeVect
       d :: Int
d = SubBasis v -> Int
forall v. FiniteDimensional v => SubBasis v -> Int
subbasisDimension (SubBasis v
forall v. FiniteDimensional v => SubBasis v
entireBasis :: SubBasis v)
       seekExtreme :: b -> [(DualVector b, b)] -> b
seekExtreme b
p₀ [] = b
p₀
       seekExtreme b
p₀ ((DualVector b
dv, b
v) : [(DualVector b, b)]
cs)
           = b -> [(DualVector b, b)] -> b
seekExtreme (b
p₀b -> b -> b
forall v. AdditiveGroup v => v -> v -> v
^+^b
vn) [(DualVector b
dw, b
w b -> b -> b
forall v. AdditiveGroup v => v -> v -> v
^-^ b
vb -> Scalar b -> b
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*((DualVector b
dvDualVector b -> b -> Scalar b
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^b
w) Scalar b -> Scalar b -> Scalar b
forall a. Fractional a => a -> a -> a
/ Scalar b
lv)) | (DualVector b
dw, b
w) <- [(DualVector b, b)]
cs]
        where vn :: b
vn = b
v b -> Scalar b -> b
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^* ((Scalar b
1 Scalar b -> Scalar b -> Scalar b
forall a. Num a => a -> a -> a
- DualVector b
dvDualVector b -> b -> Scalar b
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^b
p₀) Scalar b -> Scalar b -> Scalar b
forall a. Fractional a => a -> a -> a
/ Scalar b
lv)
              lv :: Scalar b
lv = DualVector b
dvDualVector b -> b -> Scalar b
forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^b
v

deinterlacions :: SimpleSpace a => [(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions :: [(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions [(DualVector a, a)]
l = [(DualVector a, a)]
l [(DualVector a, a)]
-> [[(DualVector a, a)]] -> [[(DualVector a, a)]]
forall a. a -> [a] -> [a]
: [(DualVector a, a)] -> [[(DualVector a, a)]]
forall a.
SimpleSpace a =>
[(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions ([(DualVector a, a)]
e [(DualVector a, a)] -> [(DualVector a, a)] -> [(DualVector a, a)]
forall a. [a] -> [a] -> [a]
++ ((DualVector a, a) -> (DualVector a, a))
-> [(DualVector a, a)] -> [(DualVector a, a)]
forall a b. (a -> b) -> [a] -> [b]
map (DualVector a, a) -> (DualVector a, a)
forall v. AdditiveGroup v => v -> v
negateV [(DualVector a, a)]
o)
 where ([(DualVector a, a)]
e,[(DualVector a, a)]
o) = [(DualVector a, a)] -> ([(DualVector a, a)], [(DualVector a, a)])
forall a. [a] -> ([a], [a])
deinterlace [(DualVector a, a)]
l
       deinterlace :: [a] -> ([a], [a])
deinterlace (a
a:a
b:[a]
xs) = (a
aa -> [a] -> [a]
forall a. a -> [a] -> [a]
:)([a] -> [a]) -> ([a] -> [a]) -> ([a], [a]) -> ([a], [a])
forall (a :: * -> * -> *) b b' c c'.
(Morphism a, ObjectPair a b b', ObjectPair a c c') =>
a b c -> a b' c' -> a (b, b') (c, c')
***(a
ba -> [a] -> [a]
forall a. a -> [a] -> [a]
:) (([a], [a]) -> ([a], [a])) -> ([a], [a]) -> ([a], [a])
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [a] -> ([a], [a])
deinterlace [a]
xs
       deinterlace [a]
xs = ([],[a]
xs)
       
-- | Simple wrapper of 'linearRegression'.
linearRegressionW :: ∀ s x m y
    . ( LinearSpace x, SimpleSpace y, SimpleSpace m
      , Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s )
         => Norm y -> (x -> (m +> y)) -> [(x,y)] -> m
linearRegressionW :: Norm y -> (x -> m +> y) -> [(x, y)] -> m
linearRegressionW Norm y
σy x -> m +> y
modelMap = LinearRegressionResult x y m -> m
forall x y m. LinearRegressionResult x y m -> m
linearFit_bestModel
                                   (LinearRegressionResult x y m -> m)
-> ([(x, y)] -> LinearRegressionResult x y m) -> [(x, y)] -> m
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (x -> m +> y) -> [(x, (y, Norm y))] -> LinearRegressionResult x y m
forall s x m y.
(LinearSpace x, SimpleSpace y, SimpleSpace m, Scalar x ~ s,
 Scalar y ~ s, Scalar m ~ s, RealFrac' s) =>
(x -> m +> y) -> [(x, (y, Norm y))] -> LinearRegressionResult x y m
linearRegression x -> m +> y
modelMap ([(x, (y, Norm y))] -> LinearRegressionResult x y m)
-> ([(x, y)] -> [(x, (y, Norm y))])
-> [(x, y)]
-> LinearRegressionResult x y m
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. ((x, y) -> (x, (y, Norm y))) -> [(x, y)] -> [(x, (y, Norm y))]
forall a b. (a -> b) -> [a] -> [b]
map ((y -> (y, Norm y)) -> (x, y) -> (x, (y, Norm y))
forall (a :: * -> * -> *) d b c.
(Morphism a, ObjectPair a d b, ObjectPair a d c) =>
a b c -> a (d, b) (d, c)
second (,Norm y
σy))

data LinearRegressionResult x y m = LinearRegressionResult {
          LinearRegressionResult x y m -> Scalar m
linearFit_χν² :: Scalar m 
           -- ^ How well the data uncertainties match the deviations from the model's
           --   synthetic data.
           -- @
           -- χν² = 1/ν · ∑ δy² / σy²
           -- @
           --   Where @ν@ is the number of degrees of freedom (data values minus model
           --   parameters), @δy = m x - yd@ is the deviation from given data to
           --   the data the model would predict (for each sample point), and @σy@ is
           --   the a-priori measurement uncertainty of the data points. 
           -- 
           --   Values @χν²>1@ indicate that the data could not be described satisfyingly;
           --   @χν²≪1@ suggests overfitting or that the data uncertainties have
           --   been postulated too high.
           -- 
           -- <http://adsabs.harvard.edu/abs/1997ieas.book.....T>
           -- 
           --   If the model is exactly determined or even underdetermined (i.e. @ν≤0@)
           --   then @χν²@ is undefined.
        , LinearRegressionResult x y m -> m
linearFit_bestModel :: m
           -- ^ The model that best corresponds to the data, in a least-squares
           --   sense WRT the supplied norm on the data points. In other words,
           --   this is the model that minimises @∑ δy² / σy²@.
        , LinearRegressionResult x y m -> Norm m
linearFit_modelUncertainty :: Norm m
        }

linearRegressionWVar :: ∀ s x m y
    . ( LinearSpace x, FiniteDimensional y, SimpleSpace m
      , Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s )
         => (x -> (m +> y)) -> [(x, (y, Norm y))] -> (m, [DualVector m])
linearRegressionWVar :: (x -> m +> y) -> [(x, (y, Norm y))] -> (m, [DualVector m])
linearRegressionWVar = case Bool
True of Bool
False -> (x -> m +> y) -> [(x, (y, Norm y))] -> (m, [DualVector m])
forall a. HasCallStack => a
undefined

linearRegression :: ∀ s x m y
    . ( LinearSpace x, SimpleSpace y, SimpleSpace m
      , Scalar x ~ s, Scalar y ~ s, Scalar m ~ s, RealFrac' s )
         => (x -> (m +> y)) -> [(x, (y, Norm y))] -> LinearRegressionResult x y m
linearRegression :: (x -> m +> y) -> [(x, (y, Norm y))] -> LinearRegressionResult x y m
linearRegression = (DualSpaceWitness y, DualSpaceWitness m)
-> (x -> m +> y)
-> [(x, (y, Norm y))]
-> LinearRegressionResult x y m
lrw (DualSpaceWitness y
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness, DualSpaceWitness m
forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness)
 where lrw :: (DualSpaceWitness y, DualSpaceWitness m)
                -> (x -> (m +> y)) -> [(x, (y, Norm y))] -> LinearRegressionResult x y m
       lrw :: (DualSpaceWitness y, DualSpaceWitness m)
-> (x -> m +> y)
-> [(x, (y, Norm y))]
-> LinearRegressionResult x y m
lrw (DualSpaceWitness y
DualSpaceWitness, DualSpaceWitness m
DualSpaceWitness) x -> m +> y
modelMap [(x, (y, Norm y))]
dataxy
         = Scalar m -> m -> Norm m -> LinearRegressionResult x y m
forall x y m.
Scalar m -> m -> Norm m -> LinearRegressionResult x y m
LinearRegressionResult (s
χ²s -> s -> s
forall a. Fractional a => a -> a -> a
/Int -> s
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ν) m
leastSquareSol Norm m
σm
        where leastSquareSol :: m
leastSquareSol = ((m -> DualVector m) -> LinearMap s m (DualVector m)
forall (f :: * -> * -> *) s u v.
(EnhancedCat f (LinearFunction s), LinearSpace u, TensorSpace v,
 Scalar u ~ s, Scalar v ~ s, Object f u, Object f v) =>
(u -> v) -> f u v
lfun ((m -> DualVector m) -> LinearMap s m (DualVector m))
-> (m -> DualVector m) -> LinearMap s m (DualVector m)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [DualVector y] -> DualVector m
forward' ([DualVector y] -> DualVector m)
-> (m -> [DualVector y]) -> m -> DualVector m
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. ((x, (y, Norm y)) -> y -> DualVector y)
-> [(x, (y, Norm y))] -> [y] -> [DualVector y]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Norm y -> y -> DualVector y
forall v. LSpace v => Norm v -> v -> DualVector v
(<$|) (Norm y -> y -> DualVector y)
-> ((x, (y, Norm y)) -> Norm y)
-> (x, (y, Norm y))
-> y
-> DualVector y
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (y, Norm y) -> Norm y
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd ((y, Norm y) -> Norm y)
-> ((x, (y, Norm y)) -> (y, Norm y)) -> (x, (y, Norm y)) -> Norm y
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (x, (y, Norm y)) -> (y, Norm y)
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd) [(x, (y, Norm y))]
dataxy
                                          ([y] -> [DualVector y]) -> (m -> [y]) -> m -> [DualVector y]
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. m -> [y]
forward)
                                 LinearMap (Scalar m) m (DualVector m) -> DualVector m -> m
forall u v.
(SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) =>
(u +> v) -> v -> u
\$ [DualVector y] -> DualVector m
forward' [Norm y
σyNorm y -> y -> DualVector y
forall v. LSpace v => Norm v -> v -> DualVector v
<$|y
y | (x
_,(y
y,Norm y
σy)) <- [(x, (y, Norm y))]
dataxy]
              χ² :: s
χ² = [s] -> s
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Norm y -> y -> Scalar y
forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm y
σy y
δy | (x
x, (y
yd, Norm y
σy)) <- [(x, (y, Norm y))]
dataxy
                                     , let δy :: y
δy = y
yd y -> y -> y
forall v. AdditiveGroup v => v -> v -> v
^-^ (x -> m +> y
modelMap x
x LinearMap s m y -> m -> y
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ m
leastSquareSol) ]
              ν :: Int
ν = [(x, (y, Norm y))] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(x, (y, Norm y))]
dataxy Int -> Int -> Int
forall a. Num a => a -> a -> a
* SubBasis y -> Int
forall v. FiniteDimensional v => SubBasis v -> Int
subbasisDimension (SubBasis y
forall v. FiniteDimensional v => SubBasis v
entireBasis :: SubBasis y)
                  Int -> Int -> Int
forall a. Num a => a -> a -> a
- SubBasis m -> Int
forall v. FiniteDimensional v => SubBasis v -> Int
subbasisDimension (SubBasis m
forall v. FiniteDimensional v => SubBasis v
entireBasis :: SubBasis m)
              forward :: m -> [y]
              forward :: m -> [y]
forward m
m = [x -> m +> y
modelMap x
x LinearMap s m y -> m -> y
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ m
m | (x
x,(y, Norm y)
_)<-[(x, (y, Norm y))]
dataxy]
              forward' :: [DualVector y] -> DualVector m
              forward' :: [DualVector y] -> DualVector m
forward' = [DualVector m] -> DualVector m
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV ([DualVector m] -> DualVector m)
-> ([DualVector y] -> [DualVector m])
-> [DualVector y]
-> DualVector m
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. ((LinearMap s m y, LinearMap s (DualVector y) (DualVector m))
 -> DualVector y -> DualVector m)
-> [(LinearMap s m y, LinearMap s (DualVector y) (DualVector m))]
-> [DualVector y]
-> [DualVector m]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (LinearMap s (DualVector y) (DualVector m)
-> DualVector y -> DualVector m
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
($) (LinearMap s (DualVector y) (DualVector m)
 -> DualVector y -> DualVector m)
-> ((LinearMap s m y, LinearMap s (DualVector y) (DualVector m))
    -> LinearMap s (DualVector y) (DualVector m))
-> (LinearMap s m y, LinearMap s (DualVector y) (DualVector m))
-> DualVector y
-> DualVector m
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (LinearMap s m y, LinearMap s (DualVector y) (DualVector m))
-> LinearMap s (DualVector y) (DualVector m)
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd) [(LinearMap s m y, LinearMap s (DualVector y) (DualVector m))]
[(m +> y, DualVector y +> DualVector m)]
modelGens
              modelGens :: [(m +> y, DualVector y +> DualVector m)]
              modelGens :: [(m +> y, DualVector y +> DualVector m)]
modelGens = ((LinearMap (Scalar m) m (DualVector (DualVector y))
-> LinearMap (Scalar m) m (DualVector (DualVector y))
forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id(LinearMap (Scalar m) m (DualVector (DualVector y))
 -> LinearMap (Scalar m) m (DualVector (DualVector y)))
-> (LinearMap (Scalar m) m (DualVector (DualVector y))
    -> DualVector y +> DualVector m)
-> LinearMap (Scalar m) m (DualVector (DualVector y))
-> (LinearMap (Scalar m) m (DualVector (DualVector y)),
    DualVector y +> DualVector m)
forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&LinearFunction
  (Scalar (DualVector y))
  (LinearMap (Scalar m) m (DualVector (DualVector y)))
  (DualVector y +> DualVector m)
-> LinearMap (Scalar m) m (DualVector (DualVector y))
-> DualVector y +> DualVector m
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr LinearFunction
  (Scalar (DualVector y))
  (LinearMap (Scalar m) m (DualVector (DualVector y)))
  (DualVector y +> DualVector m)
forall v w.
(LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) =>
(v +> DualVector w) -+> (w +> DualVector v)
adjoint) (LinearMap (Scalar m) m (DualVector (DualVector y))
 -> (LinearMap (Scalar m) m (DualVector (DualVector y)),
     DualVector y +> DualVector m))
-> ((x, (y, Norm y))
    -> LinearMap (Scalar m) m (DualVector (DualVector y)))
-> (x, (y, Norm y))
-> (LinearMap (Scalar m) m (DualVector (DualVector y)),
    DualVector y +> DualVector m)
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. x -> m +> y
x -> LinearMap (Scalar m) m (DualVector (DualVector y))
modelMap (x -> LinearMap (Scalar m) m (DualVector (DualVector y)))
-> ((x, (y, Norm y)) -> x)
-> (x, (y, Norm y))
-> LinearMap (Scalar m) m (DualVector (DualVector y))
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (x, (y, Norm y)) -> x
forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst)((x, (y, Norm y))
 -> (LinearMap (Scalar m) m (DualVector (DualVector y)),
     DualVector y +> DualVector m))
-> [(x, (y, Norm y))]
-> [(LinearMap (Scalar m) m (DualVector (DualVector y)),
     DualVector y +> DualVector m)]
forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[(x, (y, Norm y))]
dataxy
              σm :: Norm m
              σm :: Norm m
σm = [Norm m] -> Norm m
forall a. Monoid a => [a] -> a
mconcat [ LinearFunction (Scalar (DualVector m)) m (DualVector m) -> Norm m
forall v. (v -+> DualVector v) -> Norm v
Norm (LinearFunction (Scalar (DualVector m)) m (DualVector m) -> Norm m)
-> (LinearMap s m (DualVector m)
    -> LinearFunction (Scalar (DualVector m)) m (DualVector m))
-> LinearMap s m (DualVector m)
-> Norm m
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. LinearMap s m (DualVector m)
-> LinearFunction (Scalar (DualVector m)) m (DualVector m)
forall (a :: * -> * -> *) (k :: * -> * -> *) b c.
(EnhancedCat a k, Object k b, Object k c, Object a b,
 Object a c) =>
k b c -> a b c
arr (LinearMap s m (DualVector m) -> Norm m)
-> LinearMap s m (DualVector m) -> Norm m
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearMap s (DualVector y) (DualVector m)
m LinearMap s (DualVector y) (DualVector m)
-> LinearMap s m (DualVector y) -> LinearMap s m (DualVector m)
forall κ (k :: κ -> κ -> *) (a :: κ) (b :: κ) (c :: κ).
(Category k, Object k a, Object k b, Object k c) =>
k b c -> k a b -> k a c
. (LinearFunction s y (DualVector y)
-> LinearFunction
     s (LinearMap s m y) (LinearMap s m (DualVector y))
forall (f :: * -> *) (r :: * -> * -> *) (t :: * -> * -> *) a b.
(Functor f r t, Object r a, Object t (f a), Object r b,
 Object t (f b)) =>
r a b -> t (f a) (f b)
fmap LinearFunction s y (DualVector y)
y -+> DualVector y
ny LinearFunction s (LinearMap s m y) (LinearMap s m (DualVector y))
-> LinearMap s m y -> LinearMap s m (DualVector y)
forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearMap s m y
m')
                           | ((x
_,(y
_,Norm y -+> DualVector y
ny)), (LinearMap s m y
m',LinearMap s (DualVector y) (DualVector m)
m)) <- [(x, (y, Norm y))]
-> [(LinearMap s m y, LinearMap s (DualVector y) (DualVector m))]
-> [((x, (y, Norm y)),
     (LinearMap s m y, LinearMap s (DualVector y) (DualVector m)))]
forall a b. [a] -> [b] -> [(a, b)]
zip [(x, (y, Norm y))]
dataxy [(LinearMap s m y, LinearMap s (DualVector y) (DualVector m))]
[(m +> y, DualVector y +> DualVector m)]
modelGens ]
                  
-- $classesExplanation
-- To facilitate this library's abstract, base-less approach, the
-- spaces/types need to be instances of multiple typeclasses which
-- are rather tedious to write instances for. Due to the use of type
-- families, GHC's various deriving mechanisms cannot help much with this
-- either.
--
-- In the common case that your spaces /do/ have a canonical basis, the
-- instances can be generated automatically with one of the Template Haskell
-- macros 'makeLinearSpaceFromBasis' or 'makeFiniteDimensionalFromBasis'.