From b02cdc9e437cf00076a70565a1a804444d6489ea Mon Sep 17 00:00:00 2001 From: Geod24 Date: Sat, 11 Jun 2022 17:44:43 +0200 Subject: [PATCH 1/3] Trivial: Extend gitignore --- .gitignore | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index dbaf50c..4319c4c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ .dub dub.selections.json __* -*.exe \ No newline at end of file +*.exe +*.[oa] +*-test-library From ee965e8b5b73ce926832980cfa2242948f579bd1 Mon Sep 17 00:00:00 2001 From: Geod24 Date: Wed, 15 Jun 2022 07:42:40 +0200 Subject: [PATCH 2/3] Trivial: Make an assert error message readable Before naiveAnd and fastAns were not separated, leading to results such as 10.864943 instead of '1 - 0.864943'. --- source/dstats/tests.d | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/dstats/tests.d b/source/dstats/tests.d index 5c0528b..c61a7ee 100644 --- a/source/dstats/tests.d +++ b/source/dstats/tests.d @@ -2950,7 +2950,7 @@ unittest { c[1][1] = uniform(0U, 51U); double naiveAns = naive(c); double fastAns = fisherExact(c); - assert(approxEqual(naiveAns, fastAns), text(c, naiveAns, fastAns)); + assert(approxEqual(naiveAns, fastAns), text(c, " - ", naiveAns, " - ", fastAns)); } auto res = fisherExact([[19000, 80000], [20000, 90000]]); From 724cd066089c9e3443105a6968cc0d4a832a3740 Mon Sep 17 00:00:00 2001 From: Geod24 Date: Wed, 15 Jun 2022 13:43:03 +0200 Subject: [PATCH 3/3] Trivial: Remove trailing \r --- source/dstats/base.d | 518 +++++++++++++++++++++---------------------- 1 file changed, 259 insertions(+), 259 deletions(-) diff --git a/source/dstats/base.d b/source/dstats/base.d index 783771a..79f2c77 100644 --- a/source/dstats/base.d +++ b/source/dstats/base.d @@ -1493,265 +1493,265 @@ unittest { assert(rng.empty); myFile.close(); } -} - -// Used for ILP optimizations. -package template StaticIota(size_t upTo) { - static if(upTo == 0) { - alias TypeTuple!() StaticIota; - } else { - alias TypeTuple!(StaticIota!(upTo - 1), upTo - 1) StaticIota; - } +} + +// Used for ILP optimizations. +package template StaticIota(size_t upTo) { + static if(upTo == 0) { + alias TypeTuple!() StaticIota; + } else { + alias TypeTuple!(StaticIota!(upTo - 1), upTo - 1) StaticIota; + } } package: - -// If SciD is available, it is used for matrix storage and operations -// such as Cholesky decomposition. Otherwise, some old, crappy routines -// that were around since before SciD are used. -version(scid) { - import scid.matvec, scid.linalg, scid.exception; - - alias ExternalMatrixView!double DoubleMatrix; - - DoubleMatrix doubleMatrix( - size_t nRows, - size_t nColumns, - RegionAllocator alloc - ) { - return DoubleMatrix(nRows, nColumns, alloc); - } - - void choleskySolve(DoubleMatrix a, double[] b, double[] x) { - choleskyDestructive(a); - auto ans = externalVectorView(x); - auto vec = externalVectorView(b); - scid.linalg.choleskySolve(a, vec, ans); - } - - void invert(DoubleMatrix from, DoubleMatrix to) { - try { - to[] = inv(from); - } catch(Exception) { - foreach(i; 0..to.rows) foreach(j; 0..to.columns) { - to[i, j] = double.nan; - } - } - } - -} else { - version = noscid; - - struct DoubleMatrix { - double[][] arrayOfArrays; - //alias arrayOfArraysFun this; - - const pure nothrow @property { - size_t rows() { - return arrayOfArrays.length; - } - - size_t columns() { - return (arrayOfArrays.length) ? - (arrayOfArrays[0].length) : 0; - } - } - - ref double opIndex(size_t i, size_t j) { - return arrayOfArrays[i][j]; - } - } - - DoubleMatrix doubleMatrix( - size_t nRows, - size_t nColumns, - RegionAllocator alloc - ) { - return DoubleMatrix( - alloc.uninitializedArray!(double[][])(nRows, nColumns) - ); - } - - void invert(DoubleMatrix from, DoubleMatrix to) { - invert(from.arrayOfArrays, to.arrayOfArrays); - } - - // Uses Gauss-Jordan elim. w/ row pivoting to invert from. Stores the results - // in to and leaves from in an undefined state. - package void invert(double[][] from, double[][] to) { - // Normalize. - foreach(i, row; from) { - double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..from.length])); - row[] *= absMax; - to[i][] = 0; - to[i][i] = absMax; - } - - foreach(col; 0..from.length) { - size_t bestRow; - double biggest = 0; - foreach(row; col..from.length) { - if(abs(from[row][col]) > biggest) { - bestRow = row; - biggest = abs(from[row][col]); - } - } - - swap(from[col], from[bestRow]); - swap(to[col], to[bestRow]); - immutable pivotFactor = from[col][col]; - - foreach(row; 0..from.length) if(row != col) { - immutable ratio = from[row][col] / pivotFactor; - - // If you're ever looking to optimize this code, good luck. The - // bottleneck is almost ENTIRELY this one line: - from[row][] -= from[col][] * ratio; - to[row][] -= to[col][] * ratio; - } - } - - foreach(i; 0..from.length) { - immutable diagVal = from[i][i]; - from[i][] /= diagVal; - to[i][] /= diagVal; - } - } - - unittest { - auto mat = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]]; - auto toMat = [new double[3], new double[3], new double[3]]; - invert(mat, toMat); - assert(approxEqual(toMat[0], [-0.625, -0.25, 0.375])); - assert(approxEqual(toMat[1], [-0.25, 0.5, -0.25])); - assert(approxEqual(toMat[2], [0.708333, -0.25, 0.041667])); - } - - void solve(DoubleMatrix mat, double[] vec) { - solve(mat.arrayOfArrays, vec); - } - - // Solve a system of linear equations mat * b = vec for b. The result is - // stored in vec, and mat is left in an undefined state. Uses Gaussian - // elimination w/ row pivoting. Roughly 4x faster than using inversion to - // solve the system, and uses roughly half the memory. - void solve(double[][] mat, double[] vec) { - // Normalize. - foreach(i, row; mat) { - double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..mat.length])); - row[] *= absMax; - vec[i] *= absMax; - } - - foreach(col; 0..mat.length - 1) { - size_t bestRow; - double biggest = 0; - foreach(row; col..mat.length) { - if(abs(mat[row][col]) > biggest) { - bestRow = row; - biggest = abs(mat[row][col]); - } - } - - swap(mat[col], mat[bestRow]); - swap(vec[col], vec[bestRow]); - immutable pivotFactor = mat[col][col]; - - foreach(row; col + 1..mat.length) { - immutable ratio = mat[row][col] / pivotFactor; - - // If you're ever looking to optimize this code, good luck. The - // bottleneck is almost ENTIRELY this one line: - mat[row][col..$] -= mat[col][col..$] * ratio; - vec[row] -= vec[col] * ratio; - } - } - - foreach(i; 0..mat.length) { - double diagVal = mat[i][i]; - mat[i][] /= diagVal; - vec[i] /= diagVal; - } - - // Do back substitution. - for(size_t row = mat.length - 1; row != size_t.max; row--) { - auto v1 = vec[row + 1..$]; - auto v2 = mat[row][$ - v1.length..$]; - vec[row] -= dotProduct(v1, v2); - } - } - - unittest { - auto mat = [[2.0, 1, -1], [-3.0, -1, 2], [-2.0, 1, 2]]; - auto vec = [8.0, -11, -3]; - solve(mat, vec); - assert(approxEqual(vec, [2, 3, -1])); - - auto mat2 = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]]; - auto vec2 = [8.0, 6, 7]; - solve(mat2, vec2); - assert(approxEqual(vec2, [-3.875, -0.75, 4.45833])); - } - - // Cholesky decomposition functions adapted from Don Clugston's MathExtra - // lib, used w/ permission. - void choleskyDecompose(double[][] a, double[] diag) { - immutable N = diag.length; - - foreach(i; 0..N) { - const ai = a[i]; - double sum = ai[i]; - - for(sizediff_t k = i - 1; k >= 0; --k) { - sum -= ai[k] * ai[k]; - } - - if (sum > 0.0) { - diag[i] = sqrt(sum); - - foreach(j; i + 1..N) { - sum = ai[j] - dotProduct(ai[0..i], a[j][0..i]); - a[j][i] = sum / diag[i]; - } - } else { - // not positive definite (could be caused by rounding errors) - diag[i] = 0; - // make this whole row zero so they have no further effect - foreach(j; i + 1..N) a[j][i] = 0; - } - } - } - - void choleskySolve(double[][] a, double[] diag, double[] b, double[] x) { - immutable N = x.length; - - foreach(i; 0..N) { - const ai = a[i]; - - if(diag[i] > 0) { - double sum = b[i]; - sum -= dotProduct(ai[0..i], x[0..i]); - x[i] = sum / diag[i]; - } else x[i] = 0; // skip pos definite rows - } - - for(sizediff_t i = N - 1; i >= 0; --i) { - if (diag[i] > 0) { - double sum = x[i]; - for(sizediff_t k = i + 1; k < N; ++k) sum -= a[k][i] * x[k]; - x[i] = sum / diag[i]; - } else x[i] = 0; // skip pos definite rows - } - // Convert failed columns in solution to NANs if required. - foreach(i; 0..N) { - if(diag[i].isNaN() || diag[i] <= 0) x[i] = double.nan; - } - } - - void choleskySolve(DoubleMatrix a, double[] b, double[] x) { - auto alloc = newRegionAllocator(); - auto diag = alloc.uninitializedArray!(double[])(x.length); - choleskyDecompose(a.arrayOfArrays, diag); - choleskySolve(a.arrayOfArrays, diag, b, x); - } -} + +// If SciD is available, it is used for matrix storage and operations +// such as Cholesky decomposition. Otherwise, some old, crappy routines +// that were around since before SciD are used. +version(scid) { + import scid.matvec, scid.linalg, scid.exception; + + alias ExternalMatrixView!double DoubleMatrix; + + DoubleMatrix doubleMatrix( + size_t nRows, + size_t nColumns, + RegionAllocator alloc + ) { + return DoubleMatrix(nRows, nColumns, alloc); + } + + void choleskySolve(DoubleMatrix a, double[] b, double[] x) { + choleskyDestructive(a); + auto ans = externalVectorView(x); + auto vec = externalVectorView(b); + scid.linalg.choleskySolve(a, vec, ans); + } + + void invert(DoubleMatrix from, DoubleMatrix to) { + try { + to[] = inv(from); + } catch(Exception) { + foreach(i; 0..to.rows) foreach(j; 0..to.columns) { + to[i, j] = double.nan; + } + } + } + +} else { + version = noscid; + + struct DoubleMatrix { + double[][] arrayOfArrays; + //alias arrayOfArraysFun this; + + const pure nothrow @property { + size_t rows() { + return arrayOfArrays.length; + } + + size_t columns() { + return (arrayOfArrays.length) ? + (arrayOfArrays[0].length) : 0; + } + } + + ref double opIndex(size_t i, size_t j) { + return arrayOfArrays[i][j]; + } + } + + DoubleMatrix doubleMatrix( + size_t nRows, + size_t nColumns, + RegionAllocator alloc + ) { + return DoubleMatrix( + alloc.uninitializedArray!(double[][])(nRows, nColumns) + ); + } + + void invert(DoubleMatrix from, DoubleMatrix to) { + invert(from.arrayOfArrays, to.arrayOfArrays); + } + + // Uses Gauss-Jordan elim. w/ row pivoting to invert from. Stores the results + // in to and leaves from in an undefined state. + package void invert(double[][] from, double[][] to) { + // Normalize. + foreach(i, row; from) { + double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..from.length])); + row[] *= absMax; + to[i][] = 0; + to[i][i] = absMax; + } + + foreach(col; 0..from.length) { + size_t bestRow; + double biggest = 0; + foreach(row; col..from.length) { + if(abs(from[row][col]) > biggest) { + bestRow = row; + biggest = abs(from[row][col]); + } + } + + swap(from[col], from[bestRow]); + swap(to[col], to[bestRow]); + immutable pivotFactor = from[col][col]; + + foreach(row; 0..from.length) if(row != col) { + immutable ratio = from[row][col] / pivotFactor; + + // If you're ever looking to optimize this code, good luck. The + // bottleneck is almost ENTIRELY this one line: + from[row][] -= from[col][] * ratio; + to[row][] -= to[col][] * ratio; + } + } + + foreach(i; 0..from.length) { + immutable diagVal = from[i][i]; + from[i][] /= diagVal; + to[i][] /= diagVal; + } + } + + unittest { + auto mat = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]]; + auto toMat = [new double[3], new double[3], new double[3]]; + invert(mat, toMat); + assert(approxEqual(toMat[0], [-0.625, -0.25, 0.375])); + assert(approxEqual(toMat[1], [-0.25, 0.5, -0.25])); + assert(approxEqual(toMat[2], [0.708333, -0.25, 0.041667])); + } + + void solve(DoubleMatrix mat, double[] vec) { + solve(mat.arrayOfArrays, vec); + } + + // Solve a system of linear equations mat * b = vec for b. The result is + // stored in vec, and mat is left in an undefined state. Uses Gaussian + // elimination w/ row pivoting. Roughly 4x faster than using inversion to + // solve the system, and uses roughly half the memory. + void solve(double[][] mat, double[] vec) { + // Normalize. + foreach(i, row; mat) { + double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..mat.length])); + row[] *= absMax; + vec[i] *= absMax; + } + + foreach(col; 0..mat.length - 1) { + size_t bestRow; + double biggest = 0; + foreach(row; col..mat.length) { + if(abs(mat[row][col]) > biggest) { + bestRow = row; + biggest = abs(mat[row][col]); + } + } + + swap(mat[col], mat[bestRow]); + swap(vec[col], vec[bestRow]); + immutable pivotFactor = mat[col][col]; + + foreach(row; col + 1..mat.length) { + immutable ratio = mat[row][col] / pivotFactor; + + // If you're ever looking to optimize this code, good luck. The + // bottleneck is almost ENTIRELY this one line: + mat[row][col..$] -= mat[col][col..$] * ratio; + vec[row] -= vec[col] * ratio; + } + } + + foreach(i; 0..mat.length) { + double diagVal = mat[i][i]; + mat[i][] /= diagVal; + vec[i] /= diagVal; + } + + // Do back substitution. + for(size_t row = mat.length - 1; row != size_t.max; row--) { + auto v1 = vec[row + 1..$]; + auto v2 = mat[row][$ - v1.length..$]; + vec[row] -= dotProduct(v1, v2); + } + } + + unittest { + auto mat = [[2.0, 1, -1], [-3.0, -1, 2], [-2.0, 1, 2]]; + auto vec = [8.0, -11, -3]; + solve(mat, vec); + assert(approxEqual(vec, [2, 3, -1])); + + auto mat2 = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]]; + auto vec2 = [8.0, 6, 7]; + solve(mat2, vec2); + assert(approxEqual(vec2, [-3.875, -0.75, 4.45833])); + } + + // Cholesky decomposition functions adapted from Don Clugston's MathExtra + // lib, used w/ permission. + void choleskyDecompose(double[][] a, double[] diag) { + immutable N = diag.length; + + foreach(i; 0..N) { + const ai = a[i]; + double sum = ai[i]; + + for(sizediff_t k = i - 1; k >= 0; --k) { + sum -= ai[k] * ai[k]; + } + + if (sum > 0.0) { + diag[i] = sqrt(sum); + + foreach(j; i + 1..N) { + sum = ai[j] - dotProduct(ai[0..i], a[j][0..i]); + a[j][i] = sum / diag[i]; + } + } else { + // not positive definite (could be caused by rounding errors) + diag[i] = 0; + // make this whole row zero so they have no further effect + foreach(j; i + 1..N) a[j][i] = 0; + } + } + } + + void choleskySolve(double[][] a, double[] diag, double[] b, double[] x) { + immutable N = x.length; + + foreach(i; 0..N) { + const ai = a[i]; + + if(diag[i] > 0) { + double sum = b[i]; + sum -= dotProduct(ai[0..i], x[0..i]); + x[i] = sum / diag[i]; + } else x[i] = 0; // skip pos definite rows + } + + for(sizediff_t i = N - 1; i >= 0; --i) { + if (diag[i] > 0) { + double sum = x[i]; + for(sizediff_t k = i + 1; k < N; ++k) sum -= a[k][i] * x[k]; + x[i] = sum / diag[i]; + } else x[i] = 0; // skip pos definite rows + } + // Convert failed columns in solution to NANs if required. + foreach(i; 0..N) { + if(diag[i].isNaN() || diag[i] <= 0) x[i] = double.nan; + } + } + + void choleskySolve(DoubleMatrix a, double[] b, double[] x) { + auto alloc = newRegionAllocator(); + auto diag = alloc.uninitializedArray!(double[])(x.length); + choleskyDecompose(a.arrayOfArrays, diag); + choleskySolve(a.arrayOfArrays, diag, b, x); + } +}