github.com/rmera/gochem@v0.7.1/matrixhelp.go (about) 1 /* 2 * matrixhelp.go, part of gochem. 3 * 4 * 5 * Copyright 2012 Raul Mera <rmera{at}chemDOThelsinkiDOTfi> 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU Lesser General Public License as 9 * published by the Free Software Foundation; either version 2.1 of the 10 * License, or (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General 18 * Public License along with this program. If not, see 19 * <http://www.gnu.org/licenses/>. 20 * 21 * 22 * Gochem is developed at the laboratory for instruction in Swedish, Department of Chemistry, 23 * University of Helsinki, Finland. 24 * 25 * 26 */ 27 28 package chem 29 30 //A munch of unexported mathematical functions, most of them just for convenience. 31 32 import ( 33 "fmt" 34 "math" 35 36 "github.com/rmera/gochem/v3" 37 "gonum.org/v1/gonum/mat" 38 ) 39 40 const appzero float64 = 0.000000000001 //used to correct floating point 41 //errors. Everything equal or less than this is considered zero. This probably sucks. 42 43 //Computes the inverse matrix of F and puts it in target. If target is nil, it alloactes 44 //a new one. Returns target. 45 func gnInverse(F, target *v3.Matrix) (*v3.Matrix, error) { 46 if target == nil { 47 r := F.NVecs() 48 target = v3.Zeros(r) 49 } 50 err := target.Dense.Inverse(F.Dense) 51 if err != nil { 52 err = CError{err.Error(), []string{"mat.Inverse", "gnInverse"}} 53 } 54 return target, err 55 } 56 57 //cross Takes 2 3-len column or row vectors and returns a column or a row 58 //vector, respectively, with the Cross product of them. 59 //should panic 60 func cross(a, b *v3.Matrix) *v3.Matrix { 61 c := v3.Zeros(1) 62 c.Cross(a, b) 63 return c 64 } 65 66 //invSqrt return the inverse of the square root of val, or zero if 67 //|val|<appzero. Returns -1 if val is negative 68 func invSqrt(val float64) float64 { 69 if math.Abs(val) <= appzero { //Not sure if need the 70 return 0 71 } else if val < 0 { //negative 72 panic("attempted to get the square root of a negative number") 73 } 74 return 1.0 / math.Sqrt(val) 75 } 76 77 //Returns and empty, but not nil, Dense. It barely allocates memory 78 func emptyDense() (*mat.Dense, error) { 79 a := make([]float64, 0, 0) 80 return mat.NewDense(0, 0, a), nil 81 82 } 83 84 //Returns an zero-filled Dense with the given dimensions 85 //It is to be substituted by the Gonum function. 86 func gnZeros(r, c int) *mat.Dense { 87 f := make([]float64, r*c, r*c) 88 return mat.NewDense(r, c, f) 89 90 } 91 92 //Returns an identity matrix spanning span cols and rows 93 func gnEye(span int) *mat.Dense { 94 A := gnZeros(span, span) 95 for i := 0; i < span; i++ { 96 A.Set(i, i, 1.0) 97 } 98 return A 99 } 100 101 //returns a rows,cols matrix filled with gnOnes. 102 func gnOnes(rows, cols int) *mat.Dense { 103 gnOnes := gnZeros(rows, cols) 104 for i := 0; i < rows; i++ { 105 for j := 0; j < cols; j++ { 106 gnOnes.Set(i, j, 1) 107 } 108 } 109 return gnOnes 110 } 111 112 //The 2 following functions may not even be used. 113 func gnMul(A, B mat.Matrix) *mat.Dense { 114 ar, _ := A.Dims() 115 _, bc := B.Dims() 116 C := gnZeros(ar, bc) 117 C.Mul(A, B) 118 return C 119 } 120 121 func gnCopy(A mat.Matrix) *mat.Dense { 122 r, c := A.Dims() 123 B := gnZeros(r, c) 124 B.Copy(A) 125 return B 126 } 127 128 func gnT(A mat.Matrix) *mat.Dense { 129 r, c := A.Dims() 130 B := gnZeros(c, r) 131 B.Copy(A.T()) 132 return B 133 } 134 135 //This is a temporal function. It returns the determinant of a 3x3 matrix. Panics if the matrix is not 3x3 136 func det(A mat.Matrix) float64 { 137 r, c := A.Dims() 138 if r != 3 || c != 3 { 139 panic("Determinants are for now only available for 3x3 matrices") 140 } 141 return (A.At(0, 0)*(A.At(1, 1)*A.At(2, 2)-A.At(2, 1)*A.At(1, 2)) - A.At(1, 0)*(A.At(0, 1)*A.At(2, 2)-A.At(2, 1)*A.At(0, 2)) + A.At(2, 0)*(A.At(0, 1)*A.At(1, 2)-A.At(1, 1)*A.At(0, 2))) 142 } 143 144 /**These are from the current proposal for gonum, by Dan Kortschak. It will be taken out 145 * from here when gonum is implemented. The gn prefix is appended to the names to make them 146 * unimported and to allow easy use of search/replace to add the "num" prefix when I change to 147 * gonum.**/ 148 149 // A gnPanicker is a function that may panic. 150 type gnPanicker func() 151 152 // Maybe will recover a panic with a type matrix.Error or a VecError from fn, and return this error. 153 // Any other error is re-panicked. It is a small modification 154 //Maybe this funciton should be exported. 155 func gnMaybe(fn gnPanicker) error { 156 var err error 157 defer func() { 158 if r := recover(); r != nil { 159 switch e := r.(type) { 160 case v3.Error: 161 err = e 162 case mat.Error: 163 err = CError{fmt.Sprintf("goChem: Error in gonum function: %s", e), []string{"gnMaybe"}} 164 default: 165 panic(r) 166 } 167 } 168 }() 169 fn() 170 return err 171 } 172 173 //Puts A**exp on the receiver, in a pretty naive way. 174 func pow(A mat.Matrix, F *mat.Dense, exp float64) { 175 ar, ac := A.Dims() 176 fr, fc := F.Dims() 177 if ar != fr || ac != fc { 178 panic(mat.ErrShape) 179 } 180 for i := 0; i < ar; i++ { 181 for j := 0; j < ac; j++ { 182 F.Set(i, j, math.Pow(A.At(i, j), exp)) 183 } 184 185 } 186 }