1 /***********************************************************************
2  * Copyright (c) 2013-2024 Commonwealth Computer Research, Inc.
3  * All rights reserved. This program and the accompanying materials
4  * are made available under the terms of the Apache License, Version 2.0
5  * which accompanies this distribution and is available at
6  * http://www.opensource.org/licenses/apache2.0.php.
7  ***********************************************************************/
8 
9 package org.locationtech.geomesa.utils.stats
10 
11 import org.ejml.data.DMatrixRMaj
12 import org.ejml.dense.row.CommonOps_DDRM
13 import org.ejml.simple.SimpleMatrix
14 
15 object SimpleMatrixUtils {
16 
17   private implicit def toDMatrixRMaj(sm: SimpleMatrix): DMatrixRMaj = sm.getMatrix[DMatrixRMaj]
18 
19   implicit class SimpleMatrixOps(a: SimpleMatrix) {
20 
21     def +(b: Double): SimpleMatrix = a.plus(b)
22     def +(b: SimpleMatrix): SimpleMatrix = a.plus(b)
23 
24     def +=(b: Double): Unit = CommonOps_DDRM.add[DMatrixRMaj](a, b, a)
25     def +=(b: SimpleMatrix): Unit = CommonOps_DDRM.add[DMatrixRMaj](a, b, a)
26 
27     def -(b: Double): SimpleMatrix = a.minus(b)
28     def -(b: SimpleMatrix): SimpleMatrix = a.minus(b)
29 
30     def -=(b: Double): Unit = CommonOps_DDRM.subtract[DMatrixRMaj](a, b, a)
31     def -=(b: SimpleMatrix): Unit = CommonOps_DDRM.subtract[DMatrixRMaj](a, b, a)
32 
33     def *(b: Double): SimpleMatrix = a.scale(b)
34     def *(b: SimpleMatrix): SimpleMatrix = a.elementMult(b)
35 
36     def *=(b: Double): Unit = CommonOps_DDRM.scale(b, a, a)
37     def *=(b: SimpleMatrix): Unit = CommonOps_DDRM.elementMult[DMatrixRMaj](a, b, a)
38 
39     def /(b: Double): SimpleMatrix = a.divide(b)
40     def /(b: SimpleMatrix): SimpleMatrix = a.elementDiv(b)
41 
42     def /=(b: Double): Unit = CommonOps_DDRM.divide[DMatrixRMaj](a, b, a)
43     def /=(b: SimpleMatrix): Unit = CommonOps_DDRM.elementDiv[DMatrixRMaj](a, b, a)
44 
45     def **(b: Double): SimpleMatrix = a.elementPower(b)
46     def **(b: SimpleMatrix): SimpleMatrix = a.elementPower(b)
47 
48     def **=(b: Double): Unit = CommonOps_DDRM.elementPower[DMatrixRMaj](a, b, a)
49     def **=(b: SimpleMatrix): Unit = CommonOps_DDRM.elementPower[DMatrixRMaj](a, b, a)
50 
51     def diag(v: Double): SimpleMatrix = {
52       val m = new SimpleMatrix(a)
53       (0 until Math.min(m.getNumRows, m.getNumCols)).foreach(i => m.set(i, i, v))
54       m
55     }
56 
57     def |*|(b: SimpleMatrix): SimpleMatrix = a.mult(b)
58 
59     def T: SimpleMatrix = a.transpose
60 
61     def isIdenticalWithinTolerances(b: SimpleMatrix, rel_tol: Double = 1e-9, abs_tol: Double = 1e-15): Boolean = {
62       import java.lang.{Double => jDouble}
63       if (a.numRows != b.numRows || a.numCols != b.numCols) {
64         return false
65       }
66       require(rel_tol >= 0 && abs_tol >=0, "Tolerance must be greater than or equal to zero.")
67       val length = a.getNumElements
68       var i = 0
69 
70       while (i < length) {
71         val va = a.get(i)
72         val vb = b.get(i)
73         val va_nan = va != va     /* quick NaN test */
74         val vb_nan = vb != vb     /* quick NaN test */
75         if (va_nan || vb_nan) {   /* if either NaN */
76           if (va_nan != vb_nan) { /* then both should be NaN */
77             return false
78           }
79         } else {
80           val va_inf = jDouble.isInfinite(va)
81           val vb_inf = jDouble.isInfinite(vb)
82           if (va_inf || vb_inf) {
83             if (va != vb) {
84               return false
85             }
86           } else {
87             /* check absolute tolerance, important for low magnitude values (0) */
88             val diff = Math.abs(va - vb)
89             if (diff > abs_tol) {
90               val va_abs = Math.abs(va)
91               val vb_abs = Math.abs(vb)
92               /* check relative tolerance, important for high magnitude values */
93               if (diff > Math.max(va_abs, vb_abs) * rel_tol) {
94                 return false
95               }
96             }
97           }
98         }
99         i += 1
100       }
101       true
102     }
103   }
104 
105   implicit class DoubleOps(a: Double) {
106     def +(b: SimpleMatrix): SimpleMatrix = b.plus(a)
107     def -(b: SimpleMatrix): SimpleMatrix = {
108       val c = new SimpleMatrix(b.getNumRows, b.getNumCols)
109       CommonOps_DDRM.subtract[DMatrixRMaj](a, b, c)
110       c
111     }
112     def *(b: SimpleMatrix): SimpleMatrix = b.scale(a)
113     def /(b: SimpleMatrix): SimpleMatrix = {
114       val c = new SimpleMatrix(b.getNumRows, b.getNumCols)
115       CommonOps_DDRM.divide[DMatrixRMaj](a, b, c)
116       c
117     }
118     def **(b: SimpleMatrix): SimpleMatrix = {
119       val c = new SimpleMatrix(b.getNumRows, b.getNumCols)
120       CommonOps_DDRM.elementPower[DMatrixRMaj](a, b, c)
121       c
122     }
123   }
124 }