-- |
-- 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 TypeApplications     #-}
{-# 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
            , DimensionAware (..)
            , Dimensional (..)
            , StaticDimensional (..)
            , Dimension
            , 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, sRiesz, coRiesz, showsPrecAsRiesz, (.<)
            -- ** Standard decompositions
            , TensorDecomposable(..), RieszDecomposable(..)
            , tensorDecomposeShowsPrec, rieszDecomposeShowsPrec
            -- ** Constraint synonyms
            , HilbertSpace, SimpleSpace, RealSpace
            , Num'(..)
            , Fractional'
            , RealFrac', RealFloat', LinearShowable
            -- ** Coercions
            , VSCCoercion(..)
            -- ** 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.DimensionAware
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-+|> :: forall (f :: * -> * -> *) s u v.
(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
-+|>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 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
. forall s v w. (v -> w) -> LinearFunction s v w
LinearFunction forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ (v
vforall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> 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
duforall 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 :: forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \[DualVector v]
dvs -> forall v. (v -+> DualVector v) -> Norm v
Norm 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
. forall s v w. (v -> w) -> LinearFunction s v w
LinearFunction 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 [DualVector v
dv forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^* (DualVector v
dvforall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v) | DualVector v
dv <- [DualVector v]
dvs]

spanVariance :: ∀ v . LSpace v => [v] -> Variance v
spanVariance :: forall v. LSpace v => [v] -> Variance v
spanVariance = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness -> 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 :: forall v. SimpleSpace v => Norm v -> [v] -> Norm v
relaxNorm = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \Norm v
me [v]
vs -> let vs' :: [v]
vs' = forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem' Norm v
me
                     in forall v. SimpleSpace v => Norm v -> Variance v
dualNorm 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
. forall v. LSpace v => [v] -> Variance v
spanVariance forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [v]
vs' 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 :: forall v. LSpace v => Scalar v -> Norm v -> Norm v
scaleNorm = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness -> \Scalar v
μ (Norm v -+> DualVector v
n) -> forall v. (v -+> DualVector v) -> Norm v
Norm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Scalar v
μforall a. Num a => a -> Int -> a
^Int
2 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 {
    forall v. 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 ∀ v . LSpace v => Semigroup (Norm v) where
  Norm v -+> DualVector v
m <> :: Norm v -> Norm v -> Norm v
<> Norm v -+> DualVector v
n = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness @v of
    DualSpaceWitness v
DualSpaceWitness -> forall v. (v -+> DualVector v) -> Norm v
Norm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v -+> DualVector v
mforall v. AdditiveGroup v => v -> v -> v
^+^v -+> DualVector v
n
-- | @mempty|$|v ≡ 0@
instance ∀ v . LSpace v => Monoid (Seminorm v) where
  mempty :: Seminorm v
mempty = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness @v of
    DualSpaceWitness v
DualSpaceWitness -> forall v. (v -+> DualVector v) -> Norm v
Norm forall v. AdditiveGroup v => v
zeroV
  mappend :: Seminorm v -> Seminorm v -> Seminorm v
mappend = 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 :: forall v. HilbertSpace v => Norm v
euclideanNorm = forall v. (v -+> DualVector v) -> Norm v
Norm 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 :: forall v. FiniteDimensional v => Norm v
adhocNorm = forall v. (v -+> DualVector v) -> Norm v
Norm 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 :: forall v. SimpleSpace v => Norm v -> Variance v
dualNorm = forall v. LSpace v => [v] -> Variance v
spanVariance 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
. 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' :: forall v. SimpleSpace v => Variance v -> Norm v
dualNorm' = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
     DualSpaceWitness v
DualSpaceWitness -> forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm 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
. 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 :: forall v w.
(LSpace v, LSpace w, Scalar v ~ Scalar w) =>
(v +> w) -> Norm w -> Norm v
transformNorm = case ( forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v
                     , forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness w ) of
    (DualSpaceWitness v
DualSpaceWitness, DualSpaceWitness w
DualSpaceWitness)
            -> \v +> w
f (Norm w -+> DualVector w
m) -> forall v. (v -+> DualVector v) -> Norm v
Norm 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
. 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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ (forall v w.
(LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) =>
(v +> DualVector w) -+> (w +> DualVector v)
adjoint forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
f) 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
. (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 w -+> DualVector w
m forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
f)

transformVariance :: ∀ v w . (LSpace v, LSpace w, Scalar v~Scalar w)
                        => (v+>w) -> Variance v -> Variance w
transformVariance :: forall v w.
(LSpace v, LSpace w, Scalar v ~ Scalar w) =>
(v +> w) -> Variance v -> Variance w
transformVariance = case ( forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v
                     , 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) -> forall v. (v -+> DualVector v) -> Norm v
Norm 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
. 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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
f 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
. (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 DualVector v -+> DualVector (DualVector v)
m forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall v w.
(LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) =>
(v +> DualVector w) -+> (w +> DualVector v)
adjoint forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v +> w
f)

infixl 6 ^%
(^%) :: (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
v
v ^% :: forall v. (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
^% Norm v -+> DualVector v
m = v
v forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ forall a. Floating a => a -> a
sqrt ((v -+> DualVector v
mforall s v w. LinearFunction s v w -> v -> w
-+$>v
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 :: forall s. RealFrac' s => Norm s -> Maybe s
findNormalLength (Norm s -+> DualVector s
m) = case ( forall s. Num' s => ClosedScalarWitness s
closedScalarWitness :: ClosedScalarWitness s
                                 , s -+> DualVector s
mforall s v w. LinearFunction s v w -> v -> w
-+$>s
1 ) of
   (ClosedScalarWitness s
ClosedScalarWitness, DualVector s
o) | DualVector s
o forall a. Ord a => a -> a -> Bool
> DualVector s
0      -> forall a. a -> Maybe a
Just 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
. forall a. Floating a => a -> a
sqrt forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. Fractional a => a -> a
recip DualVector s
o
                            | Bool
otherwise  -> 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 :: forall s. RealFrac' s => Norm s -> s
normalLength (Norm s -+> DualVector s
m) = case ( forall s. Num' s => ClosedScalarWitness s
closedScalarWitness :: ClosedScalarWitness s
                             , s -+> DualVector s
mforall s v w. LinearFunction s v w -> v -> w
-+$>s
1 ) of
   (ClosedScalarWitness s
ClosedScalarWitness, DualVector s
o) | DualVector s
o forall a. Ord a => a -> a -> Bool
>= DualVector s
0     -> forall a. Floating a => a -> a
sqrt forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. Fractional a => a -> a
recip DualVector s
o
                            | DualVector s
o forall a. Ord a => a -> a -> Bool
< DualVector s
0      -> forall a. HasCallStack => [Char] -> a
error [Char]
"Norm fails to be positive semidefinite."
                            | Bool
otherwise  -> 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 <$| :: forall v. LSpace v => Norm v -> v -> DualVector v
<$| v
v = v -+> DualVector v
mforall 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 :: forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq (Norm v -+> DualVector v
m) v
v = (v -+> DualVector v
mforall s v w. LinearFunction s v w -> v -> w
-+$>v
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
|$| :: forall v.
(LSpace v, Floating (Scalar v)) =>
Seminorm v -> v -> Scalar v
(|$|) Seminorm v
m = forall a. Floating a => a -> a
sqrt 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
. forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Seminorm v
m

infixl 1 |&>
-- | Flipped, “ket” version of '<$|'.
-- 
-- @
-- v '<.>^' (w |&> 'euclideanNorm')  ≡  v '<.>' w
-- @
(|&>) :: ∀ v . LSpace v => DualVector v -> Variance v -> v
DualVector v
dv |&> :: forall v. LSpace v => DualVector v -> Variance v -> v
|&> Norm DualVector v -+> DualVector (DualVector v)
m = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness @v of
   DualSpaceWitness v
DualSpaceWitness -> DualVector v -+> DualVector (DualVector v)
mforall 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 :: forall v. LSpace v => Norm v -> Norm v
densifyNorm = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \(Norm LinearFunction (Scalar (DualVector v)) v (DualVector v)
m) -> forall v. (v -+> DualVector v) -> Norm v
Norm 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
. 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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall v w.
(LinearSpace v, TensorSpace w, Scalar v ~ Scalar w) =>
(v -+> w) -+> (v +> w)
sampleLinearFunction forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearFunction (Scalar (DualVector v)) 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 :: forall v. LinearSpace v => Norm v -> Maybe (Norm v)
wellDefinedNorm = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
    DualSpaceWitness v
DualSpaceWitness
        -> \(Norm v -+> DualVector v
m) -> forall v. (v -+> DualVector v) -> Norm v
Norm forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall v. TensorSpace v => v -> Maybe v
wellDefinedVector v -+> DualVector v
m

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

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

orthogonalComplementProj' :: LSpace v => [(v, DualVector v)] -> (v-+>v)
orthogonalComplementProj' :: forall v. LSpace v => [(v, DualVector v)] -> v -+> v
orthogonalComplementProj' [(v, DualVector v)]
ws = forall s v w. (v -> w) -> LinearFunction s v w
LinearFunction forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ \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 forall v. AdditiveGroup v => v -> v -> v
^-^ v
wforall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*(DualVector v
dwforall 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 :: forall v. LSpace v => Norm v -> [v] -> v -+> v
orthogonalComplementProj (Norm v -+> DualVector v
m)
      = forall v. LSpace v => [(v, DualVector v)] -> v -+> v
orthogonalComplementProj' 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
. forall a b. (a -> b) -> [a] -> [b]
map (forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id 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
mforall 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 {
      forall v. Eigenvector v -> Scalar v
ev_Eigenvalue :: Scalar v -- ^ The estimated eigenvalue @λ@.
    , forall v. Eigenvector v -> v
ev_Eigenvector :: v       -- ^ Normalised vector @v@ that gets mapped to a multiple, namely:
    , forall v. Eigenvector v -> v
ev_FunctionApplied :: v   -- ^ @f $ v ≡ λ *^ v @.
    , forall v. Eigenvector v -> v
ev_Deviation :: v         -- ^ Deviation of @v@ to @(f$v)^/λ@. Ideally, this would of course be equal.
    , forall v. 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 :: forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> Scalar v -> (v -+> v) -> [v] -> [[Eigenvector v]]
constructEigenSystem Norm v
me Scalar v
ε₀ v -+> v
f = forall a. (a -> a) -> a -> [a]
iterate (
                                             forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$
                                               forall a. Num a => a -> a
negate 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
. forall a. Num a => a -> a
abs 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
. forall v. Eigenvector v -> Scalar v
ev_Eigenvalue)
                                           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
. forall a b. (a -> b) -> [a] -> [b]
map v -> Eigenvector v
asEV
                                           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
. forall v.
(LSpace v, RealFloat (Scalar v)) =>
Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily (Scalar v
1forall a. Fractional a => a -> a -> a
/Scalar v
4) Norm v
me
                                           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)
                                         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
. forall a b. (a -> b) -> [a] -> [b]
map (v -> Eigenvector v
asEV 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
. (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
εforall a. Ord a => a -> a -> Bool
>Scalar v
ε₀       = case [Eigenvector v] -> [v]
newSys [Eigenvector v]
evs of
                         []     -> [v
fvforall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v
λ, v
dvforall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/forall a. Floating a => a -> a
sqrt(Scalar v
εforall a. Num a => a -> a -> a
+Scalar v
ε₀)]
                         v
vn:[v]
vns -> v
fvforall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v
λ forall a. a -> [a] -> [a]
: v
vn forall a. a -> [a] -> [a]
: v
dvforall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/forall a. Floating a => a -> a
sqrt(Scalar v
εforall a. Num a => a -> a -> a
+Scalar v
ε₀) forall a. a -> [a] -> [a]
: [v]
vns
         | Scalar v
εforall a. Ord a => a -> a -> Bool
>=Scalar v
0       = 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 = 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'forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
fv
              λ :: Scalar v
λ = DualVector v
fv'forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^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 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
meforall v. LSpace v => Norm v -> v -> DualVector v
<$|v
fv
              dv :: v
dv | Scalar v
λ²forall a. Ord a => a -> a -> Bool
>Scalar v
0       = v
v forall v. AdditiveGroup v => v -> v -> v
^-^ v
fvforall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*(Scalar v
λforall a. Fractional a => a -> a -> a
/Scalar v
λ²) -- for stability reasons
                 | Bool
otherwise  = forall v. AdditiveGroup v => v
zeroV


finishEigenSystem :: ∀ v . (LSpace v, RealFloat (Scalar v))
                      => Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem :: forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem Norm v
me = [Eigenvector v] -> [Eigenvector v]
go 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
. forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. Num a => a -> a
negate 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
. 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
λ₀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₀forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₀ forall v. AdditiveGroup v => v -> v -> v
^+^ v
v₁forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₁
              fv₀' :: v
fv₀' = v
fv₀forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₀ forall v. AdditiveGroup v => v -> v -> v
^+^ v
fv₁forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₀₁
              
              v₁' :: v
v₁' = v
v₀forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₀ forall v. AdditiveGroup v => v -> v -> v
^+^ v
v₁forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₁
              fv₁' :: v
fv₁' = v
fv₀forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₀ forall v. AdditiveGroup v => v -> v -> v
^+^ v
fv₁forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μ₁₁
              
              fShift₁v₀ :: v
fShift₁v₀ = v
fv₀ forall v. AdditiveGroup v => v -> v -> v
^-^ Scalar v
λ₁forall v. VectorSpace v => Scalar v -> v -> v
*^v
v₀
              
              (Scalar v
μ₀₀,Scalar v
μ₀₁) = forall {b}. Floating b => (b, b) -> (b, b)
normalised ( Scalar v
λ₀ forall a. Num a => a -> a -> a
- Scalar v
λ₁
                                     , (Norm v
me forall v. LSpace v => Norm v -> v -> DualVector v
<$| v
fShift₁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'' forall a. [a] -> [a] -> [a]
++ [Eigenvector v]
upper'
        where l :: Int
l = forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs
              lChunk :: Int
lChunk = Int
lforall a. Integral a => a -> a -> a
`quot`Int
3
              ([Eigenvector v]
loEvs, ([Eigenvector v]
midEvs, [Eigenvector v]
hiEvs)) = forall (a :: * -> * -> *) d b c.
(Morphism a, ObjectPair a d b, ObjectPair a d c) =>
a b c -> a (d, b) (d, c)
second (forall a. Int -> [a] -> ([a], [a])
splitAt forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Int
l forall a. Num a => a -> a -> a
- Int
2forall a. Num a => a -> a -> a
*Int
lChunk)
                                                    forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. Int -> [a] -> ([a], [a])
splitAt Int
lChunk [Eigenvector v]
evs
              ([Eigenvector v]
lo',[Eigenvector v]
hi') = forall a. Int -> [a] -> ([a], [a])
splitAt Int
lChunk 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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [Eigenvector v]
loEvsforall a. [a] -> [a] -> [a]
++[Eigenvector v]
hiEvs
              ([Eigenvector v]
lo'',[Eigenvector v]
mid') = forall a. Int -> [a] -> ([a], [a])
splitAt Int
lChunk 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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [Eigenvector v]
lo'forall a. [a] -> [a] -> [a]
++[Eigenvector v]
midEvs
              upper' :: [Eigenvector v]
upper'  = [Eigenvector v] -> [Eigenvector v]
go forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [Eigenvector v]
mid'forall a. [a] -> [a] -> [a]
++[Eigenvector v]
hi'
       
       asEV :: v -> v -> Eigenvector v
asEV v
v v
fv = 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
meforall v. LSpace v => Norm v -> v -> DualVector v
<$|v
v)forall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
fv
              dv :: v
dv = v
vforall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
λ forall v. AdditiveGroup v => v -> v -> v
^-^ v
fv
              ε :: Scalar v
ε = forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm v
me v
dv forall a. Fractional a => a -> a -> a
/ Scalar v
λforall a. Num a => a -> Int -> a
^Int
2
       
       normalised :: (b, b) -> (b, b)
normalised (b
x,b
y) = (b
xforall a. Fractional a => a -> a -> a
/b
r, b
yforall a. Fractional a => a -> a -> a
/b
r)
        where r :: b
r = forall a. Floating a => a -> a
sqrt forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ b
xforall a. Num a => a -> Int -> a
^Int
2 forall a. Num a => a -> a -> a
+ b
yforall 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 :: forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
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)
         | forall v. LSpace v => Seminorm v -> v -> Scalar v
normSq Norm v
me v
vPerp forall a. Ord a => a -> a -> Bool
> Scalar v
fpε  = case [[Eigenvector v]]
evss of
             [Eigenvector v]
evs':[[Eigenvector v]]
_ | forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs' forall a. Ord a => a -> a -> Bool
> Int
oldDim
               -> [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go (v
vforall a. a -> [a] -> [a]
:[v]
vs) (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Eigenvector v]
evs) [[Eigenvector v]]
evss
             [[Eigenvector v]]
_ -> let evss' :: [[Eigenvector v]]
evss' = forall a. [a] -> [a]
tail 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
. forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> Scalar v -> (v -+> v) -> [v] -> [[Eigenvector v]]
constructEigenSystem Norm v
me Scalar v
fpε (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)
                                forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall v. Eigenvector v -> v
ev_Eigenvector (forall a. [a] -> a
head forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [[Eigenvector v]]
evssforall a. [a] -> [a] -> [a]
++[[Eigenvector v]
evs]) forall a. [a] -> [a] -> [a]
++ [v
vPerp]
                  in [v] -> Int -> [[Eigenvector v]] -> [Eigenvector v]
go [v]
vs (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]
evsforall a. a -> [a] -> [a]
:[[Eigenvector v]]
evss)
        where vPerp :: v
vPerp = forall v. LSpace v => Norm v -> [v] -> v -+> v
orthogonalComplementProj Norm v
me (forall v. Eigenvector v -> v
ev_Eigenvectorforall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[Eigenvector v]
evs) forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ v
v
       fBas :: [v]
fBas = (forall v. (LSpace v, Floating (Scalar v)) => v -> Norm v -> v
^%Norm v
me) forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd (forall v w.
(FiniteDimensional v, LSpace w, Scalar w ~ Scalar v) =>
(v +> w) -> (SubBasis v, DList w)
decomposeLinMap forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
id) []
       fpε :: Scalar v
fpε = forall a. IEEE a => a
epsilon 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 :: forall v.
(FiniteDimensional v, HilbertSpace v, IEEE (Scalar v)) =>
(v +> v) -> [(Scalar v, v)]
eigen v +> v
f = forall a b. (a -> b) -> [a] -> [b]
map (forall v. Eigenvector v -> Scalar v
ev_Eigenvalue forall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&& forall v. Eigenvector v -> v
ev_Eigenvector)
   forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem forall v. HilbertSpace v => Norm v
euclideanNorm) (forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem forall v. HilbertSpace v => Norm v
euclideanNorm v +> v
f) forall a. [a] -> Int -> a
!! Int
2


-- | Approximation of the determinant.
roughDet :: (FiniteDimensional v, IEEE (Scalar v)) => (v+>v) -> Scalar v
roughDet :: forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
(v +> v) -> Scalar v
roughDet = forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem forall v. FiniteDimensional v => Norm v
adhocNorm 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
>>> forall a b. (a -> b) -> [a] -> [b]
map forall v. Eigenvector v -> Scalar v
ev_Eigenvalue 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
>>> forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product


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


normSpanningSystem :: SimpleSpace v
               => Seminorm v -> [DualVector v]
normSpanningSystem :: forall v. SimpleSpace v => Seminorm v -> [DualVector v]
normSpanningSystem = forall a b. (a -> b) -> [a] -> [b]
map forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd 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
. forall v. SimpleSpace v => Seminorm v -> [(v, DualVector v)]
normSpanningSystems

normSpanningSystems :: SimpleSpace v
               => Seminorm v -> [(v, DualVector v)]
normSpanningSystems :: forall v. SimpleSpace v => Seminorm v -> [(v, DualVector v)]
normSpanningSystems me :: Seminorm v
me@(Norm v -+> DualVector v
m)
     = forall a. [Maybe a] -> [a]
catMaybes 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
. forall a b. (a -> b) -> [a] -> [b]
map (\(v
v,Maybe (DualVector v)
d)->(v
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) 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
. forall v.
(SemiInner v, RealFrac' (Scalar v)) =>
Scalar v -> [(v, DualVector v)] -> [(v, Maybe (DualVector v))]
orthonormaliseDuals Scalar v
0
         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
. forall a b. (a -> b) -> [a] -> [b]
map (forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
idforall (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
mforall s v w. LinearFunction s v w -> v -> w
-+$>)) forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem' Seminorm v
me

normSpanningSystem' :: (FiniteDimensional v, IEEE (Scalar v))
               => Seminorm v -> [v]
normSpanningSystem' :: forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Seminorm v -> [v]
normSpanningSystem' Seminorm v
me = forall v.
(LSpace v, RealFloat (Scalar v)) =>
Scalar v -> Norm v -> [v] -> [v]
orthonormaliseFussily Scalar v
0 Seminorm v
me forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall v. FiniteDimensional v => SubBasis v -> [v]
enumerateSubBasis 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 :: forall v. SimpleSpace v => Variance v -> [v]
varianceSpanningSystem = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
                           DualSpaceWitness v
DualSpaceWitness -> 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 :: forall v.
SimpleSpace v =>
Norm v -> Norm v -> [(DualVector v, Scalar v)]
sharedNormSpanningSystem nn :: Norm v
nn@(Norm v -+> DualVector v
n) Norm v
nm
      = 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
nforall s v w. LinearFunction s v w -> v -> w
-+$>) forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' Int
0 (Norm v
nn, 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' :: forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' = DualSpaceWitness v
-> Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
snss 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)
           = forall {a}.
(Ord (Scalar a), Floating (Scalar a)) =>
Eigenvector a -> [(a, Scalar a)]
sep 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)
=<< forall a. (a -> a) -> a -> [a]
iterate (forall v.
(LSpace v, RealFloat (Scalar v)) =>
Norm v -> [Eigenvector v] -> [Eigenvector v]
finishEigenSystem Norm v
nn)
                        (forall v.
(FiniteDimensional v, IEEE (Scalar v)) =>
Norm v -> (v +> v) -> [Eigenvector v]
roughEigenSystem Norm v
nn forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ 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 DualVector v -+> DualVector (DualVector v)
n' 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
. 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 -+> DualVector v
m) 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
λforall a. Ord a => a -> a -> Bool
>=Scalar a
0       = [(a
v, 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 :: forall v.
SimpleSpace v =>
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 forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness
               forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' Int
1 (Seminorm v
combined, forall v. SimpleSpace v => Norm v -> Variance v
dualNorm Seminorm v
combined) Seminorm v
nn
 where combined :: Seminorm v
combined = forall v. LSpace v => Norm v -> Norm v
densifyNorm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Seminorm v
nnforall 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
μnforall a. Num a => a -> Int -> a
^Int
2 forall a. Ord a => a -> a -> Bool
> forall a. IEEE a => a
epsilon  = (DualVector v
v'forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*Scalar v
μn, forall a. a -> Maybe a
Just forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall a. Floating a => a -> a
sqrt (forall a. Ord a => a -> a -> a
max Scalar v
0 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Scalar v
1 forall a. Num a => a -> a -> a
- Scalar v
μnforall a. Num a => a -> Int -> a
^Int
2)forall a. Fractional a => a -> a -> a
/Scalar v
μn)
           | Bool
otherwise       = (DualVector v
v', forall a. Maybe a
Nothing)
        where v' :: DualVector v
v' = Seminorm v
combinedforall 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' :: forall v. SimpleSpace v => Seminorm v -> Seminorm v -> [v]
sharedSeminormSpanningSystem' Seminorm v
nn Seminorm v
nm
         = forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst forall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$> forall v.
SimpleSpace v =>
Int -> (Norm v, Variance v) -> Norm v -> [(v, Scalar v)]
sharedNormSpanningSystem' Int
1 (Seminorm v
combined, forall v. SimpleSpace v => Norm v -> Variance v
dualNorm Seminorm v
combined) Seminorm v
nn
 where combined :: Seminorm v
combined = forall v. LSpace v => Norm v -> Norm v
densifyNorm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ Seminorm v
nnforall 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 :: forall u v.
(SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) =>
Variance (u, v) -> u +> v
dependence = case ( forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness u
                  , 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) -> 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 ( forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd 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, v) -+> DualVector (DualVector (u, v))
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
. (forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
idforall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&forall v. AdditiveGroup v => v
zeroV) )
              forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall u v.
(SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) =>
(u +> v) -> v +> u
pseudoInverse (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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst 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, v) -+> DualVector (DualVector (u, v))
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
. (forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
idforall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&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 :: forall u v.
(SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) =>
Norm (u, v) -> (Norm u, Norm v)
summandSpaceNorms = case ( forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness u
                         , 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 = forall v. SimpleSpace v => Seminorm v -> [DualVector v]
normSpanningSystem Norm (u, v)
nuv
                   in ( forall v. LSpace v => Norm v -> Norm v
densifyNorm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm (forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fstforall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[DualVector (u, v)]
spanSys)
                      , forall v. LSpace v => Norm v -> Norm v
densifyNorm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm (forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
sndforall (f :: * -> *) (r :: * -> * -> *) a b.
(Functor f r (->), Object r a, Object r b) =>
r a b -> f a -> f b
<$>[DualVector (u, v)]
spanSys) )

sumSubspaceNorms :: ∀ u v . (LSpace u, LSpace v, Scalar u~Scalar v)
                         => Norm u -> Norm v -> Norm (u,v)
sumSubspaceNorms :: forall u v.
(LSpace u, LSpace v, Scalar u ~ Scalar v) =>
Norm u -> Norm v -> Norm (u, v)
sumSubspaceNorms = case ( forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness u
                         , 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) -> forall v. (v -+> DualVector v) -> Norm v
Norm forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ u -+> DualVector u
nu 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')
*** 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
pforall a. Ord a => a -> a -> Bool
>Int
9) forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ ([Char]
"spanNorm "forall a. [a] -> [a] -> [a]
++) 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
. forall a. Show a => a -> ShowS
shows (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 :: forall v. SimpleSpace v => [v] -> [DualVector v]
convexPolytopeHull [v]
vs = case forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness :: DualSpaceWitness v of
         DualSpaceWitness v
DualSpaceWitness
             -> [DualVector v
dvforall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/Scalar v
η | (DualVector v
dv,Scalar v
η) <- [(DualVector v, Scalar v)]
candidates, forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((forall a. Ord a => a -> a -> Bool
<=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
. (DualVector v
dvforall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^)) [v]
vs]
 where vrv :: Variance v
vrv = forall v. LSpace v => [v] -> Variance v
spanVariance [v]
vs
       nmv :: Norm v
nmv = 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
dvforall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v) | v
v <- [v]
vs
                                   , let dv :: DualVector v
dv = Norm v
nmvforall v. LSpace v => Norm v -> v -> DualVector v
<$|v
v ]

symmetricPolytopeOuterVertices :: ∀ v . SimpleSpace v => [DualVector v] -> [v]
symmetricPolytopeOuterVertices :: forall v. SimpleSpace v => [DualVector v] -> [v]
symmetricPolytopeOuterVertices [DualVector v]
dvs
         = [ forall {v}.
(LinearSpace v, Fractional (Scalar v)) =>
v -> [(DualVector v, v)] -> v
seekExtreme 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 = forall v. LSpace v => [DualVector v] -> Seminorm v
spanNorm [DualVector v]
dvs
       vrv :: Variance v
vrv = 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
dvforall v. LSpace v => DualVector v -> Variance v -> v
|&>Variance v
vrv ]
       ([[(DualVector v, v)]]
candidates, [(DualVector v, v)]
_) = forall a. Int -> Int -> [a] -> ([[a]], [a])
multiSplit Int
d (Int
2forall a. Num a => a -> a -> a
*Int
d) 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
. forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat 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
. forall a.
SimpleSpace a =>
[(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [(DualVector v, v)]
withSomeVect
       d :: Int
d = forall v. FiniteDimensional v => SubBasis v -> Int
subbasisDimension (forall v. FiniteDimensional v => SubBasis v
entireBasis :: SubBasis v)
       seekExtreme :: v -> [(DualVector v, v)] -> v
seekExtreme v
p₀ [] = v
p₀
       seekExtreme v
p₀ ((DualVector v
dv, v
v) : [(DualVector v, v)]
cs)
           = v -> [(DualVector v, v)] -> v
seekExtreme (v
p₀forall v. AdditiveGroup v => v -> v -> v
^+^v
vn) [(DualVector v
dw, v
w forall v. AdditiveGroup v => v -> v -> v
^-^ v
vforall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^*((DualVector v
dvforall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
w) forall a. Fractional a => a -> a -> a
/ Scalar v
lv)) | (DualVector v
dw, v
w) <- [(DualVector v, v)]
cs]
        where vn :: v
vn = v
v forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^* ((Scalar v
1 forall a. Num a => a -> a -> a
- DualVector v
dvforall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
p₀) forall a. Fractional a => a -> a -> a
/ Scalar v
lv)
              lv :: Scalar v
lv = DualVector v
dvforall v. LinearSpace v => DualVector v -> v -> Scalar v
<.>^v
v

deinterlacions :: SimpleSpace a => [(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions :: forall a.
SimpleSpace a =>
[(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions [(DualVector a, a)]
l = [(DualVector a, a)]
l forall a. a -> [a] -> [a]
: forall a.
SimpleSpace a =>
[(DualVector a, a)] -> [[(DualVector a, a)]]
deinterlacions ([(DualVector a, a)]
e forall a. [a] -> [a] -> [a]
++ forall a b. (a -> b) -> [a] -> [b]
map forall v. AdditiveGroup v => v -> v
negateV [(DualVector a, a)]
o)
 where ([(DualVector a, a)]
e,[(DualVector a, a)]
o) = forall {a}. [a] -> ([a], [a])
deinterlace [(DualVector a, a)]
l
       deinterlace :: [a] -> ([a], [a])
deinterlace (a
a:a
b:[a]
xs) = (a
aforall 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
bforall 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 :: forall 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
σy x -> m +> y
modelMap = forall x y m. LinearRegressionResult x y m -> m
linearFit_bestModel
                                   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
. 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 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
. forall a b. (a -> b) -> [a] -> [b]
map (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 {
          forall x y m. 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.
        , forall x y m. 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²@.
        , forall x y m. 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 :: forall 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 = case Bool
True of Bool
False -> 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 :: 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 = (DualSpaceWitness y, DualSpaceWitness m)
-> (x -> m +> y)
-> [(x, (y, Norm y))]
-> LinearRegressionResult x y m
lrw (forall v. LinearSpace v => DualSpaceWitness v
dualSpaceWitness, 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
         = forall x y m.
Scalar m -> m -> Norm m -> LinearRegressionResult x y m
LinearRegressionResult (s
χ²forall a. Fractional a => a -> a -> a
/forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ν) m
leastSquareSol Norm m
σm
        where leastSquareSol :: m
leastSquareSol = (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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ [DualVector y] -> DualVector m
forward' 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
. forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall v. LSpace v => Norm 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
. forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd 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
. forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd) [(x, (y, Norm y))]
dataxy
                                          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)
                                 forall u v.
(SimpleSpace u, SimpleSpace v, Scalar u ~ Scalar v) =>
(u +> v) -> v -> u
\$ [DualVector y] -> DualVector m
forward' [Norm y
σyforall v. LSpace v => Norm v -> v -> DualVector v
<$|y
y | (x
_,(y
y,Norm y
σy)) <- [(x, (y, Norm y))]
dataxy]
              χ² :: s
χ² = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [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 forall v. AdditiveGroup v => v -> v -> v
^-^ (x -> m +> y
modelMap x
x forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ m
leastSquareSol) ]
              ν :: Int
ν = forall (t :: * -> *) a. Foldable t => t a -> Int
length [(x, (y, Norm y))]
dataxy forall a. Num a => a -> a -> a
* forall v. FiniteDimensional v => SubBasis v -> Int
subbasisDimension (forall v. FiniteDimensional v => SubBasis v
entireBasis :: SubBasis y)
                  forall a. Num a => a -> a -> a
- forall v. FiniteDimensional v => SubBasis v -> Int
subbasisDimension (forall v. FiniteDimensional v => SubBasis v
entireBasis :: SubBasis m)
              forward :: m -> [y]
              forward :: m -> [y]
forward m
m = [x -> m +> y
modelMap x
x 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' = forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV 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
. forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
($) 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
. forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) y
snd) [(m +> y, DualVector y +> DualVector m)]
modelGens
              modelGens :: [(m +> y, DualVector y +> DualVector m)]
              modelGens :: [(m +> y, DualVector y +> DualVector m)]
modelGens = ((forall κ (k :: κ -> κ -> *) (a :: κ).
(Category k, Object k a) =>
k a a
idforall (a :: * -> * -> *) b c c'.
(PreArrow a, Object a b, ObjectPair a c c') =>
a b c -> a b c' -> a b (c, c')
&&&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 forall v w.
(LinearSpace v, LinearSpace w, Scalar v ~ Scalar w) =>
(v +> DualVector w) -+> (w +> DualVector v)
adjoint) 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
modelMap 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
. forall (a :: * -> * -> *) x y.
(PreArrow a, ObjectPair a x y) =>
a (x, y) x
fst)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 = forall a. Monoid a => [a] -> a
mconcat [ forall v. (v -+> DualVector v) -> Norm v
Norm 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
. 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 forall (f :: * -> * -> *) a b.
(Function f, Object f a, Object f b) =>
f a b -> a -> b
$ LinearMap s (DualVector y) (DualVector m)
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
. (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 y -+> DualVector y
ny 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)) <- forall a b. [a] -> [b] -> [(a, b)]
zip [(x, (y, Norm y))]
dataxy [(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'.