Skip to content

Add mahalanobisDistance #107

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add mahalanobisDistance #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 23 additions & 0 deletions LICENSE_BOOST
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
Boost Software License - Version 1.0 - August 17th, 2003

Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:

The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
65 changes: 23 additions & 42 deletions dub.sdl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,109 +60,90 @@ configuration "dips" {
buildType "unittest-blas" {
buildOptions "unittests" "debugMode" "debugInfo"
versions "mir_stat_test_blas"
dflags "-lowmem"
dflags "-lowmem" "-preview=dip1000"
}
buildType "unittest-cov-blas" {
buildOptions "unittests" "coverage" "debugMode" "debugInfo"
versions "mir_stat_test_blas"
dflags "-lowmem"
}
buildType "unittest-blas-perf" {
buildOptions "unittests" "releaseMode" "optimize" "inline"
dflags "-mcpu=native" "-preview=dip1000"
versions "mir_stat_test_blas_mahalanobis_performance"
}
configuration "openblas" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "openblas"
}

configuration "threelib" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "threelib"
}

configuration "cblas" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "cblas"
}

configuration "blas" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "blas"
}

configuration "lapack" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "lapack"
}

configuration "mkl-sequential" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-sequential"
}

configuration "mkl-sequential-ilp" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-sequential-ilp"
}

configuration "mkl-tbb-thread" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-tbb-thread"
}

configuration "mkl-tbb-thread-ilp" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-tbb-thread-ilp"
}

configuration "mkl-sequential-dll" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-sequential-dll"
}

configuration "mkl-sequential-ilp-dll" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-sequential-ilp-dll"
}

configuration "mkl-tbb-thread-dll" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-tbb-thread-dll"
}

configuration "mkl-tbb-thread-ilp-dll" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-tbb-thread-ilp-dll"
}

configuration "zerolib" {
dependency "mir-lapack" version=">=1.2.9"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "zerolib"
}

configuration "unittest-openblas" {
dependency "mir-lapack" version=">=1.2.9"
subConfiguration "mir-lapack" "openblas"
}

configuration "unittest-threelib" {
dependency "mir-lapack" version=">=1.2.9"
subConfiguration "mir-lapack" "threelib"
}

configuration "unittest-cblas" {
dependency "mir-lapack" version=">=1.2.9"
subConfiguration "mir-lapack" "cblas"
}

configuration "unittest-blas" {
dependency "mir-lapack" version=">=1.2.9"
subConfiguration "mir-lapack" "blas"
}

configuration "unittest-lapack" {
dependency "mir-lapack" version=">=1.2.9"
subConfiguration "mir-lapack" "lapack"
}

