Skip to content
Open
Show file tree
Hide file tree
Changes from 11 commits
Commits
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
3 changes: 2 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
.dub
dub.selections.json
__*
*.exe
*.exe
dstats-test-library
20 changes: 10 additions & 10 deletions source/dstats/alloc.d
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ public:
void popFront()
in {
assert(!empty);
} body {
} do {
this._length--;
if(next is null) {
do {
Expand Down Expand Up @@ -149,14 +149,14 @@ public:
@property ref Unqual!(K) front()
in {
assert(!empty);
} body {
} do {
return *frontElem;
}
} else {
@property Unqual!(K) front()
in {
assert(!empty);
} body {
} do {
return *frontElem;
}
}
Expand Down Expand Up @@ -911,7 +911,7 @@ struct AVLNodeBitwise(T) {
void left(typeof(this)* newLeft) nothrow @property
in {
assert((cast(size_t) newLeft & mask) == 0);
} body {
} do {
_left &= mask;
_left |= cast(size_t) newLeft;
assert(left is newLeft);
Expand All @@ -928,7 +928,7 @@ struct AVLNodeBitwise(T) {
void right(typeof(this)* newRight) nothrow @property
in {
assert((cast(size_t) newRight & mask) == 0);
} body {
} do {
_right &= mask;
_right |= cast(size_t) newRight;
assert(right is newRight);
Expand Down Expand Up @@ -1054,7 +1054,7 @@ private:
assert(node.left !is null);
assert( abs(node.balance) <= 2);

} body {
} do {
Node* newHead = node.left;
node.left = newHead.right;
newHead.right = node;
Expand All @@ -1070,7 +1070,7 @@ private:
in {
assert(node.right !is null);
assert( abs(node.balance) <= 2);
} body {
} do {
Node* newHead = node.right;
node.right = newHead.left;
newHead.left = node;
Expand All @@ -1087,7 +1087,7 @@ private:
assert(node is null || abs(node.balance) <= 2);
} out(ret) {
assert( abs(ret.balance) < 2);
} body {
} do {
if(node is null) {
return null;
}
Expand Down Expand Up @@ -1144,7 +1144,7 @@ private:
in {
assert(freeList);
assert(*freeList);
} body {
} do {
auto ret = *freeList;
*freeList = ret.right;
return ret;
Expand All @@ -1153,7 +1153,7 @@ private:
Node* newNode(T payload)
in {
assert(freeList, "Uninitialized StackTree!(" ~ T.stringof ~ ")");
} body {
} do {
Node* ret;
if(*freeList !is null) {
ret = popFreeList();
Expand Down
22 changes: 15 additions & 7 deletions source/dstats/base.d
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,14 @@ version(unittest)
{
static assert(nDims!(uint[]) == 1);
static assert(nDims!(float[][]) == 2);

static if (__VERSION__ < 2096)
alias approxEqual = std.math.approxEqual;
else
bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff = 1e-2, V maxAbsDiff = 1e-5)
{
return std.math.isClose(lhs, rhs, maxRelDiff, maxAbsDiff); // mimic old sloppy approxEqual for now
}
}

import std.string : strip;
Expand Down Expand Up @@ -530,7 +538,7 @@ private void averageTies(T, U)(T sortedInput, U[] ranks, size_t[] perms)
in {
assert(sortedInput.length == ranks.length);
assert(ranks.length == perms.length);
} body {
} do {
size_t tieCount = 1;
foreach(i; 1..ranks.length) {
if(sortedInput[i] == sortedInput[i - 1]) {
Expand Down Expand Up @@ -811,7 +819,7 @@ unittest {
// Values worked out by hand on paper. If you don't believe me, work
// them out yourself.
assert(auroc([4,5,6], [1,2,3]) == 1);
assert(approxEqual(auroc([8,6,7,5,3,0,9], [3,6,2,4,3,6]), 0.6904762));
assert(approxEqual(auroc([8,6,7,5,3,0,9], [3,6,2,4,3,6]), 0.690476190));
assert(approxEqual(auroc([2,7,1,8,2,8,1,8], [3,1,4,1,5,9,2,6]), 0.546875));
}

Expand Down Expand Up @@ -856,7 +864,7 @@ unittest {
double logNcomb(ulong n, ulong k)
in {
assert(k <= n);
} body {
} do {
if(n < k) return -double.infinity;
//Extra parentheses increase numerical accuracy.
return logFactorial(n) - (logFactorial(k) + logFactorial(n - k));
Expand Down Expand Up @@ -1220,7 +1228,7 @@ public:
this(uint n, uint r)
in {
assert(n >= r);
} body {
} do {
if(r > 0) {
pos = (seq(0U, r)).ptr;
pos[r - 1]--;
Expand Down Expand Up @@ -1487,7 +1495,7 @@ unittest {
rng.popFront;
assert(approxEqual(rng.front, 8.67));
rng.popFront;
assert(approxEqual(rng.front, 362435));
assert(approxEqual(rng.front, 362436));
assert(!rng.empty);
rng.popFront;
assert(rng.empty);
Expand Down Expand Up @@ -1624,7 +1632,7 @@ version(scid) {
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]));
assert(approxEqual(toMat[2], [0.70833333333, -0.25000000000, 0.04166666667]));
}

void solve(DoubleMatrix mat, double[] vec) {
Expand Down Expand Up @@ -1690,7 +1698,7 @@ version(scid) {
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]));
assert(approxEqual(vec2, [-3.875, -0.75, 4.45833333333]));
}

// Cholesky decomposition functions adapted from Don Clugston's MathExtra
Expand Down
107 changes: 57 additions & 50 deletions source/dstats/cor.d
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,14 @@ version(unittest) {
{
gen.seed(unpredictableSeed);
}

static if (__VERSION__ < 2096)
alias approxEqual = std.math.approxEqual;
else
bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff = 1e-2, V maxAbsDiff = 1e-5)
{
return std.math.isClose(lhs, rhs, maxRelDiff, maxAbsDiff); // mimic old sloppy approxEqual for now
}
}

/**Convenience function for calculating Pearson correlation.
Expand Down Expand Up @@ -136,7 +144,7 @@ if(doubleInput!(T) && doubleInput!(U)) {
unittest {
assert(approxEqual(pearsonCor([1,2,3,4,5][], [1,2,3,4,5][]).cor, 1));
assert(approxEqual(pearsonCor([1,2,3,4,5][], [10.0, 8.0, 6.0, 4.0, 2.0][]).cor, -1));
assert(approxEqual(pearsonCor([2, 4, 1, 6, 19][], [4, 5, 1, 3, 2][]).cor, -.2382314));
assert(approxEqual(pearsonCor([2, 4, 1, 6, 19][], [4, 5, 1, 3, 2][]).cor, -0.23823144500));

// Make sure everything works with lowest common denominator range type.
static struct Count {
Expand Down Expand Up @@ -169,7 +177,7 @@ unittest {
}

assert(approxEqual(pearsonCor([1,2,3,4,5,6,7,8,9,10][],
[8,6,7,5,3,0,9,3,6,2][]).cor, -0.4190758));
[8,6,7,5,3,0,9,3,6,2][]).cor, -0.41907583841));

foreach(iter; 0..1000) {
// Make sure results for the ILP-optimized and non-optimized versions
Expand Down Expand Up @@ -358,18 +366,18 @@ is(typeof(input2.front < input2.front) == bool)) {

unittest {
//Test against a few known values.
assert(approxEqual(spearmanCor([1,2,3,4,5,6].dup, [3,1,2,5,4,6].dup), 0.77143));
assert(approxEqual(spearmanCor([3,1,2,5,4,6].dup, [1,2,3,4,5,6].dup ), 0.77143));
assert(approxEqual(spearmanCor([1,2,3,4,5,6].dup, [3,1,2,5,4,6].dup), 0.77142857143));
assert(approxEqual(spearmanCor([3,1,2,5,4,6].dup, [1,2,3,4,5,6].dup ), 0.77142857143));
assert(approxEqual(spearmanCor([3,6,7,35,75].dup, [1,63,53,67,3].dup), 0.3));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [3,6,7,35,75].dup), 0.3));
assert(approxEqual(spearmanCor([1.5,6.3,7.8,4.2,1.5].dup, [1,63,53,67,3].dup), .56429));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [1.5,6.3,7.8,4.2,1.5].dup), .56429));
assert(approxEqual(spearmanCor([1.5,6.3,7.8,7.8,1.5].dup, [1,63,53,67,3].dup), .79057));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [1.5,6.3,7.8,7.8,1.5].dup), .79057));
assert(approxEqual(spearmanCor([1.5,6.3,7.8,6.3,1.5].dup, [1,63,53,67,3].dup), .63246));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [1.5,6.3,7.8,6.3,1.5].dup), .63246));
assert(approxEqual(spearmanCor([3,4,1,5,2,1,6,4].dup, [1,3,2,6,4,2,6,7].dup), .6829268));
assert(approxEqual(spearmanCor([1,3,2,6,4,2,6,7].dup, [3,4,1,5,2,1,6,4].dup), .6829268));
assert(approxEqual(spearmanCor([1.5,6.3,7.8,4.2,1.5].dup, [1,63,53,67,3].dup), 0.56428809365));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [1.5,6.3,7.8,4.2,1.5].dup), 0.56428809365));
assert(approxEqual(spearmanCor([1.5,6.3,7.8,7.8,1.5].dup, [1,63,53,67,3].dup), 0.79056941504));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [1.5,6.3,7.8,7.8,1.5].dup), 0.79056941504));
assert(approxEqual(spearmanCor([1.5,6.3,7.8,6.3,1.5].dup, [1,63,53,67,3].dup), 0.63245553203));
assert(approxEqual(spearmanCor([1,63,53,67,3].dup, [1.5,6.3,7.8,6.3,1.5].dup), 0.63245553203));
assert(approxEqual(spearmanCor([3,4,1,5,2,1,6,4].dup, [1,3,2,6,4,2,6,7].dup), 0.68292682927));
assert(approxEqual(spearmanCor([1,3,2,6,4,2,6,7].dup, [3,4,1,5,2,1,6,4].dup), 0.68292682927));
uint[] one = new uint[1000], two = new uint[1000];
foreach(i; 0..100) { //Further sanity checks for things like commutativity.
size_t lowerBound = uniform(0, one.length);
Expand Down Expand Up @@ -589,7 +597,7 @@ package KendallLowLevel kendallCorDestructiveLowLevelImpl
(R1, R2)(R1 input1, R2 input2, bool needTies)
in {
assert(input1.length == input2.length);
} body {
} do {
static ulong getMs(V)(V data) { //Assumes data is sorted.
ulong Ms = 0, tieCount = 0;
foreach(i; 1..data.length) {
Expand Down Expand Up @@ -707,7 +715,7 @@ in {
// implementation exists in this module for large N, but when N gets this
// large it's not even correct due to overflow errors.
assert(input1.length < 1 << 15);
} body {
} do {
int m1 = 0, m2 = 0;
int s = 0;

Expand Down Expand Up @@ -759,9 +767,9 @@ in {

unittest {
//Test against known values.
assert(approxEqual(kendallCor([1,2,3,4,5].dup, [3,1,7,4,3].dup), 0.1054093));
assert(approxEqual(kendallCor([1,2,3,4,5].dup, [3,1,7,4,3].dup), 0.10540925534));
assert(approxEqual(kendallCor([3,6,7,35,75].dup,[1,63,53,67,3].dup), 0.2));
assert(approxEqual(kendallCor([1.5,6.3,7.8,4.2,1.5].dup, [1,63,53,67,3].dup), .3162287));
assert(approxEqual(kendallCor([1.5,6.3,7.8,4.2,1.5].dup, [1,63,53,67,3].dup), 0.31622776602));

static void doKendallTest(T)() {
T[] one = new T[1000], two = new T[1000];
Expand Down Expand Up @@ -929,11 +937,11 @@ unittest {
uint[] consumerFear = [1, 2, 3, 4, 5, 6, 7];
double partialCor =
partial!pearsonCor(stock1Price, stock2Price, [economicHealth, consumerFear][]);
assert(approxEqual(partialCor, -0.857818));
assert(approxEqual(partialCor, -0.85781752401));

double spearmanPartial =
partial!spearmanCor(stock1Price, stock2Price, economicHealth, consumerFear);
assert(approxEqual(spearmanPartial, -0.7252));
assert(approxEqual(spearmanPartial, -0.72521216681));
}

private __gshared TaskPool emptyPool;
Expand Down Expand Up @@ -1225,8 +1233,8 @@ private void pearsonSpearmanCov(bool makeNewMatrix, RoR, Matrix)
break;
case CorCovType.spearman:
foreach(row; pool.parallel(normalized)) {
auto alloc = newRegionAllocator();
auto buf = alloc.uninitializedArray!(double[])(row.length);
auto alloc2 = newRegionAllocator();
auto buf = alloc2.uninitializedArray!(double[])(row.length);
rank(row, buf);

// Need to find mean, stdev separately for every row b/c
Expand Down Expand Up @@ -1259,7 +1267,7 @@ private void dotMatrix(Matrix)(
}

assert(ret.rows == rows.length);
} body {
} do {
// HACK: Before the multithreaded portion of this algorithm
// starts, make sure that there's no need to unshare ret if it's
// using ref-counted COW semantics.
Expand Down Expand Up @@ -1298,39 +1306,38 @@ unittest {

// Values from R.

alias approxEqual ae; // Save typing.
assert(ae(pearsonRoR[0][0], 1));
assert(ae(pearsonRoR[1][1], 1));
assert(ae(pearsonRoR[2][2], 1));
assert(ae(pearsonRoR[1][0], 0.3077935));
assert(ae(pearsonRoR[2][0], -0.9393364));
assert(ae(pearsonRoR[2][1], -0.6103679));

assert(ae(spearmanRoR[0][0], 1));
assert(ae(spearmanRoR[1][1], 1));
assert(ae(spearmanRoR[2][2], 1));
assert(ae(spearmanRoR[1][0], 0.3162278));
assert(ae(spearmanRoR[2][0], -0.9486833));
assert(ae(spearmanRoR[2][1], -0.5));

assert(ae(kendallRoR[0][0], 1));
assert(ae(kendallRoR[1][1], 1));
assert(ae(kendallRoR[2][2], 1));
assert(ae(kendallRoR[1][0], 0.1825742));
assert(ae(kendallRoR[2][0], -0.9128709));
assert(ae(kendallRoR[2][1], -0.4));

assert(ae(covRoR[0][0], 1.66666));
assert(ae(covRoR[1][1], 14.25));
assert(ae(covRoR[2][2], 4.25));
assert(ae(covRoR[1][0], 1.5));
assert(ae(covRoR[2][0], -2.5));
assert(ae(covRoR[2][1], -4.75));
assert(approxEqual(pearsonRoR[0][0], 1));
assert(approxEqual(pearsonRoR[1][1], 1));
assert(approxEqual(pearsonRoR[2][2], 1));
assert(approxEqual(pearsonRoR[1][0], 0.30779350563));
assert(approxEqual(pearsonRoR[2][0], -0.93933643663));
assert(approxEqual(pearsonRoR[2][1], -0.61036793789));

assert(approxEqual(spearmanRoR[0][0], 1));
assert(approxEqual(spearmanRoR[1][1], 1));
assert(approxEqual(spearmanRoR[2][2], 1));
assert(approxEqual(spearmanRoR[1][0], 0.31622776602));
assert(approxEqual(spearmanRoR[2][0], -0.94868329805));
assert(approxEqual(spearmanRoR[2][1], -0.5));

assert(approxEqual(kendallRoR[0][0], 1));
assert(approxEqual(kendallRoR[1][1], 1));
assert(approxEqual(kendallRoR[2][2], 1));
assert(approxEqual(kendallRoR[1][0], 0.18257418584));
assert(approxEqual(kendallRoR[2][0], -0.91287092918));
assert(approxEqual(kendallRoR[2][1], -0.4));

assert(approxEqual(covRoR[0][0], 1.66666666667));
assert(approxEqual(covRoR[1][1], 14.25));
assert(approxEqual(covRoR[2][2], 4.25));
assert(approxEqual(covRoR[1][0], 1.5));
assert(approxEqual(covRoR[2][0], -2.5));
assert(approxEqual(covRoR[2][1], -4.75));

version(scid) {
static bool test(double[][] a, SymmetricMatrix!double b) {
foreach(i; 0..3) foreach(j; 0..i + 1) {
if(!ae(a[i][j], b[i, j])) return false;
if(!approxEqual(a[i][j], b[i, j])) return false;
}

return true;
Expand Down
Loading