factor: trim using lists

[factor.git] / basis / math / matrices / matrices.factor

diff --git a/basis/math/matrices/matrices.factor b/basis/math/matrices/matrices.factor
index 7af1b77210c76dd12cd964e7481a279aa382dcf1..eeddbc06cd0222434eba17da402d2d74b0a8ceb3 100644 (file)
--- a/basis/math/matrices/matrices.factor
+++ b/basis/math/matrices/matrices.factor
@@ -1,11 +1,8 @@
-! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
+! Copyright (C) 2005, 2010, 2018, 2020 Slava Pestov, Joe Groff, and Cat Stevens.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: accessors arrays classes.singleton columns combinators
-combinators.short-circuit combinators.smart formatting fry
-grouping kernel locals math math.bits math.functions math.order
-math.private math.ranges math.statistics math.vectors
-math.vectors.private sequences sequences.deep sequences.private
-slots.private summary ;
+USING: arrays combinators combinators.short-circuit kernel math
+math.functions math.order math.private math.vectors ranges
+sequences sequences.deep sequences.private slots.private ;
 IN: math.matrices
 
 ! defined here because of issue #1943
@@ -190,7 +187,7 @@ M: matrix <square-cols>
 
 <PRIVATE ! implementation details of <lower-matrix> and <upper-matrix>
 : dimension-range ( matrix -- dim range )
-    dimension [ <coordinate-matrix> ] [ first [1,b] ] bi ;
+    dimension [ <coordinate-matrix> ] [ first [1..b] ] bi ;
 
 : upper-matrix-indices ( matrix -- matrix' )
     dimension-range <reversed> [ tail-slice* >array ] 2map concat ;
@@ -214,6 +211,9 @@ DEFER: matrix-set-nths
 : matrix-map ( matrix quot: ( ... elt -- ... elt' ) -- matrix' )
     '[ _ map ] map ; inline
 
+: matrix-map-index ( matrix quot: ( ... elt i j -- ... elt' ) -- matrix' )
+    '[ [ swap @ ] curry map-index ] map-index ; inline
+
 : column-map ( matrix quot: ( ... col -- ... col' ) -- matrix' )
     [ transpose ] dip map transpose ; inline
 
@@ -248,18 +248,70 @@ DEFER: matrix-set-nths
 : m*  ( m1 m2 -- m ) [ v* ] 2map ;
 : m/  ( m1 m2 -- m ) [ v/ ] 2map ;
 
-: v.m ( v m -- p ) flip [ vdot ] with map ;
-: m.v ( m v -- p ) [ vdot ] curry map ;
-: m. ( m m -- m ) flip [ swap m.v ] curry map ;
+: vdotm ( v m -- p ) flip [ vdot ] with map ;
+: mdotv ( m v -- p ) [ vdot ] curry map ;
+: mdot ( m m -- m ) flip [ swap mdotv ] curry map ;
 
 : m~  ( m1 m2 epsilon -- ? ) [ v~ ] curry 2all? ;
 
 : mmin ( m -- n ) [ 1/0. ] dip [ [ min ] each ] each ;
 : mmax ( m -- n ) [ -1/0. ] dip [ [ max ] each ] each ;
-: mnorm ( m -- m' ) dup mmax abs m/n ;
-: m-infinity-norm ( m -- n ) [ [ abs ] map-sum ] map supremum ;
-: m-1norm ( m -- n ) flip m-infinity-norm ;
-: frobenius-norm ( m -- n ) [ [ sq ] map-sum ] map-sum sqrt ;
+
+: matrix-l-infinity-norm ( m -- n )
+    dup zero-matrix? [ drop 0 ] [
+        [ [ abs ] map-sum ] map supremum
+    ] if ; inline foldable
+
+: matrix-l1-norm ( m -- n )
+    dup zero-matrix? [ drop 0 ] [
+        flip matrix-l-infinity-norm
+    ] if ; inline foldable
+
+: matrix-l2-norm ( m -- n )
+    dup zero-matrix? [ drop 0 ] [
+        [ [ sq ] map-sum ] map-sum sqrt
+    ] if ; inline foldable
+
+! XXX: M: zero-matrix l1-norm drop 0 ; inline
+! XXX: M: matrix l1-norm matrix-l1-norm ; inline
+
+! XXX: M: zero-matrix l2-norm drop 0 ; inline
+! XXX: M: matrix l2-norm matrix-l2-norm ; inline
+
+! XXX: M: zero-matrix l-infinity-norm drop 0 ; inline
+! XXX: M: matrix l-infinity-norm matrix-l-infinity-norm ; inline
+
+ALIAS: frobenius-norm matrix-l2-norm
+ALIAS: hilbert-schmidt-norm matrix-l2-norm
+
+:: matrix-p-q-norm ( m p q -- n )
+    m dup zero-matrix? [ drop 0 ] [
+        [ [ sq ] map-sum q p / ^ ] map-sum q recip ^
+    ] if ; inline foldable
+
+: matrix-p-norm-entrywise ( m p -- n )
+    [ flatten1 V{ } like ] dip p-norm-default ; inline
+
+! XXX: M: zero-matrix p-norm-default 2drop 0 ; inline
+! XXX: M: matrix p-norm-default matrix-p-norm-entrywise ; inline
+
+: matrix-p-norm ( m p -- n )
+    over zero-matrix? [ 2drop 0 ] [
+        {
+            { [ dup 1 number= ] [ drop matrix-l1-norm ] }
+            { [ dup 2 number= ] [ drop matrix-l2-norm ] }
+            { [ dup fp-infinity? ] [ drop matrix-l-infinity-norm ] }
+            [ matrix-p-norm-entrywise ]
+        } cond
+    ] if ; inline foldable
+
+! XXX: M: zero-matrix p-norm 2drop 0 ; inline
+! XXX: M: matrix p-norm matrix-p-norm ; inline
+
+: matrix-normalize ( m -- m' )
+    dup zero-matrix? [
+        dup mabs mmax m/n
+    ] unless ; inline foldable
 
 ! well-defined for square matrices; but works on nonsquare too
 : main-diagonal ( matrix -- seq )