diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,6 @@ + +0.1 +=== + +- Harpie.Fixed orphans + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2016, Tony Day + +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 Tony Day 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. diff --git a/harpie-numhask.cabal b/harpie-numhask.cabal new file mode 100644 --- /dev/null +++ b/harpie-numhask.cabal @@ -0,0 +1,69 @@ +cabal-version: 3.0 +name: harpie-numhask +version: 0.1.0.0 +license: BSD-3-Clause +license-file: LICENSE +copyright: Tony Day (c) 2024 +category: math +author: Tony Day +maintainer: tonyday567@gmail.com +homepage: https://github.com/tonyday567/harpie-numhask#readme +bug-reports: https://github.com/tonyday567/harpie-numhask/issues +synopsis: numhask shim for harpie. +description: + This package provides orphans and interface into the [numhask](https://hackage.haskell.org/package/numhask) API for [harpie](https://hackage.haskell.org/package/harpie). + + == Usage + + >>> {-# LANGUAGE RebindableSyntax #-} + >>> import NumHask.Prelude + >>> import Harpie.Fixed as F + >>> import Harpie.Numhask qualified as N + +build-type: Simple +tested-with: + , GHC == 9.10.1 + , GHC == 9.6.5 + , GHC == 9.8.2 +extra-doc-files: ChangeLog.md + +source-repository head + type: git + location: https://github.com/tonyday567/harpie-numhask + +common ghc-options-stanza + ghc-options: + -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wpartial-fields + -Wredundant-constraints + +common ghc2021-stanza + default-language: GHC2021 + +library + import: ghc-options-stanza + import: ghc2021-stanza + hs-source-dirs: src + build-depends: + , adjunctions >=4.4.2 && <4.5 + , base >=4.14 && <5 + , first-class-families >=0.8 && <0.9 + , harpie >=0.1.0 && <0.2 + , numhask >=0.12 && <0.13 + exposed-modules: Harpie.NumHask + +test-suite doctests + import: ghc2021-stanza + main-is: doctests.hs + hs-source-dirs: test + build-depends: + , base >=4.14 && <5 + , doctest-parallel >=0.3 && <0.4 + , harpie-numhask + default-extensions: RebindableSyntax + ghc-options: -threaded + type: exitcode-stdio-1.0 diff --git a/src/Harpie/NumHask.hs b/src/Harpie/NumHask.hs new file mode 100644 --- /dev/null +++ b/src/Harpie/NumHask.hs @@ -0,0 +1,284 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE TypeFamilies #-} +{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-} +{-# OPTIONS_GHC -Wno-orphans #-} +{-# OPTIONS_GHC -Wno-redundant-constraints #-} + +-- | numhask orphan instances and shim for harpie. +module Harpie.NumHask + ( -- * Usage + -- $usage + ident, + undiag, + mult, + invtri, + inverse, + chol, + ) +where + +import Data.Functor.Rep +import Fcf qualified +import GHC.TypeNats +import Harpie.Fixed as F hiding (chol, ident, inverse, invtri, mult, undiag) +import Harpie.Shape (DeleteDims, Eval, GetDims, KnownNats, Rank, type (++)) +import Harpie.Shape qualified as S +import NumHask.Prelude as P hiding (Min, cycle, diff, drop, empty, find, length, repeat, sequence, take, zipWith) + +-- $setup +-- +-- >>> :m -Prelude +-- >>> :set -XDataKinds +-- >>> :set -XRebindableSyntax +-- >>> import NumHask.Prelude hiding (cycle, repeat, empty, diff, take, drop, zipWith, find) +-- >>> import Harpie.Fixed qualified as F +-- >>> import Harpie.Fixed (Array, array, range, shape, toDynamic) +-- >>> import Harpie.NumHask +-- >>> import Prettyprinter hiding (dot,fill) +-- +-- >>> s = 1 :: Array '[] Int +-- >>> s +-- [1] +-- >>> shape s +-- [] +-- >>> pretty s +-- 1 +-- >>> let v = range @'[3] +-- >>> pretty v +-- [0,1,2] +-- >>> let m = range @[2,3] +-- >>> pretty m +-- [[0,1,2], +-- [3,4,5]] +-- >>> a = range @[2,3,4] +-- >>> a +-- [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] +-- >>> pretty a +-- [[[0,1,2,3], +-- [4,5,6,7], +-- [8,9,10,11]], +-- [[12,13,14,15], +-- [16,17,18,19], +-- [20,21,22,23]]] +-- >>> toDynamic a +-- UnsafeArray [2,3,4] [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] + +-- * NumHask heirarchy + +instance + ( Additive a, + KnownNats s + ) => + Additive (Array s a) + where + (+) = liftR2 (+) + + zero = pureRep zero + +instance + ( Subtractive a, + KnownNats s + ) => + Subtractive (Array s a) + where + negate = fmapRep negate + +instance + (Multiplicative a) => + MultiplicativeAction (Array s a) + where + type Scalar (Array s a) = a + (|*) r s = fmap (s *) r + +instance (Additive a) => AdditiveAction (Array s a) where + type AdditiveScalar (Array s a) = a + (|+) r s = fmap (s +) r + +instance + (Subtractive a) => + SubtractiveAction (Array s a) + where + (|-) r s = fmap (\x -> x - s) r + +instance + (Divisive a) => + DivisiveAction (Array s a) + where + (|/) r s = fmap (/ s) r + +instance (KnownNats s, JoinSemiLattice a) => JoinSemiLattice (Array s a) where + (\/) = liftR2 (\/) + +instance (KnownNats s, MeetSemiLattice a) => MeetSemiLattice (Array s a) where + (/\) = liftR2 (/\) + +instance (KnownNats s, Subtractive a, Epsilon a) => Epsilon (Array s a) where + epsilon = konst epsilon + +instance (FromInteger a) => FromInteger (Array ('[] :: [Nat]) a) where + fromInteger x = toScalar (fromInteger x) + +instance (FromRational a) => FromRational (Array ('[] :: [Nat]) a) where + fromRational x = toScalar (fromRational x) + +-- | The identity array. +-- +-- >>> pretty $ ident @[3,3] +-- [[1,0,0], +-- [0,1,0], +-- [0,0,1]] +ident :: (KnownNats s, Additive a, Multiplicative a) => Array s a +ident = tabulate (bool zero one . S.isDiag . S.fromFins) + +-- | Expand the array to form a diagonal array +-- +-- >>> pretty $ undiag (range @'[3]) +-- [[0,0,0], +-- [0,1,0], +-- [0,0,2]] +undiag :: + forall s' a s. + ( KnownNats s, + KnownNats s', + s' ~ Eval ((++) s s), + Additive a + ) => + Array s a -> + Array s' a +undiag a = tabulate (\xs -> bool zero (index a (S.UnsafeFins $ pure $ S.getDim 0 (S.fromFins xs))) (S.isDiag (S.fromFins xs))) + +-- | Array multiplication. +-- +-- matrix multiplication +-- +-- >>> pretty $ mult m (F.transpose m) +-- [[5,14], +-- [14,50]] +-- +-- inner product +-- +-- >>> pretty $ mult v v +-- 5 +-- +-- matrix-vector multiplication +-- +-- >>> pretty $ mult v (F.transpose m) +-- [5,14] +-- +-- >>> pretty $ mult m v +-- [5,14] +mult :: + forall a ds0 ds1 s0 s1 so0 so1 st si. + ( Ring a, + KnownNats s0, + KnownNats s1, + KnownNats ds0, + KnownNats ds1, + KnownNats so0, + KnownNats so1, + KnownNats st, + KnownNats si, + so0 ~ Eval (DeleteDims ds0 s0), + so1 ~ Eval (DeleteDims ds1 s1), + si ~ Eval (GetDims ds0 s0), + si ~ Eval (GetDims ds1 s1), + st ~ Eval ((++) so0 so1), + ds0 ~ '[Eval ((Fcf.-) (Eval (Rank s0)) 1)], + ds1 ~ '[0] + ) => + Array s0 a -> + Array s1 a -> + Array st a +mult = dot sum (*) + +instance + ( Multiplicative a, + P.Distributive a, + Subtractive a, + KnownNat m + ) => + Multiplicative (Matrix m m a) + where + (*) = mult + + one = ident + +instance + ( Multiplicative a, + P.Distributive a, + Subtractive a, + Eq a, + ExpField a, + KnownNat m + ) => + Divisive (Matrix m m a) + where + recip a = invtri (transpose (chol a)) * invtri (chol a) + +-- | Inverse of a square matrix. +-- +-- > mult (inverse a) a == a +-- +-- >>> e = array @[3,3] @Double [4,12,-16,12,37,-43,-16,-43,98] +-- >>> pretty (inverse e) +-- [[49.36111111111111,-13.555555555555554,2.1111111111111107], +-- [-13.555555555555554,3.7777777777777772,-0.5555555555555555], +-- [2.1111111111111107,-0.5555555555555555,0.1111111111111111]] +inverse :: (Eq a, ExpField a, KnownNat m) => Matrix m m a -> Matrix m m a +inverse a = mult (invtri (transpose (chol a))) (invtri (chol a)) + +-- | [Inversion of a Triangular Matrix](https://math.stackexchange.com/questions/1003801/inverse-of-an-invertible-upper-triangular-matrix-of-order-3) +-- +-- >>> t = array @[3,3] @Double [1,0,1,0,1,2,0,0,1] +-- >>> pretty (invtri t) +-- [[1.0,0.0,-1.0], +-- [0.0,1.0,-2.0], +-- [0.0,0.0,1.0]] +-- +-- > ident == mult t (invtri t) +-- True +invtri :: forall a n. (KnownNat n, ExpField a, Eq a) => Matrix n n a -> Matrix n n a +invtri a = sum (fmap (l ^) (iota @n)) * ti + where + ti = undiag (fmap recip (diag a)) + tl = a - undiag (diag a) + l = negate (ti * tl) + +-- | cholesky decomposition +-- +-- Uses the algorithm. +-- +-- >>> e = array @[3,3] @Double [4,12,-16,12,37,-43,-16,-43,98] +-- >>> pretty (chol e) +-- [[2.0,0.0,0.0], +-- [6.0,1.0,0.0], +-- [-8.0,5.0,3.0]] +-- >>> mult (chol e) (F.transpose (chol e)) == e +-- True +chol :: (KnownNat m, ExpField a) => Matrix m m a -> Matrix m m a +chol a = + let l = + unsafeTabulate + ( \[i, j] -> + bool + ( one + / unsafeIndex l [j, j] + * ( unsafeIndex a [i, j] + - sum + ( (\k -> unsafeIndex l [i, k] * unsafeIndex l [j, k]) + <$> ([0 .. (j - 1)] :: [Int]) + ) + ) + ) + ( sqrt + ( unsafeIndex a [i, i] + - sum + ( (\k -> unsafeIndex l [j, k] ^ (2 :: Int)) + <$> ([0 .. (j - 1)] :: [Int]) + ) + ) + ) + (i == j) + ) + in l diff --git a/test/doctests.hs b/test/doctests.hs new file mode 100644 --- /dev/null +++ b/test/doctests.hs @@ -0,0 +1,8 @@ +module Main where + +import System.Environment (getArgs) +import Test.DocTest (mainFromCabal) +import Prelude (IO, (=<<)) + +main :: IO () +main = mainFromCabal "harpie-numhask" =<< getArgs