configuration "unittest-mkl-sequential" {
dependency "mir-lapack" version=">=1.2.9"
configuration "unittest-blas-perf-mkl-sequential" {
dependency "mir-random" version=">=2.2.19"
dependency "mir-lapack" version=">=1.2.12"
subConfiguration "mir-lapack" "mkl-sequential"
}
6 changes: 5 additions & 1 deletion meson.build
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,14 @@ endif
has_cpp_headers = false

sources_list = [
'mir/math/internal/benchmark',
'mir/math/internal/fp_powi',
'mir/math/internal/powi',
'mir/math/internal/linearAlgebra',
'mir/math/internal/log_beta',
'mir/math/internal/log_binomial',
'mir/math/internal/log1p',
'mir/math/internal/lubeck2',
'mir/math/internal/powi',
'mir/math/internal/xlogy',
'mir/stat/descriptive/multivariate',
'mir/stat/descriptive/univariate',
Expand Down
122 changes: 107 additions & 15 deletions source/mir/math/internal/benchmark.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,36 +6,37 @@ import std.traits: isMutable;
package(mir)
template benchmarkValues(fun...)
{
Duration[fun.length] benchmarkValues(T)(size_t n, out T[fun.length] values)
Duration[fun.length] benchmarkValues(T)(size_t numberSimulations, out T[fun.length] values)
{
import std.datetime.stopwatch: StopWatch, AutoStart;
import std.datetime.stopwatch: StopWatch, AutoStart;
Duration[fun.length] result;
auto sw = StopWatch(AutoStart.yes);

foreach (i, unused; fun) {
values[i] = 0;
sw.reset();
foreach (size_t j; 1 .. n) {
foreach (size_t j; 1 .. numberSimulations) {
values[i] += fun[i]();
}
result[i] = sw.peek();
values[i] /= n;
values[i] /= numberSimulations;
}

return result;
}
}

// handles 1-dimension case
package(mir)
template benchmarkRandom(fun...)
{
Duration[fun.length] benchmarkRandom(T)(size_t n, size_t m, out T[fun.length] values)
Duration[fun.length] benchmarkRandom(T)(size_t numberSimulations, size_t m, out T[fun.length] values)
if (isMutable!T)
{
import mir.ndslice.allocation: stdcFreeSlice, stdcUninitSlice;
import mir.random.engine: Random, threadLocalPtr;
import mir.random.variable: NormalVariable;
import std.datetime.stopwatch: StopWatch, AutoStart;
import std.datetime.stopwatch: StopWatch, AutoStart;

Random* gen = threadLocalPtr!Random;
auto rv = NormalVariable!T(0, 1);
Expand All @@ -47,31 +48,32 @@ template benchmarkRandom(fun...)
foreach (i, unused; fun) {
values[i] = 0;
sw.reset();
foreach (size_t j; 1 .. n) {
foreach (size_t j; 1 .. numberSimulations) {
sw.stop();
foreach (ref e; r)
e = rv(gen);
sw.start();
values[i] += fun[i](r);
}
result[i] = sw.peek();
values[i] /= n;
values[i] /= numberSimulations;
}
r.stdcFreeSlice;
r.stdcFreeSlice;
return result;
}
}

// handles 2-dimension case
package(mir)
template benchmarkRandom2(fun...)
{
Duration[fun.length] benchmarkRandom2(T)(size_t n, size_t m, out T[fun.length] values)
Duration[fun.length] benchmarkRandom2(T)(size_t numberSimulations, size_t m, out T[fun.length] values)
if (isMutable!T)
{
import mir.ndslice.allocation: stdcFreeSlice, stdcUninitSlice;
import mir.random.engine: Random, threadLocalPtr;
import mir.random.variable: NormalVariable;
import std.datetime.stopwatch: StopWatch, AutoStart;
import std.datetime.stopwatch: StopWatch, AutoStart;

Random* gen = threadLocalPtr!Random;
auto rv = NormalVariable!T(0, 1);
Expand All @@ -84,7 +86,7 @@ template benchmarkRandom2(fun...)
foreach (i, unused; fun) {
values[i] = 0;
sw.reset();
foreach (size_t j; 1 .. n) {
foreach (size_t j; 1 .. numberSimulations) {
sw.stop();
foreach (size_t k; 0 .. m) {
r1[k] = rv(gen);
Expand All @@ -94,10 +96,100 @@ template benchmarkRandom2(fun...)
values[i] += fun[i](r1, r2);
}
result[i] = sw.peek();
values[i] /= n;
values[i] /= numberSimulations;
}
r1.stdcFreeSlice;
r2.stdcFreeSlice;
r1.stdcFreeSlice;
r2.stdcFreeSlice;
return result;
}
}

// handles N-dimension case with additional inputs for mahalanobis
package(mir)
template benchmarkRandomN_mahal(bool[] usePrecision, bool[] useCholesky, fun...)
{
import mir.ndslice;
import mir.ndslice.slice: Slice;
import mir.rc.array: RCI;

Duration[fun.length] benchmarkRandomN_mahal(T)(size_t numberSimulations, size_t m, size_t n, Slice!(RCI!T, 1) mu, Slice!(RCI!T, 2) sigma, ref Slice!(RCI!T, 2) values)
if (isMutable!T)
{
import mir.math.internal.lubeck2: mlinverse;
import mir.math.internal.linearAlgebra: choleskyDecompose;
import mir.ndslice.allocation: stdcFreeSlice, stdcUninitSlice, mininitRcslice;
import mir.random.engine: Random, threadLocalPtr;
import mir.random.ndvariable: multivariateNormalVariable;
import mir.rc.array: RCArray;
import std.datetime.stopwatch: StopWatch, AutoStart;

auto muScope = mu.lightScope;
auto sigmaScope = sigma.lightScope;
auto invSigma = mlinverse(sigma);
auto invSigmaScope = invSigma.lightScope;
auto cholesky = choleskyDecompose(sigma);
auto choleskyScope = cholesky.lightScope;
auto invCholesky = choleskyDecompose(invSigma);
auto invCholeskyScope = invCholesky.lightScope;

Random* gen = threadLocalPtr!Random;
auto sigmaCopy = stdcUninitSlice!T(n, n);
foreach (size_t i; 0 .. n) {
foreach (size_t j; i .. n) {
sigmaCopy[i, j] = sigma[i, j];
if (i != j) {
sigmaCopy[j, i] = sigmaCopy[i, j];
}
}
}
auto rv = multivariateNormalVariable!T(muScope, sigmaCopy, false);

Duration[fun.length] result;
auto r = stdcUninitSlice!T(m, n);
auto temp1 = mininitRcslice!T(m);
auto temp2 = RCArray!T(n);
auto sw = StopWatch(AutoStart.yes);
// Make sure initially zeros
foreach (size_t i; 0 .. fun.length) {
foreach (size_t k; 0 .. m) {
values[i, k] = 0.0;
}
}
foreach (i, unused; fun) {
sw.reset();
foreach (size_t j; 1 .. numberSimulations) {
sw.stop();
foreach (size_t k; 0 .. m) {
rv(gen, temp2[]);
foreach (size_t l; 0 .. n) {
r[k, l] -= temp2[l] - muScope[l];
}
}
sw.start();
static if (usePrecision[i] == true) {
static if (useCholesky[i] == false) {
temp1 = fun[i](r, invSigmaScope);
} else {
temp1 = fun[i](r, invCholeskyScope);
}
} else {
static if (useCholesky[i] == false) {
temp1 = fun[i](r, sigmaScope);
} else {
temp1 = fun[i](r, choleskyScope);
}
}
foreach (size_t k; 0 .. m) {
values[i, k] += temp1[k];
}
}
result[i] = sw.peek();
foreach (size_t k; 0 .. m) {
values[i, k] /= numberSimulations;
}
}
r.stdcFreeSlice;
sigmaCopy.stdcFreeSlice;
return result;
}
}
Loading