hmatrix-backprop (empty) → 0.1.0.0
raw patch · 8 files changed
+1816/−0 lines, 8 filesdep +ANumdep +backpropdep +basesetup-changed
Dependencies added: ANum, backprop, base, finite-typelits, ghc-typelits-knownnat, ghc-typelits-natnormalise, hedgehog, hmatrix, hmatrix-backprop, hmatrix-vector-sized, microlens, microlens-platform, vector, vector-sized
Files
- CHANGELOG.md +11/−0
- LICENSE +30/−0
- README.md +49/−0
- Setup.hs +2/−0
- hmatrix-backprop.cabal +76/−0
- src/Numeric/LinearAlgebra/Static/Backprop.hs +1212/−0
- test/Nudge.hs +177/−0
- test/Spec.hs +259/−0
+ CHANGELOG.md view
@@ -0,0 +1,11 @@+Changelog+=========++Version 0.1.0.0+---------------++*Feb 10, 2018*++<https://github.com/mstksg/hmatrix-backprop/releases/tag/v0.1.0.0>++* Initial release
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Justin Le (c) 2018++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Justin Le nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,49 @@+hmatrix-backprop+================++[](https://hackage.haskell.org/package/hmatrix-backprop)+[](https://travis-ci.org/mstksg/hmatrix-backprop)++*[hmatrix][]* operations lifted for *[backprop][]*.++[hmatrix]: http://hackage.haskell.org/package/hmatrix+[backprop]: http://hackage.haskell.org/package/backprop++Meant to act as a drop-in replacement to the API of+[Numeric.LinearAlgebra.Static][static]. Just change your imports, and your+functions are automatically backpropagatable. Useful types are all+re-exported.++[static]: https://hackage.haskell.org/package/hmatrix-0.18.2.0/docs/Numeric-LinearAlgebra-Static.html++Formulas for gradients come from the following papers:++* <https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf>+* <http://www.dtic.mil/dtic/tr/fulltext/u2/624426.pdf>+* <http://www.cs.cmu.edu/~zkolter/course/15-884/linalg-review.pdf>+* <https://arxiv.org/abs/1602.07527>++Some functions are not yet implemented! See module documentation for details.+PR's definitely appreciated :)++Tests+-----++Currently numeric tests are implemented as property tests using hedgehog, but+it is possible that the answers might differ from the true values by an amount+undetectable by property tests.++All functions currently are tested except for the higher-order functions.++They are tested by "nudging" components of inputs and checking if the change in+the function outputs match what is expected from the backpropagated gradient.++TODO+----++Apart from the exact API of hmatrix, it'd be nice to have:++1. Statically sized convolutions. Should probably add this to hmatrix instead+ first, though.++
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ hmatrix-backprop.cabal view
@@ -0,0 +1,76 @@+-- This file has been generated from package.yaml by hpack version 0.20.0.+--+-- see: https://github.com/sol/hpack+--+-- hash: 0d60d828b601d4eb24c4978c049e961589bc85fcabfd9c4c71ee0103caf1477e++name: hmatrix-backprop+version: 0.1.0.0+synopsis: hmatrix operations lifted for backprop+description: hmatrix operations lifted for backprop.+ .+ Meant to act as a drop-in replacement to the API of+ Numeric.LinearAlgebra.Static. Just change your imports, and your+ functions are automatically backpropagatable.+ .+ See README on Github at <https://github.com/mstksg/hmatrix-backprop#readme>+category: Math+homepage: https://github.com/mstksg/hmatrix-backprop#readme+bug-reports: https://github.com/mstksg/hmatrix-backprop/issues+author: Justin Le+maintainer: justin@jle.im+copyright: (c) Justin Le 2018+license: BSD3+license-file: LICENSE+build-type: Simple+cabal-version: >= 1.10++extra-source-files:+ CHANGELOG.md+ README.md++source-repository head+ type: git+ location: https://github.com/mstksg/hmatrix-backprop++library+ hs-source-dirs:+ src+ ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wredundant-constraints+ build-depends:+ ANum >=0.2+ , backprop >=0.1.2+ , base >=4.7 && <5+ , ghc-typelits-knownnat+ , ghc-typelits-natnormalise+ , hmatrix >=0.18+ , hmatrix-vector-sized >=0.1+ , microlens+ , vector+ , vector-sized >=0.6+ exposed-modules:+ Numeric.LinearAlgebra.Static.Backprop+ other-modules:+ Paths_hmatrix_backprop+ default-language: Haskell2010++test-suite hmatrix-backprop-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ hs-source-dirs:+ test+ ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ backprop >=0.1.2+ , base >=4.7 && <5+ , finite-typelits+ , hedgehog+ , hmatrix >=0.18+ , hmatrix-backprop+ , hmatrix-vector-sized >=0.1+ , microlens+ , microlens-platform+ , vector-sized >=0.6+ other-modules:+ Nudge+ default-language: Haskell2010
+ src/Numeric/LinearAlgebra/Static/Backprop.hs view
@@ -0,0 +1,1212 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}+{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}++-- |+-- Module : Numeric.LinearAlgebra.Static.Backprop+-- Copyright : (c) Justin Le 2018+-- License : BSD3+--+-- Maintainer : justin@jle.im+-- Stability : experimental+-- Portability : non-portable+--+-- A wrapper over "Numeric.LinearAlgebra.Static" (type-safe vector and+-- matrix operations based on blas/lapack) that allows its operations to+-- work with <https://hackage.haskell.org/package/backprop backprop>.+--+-- In short, these functions are "lifted" to work with 'BVar's.+--+-- Using 'evalBP' will run the original operation:+--+-- @+-- 'evalBP' :: (forall s. 'Reifies' s 'W'. 'BVar' s a -> 'BVar' s b) -> a -> b+-- @+--+-- But using 'gradBP' or 'backprop' will give you the gradient:+--+-- @+-- 'gradBP' :: (forall s. 'Reifies' s 'W'. 'BVar' s a -> 'BVar' s b) -> a -> a+-- @+--+-- These can act as a drop-in replacement to the API of+-- "Numeric.LinearAlgebra.Static". Just change your imports, and your+-- functions are automatically backpropagatable. Useful types are all+-- re-exported.+--+-- Formulas for gradients come from the following papers:+--+-- * https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf+-- * http://www.dtic.mil/dtic/tr/fulltext/u2/624426.pdf+-- * http://www.cs.cmu.edu/~zkolter/course/15-884/linalg-review.pdf+-- * https://arxiv.org/abs/1602.07527+--+-- Some functions are notably unlifted:+--+-- * 'H.svd': I can't find any resources that allow you to backpropagate+-- if the U and V matrices are used! If you find one, let me know, or+-- feel free to submit a PR! Because of this, Currently only a version+-- that exports only the singular values is exported.+-- * 'H.svdTall', 'H.svdFlat': Not sure where to start for these+-- * 'qr': Same story.+-- https://github.com/tensorflow/tensorflow/issues/6504 might yield+-- a clue?+-- * 'H.her': No 'Num' instance for 'H.Her' makes this impossible at+-- the moment with the current backprop API+-- * 'H.exmp': Definitely possible, but I haven't dug deep enough to+-- figure it out yet! There is a description here+-- https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf but it+-- requires some things I am not familiar with yet. Feel free to+-- submit a PR!+-- * 'H.sqrtm': Also likely possible. Maybe try to translate+-- http://people.cs.umass.edu/~smaji/projects/matrix-sqrt/ ? PRs+-- welcomed!+-- * 'H.linSolve': Haven't figured out where to start!+-- * 'H.</>': Same story+-- * Functions returning existential types, like 'H.withNullSpace',+-- 'H.withOrth', 'H.withRows', etc.; not quite sure what the best way+-- to handle these are at the moment.+-- * 'H.withRows' and 'H.withColumns' made "type-safe", without+-- existential types, with 'fromRows' and 'fromColumns'.+--+-- Some other notes:+--+-- * Added 'sumElements', as well, for convenience.++module Numeric.LinearAlgebra.Static.Backprop (+ -- * Vector+ H.R+ , H.ℝ+ , vec2+ , vec3+ , vec4+ , (&)+ , (#)+ , split+ , headTail+ , vector+ , linspace+ , H.range+ , H.dim+ -- * Matrix+ , H.L+ , H.Sq+ , row+ , col+ , (|||)+ , (===)+ , splitRows+ , splitCols+ , unrow+ , uncol+ , tr+ , H.eye+ , diag+ , matrix+ -- * Complex+ , H.ℂ+ , H.C+ , H.M+ , H.𝑖+ -- * Products+ , (<>)+ , (#>)+ , (<.>)+ -- * Factorizations+ , svd+ , svd_+ , H.Eigen+ , eigensystem+ , eigenvalues+ , chol+ -- * Norms+ , H.Normed+ , norm_0+ , norm_1V+ , norm_1M+ , norm_2V+ , norm_2M+ , norm_InfV+ , norm_InfM+ -- * Misc+ , mean+ , meanCov+ , meanL+ , cov+ , H.Disp(..)+ -- ** Domain+ , H.Domain+ , mul+ , app+ , dot+ , cross+ , diagR+ , dvmap+ , dvmap'+ , dmmap+ , dmmap'+ , outer+ , zipWithVector+ , zipWithVector'+ , det+ , invlndet+ , lndet+ , inv+ -- ** Conversions+ , toRows+ , toColumns+ , fromRows+ , fromColumns+ -- ** Misc Operations+ , konst+ , sumElements+ , extractV+ , extractM+ , create+ , H.Diag+ , takeDiag+ , H.Sym+ , sym+ , mTm+ , unSym+ , (<·>)+ ) where++import Data.ANum+import Data.Maybe+import Data.Proxy+import Foreign.Storable+import GHC.TypeLits+import Lens.Micro hiding ((&))+import Numeric.Backprop+import Numeric.Backprop.Op+import Numeric.Backprop.Tuple+import Unsafe.Coerce+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Sized as SVG+import qualified Data.Vector.Sized as SV+import qualified Data.Vector.Storable.Sized as SVS+import qualified Numeric.LinearAlgebra as HU+import qualified Numeric.LinearAlgebra.Devel as HU+import qualified Numeric.LinearAlgebra.Static as H+import qualified Numeric.LinearAlgebra.Static.Vector as H++#if MIN_VERSION_base(4,11,0)+import Prelude hiding ((<>))+#endif++vec2+ :: Reifies s W+ => BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s (H.R 2)+vec2 = liftOp2 $ opIsoN+ (\(x ::< y ::< Ø) -> H.vec2 x y )+ (\(HU.toList.H.extract->[dx,dy]) -> dx ::< dy ::< Ø)+{-# INLINE vec2 #-}++vec3+ :: Reifies s W+ => BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s (H.R 3)+vec3 = liftOp3 $ opIsoN+ (\(x ::< y ::< z ::< Ø) -> H.vec3 x y z )+ (\(HU.toList.H.extract->[dx,dy,dz]) -> dx ::< dy ::< dz ::< Ø)+{-# INLINE vec3 #-}++vec4+ :: Reifies s W+ => BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s (H.R 4)+vec4 vX vY vZ vW = liftOp o (vX :< vY :< vZ :< vW :< Ø)+ where+ o :: Op '[H.ℝ, H.ℝ, H.ℝ, H.ℝ] (H.R 4)+ o = opIsoN+ (\(x ::< y ::< z ::< w ::< Ø) -> H.vec4 x y z w )+ (\(HU.toList.H.extract->[dx,dy,dz,dw]) -> dx ::< dy ::< dz ::< dw ::< Ø)+ {-# INLINE o #-}+{-# INLINE vec4 #-}++(&) :: (Reifies s W, KnownNat n, 1 <= n, KnownNat (n + 1))+ => BVar s (H.R n)+ -> BVar s H.ℝ+ -> BVar s (H.R (n + 1))+(&) = liftOp2 $ opIsoN+ (\(xs ::< y ::< Ø) -> xs H.& y )+ (\(H.split->(dxs,dy)) -> dxs ::< fst (H.headTail dy) ::< Ø)+infixl 4 &+{-# INLINE (&) #-}++(#) :: (Reifies s W, KnownNat n, KnownNat m)+ => BVar s (H.R n)+ -> BVar s (H.R m)+ -> BVar s (H.R (n + m))+(#) = liftOp2 $ opIsoN+ (\(x ::< y ::< Ø) -> x H.# y )+ (\(H.split->(dX,dY)) -> dX ::< dY ::< Ø)+infixl 4 #+{-# INLINE (#) #-}++split+ :: forall p n s. (Reifies s W, KnownNat p, KnownNat n, p <= n)+ => BVar s (H.R n)+ -> (BVar s (H.R p), BVar s (H.R (n - p)))+split v = (t ^^. _1, t ^^. _2) -- should we just return the T2 ?+ where+ t = liftOp1 (opIso (tupT2 . H.split)+ (uncurryT2 (H.#))+ ) v+ {-# NOINLINE t #-}+{-# INLINE split #-}++headTail+ :: (Reifies s W, KnownNat n, 1 <= n)+ => BVar s (H.R n)+ -> (BVar s H.ℝ, BVar s (H.R (n - 1)))+headTail v = (t ^^. _1, t ^^. _2)+ where+ t = liftOp1 (opIso (tupT2 . H.headTail)+ (\(T2 d dx) -> (H.konst d :: H.R 1) H.# dx)+ ) v+ {-# NOINLINE t #-}+{-# INLINE headTail #-}++-- | Potentially extremely bad for anything but short lists!!!+vector+ :: forall n s. (Reifies s W, KnownNat n)+ => SV.Vector n (BVar s H.ℝ)+ -> BVar s (H.R n)+vector vs =+ liftOp1 (opIso (H.vecR . SVG.convert) (SVG.convert . H.rVec))+ (collectVar vs)+{-# INLINE vector #-}++linspace+ :: forall n s. (Reifies s W, KnownNat n)+ => BVar s H.ℝ+ -> BVar s H.ℝ+ -> BVar s (H.R n)+linspace = liftOp2 . op2 $ \l u ->+ ( H.linspace (l, u)+ , \d -> let n1 = fromInteger $ natVal (Proxy @n) - 1+ dDot = ((H.range - 1) H.<.> d) / n1+ dSum = HU.sumElements . H.extract $ d+ in (dSum - dDot, dDot)+ )+{-# INLINE linspace #-}++row :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s (H.L 1 n)+row = liftOp1 $ opIso H.row H.unrow+{-# INLINE row #-}++col :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s (H.L n 1)+col = liftOp1 $ opIso H.col H.uncol+{-# INLINE col #-}++(|||) :: (Reifies s W, KnownNat c, KnownNat r1, KnownNat (r1 + r2))+ => BVar s (H.L c r1)+ -> BVar s (H.L c r2)+ -> BVar s (H.L c (r1 + r2))+(|||) = liftOp2 $ opIsoN+ (\(x ::< y ::< Ø) -> x H.||| y )+ (\(H.splitCols->(dX,dY)) -> dX ::< dY ::< Ø)+infixl 3 |||+{-# INLINE (|||) #-}++(===) :: (Reifies s W, KnownNat c, KnownNat r1, KnownNat (r1 + r2))+ => BVar s (H.L r1 c)+ -> BVar s (H.L r2 c)+ -> BVar s (H.L (r1 + r2) c)+(===) = liftOp2 $ opIsoN+ (\(x ::< y ::< Ø) -> x H.=== y )+ (\(H.splitRows->(dX,dY)) -> dX ::< dY ::< Ø)+infixl 2 ===+{-# INLINE (===) #-}++splitRows+ :: forall p m n s. (Reifies s W, KnownNat p, KnownNat m, KnownNat n, p <= m)+ => BVar s (H.L m n)+ -> (BVar s (H.L p n), BVar s (H.L (m - p) n))+splitRows v = (t ^^. _1, t ^^. _2)+ where+ t = liftOp1 (opIso (tupT2 . H.splitRows)+ (\(T2 dx dy) -> dx H.=== dy)+ ) v+ {-# NOINLINE t #-}+{-# INLINE splitRows #-}++splitCols+ :: forall p m n s. (Reifies s W, KnownNat p, KnownNat m, KnownNat n, KnownNat (n - p), p <= n)+ => BVar s (H.L m n)+ -> (BVar s (H.L m p), BVar s (H.L m (n - p)))+splitCols v = (t ^^. _1, t ^^. _2)+ where+ t = liftOp1 (opIso (tupT2 . H.splitCols)+ (uncurryT2 (H.|||))+ ) v+ {-# NOINLINE t #-}+{-# INLINE splitCols #-}++unrow+ :: (Reifies s W, KnownNat n)+ => BVar s (H.L 1 n)+ -> BVar s (H.R n)+unrow = liftOp1 $ opIso H.unrow H.row+{-# INLINE unrow #-}++uncol+ :: (Reifies s W, KnownNat n)+ => BVar s (H.L n 1)+ -> BVar s (H.R n)+uncol = liftOp1 $ opIso H.uncol H.col+{-# INLINE uncol #-}++tr :: (Reifies s W, HU.Transposable m mt, HU.Transposable mt m, Num m, Num mt)+ => BVar s m+ -> BVar s mt+tr = liftOp1 $ opIso H.tr H.tr+{-# INLINE tr #-}++diag+ :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s (H.Sq n)+diag = liftOp1 . op1 $ \x -> (H.diag x, H.takeDiag)+{-# INLINE diag #-}++-- | Potentially extremely bad for anything but short lists!!!+matrix+ :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => [BVar s H.ℝ]+ -> BVar s (H.L m n)+matrix vs = case SV.fromList @(m * n) vs of+ Nothing -> error "matrix: invalid number of elements"+ Just vs' ->+ liftOp1 (opIso (fromJust . H.create . HU.reshape n . VG.convert . SV.fromSized)+ (SV.concatMap (SVG.convert . H.rVec) . H.lRows)+ )+ (collectVar vs')+ where+ n = fromInteger $ natVal (Proxy @n)+{-# INLINE matrix #-}++-- | Matrix product+(<>)+ :: (Reifies s W, KnownNat m, KnownNat k, KnownNat n)+ => BVar s (H.L m k)+ -> BVar s (H.L k n)+ -> BVar s (H.L m n)+(<>) = mul+infixr 8 <>+{-# INLINE (<>) #-}++-- | Matrix-vector product+(#>)+ :: (Reifies s W, KnownNat m, KnownNat n)+ => BVar s (H.L m n)+ -> BVar s (H.R n)+ -> BVar s (H.R m)+(#>) = app+infixr 8 #>+{-# INLINE (#>) #-}++-- | Dot product+(<.>)+ :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s (H.R n)+ -> BVar s H.ℝ+(<.>) = dot+infixr 8 <.>+{-# INLINE (<.>) #-}++-- | Can only get the singular values, for now. Let me know if you find an+-- algorithm that can compute the gradients based on differentials for the+-- other matricies!+--+svd :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => BVar s (H.L m n)+ -> BVar s (H.R n)+svd = liftOp1 . op1 $ \x ->+ let (u, σ, v) = H.svd x+ in ( σ+ , \(dΣ :: H.R n) -> (u H.<> H.diagR 0 dΣ) H.<> H.tr v+ -- must manually associate because of bug in diagR in+ -- hmatrix-0.18.2.0+ )+{-# INLINE svd #-}++-- | Version of 'svd' that returns the full SVD, but if you attempt to find+-- the gradient, it will fail at runtime if you ever use U or V.+svd_+ :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => BVar s (H.L m n)+ -> (BVar s (H.L m m), BVar s (H.R n), BVar s (H.L n n))+svd_ r = (t ^^. _1, t ^^. _2, t ^^. _3)+ where+ o :: Op '[H.L m n] (T3 (H.L m m) (H.R n) (H.L n n))+ o = op1 $ \x ->+ let (u, σ, v) = H.svd x+ in ( T3 u σ v+ , \(T3 dU dΣ dV) ->+ if H.norm_0 dU == 0 && H.norm_0 dV == 0+ then (u H.<> H.diagR 0 dΣ) H.<> H.tr v+ else error "svd_: Cannot backprop if U and V are used."+ )+ {-# INLINE o #-}+ t = liftOp1 o r+ {-# NOINLINE t #-}+{-# INLINE svd_ #-}++helpEigen :: KnownNat n => H.Sym n -> (H.R n, H.L n n, H.L n n, H.L n n)+helpEigen x = (l, v, H.inv v, H.tr v)+ where+ (l, v) = H.eigensystem x+{-# INLINE helpEigen #-}++-- | /NOTE/ The gradient is not necessarily symmetric! The gradient is not+-- meant to be retireved directly; insteadl, 'eigenvalues' is meant to be+-- used as a part of a larger computation, and the gradient as an+-- intermediate step.+eigensystem+ :: forall n s. (Reifies s W, KnownNat n)+ => BVar s (H.Sym n)+ -> (BVar s (H.R n), BVar s (H.L n n))+eigensystem u = (t ^^. _1, t ^^. _2)+ where+ o :: Op '[H.Sym n] (T2 (H.R n) (H.L n n))+ o = op1 $ \x ->+ let (l, v, vInv, vTr) = helpEigen x+ lRep = H.rowsL . SV.replicate $ l+ fMat = (1 - H.eye) * (lRep - H.tr lRep)+ in ( T2 l v+ , \(T2 dL dV) -> unsafeCoerce $+ H.tr vInv+ H.<> (H.diag dL + fMat * (vTr H.<> dV))+ H.<> vTr+ )+ {-# INLINE o #-}+ t = liftOp1 o u+ {-# NOINLINE t #-}+{-# INLINE eigensystem #-}++-- | /NOTE/ The gradient is not necessarily symmetric! The gradient is not+-- meant to be retireved directly; insteadl, 'eigenvalues' is meant to be+-- used as a part of a larger computation, and the gradient as an+-- intermediate step.+eigenvalues+ :: forall n s. (Reifies s W, KnownNat n)+ => BVar s (H.Sym n)+ -> BVar s (H.R n)+eigenvalues = liftOp1 . op1 $ \x ->+ let (l, _, vInv, vTr) = helpEigen x+ in ( l+ , \dL -> unsafeCoerce $+ H.tr vInv H.<> H.diag dL H.<> vTr+ )+{-# INLINE eigenvalues #-}++-- | Algorithm from https://arxiv.org/abs/1602.07527+--+-- The paper also suggests a potential imperative algorithm that might+-- help. Need to benchmark to see what is best.+--+-- /NOTE/ The gradient is not necessarily symmetric! The gradient is not+-- meant to be retireved directly; insteadl, 'eigenvalues' is meant to be+-- used as a part of a larger computation, and the gradient as an+-- intermediate step.+chol+ :: forall n s. (Reifies s W, KnownNat n)+ => BVar s (H.Sym n)+ -> BVar s (H.Sq n)+chol = liftOp1 . op1 $ \x ->+ let l = H.chol x+ lInv = H.inv l+ phi :: H.Sq n+ phi = H.build $ \i j -> case compare i j of+ LT -> 1+ EQ -> 0.5+ GT -> 0+ in ( l+ , \dL -> let s = H.tr lInv H.<> (phi * (H.tr l H.<> dL)) H.<> lInv+ in unsafeCoerce $ s + H.tr s - H.eye * s+ )+{-# INLINE chol #-}++-- | Number of non-zero items+norm_0+ :: (Reifies s W, H.Normed a, Num a)+ => BVar s a+ -> BVar s H.ℝ+norm_0 = liftOp1 . op1 $ \x -> (H.norm_0 x, const 0)+{-# INLINE norm_0 #-}++-- | Sum of absolute values+norm_1V+ :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s H.ℝ+norm_1V = liftOp1 . op1 $ \x -> (H.norm_1 x, (* signum x) . H.konst)+{-# INLINE norm_1V #-}++-- | Maximum 'H.norm_1' of columns+norm_1M+ :: (Reifies s W, KnownNat n, KnownNat m)+ => BVar s (H.L n m)+ -> BVar s H.ℝ+norm_1M = liftOp1 . op1 $ \x ->+ let n = H.norm_1 x+ in (n, \d -> let d' = H.konst d+ in H.colsL+ . SV.map (\c -> if H.norm_1 c == n+ then d' * signum c+ else 0+ )+ . H.lCols+ $ x+ )+{-# INLINE norm_1M #-}++-- | Square root of sum of squares+--+-- Be aware that gradient diverges when the norm is zero+norm_2V+ :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s H.ℝ+norm_2V = liftOp1 . op1 $ \x ->+ let n = H.norm_2 x+ in (n, \d -> x * H.konst (d / n))+{-# INLINE norm_2V #-}++-- | Maximum singular value+norm_2M+ :: (Reifies s W, KnownNat n, KnownNat m)+ => BVar s (H.L n m)+ -> BVar s H.ℝ+norm_2M = liftOp1 . op1 $ \x ->+ let n = H.norm_2 x+ (head.H.toColumns->u1,_,head.H.toColumns->v1) = H.svd x+ in (n, \d -> H.konst d * (u1 `H.outer` v1))+{-# INLINE norm_2M #-}++-- | Maximum absolute value+norm_InfV+ :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s H.ℝ+norm_InfV = liftOp1 . op1 $ \x ->+ let n :: H.ℝ+ n = H.norm_Inf x+ in (n, \d -> H.vecR+ . SVS.map (\e -> if abs e == n+ then signum e * d+ else 0+ )+ . H.rVec+ $ x+ )+{-# ANN norm_InfV "HLint: ignore Use camelCase" #-}+{-# INLINE norm_InfV #-}++-- | Maximum 'H.norm_1' of rows+norm_InfM+ :: (Reifies s W, KnownNat n, KnownNat m)+ => BVar s (H.L n m)+ -> BVar s H.ℝ+norm_InfM = liftOp1 . op1 $ \x ->+ let n = H.norm_Inf x+ in (n, \d -> let d' = H.konst d+ in H.rowsL+ . SV.map (\c -> if H.norm_1 c == n+ then d' * signum c+ else 0+ )+ . H.lRows+ $ x+ )+{-# ANN norm_InfM "HLint: ignore Use camelCase" #-}+{-# INLINE norm_InfM #-}++mean+ :: (Reifies s W, KnownNat n, 1 <= n)+ => BVar s (H.R n)+ -> BVar s H.ℝ+mean = liftOp1 . op1 $ \x -> (H.mean x, H.konst . (/ H.norm_0 x))+{-# INLINE mean #-}++gradCov+ :: forall m n. (KnownNat m, KnownNat n)+ => H.L m n+ -> H.R n+ -> H.Sym n+ -> H.L m n+gradCov x μ dσ = H.rowsL+ . SV.map (subtract (dDiffsSum / m))+ . H.lRows+ $ dDiffs+ where+ diffs = H.rowsL . SV.map (subtract μ) . H.lRows $ x+ dDiffs = H.konst (2/n) * (diffs H.<> H.tr (H.unSym dσ))+ dDiffsSum = sum . H.toRows $ dDiffs+ m = fromIntegral $ natVal (Proxy @m)+ n = fromIntegral $ natVal (Proxy @n)+{-# INLINE gradCov #-}++-- | Mean and covariance. If you know you only want to use one or the+-- other, use 'meanL' or 'cov'.+meanCov+ :: forall m n s. (Reifies s W, KnownNat n, KnownNat m, 1 <= m)+ => BVar s (H.L m n)+ -> (BVar s (H.R n), BVar s (H.Sym n))+meanCov v = (t ^^. _1, t ^^. _2)+ where+ m = fromInteger $ natVal (Proxy @m)+ t = ($ v) . liftOp1 . op1 $ \x ->+ let (μ, σ) = H.meanCov x+ in ( T2 μ σ+ , \(T2 dμ dσ) ->+ let gradMean = H.rowsL+ . SV.replicate+ $ (dμ / H.konst m)+ in gradMean + gradCov x μ dσ+ )+ {-# NOINLINE t #-}+{-# INLINE meanCov #-}++-- | 'meanCov', but if you know you won't use the covariance.+meanL+ :: forall m n s. (Reifies s W, KnownNat n, KnownNat m, 1 <= m)+ => BVar s (H.L m n)+ -> BVar s (H.R n)+meanL = liftOp1 . op1 $ \x ->+ ( fst (H.meanCov x)+ , H.rowsL . SV.replicate . (/ H.konst m)+ )+ where+ m = fromInteger $ natVal (Proxy @m)+{-# INLINE meanL #-}++-- | 'cov', but if you know you won't use the covariance.+cov+ :: forall m n s. (Reifies s W, KnownNat n, KnownNat m, 1 <= m)+ => BVar s (H.L m n)+ -> BVar s (H.Sym n)+cov = liftOp1 . op1 $ \x ->+ let (μ, σ) = H.meanCov x+ in (σ, gradCov x μ)+{-# INLINE cov #-}++mul :: ( Reifies s W+ , KnownNat m+ , KnownNat k+ , KnownNat n+ , H.Domain field vec mat+ , Num (mat m k)+ , Num (mat k n)+ , Num (mat m n)+ , HU.Transposable (mat m k) (mat k m)+ , HU.Transposable (mat k n) (mat n k)+ )+ => BVar s (mat m k)+ -> BVar s (mat k n)+ -> BVar s (mat m n)+mul = liftOp2 . op2 $ \x y ->+ ( x `H.mul` y+ , \d -> (d `H.mul` H.tr y, H.tr x `H.mul` d)+ )+{-# INLINE mul #-}++app :: ( Reifies s W+ , KnownNat m+ , KnownNat n+ , H.Domain field vec mat+ , Num (mat m n)+ , Num (vec n)+ , Num (vec m)+ , HU.Transposable (mat m n) (mat n m)+ )+ => BVar s (mat m n)+ -> BVar s (vec n)+ -> BVar s (vec m)+app = liftOp2 . op2 $ \xs y ->+ ( xs `H.app` y+ , \d -> (d `H.outer` y, H.tr xs `H.app` d)+ )+{-# INLINE app #-}++dot :: ( Reifies s W+ , KnownNat n+ , H.Domain field vec mat+ , H.Sized field (vec n) d+ , Num (vec n)+ )+ => BVar s (vec n)+ -> BVar s (vec n)+ -> BVar s field+dot = liftOp2 . op2 $ \x y ->+ ( x `H.dot` y+ , \d -> let d' = H.konst d+ in (d' * y, x * d')+ )+{-# INLINE dot #-}++cross+ :: ( Reifies s W+ , H.Domain field vec mat+ , Num (vec 3)+ )+ => BVar s (vec 3)+ -> BVar s (vec 3)+ -> BVar s (vec 3)+cross = liftOp2 . op2 $ \x y ->+ ( x `H.cross` y+ , \d -> (y `H.cross` d, d `H.cross` x)+ )+{-# INLINE cross #-}++-- | Create matrix with diagonal, and fill with default entries+diagR+ :: forall m n k field vec mat s.+ ( Reifies s W+ , H.Domain field vec mat+ , Num (vec k)+ , Num (mat m n)+ , KnownNat m+ , KnownNat n+ , KnownNat k+ , HU.Container HU.Vector field+ , H.Sized field (mat m n) HU.Matrix+ , H.Sized field (vec k) HU.Vector+ )+ => BVar s field -- ^ default value+ -> BVar s (vec k) -- ^ diagonal+ -> BVar s (mat m n)+diagR = liftOp2 . op2 $ \c x ->+ ( H.diagR c x+ , \d -> ( HU.sumElements . H.extract $ H.diagR 1 (0 :: vec k) * d+ , fromJust . H.create . HU.takeDiag . H.extract $ d+ )+ )+{-# INLINE diagR #-}++dvmap+ :: ( Reifies s W+ , Num (vec n)+ , Storable field+ , Storable (field, field)+ , H.Sized field (vec n) HU.Vector+ )+ => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)+ -> BVar s (vec n)+ -> BVar s (vec n)+dvmap f = liftOp1 . op1 $ \x ->+ let (y, dx) = HU.unzipVector $ VG.map (backprop f) (H.extract x)+ in ( fromJust (H.create y)+ , \d -> d * fromJust (H.create dx)+ )+{-# INLINE dvmap #-}++-- TODO: Can be made more efficient if backprop exports+-- a custom-total-derivative version++-- | A version of 'dvmap' that is less performant but is based on+-- 'H.zipWithVector' from 'H.Domain'.+dvmap'+ :: ( Reifies s W+ , KnownNat n+ , H.Domain field vec mat+ , Num (vec n)+ , Num field+ )+ => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)+ -> BVar s (vec n)+ -> BVar s (vec n)+dvmap' f = liftOp1 . op1 $ \x ->+ ( H.dvmap (evalBP f) x+ , (H.dvmap (gradBP f) x *)+ )+{-# INLINE dvmap' #-}++dmmap+ :: forall n m mat field s.+ ( Reifies s W+ , KnownNat m+ , Num (mat n m)+ , Storable (field, field)+ , H.Sized field (mat n m) HU.Matrix+ , HU.Element field+ )+ => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)+ -> BVar s (mat n m)+ -> BVar s (mat n m)+dmmap f = liftOp1 . op1 $ \x ->+ let (y', dx') = HU.unzipVector+ . VG.map (backprop f)+ . HU.flatten+ $ H.extract x+ in ( fromJust . H.create . HU.reshape m $ y'+ , \d -> (* d) . fromJust . H.create . HU.reshape m $ dx'+ )+ where+ m :: Int+ m = fromInteger $ natVal (Proxy @m)+{-# INLINE dmmap #-}++dmmap'+ :: ( Reifies s W+ , KnownNat n+ , KnownNat m+ , H.Domain field vec mat+ , Num (mat n m)+ , Num field+ )+ => (forall s'. Reifies s' W => BVar s' field -> BVar s' field)+ -> BVar s (mat n m)+ -> BVar s (mat n m)+dmmap' f = liftOp1 . op1 $ \x ->+ ( H.dmmap (evalBP f) x+ , (H.dmmap (gradBP f) x *)+ )+{-# INLINE dmmap' #-}++outer+ :: ( Reifies s W+ , KnownNat m+ , KnownNat n+ , H.Domain field vec mat+ , HU.Transposable (mat n m) (mat m n)+ , Num (vec n)+ , Num (vec m)+ , Num (mat n m)+ )+ => BVar s (vec n)+ -> BVar s (vec m)+ -> BVar s (mat n m)+outer = liftOp2 . op2 $ \x y ->+ ( x `H.outer` y+ , \d -> ( d `H.app` y+ , H.tr d `H.app` x)+ )+{-# INLINE outer #-}++zipWithVector+ :: ( Reifies s W+ , Num (vec n)+ , Storable field+ , Storable (field, field, field)+ , H.Sized field (vec n) HU.Vector+ )+ => (forall s'. Reifies s' W => BVar s' field -> BVar s' field -> BVar s' field)+ -> BVar s (vec n)+ -> BVar s (vec n)+ -> BVar s (vec n)+zipWithVector f = liftOp2 . op2 $ \(H.extract->x) (H.extract->y) ->+ let (z,dx,dy) = VG.unzip3+ $ VG.zipWith (\x' y' ->+ let (z', (dx', dy')) = backprop2 f x' y'+ in (z', dx', dy')+ ) x y+ in ( fromJust (H.create z)+ , \d -> (d * fromJust (H.create dx), d * fromJust (H.create dy))+ )+{-# INLINE zipWithVector #-}++-- | A version of 'zipWithVector' that is less performant but is based on+-- 'H.zipWithVector' from 'H.Domain'.+zipWithVector'+ :: ( Reifies s W+ , KnownNat n+ , H.Domain field vec mat+ , Num (vec n)+ , Num field+ )+ => (forall s'. Reifies s' W => BVar s' field -> BVar s' field -> BVar s' field)+ -> BVar s (vec n)+ -> BVar s (vec n)+ -> BVar s (vec n)+zipWithVector' f = liftOp2 . op2 $ \x y ->+ ( H.zipWithVector (evalBP2 f) x y+ , \d -> let dx = H.zipWithVector (\x' -> fst . gradBP2 f x') x y+ dy = H.zipWithVector (\x' -> snd . gradBP2 f x') x y+ in (d * dx, d * dy)+ )+{-# INLINE zipWithVector' #-}++det :: ( Reifies s W+ , KnownNat n+ , Num (mat n n)+ , H.Domain field vec mat+ , H.Sized field (mat n n) d+ , HU.Transposable (mat n n) (mat n n)+ )+ => BVar s (mat n n)+ -> BVar s field+det = liftOp1 . op1 $ \x ->+ let xDet = H.det x+ xInv = H.inv x+ in ( xDet, \d -> H.konst (d * xDet) * H.tr xInv )+{-# INLINE det #-}++-- | The inverse and the natural log of the determinant together. If you+-- know you don't need the inverse, it is best to use 'lndet'.+invlndet+ :: forall n mat field vec d s.+ ( Reifies s W+ , KnownNat n+ , Num (mat n n)+ , H.Domain field vec mat+ , H.Sized field (mat n n) d+ , HU.Transposable (mat n n) (mat n n)+ )+ => BVar s (mat n n)+ -> (BVar s (mat n n), (BVar s field, BVar s field))+invlndet v = (t ^^. _1, (t ^^. _2, t ^^. _3))+ where+ o :: Op '[mat n n] (T3 (mat n n) field field)+ o = op1 $ \x ->+ let (i,(ldet, s)) = H.invlndet x+ iTr = H.tr i+ in ( T3 i ldet s+ , \(T3 dI dLDet _) ->+ let gradI = - iTr `H.mul` dI `H.mul` iTr+ gradLDet = H.konst dLDet * H.tr i+ in gradI + gradLDet+ )+ {-# INLINE o #-}+ t = liftOp1 o v+ {-# NOINLINE t #-}+{-# INLINE invlndet #-}++-- | The natural log of the determinant.+lndet+ :: forall n mat field vec d s.+ ( Reifies s W+ , KnownNat n+ , Num (mat n n)+ , H.Domain field vec mat+ , H.Sized field (mat n n) d+ , HU.Transposable (mat n n) (mat n n)+ )+ => BVar s (mat n n)+ -> BVar s field+lndet = liftOp1 . op1 $ \x ->+ let (i,(ldet,_)) = H.invlndet x+ in (ldet, (* H.tr i) . H.konst)+{-# INLINE lndet #-}++inv :: ( Reifies s W+ , KnownNat n+ , Num (mat n n)+ , H.Domain field vec mat+ , HU.Transposable (mat n n) (mat n n)+ )+ => BVar s (mat n n)+ -> BVar s (mat n n)+inv = liftOp1 . op1 $ \x ->+ let xInv = H.inv x+ xInvTr = H.tr xInv+ in ( xInv, \d -> - xInvTr `H.mul` d `H.mul` xInvTr )+{-# INLINE inv #-}++toRows+ :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => BVar s (H.L m n)+ -> SV.Vector m (BVar s (H.R n))+toRows = sequenceVar . liftOp1 (opIso H.lRows H.rowsL)+{-# INLINE toRows #-}++toColumns+ :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => BVar s (H.L m n)+ -> SV.Vector n (BVar s (H.R m))+toColumns = sequenceVar . liftOp1 (opIso H.lCols H.colsL)+{-# INLINE toColumns #-}++fromRows+ :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => SV.Vector m (BVar s (H.R n))+ -> BVar s (H.L m n)+fromRows = liftOp1 (opIso H.rowsL H.lRows) . collectVar+{-# INLINE fromRows #-}++fromColumns+ :: forall m n s. (Reifies s W, KnownNat m, KnownNat n)+ => SV.Vector n (BVar s (H.R m))+ -> BVar s (H.L m n)+fromColumns = liftOp1 (opIso H.colsL H.lCols) . collectVar+{-# INLINE fromColumns #-}++konst+ :: forall t s d q. (Reifies q W, H.Sized t s d, HU.Container d t, Num s)+ => BVar q t+ -> BVar q s+konst = liftOp1 . op1 $ \x ->+ ( H.konst x+ , HU.sumElements . H.extract+ )+{-# INLINE konst #-}++sumElements+ :: forall t s d q. (Reifies q W, H.Sized t s d, HU.Container d t, Num s)+ => BVar q s+ -> BVar q t+sumElements = liftOp1 . op1 $ \x ->+ ( HU.sumElements . H.extract $ x+ , H.konst+ )+{-# INLINE sumElements #-}++-- | If there are extra items in the total derivative, they are dropped.+-- If there are missing items, they are treated as zero.+extractV+ :: forall t s q.+ ( Reifies q W+ , H.Sized t s HU.Vector+ , Num s+ , HU.Konst t Int HU.Vector+ , HU.Container HU.Vector t+ , Num (HU.Vector t)+ )+ => BVar q s+ -> BVar q (HU.Vector t)+extractV = liftOp1 . op1 $ \x ->+ let n = H.size x+ in ( H.extract x+ , \d -> let m = HU.size d+ m' = case compare n m of+ LT -> HU.subVector 0 n d+ EQ -> d+ GT -> HU.vjoin [d, HU.konst 0 (n - m)]+ in fromJust . H.create $ m'+ )+{-# INLINE extractV #-}++-- | If there are extra items in the total derivative, they are dropped.+-- If there are missing items, they are treated as zero.+extractM+ :: forall t s q.+ ( Reifies q W+ , H.Sized t s HU.Matrix+ , Num s+ , HU.Konst t (Int, Int) HU.Matrix+ , HU.Container HU.Matrix t+ , Num (HU.Matrix t)+ )+ => BVar q s+ -> BVar q (HU.Matrix t)+extractM = liftOp1 . op1 $ \x ->+ let (xI,xJ) = H.size x+ in ( H.extract x+ , \d -> let (dI,dJ) = HU.size d+ m' = case (compare xI dI, compare xJ dJ) of+ (LT, LT) -> d HU.?? (HU.Take xI, HU.Take xJ)+ (LT, EQ) -> d HU.?? (HU.Take xI, HU.All)+ (LT, GT) -> d HU.?? (HU.Take xI, HU.All)+ HU.||| HU.konst 0 (xI, xJ - dJ)+ (EQ, LT) -> d HU.?? (HU.All , HU.Take xJ)+ (EQ, EQ) -> d+ (EQ, GT) -> d HU.?? (HU.All, HU.All)+ HU.||| HU.konst 0 (xI, xJ - dJ)+ (GT, LT) -> d HU.?? (HU.All, HU.Take xJ)+ HU.=== HU.konst 0 (xI - dI, xJ)+ (GT, EQ) -> d HU.?? (HU.All, HU.All)+ HU.=== HU.konst 0 (xI - dI, xJ)+ (GT, GT) -> HU.fromBlocks+ [[d,0 ]+ ,[0,HU.konst 0 (xI - dI, xJ - dJ)]+ ]+ in fromJust . H.create $ m'+ )+{-# INLINE extractM #-}++create+ :: forall t s d q. (Reifies q W, H.Sized t s d, Num s, Num (d t))+ => BVar q (d t)+ -> Maybe (BVar q s)+create = fmap (unANum . sequenceVar) . liftOp1 $+ opIso (ANum . H.create)+ (maybe 0 H.extract . unANum )+{-# INLINE create #-}+++takeDiag+ :: ( Reifies s W+ , KnownNat n+ , H.Diag (mat n n) (vec n)+ , H.Domain field vec mat+ , Num (vec n)+ , Num (mat n n)+ , Num field+ )+ => BVar s (mat n n)+ -> BVar s (vec n)+takeDiag = liftOp1 . op1 $ \x ->+ ( H.takeDiag x+ , H.diagR 0+ )+{-# INLINE takeDiag #-}++-- |+-- \[+-- \frac{1}{2} (M + M^T)+-- \]+sym :: (Reifies s W, KnownNat n)+ => BVar s (H.Sq n)+ -> BVar s (H.Sym n)+sym = liftOp1 . op1 $ \x ->+ ( H.sym x+ , H.unSym . H.sym . H.unSym+ )+{-# INLINE sym #-}++-- |+-- \[+-- M^T M+-- \]+mTm :: (Reifies s W, KnownNat m, KnownNat n)+ => BVar s (H.L m n)+ -> BVar s (H.Sym n)+mTm = liftOp1 . op1 $ \x ->+ ( H.mTm x+ , \d -> 2 * (x H.<> H.unSym d)+ )+{-# INLINE mTm #-}++-- | Warning: the gradient is going necessarily symmetric, and so is /not/+-- meant to be used directly. Rather, it is meant to be used in the middle+-- (or at the end) of a longer computation.+unSym+ :: (Reifies s W, KnownNat n)+ => BVar s (H.Sym n)+ -> BVar s (H.Sq n)+unSym = liftOp1 (opIso H.unSym unsafeCoerce)+{-# INLINE unSym #-}++-- | Unicode synonym for '<.>>'+(<·>)+ :: (Reifies s W, KnownNat n)+ => BVar s (H.R n)+ -> BVar s (H.R n)+ -> BVar s H.ℝ+(<·>) = dot+infixr 8 <·>+{-# INLINE (<·>) #-}
+ test/Nudge.hs view
@@ -0,0 +1,177 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE TypeFamilies #-}++module Nudge where++import Control.Monad+import Data.Bifunctor+import Data.Finite+import Data.Kind+import Data.Maybe+import Data.Proxy+import GHC.TypeLits+import Hedgehog+import Lens.Micro+import Lens.Micro.Platform ()+import Numeric.Backprop+import Numeric.Backprop.Tuple+import qualified Data.Ix as Ix+import qualified Data.Vector.Sized as SV+import qualified Hedgehog.Gen as Gen+import qualified Hedgehog.Range as Range+import qualified Numeric.LinearAlgebra as HU+import qualified Numeric.LinearAlgebra.Static as H+import qualified Numeric.LinearAlgebra.Static.Backprop as B++nudge :: Double+nudge = 1e-6++eps :: Double+eps = 1e-11++class (Num c, Show c, Show (TIx c)) => Testing c where+ type TIx c :: Type+ allIx :: c -> [TIx c]+ ixLens :: TIx c -> Lens' c Double+ scalarize :: Reifies s W => BVar s c -> BVar s Double+ genTest :: Gen c++sized+ :: forall s t d. H.Sized t s d+ => Lens' s (d t)+sized f = fmap (fromJust . H.create) . f . H.extract++ixContainer+ :: forall t d. HU.Container d t+ => HU.IndexOf d+ -> Lens' (d t) t+ixContainer i = lens (`HU.atIndex` i)+ (\xs x -> HU.accum xs (\_ _ -> x) [(i, x)])++instance Testing Double where+ type TIx Double = ()+ allIx _ = [()]+ ixLens _ = id+ scalarize = abs+ genTest = Gen.filter ((> eps) . (**2)) $+ Gen.double (Range.linearFracFrom 0 (-5) 5)++instance KnownNat n => Testing (H.R n) where+ type TIx (H.R n) = Int+ allIx v = [0 .. H.size v - 1]+ ixLens i = sized . ixContainer i+ scalarize = B.norm_2V+ genTest = H.vector <$> replicateM n genTest+ where+ n = fromInteger $ natVal (Proxy @n)++instance (KnownNat n, KnownNat m) => Testing (H.L n m) where+ type TIx (H.L n m) = (Int, Int)+ allIx m = Ix.range ((0,0), bimap pred pred (H.size m))+ ixLens i = sized . ixContainer i+ scalarize = sqrt . B.sumElements . (**2)+ genTest = H.matrix <$> replicateM nm genTest+ where+ nm = fromInteger $ natVal (Proxy @n) * natVal (Proxy @m)++instance Testing (HU.Vector Double) where+ type TIx (HU.Vector Double) = Int+ allIx v = [0 .. HU.size v - 1]+ ixLens = ixContainer+ scalarize = liftOp1 . op1 $ \xs -> (HU.sumElements xs, (`HU.konst` HU.size xs))+ genTest = HU.fromList <$> replicateM 3 genTest++instance Testing (HU.Matrix Double) where+ type TIx (HU.Matrix Double) = (Int, Int)+ allIx m = Ix.range ((0,0), bimap pred pred (HU.size m))+ ixLens = ixContainer+ scalarize = liftOp1 . op1 $ \xs -> (HU.sumElements xs, (`HU.konst` HU.size xs))+ genTest = HU.fromLists <$> (replicateM 3 . replicateM 2) genTest++instance (KnownNat n, Testing a) => Testing (SV.Vector n a) where+ type TIx (SV.Vector n a) = (Finite n, TIx a)+ allIx = fst . SV.imapM (\i x -> ((fromIntegral i,) <$> allIx x , x))+ ixLens (i,j) = SV.ix i . ixLens j+ scalarize = scalarize . liftOp1 o . (^ (2 :: Int))+ where+ o :: Op '[SV.Vector n a] a+ o = op1 $ \xs -> (SV.sum xs, SV.replicate)+ genTest = SV.replicateM genTest++instance (Testing a, Testing b) => Testing (T2 a b) where+ type TIx (T2 a b) = Either (TIx a) (TIx b)+ allIx (T2 x y) = (Left <$> allIx x)+ ++ (Right <$> allIx y)+ ixLens (Left i) = _1 . ixLens i+ ixLens (Right j) = _2 . ixLens j+ scalarize t = B.norm_2V (B.vec2 (scalarize (t ^^. _1))+ (scalarize (t ^^. _2))+ )+ genTest = T2 <$> genTest <*> genTest++instance (Testing a, Testing b, Testing c, Num a, Num b, Num c) => Testing (T3 a b c) where+ type TIx (T3 a b c) = Either (TIx a) (Either (TIx b) (TIx c))+ allIx (T3 x y z) = (Left <$> allIx x)+ ++ (Right . Left <$> allIx y)+ ++ (Right . Right <$> allIx z)+ ixLens (Left i ) = _1 . ixLens i+ ixLens (Right (Left j)) = _2 . ixLens j+ ixLens (Right (Right k)) = _3 . ixLens k+ scalarize t = B.norm_2V (B.vec3 (scalarize (t ^^. _1))+ (scalarize (t ^^. _2))+ (scalarize (t ^^. _3))+ )+ genTest = T3 <$> genTest <*> genTest <*> genTest++validGrad+ :: Monad m+ => Lens' c Double+ -> c+ -> c+ -> (c -> Double)+ -> PropertyT m (Double, Double)+validGrad l x0 g f = forAll $ Gen.double (Range.constantFrom 0 (-nudge) nudge) <&> \d ->+ let x = x0 & l %~ (+d)+ old = f x0 + (g ^. l) * d+ new = f x+ in (old, new)++nudgeProp+ :: forall c d. (Testing c, Testing d)+ => (forall s. Reifies s W => BVar s c -> BVar s d)+ -> Property+nudgeProp f = property $ do+ (inp, i) <- forAll $ do+ inp <- genTest+ i <- Gen.element (allIx inp)+ return (inp, i)+ let (r,gr) = backprop (scalarize . f) inp+ when (r**2 < eps) discard+ (old, new) <- validGrad (ixLens i) inp gr (evalBP (scalarize . f))+ footnoteShow (r, gr, old, new, (old - new)**2, ((old - new)/old)**2)+ assert $ ((old - new)/old)**2 < eps++nudgeProp2+ :: forall c d e. (Testing c, Testing d, Testing e)+ => (forall s. Reifies s W => BVar s c -> BVar s d -> BVar s e)+ -> Property+nudgeProp2 f = property $ do+ (inpC, inpD, i) <- forAll $ do+ inpC <- genTest+ inpD <- genTest+ i <- Gen.element (allIx (T2 inpC inpD))+ return (inpC, inpD, i)+ let (r, gr) = second tupT2 $ backprop2 (\x -> scalarize . f x) inpC inpD+ when (r**2 < eps) discard+ (old, new) <- validGrad (ixLens i) (T2 inpC inpD) gr+ (evalBP (\t -> scalarize $ f (t ^^. _1) (t ^^. _2)))+ footnoteShow (r, gr, old, new, (old - new)**2, ((old - new)/old)**2)+ assert $ ((old - new)/old)**2 < eps+
+ test/Spec.hs view
@@ -0,0 +1,259 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeApplications #-}++import Control.Monad+import Data.Bifunctor+import Data.Maybe+import Hedgehog+import Lens.Micro+import Nudge+import Numeric.Backprop+import Numeric.Backprop.Tuple+import Numeric.LinearAlgebra.Static (L, R)+import System.Exit+import System.IO+import qualified Numeric.LinearAlgebra.Static.Backprop as B++prop_vec2 :: Property+prop_vec2 = nudgeProp2 B.vec2++prop_vec3 :: Property+prop_vec3 = nudgeProp @(T3 Double Double Double)+ (\t -> B.vec3 (t ^^. _1) (t ^^. _2) (t ^^. _3))++prop_vec4 :: Property+prop_vec4 = nudgeProp2 @(T2 Double Double) @(T2 Double Double)+ (\x y -> B.vec4 (x ^^. _1) (x ^^. _2) (y ^^. _1) (y ^^. _2))++prop_snoc :: Property+prop_snoc = nudgeProp2 @(R 3) (B.&)++prop_append :: Property+prop_append = nudgeProp2 @(R 3) @(R 2) (B.#)++prop_split1 :: Property+prop_split1 = nudgeProp @(R 3) (fst . B.split @2)++prop_split2 :: Property+prop_split2 = nudgeProp @(R 3) (snd . B.split @2)++prop_headTail1 :: Property+prop_headTail1 = nudgeProp @(R 3) (fst . B.headTail)++prop_headTail2 :: Property+prop_headTail2 = nudgeProp @(R 3) (snd . B.headTail)++prop_vector :: Property+prop_vector = nudgeProp (B.vector @3 . sequenceVar)++prop_linspace :: Property+prop_linspace = nudgeProp2 (B.linspace @3)++prop_row :: Property+prop_row = nudgeProp @(R 3) B.row++prop_col :: Property+prop_col = nudgeProp @(R 3) B.col++prop_horzcat :: Property+prop_horzcat = nudgeProp2 @(L 3 2) @(L 3 1) (B.|||)++prop_vertcat :: Property+prop_vertcat = nudgeProp2 @(L 2 3) @(L 1 3) (B.===)++prop_splitRows1 :: Property+prop_splitRows1 = nudgeProp @(L 2 3) (fst . B.splitRows @1)++prop_splitRows2 :: Property+prop_splitRows2 = nudgeProp @(L 2 3) (snd . B.splitRows @1)++prop_splitCols1 :: Property+prop_splitCols1 = nudgeProp @(L 3 2) (fst . B.splitCols @1)++prop_splitCols2 :: Property+prop_splitCols2 = nudgeProp @(L 3 2) (snd . B.splitCols @1)++prop_unrow :: Property+prop_unrow = nudgeProp @(L 1 3) B.unrow++prop_uncol :: Property+prop_uncol = nudgeProp @(L 3 1) B.uncol++prop_tr :: Property+prop_tr = nudgeProp @(L 3 2) B.tr++prop_diag :: Property+prop_diag = nudgeProp @(R 3) B.diag++prop_svd :: Property+prop_svd = nudgeProp @(L 3 2) B.svd++prop_svd_ :: Property+prop_svd_ = nudgeProp @(L 3 2) ((\(_,x,_) -> x) . B.svd_)++prop_eigensystem1 :: Property+prop_eigensystem1 = nudgeProp @(L 3 2) (fst . B.eigensystem . B.mTm)++prop_eigensystem2 :: Property+prop_eigensystem2 = nudgeProp @(L 3 2) (snd . B.eigensystem . B.mTm)++prop_eigenvalues :: Property+prop_eigenvalues = nudgeProp @(L 3 2) (B.eigenvalues . B.mTm)++prop_chol :: Property+prop_chol = nudgeProp @(L 3 2) (B.chol . B.mTm)++prop_norm_0V :: Property+prop_norm_0V = nudgeProp @(R 3) B.norm_0++prop_norm_0M :: Property+prop_norm_0M = nudgeProp @(L 3 2) B.norm_0++prop_norm_1V :: Property+prop_norm_1V = nudgeProp @(R 3) B.norm_1V++prop_norm_1M :: Property+prop_norm_1M = nudgeProp @(L 3 2) B.norm_1M++prop_norm_2V :: Property+prop_norm_2V = nudgeProp @(R 3) B.norm_2V++prop_norm_2M :: Property+prop_norm_2M = nudgeProp @(L 3 2) B.norm_2M++prop_norm_InfV :: Property+prop_norm_InfV = nudgeProp @(R 3) B.norm_InfV++prop_norm_InfM :: Property+prop_norm_InfM = nudgeProp @(L 3 2) B.norm_InfM++prop_mean :: Property+prop_mean = nudgeProp @(R 3) B.mean++prop_meanCov1 :: Property+prop_meanCov1 = nudgeProp @(L 3 2) (fst . B.meanCov)++prop_meanCov2 :: Property+prop_meanCov2 = nudgeProp @(L 3 2) (B.unSym . snd . B.meanCov)++prop_meanL :: Property+prop_meanL = nudgeProp @(L 3 2) B.meanL++prop_cov :: Property+prop_cov = nudgeProp @(L 3 2) (B.unSym . B.cov)++prop_mul :: Property+prop_mul = nudgeProp2 @(L 3 2) @(L 2 3) B.mul++prop_app :: Property+prop_app = nudgeProp2 @(L 3 2) @(R 2) B.app++prop_dot :: Property+prop_dot = nudgeProp2 @(R 3) @(R 3) B.dot++prop_cross :: Property+prop_cross = nudgeProp2 @(R 3) B.cross++-- TODO: bug in diagR?+-- prop_diagR :: Property+-- prop_diagR = nudgeProp2 genDouble (genVec @3) (B.diagR @5 @4)++-- TODO: Mappers+-- , dvmap+-- , dvmap'+-- , dmmap+-- , dmmap'++prop_outer :: Property+prop_outer = nudgeProp2 @(R 3) @(R 2) B.outer++-- TODO: Zippers+-- , zipWithVector+-- , zipWithVector'++prop_det :: Property+prop_det = nudgeProp @(L 3 3) B.det++prop_invlndet1 :: Property+prop_invlndet1 = nudgeProp @(L 3 3) (fst . B.invlndet)++prop_invlndet2 :: Property+prop_invlndet2 = nudgeProp @(L 3 3) (fst . snd . B.invlndet)++prop_invlndet3 :: Property+prop_invlndet3 = nudgeProp @(L 3 3) (snd . snd . B.invlndet)++prop_lndet :: Property+prop_lndet = nudgeProp @(L 3 3) B.lndet++-- TODO: more general invertible matrix+prop_inv :: Property+prop_inv = nudgeProp @(L 3 2) (B.inv . B.unSym . B.mTm)++prop_toRows :: Property+prop_toRows = nudgeProp @(L 3 2) (collectVar . B.toRows)++prop_toColumns :: Property+prop_toColumns = nudgeProp @(L 2 3) (collectVar . B.toColumns)++prop_fromRows :: Property+prop_fromRows = nudgeProp (B.fromRows @3 @2 . sequenceVar)++prop_fromColumns :: Property+prop_fromColumns = nudgeProp (B.fromColumns @2 @3 . sequenceVar)++prop_konstV :: Property+prop_konstV = nudgeProp (B.konst @_ @(B.R 3))++prop_konstM :: Property+prop_konstM = nudgeProp (B.konst @_ @(B.L 3 2))++prop_sumElementsV :: Property+prop_sumElementsV = nudgeProp @(R 3) B.sumElements++prop_sumElementsM :: Property+prop_sumElementsM = nudgeProp @(L 3 2) B.sumElements++prop_extractV :: Property+prop_extractV = nudgeProp (B.extractV @_ @(R 3))++prop_extractM :: Property+prop_extractM = nudgeProp (B.extractM @_ @(L 3 2))++prop_createV :: Property+prop_createV = nudgeProp (fromMaybe 0 . B.create @_ @(R 3))++prop_createM :: Property+prop_createM = nudgeProp (fromMaybe 0 . B.create @_ @(L 3 2))++prop_takeDiag :: Property+prop_takeDiag = nudgeProp @(L 3 3) B.takeDiag++prop_sym :: Property+prop_sym = nudgeProp @(L 3 3) (B.unSym . B.sym)++prop_mTm :: Property+prop_mTm = nudgeProp @(L 3 2) (B.unSym . B.mTm)++prop_unSym :: Property+prop_unSym = nudgeProp @(L 3 3) (B.unSym . B.sym)++tryGroup :: (forall a. Num a => a) -> Group -> Group+tryGroup n Group{..} =+ Group groupName+ ((map . second) (withDiscards n . withTests n)+ groupProperties+ )++main :: IO ()+main = do+ hSetBuffering stdout LineBuffering+ hSetBuffering stderr LineBuffering++ results <- checkParallel (tryGroup 100 $$(discover))++ unless results exitFailure+