From d66c1929a587bbb9abe2870dabc19319938f4570 Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Fri, 1 Oct 2021 23:57:07 -0400
Subject: [PATCH] fix tests and bugs revealed by new tests

---
 math.scad            |  9 +++++----
 tests/test_math.scad | 11 +++++++++++
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/math.scad b/math.scad
index 3a9a818..bb5c67d 100644
--- a/math.scad
+++ b/math.scad
@@ -545,10 +545,12 @@ function gaussian_rands(N=1, mean=0, cov=1, seed=undef) =
          rdata = [for (i = count(dim*N,0,2)) sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])]
     )
     dim==1 ? add_scalar(sqrt(cov)*rdata,mean) :
+    assert(is_matrix_symmetric(cov),"Supplied covariance matrix is not symmetric")
     let(
         L = cholesky(cov)
     )
-    array_group(rdata,dim)*transpose(L);
+    assert(is_def(L), "Supplied covariance matrix is not positive definite")
+    move(mean,array_group(rdata,dim)*transpose(L));
 
 
 // Function: spherical_random_points()
@@ -1092,14 +1094,13 @@ function _back_substitute(R, b, x=[]) =
 function cholesky(A) =
   assert(is_matrix(A,square=true),"A must be a square matrix")
   assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix")
-  echo(A=A,len=len(A))
   _cholesky(A,ident(len(A)), len(A));
 
-function _cholesky(A,L,n) = let(ffee=echo(insideA=A,L,n))
+function _cholesky(A,L,n) = 
     A[0][0]<0 ? undef :     // Matrix not positive definite
     len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
     let(
-        i = n+1-len(A),ff=echo(i=i,lenA=len(A))
+        i = n+1-len(A)
     )
     let(
         sqrtAii = sqrt(A[0][0]),
diff --git a/tests/test_math.scad b/tests/test_math.scad
index 50adce3..a5a5843 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -368,6 +368,17 @@ module test_gaussian_rands() {
     assert_equal(len(nums3), 1000);
     assert_equal(nums1, nums3);
     assert(nums1!=nums2);
+
+    R = [[4,2],[2,17]];
+    data = gaussian_rands(100000,[0,0],R,seed=49);
+    assert(approx(mean(data), [0,0], eps=1e-2));
+    assert(approx(transpose(data)*data/len(data), R, eps=2e-2));
+    
+    R2 = [[4,2,-1],[2,17,4],[-1,4,11]];
+    data3 = gaussian_rands(100000,[1,2,3],R2,seed=97);
+    assert(approx(mean(data3),[1,2,3], eps=1e-2));
+    cdata = move(-mean(data3),data3);    
+    assert(approx(transpose(cdata)*cdata/len(cdata),R2,eps=.1));
 }
 test_gaussian_rands();