Changeset 28

Timestamp:

04/23/09 23:30:02 (3 years ago)

Author:

dsimcha

Message:

Make dstats work with ranges.

Files:

• trunk/alloc.d (modified) (11 diffs)
• trunk/base.d (modified) (34 diffs)
• trunk/cor.d (modified) (15 diffs)
• trunk/distrib.d (modified) (16 diffs)
• trunk/gamma.d (modified) (1 diff)
• trunk/infotheory.d (modified) (9 diffs)
• trunk/random.d (modified) (19 diffs)
• trunk/sort.d (modified) (24 diffs)
• trunk/stats.d (modified) (1 diff)
• trunk/summary.d (modified) (25 diffs)
• trunk/tests.d (modified) (42 diffs)

trunk/alloc.d

@@ -1 +1 @@
-/**Allocation functions.
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -28 +61 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-
-module dstats.alloc;
-
-
+
+
-
-version(unittest) {
-
+, std.conv, dstats.sort
-    void main() {}
-}
@@ -55 +90 @@
-        const bool isReferenceType = false;
-    } else static if (Types.length == 1) {
-Mutable
+Unqual
-                           int, uint, long, ulong, float, double, real, ifloat,
-                           idouble, ireal, cfloat, cdouble, creal, char, dchar,
@@ -131 +166 @@
-}
-
+private template Appends(T, U) {
+    enum bool Appends = AppendsImpl!(T, U).ret;
+}
+
+private template AppendsImpl(T, U) {
+    T[] a;
+    U b;
+    enum bool ret = is(typeof(a ~= b));
+}
+
-///Appends to an array, deleting the old array if it has to be realloced.
-void appendDelOld(T, U)(ref T[] to, U from)
- (is(Mutable!(T) : Mutable!(U)) || is(Mutable!(T[0]) : Mutable!(U)
+(Appends!(T, U
-    auto oldPtr = to.ptr;
-    to ~= from;
@@ -203 +248 @@
-        Stack!(Block) inUse;
-        Stack!(void*) freelist;
-
-        ~this() {  // Blocks are pretty large.  Prevent false ptrs.
-            ntFree(lastAlloc.ptr);
-            while(nblocks > 1) {
-                ntFree((inUse.pop()).space);
-                nblocks--;
-            }
-            ntFree(space);
-            while(nfree > 0) {
-                ntFree(freelist.pop);
-                nfree--;
-            }
-        }
-    }
-
@@ -230 +262 @@
-
-    static void die() nothrow {
-
+.c.stdio.std
-        exit(1);
-    }
@@ -457 +489 @@
-        totalLen += array.length;
-    }
-Mutable
+Unqual
-
-    size_t offset = 0;
@@ -467 +499 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
@@ -495 +617 @@
- * are allocated, due to caching of data stored in thread-local
- * storage.*/
-nvariant
+mmutable
-    "TempAlloc.frameInit; scope(exit) TempAlloc.frameFree;";
-
@@ -509 +631 @@
-         TempAlloc(TempAlloc.alignBytes);
-     }
-
+, to!string(TempAlloc.getState.nblocks)
-     assert(TempAlloc.getState.nfree == 0);
-     foreach(i; 0..nIter) {
@@ -580 +702 @@
-        TempAlloc.free;
-    }
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

trunk/base.d

@@ -1 +1 @@
-/**Relatively low-level primitives on which to build higher-level math/stat
- * functionality.  Some are used internally, some are just things that may be
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -31 +64 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-
-module dstats.base;
-
-public import std.math, std.traits, dstats.gamma, dstats.alloc;
-
-
+ std.typecons,
+, std.range, std.bitmanip
-
-invariant real[] staticLogFacTable;
@@ -43 +77 @@
-enum : size_t {
-    staticFacTableLen = 10_000,
+}
+
+// Tests whether T is an input range whose elements can be implicitly
+// converted to reals.
+template realInput(T) {
+    enum realInput = isInputRange!(T) && is(ElementType!(T) : real);
+}
+
+// See Bugzilla 2873.  This can be removed once that's fixed.
+template hasLength(R) {
+    enum bool hasLength = is(typeof(R.init.length) : ulong) ||
+                      is(typeof(R.init.length()) : ulong);
-}
-
@@ -78 +124 @@
-}
-
+/**Converts any range to an array on the GC heap by the most efficient means
+ * available.  If it is already an array, duplicates the range.*/
+Unqual!(ElementType!(T))[] toArray(T)(T range) if(isInputRange!(T)) {
+    static if(isArray!(T)) {
+        return range.dup;
+    } else static if(hasLength!(T)) {
+        auto ret = newVoid!(Unqual!(ElementType!(T)))(range.length);
+        static if(is(typeof(ret[] = range[]))) {
+            ret[] = range[];
+        } else {
+            size_t pos = 0;
+            foreach(elem; range) {
+                ret[pos++] = elem;
+            }
+        }
+        return ret;
+    } else {
+        Unqual!(ElementType!(T))[] ret;
+        mate(range, appender(&ret));
+        return ret;
+    }
+}
+
+/**Writes the contents of an input range to an output range.
+ *
+ * Returns:  The output range.*/
+O mate(I, O)(I input, O output)
+if(isInputRange!(I) && isOutputRange!(O, ElementType!(I))) {
+    foreach(elem; input) {
+        output.put(elem);
+    }
+    return output;
+}
+
-/**Bins data into nbin equal width bins, indexed from
- * 0 to nbin - 1, with 0 being the smallest bin, etc.
- * The values returned are the counts for each bin.  Returns results on the GC
-
-
+
+
+
+
+
-in {
-    assert(data.length > 0);
-    assert(nbin > 0);
-} body {
-
+
+
-    foreach(elem; data[1..$]) {
-        if(elem > max)
@@ -94 +177 @@
-            min = elem;
-    }
-T
+E
-
-    Ret[] bins;
@@ -130 +213 @@
-/**Bins data into nbin equal width bins, indexed from
- * 0 to nbin - 1, with 0 being the smallest bin, etc.
-
-
-
+
+
+
+
+
+
+
+
+
+
+
-in {
-data.length > 0
+nbin <= Ret.max + 1
-    assert(nbin > 0);
-} body {
-
-
+
+
+
+
+
-        if(elem > max)
-            max = elem;
@@ -144 +238 @@
-            min = elem;
-    }
-T
+E
-
-    Ret[] bins;
@@ -171 +265 @@
-    double[] data = [0.0, .01, .03, .05, .11, .3, .5, .7, .89, 1];
-    auto res = bin(data, 10);
-short
+byte
-    res = bin(data, 10, Alloc.STACK);
-short
+byte
-    TempAlloc.free;
-    writeln("Passed bin unittest.");
@@ -180 +274 @@
-/**Bins data into nbin equal frequency bins, indexed from
- * 0 to nbin - 1, with 0 being the smallest bin, etc.
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-in {
-    assert(data.length > 0);
-    assert(nbin > 0);
-    assert(nbin <= data.length);
+    assert(nbin <= Ret.max + 1);
-} body {
-    Ret[] result;
@@ -198 +301 @@
-    foreach(i, ref e; perm)
-        e = i;
-data.tempdup
+tempdup(data)
-    qsort(dd, perm);
-    TempAlloc.free;
@@ -217 +320 @@
-    double[] data = [5U, 1, 3, 8, 30, 10, 7];
-    auto res = frqBin(data, 3);
-short
+byte
-    data = [3, 1, 4, 1, 5, 9, 2, 6, 5];
-    res = frqBin(data, 4, Alloc.STACK);
-short
+byte
-    data = [3U, 1, 4, 1, 5, 9, 2, 6, 5, 3, 4, 8, 9, 7, 9, 2];
-    res = frqBin(data, 4);
-short
+byte
-    TempAlloc.free;
-    writeln("Passed frqBin unittest.");
@@ -229 +332 @@
-
-/**Generates a sequence from [start..end] by increment.  Includes start,
-
+  Does so eagerly as an array.
- *
- * Examples:
@@ -254 +357 @@
-    auto s = seq(0, 5);
-    assert(s == [0, 1, 2, 3, 4]);
-stderr, 
+
-}
-
-/**Given an input array, outputs an array containing the rank from
- * [1, input.length] corresponding to each element.  Ties are dealt with by
-array
+range
- * it.  Return type is float[] by default, but if you are sure you have no ties,
- * ints can be used for efficiency, and if you need more precision when
- * averaging ties, you can use double or real.
+ *
+ * Works with any input range.
- *
- * Examples:
@@ -271 +376 @@
- * ---*/
-Ret[] rank(Ret = float, T)(const T[] input) {
-input.tempdup
+tempdup(input)
-    scope(exit) TempAlloc.free;
-    return rankSort!(Ret)(iDup);
-}
-
-array
+range
- * duping it and working on a copy.  The array returned will still be
- * identical to that returned by rank(), i.e. the rank of each element will
- * correspond to the ranks of the elements in the input array before sorting.
+ *
+ * Works with any random access range with a length property.
- *
- * Examples:
@@ -287 +394 @@
- * assert(test == [1U, 2, 3, 4, 5]);
- * ---*/
-
+
+
-    Ret[] ranks = newVoid!(Ret)(input.length);
-    rankSort!(Ret, T)(input, ranks);
@@ -294 +402 @@
-
-// Speed hack used internally.
-[]
+
-    size_t[] perms = newStack!(size_t)(input.length);
-    scope(exit) TempAlloc.free;
@@ -316 +424 @@
-    assert(rankSort(test) == [3.5f, 5f, 3.5f, 1f, 2f]);
-    assert(test == [1U,2,3,3,5]);
-stderr, 
+
-}
-
-// Used internally by rank().
-[]
+
-in {
-    assert(sortedInput.length == ranks.length);
@@ -350 +458 @@
-}
-
-
+  Works w/ any input range.
- *
- * Examples:
@@ -360 +468 @@
- * assert(frq[4] == 1);
- * ---*/
-
-
+
+
+
-    foreach(i; input) {
-        output[i]++;
@@ -377 +486 @@
-}
-
-unittest {
-    uint[int] temp=frequency([1,2,2,1,2]);
-    assert(temp[1]==2);
-    assert(temp[2]==3);
-    writefln("Passed frequency unittest.");
-}
-
-///
-int
+T
-    if (num > 0) return 1;
-    if (num < 0) return -1;
@@ -403 +505 @@
- * function is used, because caching would take up too much memory, and if
- * done lazily, would cause threading issues.*/
- 
+
-    //Input is uint, can't be less than 0, no need to check.
-    if(n < staticFacTableLen) {
@@ -440 +542 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+//
+//
+//
+//
+//
+//
+//
+//
+
+//
+//
+//
+//
+//
+//
+//
+//
-
-/**A struct that generates all possible permutations of a sequence,
@@ -463 +565 @@
- * are references to the internal permutation state.  This is dangerous, but
- * necessary for performance.  Therefore, you
-
+
+
- *
- * Examples:
@@ -488 +591 @@
-
-public:
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-     * so that the original sequence is not modified.*/
-
-
-
-
+
+
+
+
+
+
-    }
-
@@ -545 +638 @@
-}
-
+private template PermRet(T...) {
+    static if(isForwardRange!(T[0])) {
+        alias Perm!(ElementType!(T[0])) PermRet;
+    } else static if(T.length == 1) {
+        alias Perm!uint PermRet;
+    } else alias Perm!(T[0]) PermRet;
+}
+
+/**Create a Perm struct from a range or of a set of bounds.
+ *
+ * Note:  PermRet is just a template to figure out what this should return.
+ * I would use auto if not for bug 2251.
+ *
+ * Examples:
+ * ---
+ * auto p = perm([1,2,3]);
+ * auto p = perm(5);  // Permutations of integers on range [0, 5].
+ * auto p = perm(-1, 2); // Permutations of integers on range [-1, 2].
+ * ---
+ */
+PermRet!(T) perm(T...)(T stuff) {
+    alias typeof(return) rt;
+    static if(isForwardRange!(T[0])) {
+        return rt(stuff);
+    } else static if(T.length == 1) {
+        static assert(isIntegral!(T[0]));
+        return rt(seq(0U, cast(uint) stuff[0]));
+    } else {
+        return rt(seq(stuff));
+    }
+}
+
+
-unittest {
-    double[][] res;
-
-
-
+
+
-        res ~= p.dup;
-    }
-
-
-
-
-
-
+
+
+
+
+
+
+
-    assert(res.length == 6);
-    uint[][] res2;
-
-
+
-    foreach(p; perm2) {
-        res2 ~= p.dup;
-    }
-
-
-
-
-
-
+
+
+
+
+
+
+
-    assert(res2.length == 6);
-
@@ -576 +702 @@
-    // them, and they contain what they're supposed to contain, the result is
-    // correct.
-PermU
+perm
-    bool[uint[]] table;
-    foreach(p; perm3) {
@@ -585 +711 @@
-        assert(elem.dup.insertionSort == [0U, 1, 2, 3, 4, 5]);
-    }
-PermU
+perm
-    bool[uint[]] table2;
-    foreach(p; perm4) {
@@ -594 +720 @@
-        assert(elem.dup.insertionSort == [0U, 1, 2, 3, 4]);
-    }
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
@@ -602 +756 @@
- * are const references to the internal state of the Comb object.  This is
- * dangerous but necessary for performance.  If you want to save them past the
-
+
+
+
- * Examples:
- * ---
@@ -719 +875 @@
-unittest {
-    // Test indexing verison first.
-Comb!(uint)
+comb
-    uint[][] vals;
-    foreach(c; comb1) {
-        vals ~= c.dup;
-    }
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
-    assert(vals.length == 10);
-
-    // Now, test the array version.
-Comb!(uint)
+comb
-    vals = null;
-    foreach(c; comb2) {
-        vals ~= c.dup;
-    }
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
-    assert(vals.length == 10);
-
@@ -777 +935 @@
-}
-
-// This should be encapsulated within StackHash.  Workaround for bug 1629.
-struct SHNode(K, V) {
-    alias SHNode!(K, V) SomeType;
-    SomeType* next;
-    Mutable!(K) key;
-    Mutable!(V) val;
-}
-
-/* A hash table that uses TempAlloc to allocate space.  Useful for building
- * a quick symbol table in a performance-critical function that can't be
- * performing tons of heap allocations.  Intentionally lacking ddoc because
- * the design or even existence of this is still likely to change.*/
-struct StackHash(K, V) {
-private:
-    alias SHNode!(K, V) Node;
-    Node[] roots;
-    TempAlloc.State TAState;
-    TypeInfo keyTI;
-    size_t _length;
-    Node* usedSentinel;
-
-    Node* newNode(K key) {
-        Node* ret = cast(Node*) TempAlloc(Node.sizeof, TAState);
-        ret.key =  key;
-        ret.val =  V.init;
-        ret.next = null;
-        return ret;
-    }
-
-    Node* newNode(K key, V val) {
-        Node* ret = cast(Node*) TempAlloc(Node.sizeof, TAState);
-        ret.key =  key;
-        ret.val = val;
-        ret.next = null;
-        return ret;
-    }
-
-    hash_t getHash(K key) {
-        static if(is(K : long) && K.sizeof <= hash_t.sizeof) {
-            hash_t hash = cast(hash_t) key;
-        } else static if(__traits(compiles, key.toHash())) {
-            hash_t hash = key.toHash();
-        } else {
-            hash_t hash = keyTI.getHash(&key);
-        }
-        hash %= roots.length;
-        return hash;
-    }
-
-
-public:
-    this(size_t nElem) {
-        // Obviously, the caller can never mean zero, because this struct
-        // can't work at all with nElem == 0, so assume it's a mistake and fix
-        // it here.
-        if(nElem == 0)
-            nElem++;
-        TAState = TempAlloc.getState;
-        roots = newStack!(Node)(nElem, TAState);
-        usedSentinel = cast(Node*) roots.ptr;
-        foreach(ref root; roots) {
-            root.key =  K.init;
-            root.val = V.init;
-            root.next = usedSentinel;
-        }
-        keyTI = typeid(K);
-    }
-
-    ref V opIndex(K key) {
-        hash_t hash = getHash(key);
-
-        if(roots[hash].next == usedSentinel) {
-            roots[hash].key =  key;
-            roots[hash].next = null;
-            _length++;
-            return roots[hash].val;
-        } else if(roots[hash].key == key) {
-            return roots[hash].val;
-        } else {  // Collision.  Start chaining.
-            Node** next = &(roots[hash].next);
-            while(*next !is null) {
-                if((**next).key ==  key) {
-                    return (**next).val;
-                }
-                next = &((**next).next);
-            }
-            *next = newNode(key);
-            _length++;
-            return (**next).val;
-        }
-    }
-
-    V opIndexAssign(V val, K key) {
-        hash_t hash = getHash(key);
-
-        if(roots[hash].next == usedSentinel) {
-            roots[hash].key =  key;
-            roots[hash].val = val;
-            roots[hash].next = null;
-            _length++;
-            return val;
-        } else if(roots[hash].key ==  key) {
-            roots[hash].val = val;
-            return val;
-        } else {  // Collision.  Start chaining.
-            Node** next = &(roots[hash].next);
-            while(*next !is null) {
-                if((**next).key == key) {
-                    (**next).val = val;
-                    return val;
-                }
-                next = &((**next).next);
-            }
-            _length++;
-            *next = newNode(key, val);
-            return val;
-        }
-    }
-
-    V[] values() {
-        auto space = newVoid!(V)(_length);
-        return values(space);
-    }
-
-    V[] valStack() {
-        auto space = newStack!(V)(_length);
-        return values(space);
-    }
-
-    V[] values(V[] space) {
-        size_t pos;
-        foreach(r; roots) {
-            if(r.next == usedSentinel)
-                continue;
-            space[pos++] = r.val;
-            Node* next = r.next;
-            while(next !is null) {
-                space[pos++] = next.val;
-                next = next.next;
-            }
-        }
-        return space;
-    }
-
-    Mutable!(K)[] keys() const {
-        auto space = newVoid!(Mutable!(K))(_length);
-        return keys(space);
-    }
-
-    Mutable!(K)[] keyStack() const {
-        auto space = newStack!(Mutable!(K))(_length);
-        return keys(space);
-    }
-
-    Mutable!(K)[] keys(Mutable!(K)[] space) const {
-        size_t pos;
-        foreach(r; roots) {
-            if(r.next == usedSentinel)
-                continue;
-            space[pos++] = r.key;
-            const(Node)* next = r.next;
-            while(next !is null) {
-                space[pos++] = next.key;
-                next = next.next;
-            }
-        }
-        return space;
-    }
-
-    int opApply(int delegate(ref V value) dg) {
-        int res = 0;
-        outer:
-        foreach(r; roots) {
-            if(r.next == usedSentinel)
-                continue;
-            res = dg(r.val);
-            if (res) break;
-            Node* next = r.next;
-            while(next !is null) {
-                res = dg(next.val);
-                if(res) break outer;
-                next = next.next;
-            }
-        }
-        return res;
-   }
-
-   int opApply(int delegate(ref K key, ref V value) dg) {
-        int res = 0;
-        outer:
-        foreach(r; roots) {
-            if(r.next == usedSentinel)
-                continue;
-            res = dg(r.key, r.val);
-            if (res) break;
-            Node* next = r.next;
-            while(next !is null) {
-                res = dg(next.key, next.val);
-                if(res) break outer;
-                next = next.next;
-            }
-        }
-        return res;
-   }
-
-   V* opIn_r(K key) {
-        hash_t hash = getHash(key);
-
-        if(roots[hash].next == usedSentinel) {
-            return null;
-        } else if(roots[hash].key == key) {
-            return &(roots[hash].val);
-        } else {  // Collision.  Start chaining.
-            Node* next = roots[hash].next;
-            while(next !is null) {
-                if(next.key == key) {
-                    return &(next.val);
-                }
-                next = next.next;
-            }
-            return null;
-        }
-   }
-
-   void remove(K key) {
-        hash_t hash = getHash(key);
-
-        Node** next = &(roots[hash].next);
-        if(roots[hash].next == usedSentinel) {
-            return;
-        } else if(roots[hash].key == key) {
-            _length--;
-            if(roots[hash].next is null) {
-                roots[hash].key = K.init;
-                roots[hash].val = V.init;
-                roots[hash].next = usedSentinel;
-                return;
-            } else {
-                roots[hash].key = (**next).key;
-                roots[hash].val = (**next).val;
-                roots[hash].next = (**next).next;
-                return;
-            }
-        } else {  // Collision.  Start chaining.
-            while(*next !is null) {
-                if((**next).key == key) {
-                    _length--;
-                    *next = (**next).next;
-                    break;
-                }
-                next = &((**next).next);
-            }
-            return;
-        }
-   }
-
-   size_t length() {
-       return _length;
-   }
-
-   real efficiency() {
-       uint used = 0;
-       foreach(root; roots) {
-           if(root.next != usedSentinel) {
-               used++;
-           }
-       }
-       return cast(real) used / roots.length;
-   }
-}
-
-unittest {
-    alias StackHash!(string, uint) mySh;
-    mixin(newFrame);
-    auto data = mySh(2);  // Make sure we get some collisions.
-    data["foo"] = 1;
-    data["bar"] = 2;
-    data["baz"] = 3;
-    data["waldo"] = 4;
-    assert(!("foobar" in data));
-    assert(*("foo" in data) == 1);
-    assert(*("bar" in data) == 2);
-    assert(*("baz" in data) == 3);
-    assert(*("waldo" in data) == 4);
-    assert(data["foo"] == 1);
-    assert(data["bar"] == 2);
-    assert(data["baz"] == 3);
-    assert(data["waldo"] == 4);
-    assert(data.keys.sort == ["bar", "baz", "foo", "waldo"]);
-    assert(data.values.sort == [1U, 2, 3, 4]);
-    foreach(k, v; data) {
-        assert(data[k] == v);
-    }
-    foreach(v; data) {
-        assert(v > 0 && v < 5);
-    }
-
-    // Test remove.
-
-    alias StackHash!(uint, uint) mySh2;
-    auto foo = mySh2(7);
-    for(uint i = 0; i < 20; i++) {
-        foo[i] = i;
-    }
-    assert(foo.length == 20);
-    for(uint i = 0; i < 20; i += 2) {
-        foo.remove(i);
-    }
-    for(uint i = 0; i < 20; i++) {
-        if(i & 1) {
-            assert(*(i in foo) == i);
-        } else {
-            assert(!(i in foo));
-        }
-    }
-    auto vals = foo.values;
-    assert(foo.length == 10);
-    assert(vals.qsort == [1U, 3, 5, 7, 9, 11, 13, 15, 17, 19]);
-
-    writeln("Passed StackHash test.");
-}
-
-/* A hash set that uses TempAlloc to allocate space.  Useful for building
- * a quick set in a performance-critical function that can't be
- * performing tons of heap allocations.  Intentionally lacking ddoc because
- * the design or even existence of this is still likely to change.*/
-struct StackSet(K) {
-private:
-    // Choose smallest representation of the data.
-    struct Node1 {
-        Node1* next;
-        K key;
-    }
-
-    struct Node2 {
-        K key;
-        Node2* next;
-    }
-
-    static if(Node1.sizeof < Node2.sizeof) {
-        alias Node1 Node;
-    } else {
-        alias Node2 Node;
-    }
-
-    Node[] roots;
-    TempAlloc.State TAState;
-    TypeInfo keyTI;
-    size_t _length;
-    Node* usedSentinel;
-
-    Node* newNode(K key) {
-        Node* ret = cast(Node*) TempAlloc(Node.sizeof, TAState);
-        ret.key = key;
-        ret.next = null;
-        return ret;
-    }
-
-    hash_t getHash(K key) {
-        static if(is(K : long) && K.sizeof <= hash_t.sizeof) {
-            hash_t hash = cast(hash_t) key;
-        } else static if(__traits(compiles, key.toHash())) {
-            hash_t hash = key.toHash();
-        } else {
-            hash_t hash = keyTI.getHash(&key);
-        }
-        hash %= roots.length;
-        return hash;
-    }
-
-public:
-    this(size_t nElem) {
-        // Obviously, the caller can never mean zero, because this struct
-        // can't work at all with nElem == 0, so assume it's a mistake and fix
-        // it here.
-        if(nElem == 0)
-            nElem++;
-        TAState = TempAlloc.getState;
-        roots = newStack!(Node)(nElem, TAState);
-        usedSentinel = cast(Node*) roots.ptr;
-        foreach(ref root; roots) {
-            root.key = K.init;
-            root.next = usedSentinel;
-        }
-        keyTI = typeid(K);
-    }
-
-    void insert(K key) {
-        hash_t hash = getHash(key);
-
-        if(roots[hash].next == usedSentinel) {
-            roots[hash].key = key;
-            roots[hash].next = null;
-            _length++;
-            return;
-        } else if(roots[hash].key == key) {
-            return;
-        } else {  // Collision.  Start chaining.
-            Node** next = &(roots[hash].next);
-            while(*next !is null) {
-                if((**next).key == key) {
-                    return;
-                }
-                next = &((**next).next);
-            }
-            *next = newNode(key);
-            _length++;
-            return;
-        }
-    }
-
-    K[] elems() {
-        auto space = newVoid!(K)(_length);
-        return elems(space);
-    }
-
-    K[] elemStack() {
-        auto space = newStack!(K)(_length);
-        return elems(space);
-    }
-
-    K[] elems(K[] space) {
-        size_t pos;
-        foreach(r; roots) {
-            if(r.next == usedSentinel)
-                continue;
-            space[pos++] = r.key;
-            Node* next = r.next;
-            while(next !is null) {
-                space[pos++] = next.key;
-                next = next.next;
-            }
-        }
-        return space;
-    }
-
-   int opApply(int delegate(ref K key) dg) {
-        int res = 0;
-        outer:
-        foreach(r; roots) {
-            if(r.next == usedSentinel)
-                continue;
-            res = dg(r.key);
-            if (res) break;
-            Node* next = r.next;
-            while(next !is null) {
-                res = dg(next.key);
-                if(res) break outer;
-                next = next.next;
-            }
-        }
-        return res;
-   }
-
-   bool opIn_r(K key) {
-        hash_t hash = getHash(key);
-
-        if(roots[hash].next == usedSentinel) {
-            return false;
-        } else if(roots[hash].key == key) {
-            return true;
-        } else {  // Collision.  Start chaining.
-            Node* next = roots[hash].next;
-            while(next !is null) {
-                if(next.key == key) {
-                    return true;
-                }
-                next = next.next;
-            }
-            return false;
-        }
-   }
-
-   void remove(K key) {
-        hash_t hash = getHash(key);
-
-        Node** next = &(roots[hash].next);
-        if(roots[hash].next == usedSentinel) {
-            return;
-        } else if(roots[hash].key == key) {
-            _length--;
-            if(roots[hash].next is null) {
-                roots[hash].key = K.init;
-                roots[hash].next = usedSentinel;
-                return;
-            } else {
-                roots[hash].key = (**next).key;
-                roots[hash].next = (**next).next;
-                return;
-            }
-        } else {  // Collision.  Start chaining.
-            while(*next !is null) {
-                if((**next).key == key) {
-                    _length--;
-                    *next = (**next).next;
-                    break;
-                }
-                next = &((**next).next);
-            }
-            return;
-        }
-   }
-
-   size_t length() {
-       return _length;
-   }
-}
-
-unittest {
-    mixin(newFrame);
-    alias StackSet!(uint) mySS;
-    mySS set = mySS(12);
-    foreach(i; 0..20) {
-        set.insert(i);
-    }
-    assert(set.elems.qsort == seq(0U, 20U));
-    assert(set.elemStack.qsort == seq(0U, 20U));
-
-    for(uint i = 0; i < 20; i += 2) {
-        set.remove(i);
-    }
-
-    foreach(i; 0..20) {
-        if(i & 1) {
-            assert(i in set);
-        } else {
-            assert(!(i in set));
-        }
-    }
-    uint[] contents;
-    foreach(elem; set) {
-        contents ~= elem;
-    }
-    assert(contents.qsort == [1U,3,5,7,9,11,13,15,17,19]);
-    writeln("Passed StackSet test.");
-}
-
-/**Computes the intersect of two arrays, i.e. the elements that are in both
- * arrays.  Time and space complexity are O(first.length + second.length).
@@ -1343 +964 @@
-
-    size_t pos;
-
-
-
-
-
+
+
+
+
+
+
+
+
+
-    return result[0..pos];
-}
@@ -1402 +1027 @@
-        }
-    }
-
-    return result[0..resPos];
-}
@@ -1415 +1039 @@
-    foreach(i; 0..1000) {
-        foreach(ref f; first)
-gen, 
+
-        foreach(ref s; second)
-gen, 
+
-        auto hash = qsort!("a > b")(intersect(first, second));
-        auto sort = intersectSorted!("a > b")

trunk/cor.d

@@ -1 +1 @@
-/**Pearson, Spearman and Kendall correlations, covariance.
- *
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -29 +61 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-
-module dstats.cor;
-
-
+, std.range
-
-import dstats.sort, dstats.base, dstats.alloc;
@@ -49 +82 @@
-/**Pearson correlation.  When the term correlation is used unqualified, it is
- * usually referring to this quantity.  Pearson correlation assumes the input
-*/
-
-
-
-} body
+
+
+
+
+if(realInput!(T) && realInput!(U))
-    OnlinePcor corCalc;
-
-
+
+
+
+
-    }
-
@@ -63 +98 @@
-
-unittest {
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    writefln("Passed pcor unittest.");
-}
@@ -78 +134 @@
-public:
-    ///
-addElemen
+pu
-        _k++;
-        real kNeg1 = 1.0L / _k;
@@ -135 +191 @@
-
-///
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-in {
-    assert(input1.length == input2.length);
-} body {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-    if(input1.length < 2)
-        return real.nan;
@@ -171 +233 @@
-    iDup = (cast(T*) TempAlloc.malloc(largerSize * input1.length))
-           [0..input1.length];
-iDup[] = input1[]
+rangeCopy(iDup, input1)
-    qsort(iDup, perms);
-
@@ -189 +251 @@
-    // for the larger of T, U.
-    U[] iDup2 = (cast(U*) iDup.ptr)[0..input2.length];
-iDup2[] = input2[]
+rangeCopy(iDup2, input2)
-
-    qsort(iDup2, perms);
@@ -202 +264 @@
-unittest {
-    //Test against a few known values.
-, [3,1,2,5,4,6]
-, [1,2,3,4,5,6]
-, [1,63,53,67,3]
-, [3,6,7,35,75]
-, [1,63,53,67,3]
-, [1.5,6.3,7.8,4.2,1.5]
-, [1,63,53,67,3]
-, [1.5,6.3,7.8,7.8,1.5]
-, [1,63,53,67,3]
-, [1.5,6.3,7.8,6.3,1.5]
-, [1,3,2,6,4,2,6,7]
-, [3,4,1,5,2,1,6,4]
+.dup, [3,1,2,5,4,6].dup
+.dup, [1,2,3,4,5,6].dup
+.dup, [1,63,53,67,3].dup
+.dup, [3,6,7,35,75].dup
+.dup, [1,63,53,67,3].dup
+.dup, [1.5,6.3,7.8,4.2,1.5].dup
+.dup, [1,63,53,67,3].dup
+.dup, [1.5,6.3,7.8,7.8,1.5].dup
+.dup, [1,63,53,67,3].dup
+.dup, [1.5,6.3,7.8,6.3,1.5].dup
+.dup, [1,3,2,6,4,2,6,7].dup
+.dup, [3,4,1,5,2,1,6,4].dup
-    uint[] one = new uint[1000], two = new uint[1000];
-    foreach(i; 0..100) {  //Further sanity checks for things like commutativity.
-gen, 
-gen, 
+
+
-        if(lowerBound > upperBound) swap(lowerBound, upperBound);
-        foreach(ref o; one) {
-gen, 
+
-        }
-        foreach(ref o; two) {
-gen, 
+
-        }
-        real sOne =
@@ -245 +307 @@
-        assert(approxEqual(sFour, sFive) || (isnan(sFour) && isnan(sFive)));
-    }
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
@@ -255 +341 @@
-  * standard formulas, and therefore unlikely to have weird, subtle bugs.*/
-
-
+private 
-in {
-    assert(input1.length == input2.length);
@@ -294 +380 @@
- * bubble sort distance, or the number of swaps that would be needed in a
- * bubble sort to sort input2 into the same order as input1.  It is
-
-
-
+
+
+
+
+
+
+
+
-    scope(exit) TempAlloc.free;
-input2.tempdup
+tempdup(input2)
-    scope(exit) TempAlloc.free;
-    return kcorDestructive(i1d, i2d);
@@ -304 +395 @@
-
-/**Kendall's Tau O(N log N) destroys input vectors but requires
-
+
+
-real kcorDestructive(T, U)(T[] input1, U[] input2)
-in {
@@ -369 +461 @@
-unittest {
-    //Test against known values.
-, [3,1,7,4,3]
-,[1,63,53,67,3]
-, [1,63,53,67,3]
+.dup, [3,1,7,4,3].dup
+.dup,[1,63,53,67,3].dup
+.dup, [1,63,53,67,3].dup
-    uint[] one = new uint[1000], two = new uint[1000];
-    // Test complex, fast implementation against straightforward,
-    // slow implementation.
-    foreach(i; 0..100) {
-gen, 
-gen, 
+
+
-        if(lowerBound > upperBound) swap(lowerBound, upperBound);
-        foreach(ref o; one) {
-gen, 
+
-        }
-        foreach(ref o; two) {
-gen, 
+
-        }
-        real kOne =
@@ -389 +481 @@
-        real kTwo =
-             kcorOld(one[lowerBound..upperBound], two[lowerBound..upperBound]);
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    writefln("Passed kcor unittest.");
-}

trunk/distrib.d

@@ -1 +1 @@
-/**Probability distribution CDFs, PDFs/PMFs, and a few inverse CDFs.
- *
-
-
+*/
+/
- * Acknowledgements:  Some of this module was borrowed the mathstat module
- * of Don Clugston's MathExtra library.  This was done to create a
@@ -251 +251 @@
-
-unittest {
-   Random gen = Random(unpredictableSeed);
-   foreach(i; 0..1_000) {
-       // Restricted variable ranges are because, in the tails, more than one
@@ -257 +256 @@
-       // Obviously, this is one of those corner cases that nothing can be
-       // done about.
-gen, 
-gen, 
+
+
-       real pVal = poissonCDF(k, lambda);
-       assert(invPoissonCDF(pVal, lambda) == k);
@@ -377 +376 @@
-       // Obviously, this is one of those corner cases that nothing can be
-       // done about.  Using small n's, moderate p's prevents this.
-gen, 
-gen, 
-gen, 
+
+
+
-       real pVal = binomialCDF(k, n, p);
-       assert(invBinomialCDF(pVal, n, p) == k);
@@ -505 +504 @@
-
-real hyperExact(ulong x, ulong n1, ulong n2, ulong n, ulong startAt = 0) {
-nvariant
+mmutable
-        logFactorial(n) + logFactorial(n1 + n2 - n) - logFactorial(n1 + n2);
-    real sum = 0;
@@ -663 +662 @@
-
-private {
-nvariant
+mmutable
-   -0x1.758f4d969484bfdcp-7,    // -0.011400139698853582732
-   0x1.53cee17a59259dd2p-3, // 0.16592193750979583221
@@ -674 +673 @@
-];
-
-nvariant
+mmutable
-   -0x1.64b92ae791e64bb2p-7,    // -0.010886331510064192632
-   0x1.7585c7d597298286p-3, // 0.1823840725000038842
@@ -685 +684 @@
-];
-
-nvariant
+mmutable
-   0x1.20ceea49ea142f12p-13,    // 0.00013771451113809605662
-   0x1.cbe8a7267aea80bp-7,  // 0.014035302749980729871
@@ -698 +697 @@
-];
-
-nvariant
+mmutable
-   0x1.3a4ce1406cea98fap-13,    // 0.00014987006762866754669
-   0x1.f45332623335cda2p-7, // 0.015268706895221911913
@@ -711 +710 @@
-];
-
-nvariant
+mmutable
-   0x1.8c124a850116a6d8p-21,    // 7.3774056430545041787e-07
-   0x1.534abda3c2fb90bap-13,    // 0.0001617870121822776094
@@ -722 +721 @@
-];
-
-nvariant
+mmutable
-   0x1.af03f4fc0655e006p-21,    // 8.0282885006885383316e-07
-   0x1.713192048d11fb2p-13, // 0.00017604524340842589303
@@ -733 +732 @@
-];
-
-nvariant
+mmutable
-   -0x1.55da447ae3806168p-34,   // -7.7728283809481633868e-11
-   -0x1.145635641f8778a6p-24,   // -6.4339663876133447143e-08
@@ -744 +743 @@
-];
-
-nvariant
+mmutable
-   -0x1.74022dd5523e6f84p-34,   // -8.4584942637876803775e-11
-   -0x1.2cb60d61e29ee836p-24,   // -7.0014768675591937804e-08
@@ -833 +832 @@
-    // Just make sure invNormalCDF works like it should as the inverse.
-    foreach(i; 0..1000) {
-gen, 
-gen, 
-gen, 
+
+
+
-        real inv = invNormalCDF(x, mean, sd);
-        real rec = normalCDF(inv, mean, sd);
-        assert(approxEqual(x, rec));
-    }
-stderr, 
+
-}
-
@@ -919 +918 @@
-assert(approxEqual(studentsTCDF(invStudentsTCDF(0.4L, 18), 18), .4L));
-assert(approxEqual(studentsTCDF( invStudentsTCDF(0.9L, 11), 11), 0.9L));
-stderr, 
+
-}
-
@@ -1161 +1160 @@
-    uint nSkipped;
-    foreach(i; 0..1000) {
-gen, 
-gen, 
-gen, 
+
+
+
-        real pVal = negBinomCDF(k, n, p);
-

trunk/gamma.d

@@ -5 +5 @@
- * is trying to integrate the rest of this library into Tango,
- * tango.math.GammaFunction would be a drop-in replacement.
-
-
+*/
+/
- *            Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
- * Authors:   Stephen L. Moshier (original C code). Conversion to D by Don Clugston

trunk/infotheory.d

@@ -3 +3 @@
- * quantities, i.e, entropy, mutual info, etc. are output in bits.
- *
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-module dstats.infotheory;
-
-
+
+
-
-import dstats.sort, dstats.summary, dstats.base, dstats.alloc;
@@ -42 +79 @@
-}
-
-n array
+ forward range
- * treated as frequency counts of a set of discrete observations.
- *
@@ -53 +90 @@
- * ---
- */
-
-
-
-
-
+
+
+
+
+
+
+
+
-    real nNeg1 = 1.0L / n;
-    real entropy = 0;
@@ -70 +110 @@
-
-unittest {
-
+.dup
-    assert(approxEqual(uniform3, log2(3)));
-
+.dup
-    assert(approxEqual(uniform4, 2));
-
-
-
+.dup
+.dup
+.dup
-    writefln("Passed entropyCounts unittest.");
-}
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-struct ObsEnt(T...) {
-atomsOfArrayTuple!(T)
+T
-
-    hash_t toHash() {
@@ -99 +255 @@
-    }
-
-
+ref 
-        foreach(ti, elem; this.tupleof) {
-            if(elem != rhs.tupleof[ti])
@@ -106 +262 @@
-        return true;
-    }
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- *
- * Examples:
@@ -118 +292 @@
- * assert(approxEqual(entropyFoo, log2(3)));
- * int[] bar = [1, 2, 3, 1, 2, 3, 1, 2, 3];
-jointEntropyFooBar = entropy(foo, bar
-jointEntropy
+HFooBar = entropy(joint(foo, bar)
+H
- * ---
- */
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-    TempAlloc.frameInit;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
-    TempAlloc.frameFree;
-
-    return ans;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    writeln("Passed entropy unittest.");
-}
+
+/**Calculate the conditional entropy H(data | cond).*/
+real condEntropy(T, U)(T data, U cond)
+if(isInputRange!(T) && isInputRange!(U)) {
+    return entropy(joint(data, cond)) - entropy(cond);
+}
+
+unittest {
+    // This shouldn't be easy to screw up.  Just really basic.
+    int[] foo = [1,2,2,1,1];
+    int[] bar = [1,2,3,1,2];
+    assert(approxEqual(entropy(foo) - condEntropy(foo, bar),
+           mutualInfo(foo, bar)));
+    writeln("Passed condEntroy unittest.");
+}
+
-
-
-/**Calculates the mutual information of two vectors of observations.
- */
-
-
-
-
-
+
+
+
-}
-
-unittest {
-    // Values from R, but converted from base e to base 2.
-
-
-
-
-
-
+.dup
+.dup
+.dup
+.dup
+.dup
+.dup
-
-    writeln("Passed mutualInfo unittest.");
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-}
-
-unittest {
-    // Values from Matlab mi package by Hanchuan Peng.
-, [3,1,2,3,4,2,1,2]
-
+.dup, [3,1,2,3,4,2,1,2].dup
+.dup
-    assert(approxEqual(res, 0.4387));
-
+
+
-    assert(approxEqual(res, 1.3510));
-    writeln("Passed condMutualInfo unittest.");
-}
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-    alias binaryFun!(compFun) comp;
-nvariant size_t
-nvariant real
+mmutable
+mmutable
-
-    real sum = 0.0L;
-    real nSame = 1.0L;
-
-
-
+
+
+
+
-            nSame++;
-        else {
@@ -257 +449 @@
-            sum -= p * log2(p);
-        }
-pos++
+last = elem
-    }
-    // Handle last run.
@@ -272 +464 @@
-}
-
-template atomsOfArrayTuple(T...) {
-    mixin("alias TypeTuple!(Mutable!(" ~ (atomType!(T[0])).stringof ~ ")" ~
-          atomsOfArrayTupleImpl!(T[1..$]) ~ " atomsOfArrayTuple;");
-}
-
-template atomsOfArrayTupleImpl(T...) {
-    static if(T.length == 0)
-        const char[] atomsOfArrayTupleImpl = ")";
-    else const char[] atomsOfArrayTupleImpl = ", Mutable!(" ~
-                      (atomType!(T[0])).stringof ~ ")" ~
-                      atomsOfArrayTupleImpl!(T[1..$]);
-}
-
-template atomType(T : T[]) {
-    alias T atomType;
-}
-
-// Verify that there are no TempAlloc memory leaks anywhere in the code covered
-// by the unittest.  This should always be the last unittest of the module.

trunk/random.d

@@ -2 +2 @@
- * distributions.
- *
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -30 +62 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-
-
@@ -48 +81 @@
-
-///Generates a random number from normal(mean, sd) distribution.
-)(ref RGen gen, real mean = 0.0L, real sd = 1.0L
-gen, 0.0L, 1.0L)
+ = Random)(real mean = 0.0L, real sd = 1.0L, ref RGen gen = rndGen
+0.0L, 1.0L, gen);
-    return invNormalCDF(p) * sd + mean;
-}
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real[] observ = new real[100_000];
-    foreach(ref elem; observ)
-gen
+
-    auto ksRes = ksPval(observ, parametrize!(normalCDF)(0.0L, 1.0L));
-    writeln("100k samples from normal(0, 1):  K-S P-val:  ", ksRes);
@@ -70 +101 @@
-
-///Random Cauchy-distributed number.
-)(ref RGen gen, real X0 = 0, real gamma = 1
-gen, 0.0L, 1.0L)
+ = Random)(real X0 = 0, real gamma = 1, ref RGen gen = rndGen
+0.0L, 1.0L, gen);
-    return invCauchyCDF(p, X0, gamma);
-}
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real[] observ = new real[100_000];
-    foreach(ref elem; observ)
-gen
+
-    auto ksRes = ksPval(observ, parametrize!(cauchyCDF)(0.0L, 1.0L));
-    writeln("100k samples from Cauchy(0, 1):  K-S P-val:  ", ksRes);
@@ -91 +120 @@
-}
-
-)(ref RGen gen, real df
-gen, 0.0L, 1.0L
+ = Random)(real df, ref RGen gen = rndGen
+0.0L, 1.0L, gen
-    return invStudentsTCDF(pVal, df);
-}
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real[] observ = new real[10_000];
-    foreach(ref elem; observ)
-gen, 
+
-    auto ksRes = ksPval(observ, parametrize!(studentsTCDF)(5));
-    writeln("10k samples from T(5):  K-S P-val:  ", ksRes);
@@ -113 +140 @@
-
-/**Generates a random number from the Poisson distribution.*/
-)(ref RGen gen, real lambda
-gen, 0.0L, 1.0L
+ = Random)(real lambda, ref RGen gen = rndGen
+0.0L, 1.0L, gen
-    uint guess = cast(uint) max(round(
-          invNormalCDF(pVal, lambda, sqrt(lambda)) + 0.5), 0.0L);
@@ -142 +169 @@
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real lambda = 4L;
-    uint[] observ = new uint[100_000];
-    foreach(ref elem; observ)
-
-
-
-
+
+
-    writeln("\tMean Expected: ", lambda,
-            "  Observed:  ", mean(observ));
@@ -164 +187 @@
-
-///Bernoulli r.v. with probability p of equaling 1.
-)(ref RGen gen, real P = 0.5
-gen, 0.0L, 1.0L
+ = Random)(real P = 0.5, ref RGen gen = rndGen
+0.0L, 1.0L, gen
-    return cast(uint) (pVal <= P);
-}
@@ -171 +194 @@
-// bother writing a separate test for.
-
-import std.stdio;
-long gue;
-
-/**Generates a random number from the binionial distribution.*/
-)(ref RGen gen, uint n, real p
+ = Random)(uint n, real p, ref RGen gen = rndGen
-    // Generate p-value, get normal approx. as starting point, search for
-    // binomial element that matches this starting point.  Normal approx.
-    // is usually within 1 or 2 elements.
-gen, 0.0L, 1.0L
+0.0L, 1.0L, gen
-    uint guess = cast(uint) max(round(
-          invNormalCDF(pVal, n * p, sqrt(n * p * (1 - p)))) + 0.5, 0);
@@ -211 +231 @@
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real p = 0.7;
-    uint n = 10;
-    uint[] observ = new uint[100_000];
-    foreach(ref elem; observ)
-gen, 
+
-    auto ksRes = ksPval(observ, parametrize!(binomialCDF)(n, p));
-
-
+
-    writeln("\tMean Expected: ", n * p,
-            "  Observed:  ", mean(observ));
@@ -237 +254 @@
- * Bugs:  Slow for large values of N.  Fairly naive implementation.
- * Reasonably fast, though, for small values (<100) of N.*/
-)(ref RGen gen, uint N1, uint N2, uint N
+ = Random)(uint N1, uint N2, uint N, ref RGen gen = rndGen
-in {
-    assert(N <= (N1 + N2));
@@ -247 +264 @@
-    NR[1] = N1;
-    foreach(i; 0..N) {
-gen, NR[1] / (NR[0] + NR[1])
+NR[1] / (NR[0] + NR[1]), gen
-        NR[X] -= 1.0L;
-        result += X;
@@ -255 +272 @@
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    uint n1 = 20, n2 = 30, n = 10;
-    uint[] observ = new uint[100_000];
-    foreach(ref elem; observ)
-gen, 
+
-    auto ksRes = ksPval(observ, parametrize!(hypergeometricCDF)(n1, n2, n));
-
-
+
-    writeln("\tMean Expected: ", n * cast(real) n1 / (n1 + n2),
-            "  Observed:  ", mean(observ));
@@ -278 +292 @@
- * set n = 1, since geometric distribution is just a special case of neg.
- * binomial.*/
-)(ref RGen gen, uint n, real p
+ = Random)(uint n, real p, ref RGen gen = rndGen
-in {
-    assert(p >= 0 && p <= 1);
-    assert(n > 0);
-} body {
-gen, 0.0L, 1.0L
+0.0L, 1.0L, gen
-    // Normal or gamma approx, then adjust.
-    real mean = n * (1 - p) / p;
@@ -348 +362 @@
-    uint[] observ = new uint[100_000];
-    foreach(ref elem; observ)
-
-
-
-
+
+
-    writeln("\tMean Expected: ", n * (1 - p) / p,
-            "  Observed:  ", mean(observ));
@@ -365 +377 @@
-
-///
-)(ref RGen gen, real mu = 0, real b = 1
-gen, 0.0L, 1.0L
+ = Random)(real mu = 0, real b = 1, ref RGen gen = rndGen
+0.0L, 1.0L, gen
-    return invLaplaceCDF(p, mu, b);
-}
@@ -375 +387 @@
-    real[] observ = new real[100_000];
-    foreach(ref elem; observ)
-gen
+
-    auto ksRes = ksPval(observ, parametrize!(laplaceCDF)(0.0L, 1.0L));
-    writeln("100k samples from Laplace(0, 1):  K-S P-val:  ", ksRes);
@@ -389 +401 @@
- * gamma distribution, with b = 1.  However, it's a common special
- * case and this function is a lot faster than the more general rGamma.*/
-)(ref RGen gen, real lambda
-gen, 0.0L, 1.0L
+ = Random)(real lambda, ref RGen gen = rndGen
+0.0L, 1.0L, gen
-    return -log(p) / lambda;
-}
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real[] observ = new real[100_000];
-    foreach(ref elem; observ)
-gen, 
+
-    auto ksRes = ksPval(observ, parametrize!(gammaCDF)(2, 1));
-    writeln("100k samples from exponential(2):  K-S P-val:  ", ksRes);
@@ -414 +424 @@
- * this is a common special case, a separate exponential random
- * number generator is provided for speed.*/
-)(ref RGen gen, real a = 1, real b = 1
-gen, 0.0L, 1.0L
+ = Random)(real a = 1, real b = 1, ref RGen gen = rndGen
+0.0L, 1.0L, gen
-    return invGammaCDFR(p, a, b);
-}
-
-unittest {
-    Random gen;
-    gen.seed(unpredictableSeed);
-    real[] observ = new real[100_000];
-    foreach(ref elem; observ)
-gen, 
+
-    auto ksRes = ksPval(observ, parametrize!(gammaCDF)(2, 3));
-    writeln("100k samples from gamma(2, 3):  K-S P-val:  ", ksRes);

trunk/sort.d

@@ -23 +23 @@
- * ---
- *
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -51 +83 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-
-module dstats.sort;
-
-import std.traits, std.algorithm, std.math, std.functional, std.math,
-
+, std.range, std.array
-
-import dstats.alloc;
@@ -62 +95 @@
-version(unittest) {
-    import std.stdio, std.random;
-    Random gen;
-
-    void main (){
@@ -68 +100 @@
-}
-
-
+
+
-    if(input.length < 2) return;
-T
+ElementType!(T)
-    foreach(i; 1..input.length) {
-        input[i-1] = input[i];
@@ -77 +110 @@
-}
-
-
+
+
-    if(input.length < 2) return;
-T
+ElementType!(T)
-    for(size_t i = input.length - 1; i > 0; i--) {
-        input[i] = input[i-1];
-    }
-    input[0] = temp;
-}
-
-// For testing purposes for stable sorts.  Shuffles N arrays in lockstep.
-version(unittest) {
-    void randomMultiShuffle(SomeRandomGen, T...)(ref SomeRandomGen r, T array) {
-        foreach (i; 0 .. array[0].length) {
-            // generate a random number i .. n
-            immutable which = i + uniform!(size_t)(r, 0u, array[0].length - i);
-            foreach(ti, element; array) {
-                swap(element[i], element[which]);
-            }
-        }
-    }
-}
-
@@ -242 +263 @@
-
-unittest {
-    gen.seed(unpredictableSeed);
-    {  // Test integer.
-        uint[] test = new uint[1_000];
-        foreach(ref e; test) {
-gen, 
+
-        }
-        auto test2 = test.dup;
-        foreach(i; 0..1_000) {
-MultiShuffle(gen, test, test2
-gen, 
+Shuffle(zip(test, test2)
+
-            qsort(test[0..len], test2[0..len]);
-            assert(isSorted(test[0..len]));
@@ -260 +280 @@
-        double[] test = new double[1_000];
-        foreach(ref e; test) {
-gen, 
+
-        }
-        auto test2 = test.dup;
-        foreach(i; 0..1_000) {
-MultiShuffle(gen, test, test2
-gen, 
+Shuffle(zip(test, test2)
+
-            qsort!("a > b")(test[0..len], test2[0..len]);
-            assert(isSorted!("a > b")(test[0..len]));
@@ -330 +350 @@
-    uint[] temp1 = new uint[1_000], temp2 = new uint[1_000];
-    foreach(ref e; test) {
-gen, 
+
-    }
-    foreach(i; 0..100) {
@@ -337 +357 @@
-            e = j;
-        }
-MultiShuffle(gen, 
-gen, 
+Shuffle(
+
-        // Testing bubble sort distance against bubble sort,
-        // since bubble sort distance computed by bubble sort
@@ -361 +381 @@
-            e = j;
-        }
-MultiShuffle(gen, 
-gen, 
+Shuffle(
+
-        if(i & 1)  // Test both temp and non-temp branches.
-            mergeSort(test[0..len], stability[0..len]);
@@ -378 +398 @@
-        double[] testL = new double[1_000];
-        foreach(ref e; testL) {
-gen, 
+
-        }
-        auto testL2 = testL.dup;
-        foreach(i; 0..1_000) {
-MultiShuffle(gen, testL, testL2
-gen, 
+Shuffle(zip(testL, testL2)
+
-            mergeSort!("a > b")(testL[0..len], testL2[0..len]);
-            assert(isSorted!("a > b")(testL[0..len]));
@@ -562 +582 @@
-    uint[] test = new uint[1_000], stability = new uint[1_000];
-    foreach(ref e; test) {
-gen, 
+
-    }
-    uint[] test2 = test.dup;
@@ -569 +589 @@
-            e = j;
-        }
-MultiShuffle(gen, test, test2
-gen, 
+Shuffle(zip(test, test2)
+
-        mergeSortInPlace(test[0..len], test2[0..len], stability[0..len]);
-        assert(isSorted(test[0..len]));
@@ -644 +664 @@
-
-    foreach(array; data) {
-rotate(array[half1..middle + half2], array.ptr + middle
+bringToFront(array[half1..middle], array[middle..middle + half2]
-    }
-    size_t newMiddle = half1 + half2;
@@ -688 +708 @@
-
-unittest {
-    gen.seed(unpredictableSeed);
-    uint[] test = new uint[1_000];
-    foreach(ref e; test) {
-gen, 
+
-    }
-    auto test2 = test.dup;
-    foreach(i; 0..1_000) {
-MultiShuffle(gen, test, test2
-gen, 
+Shuffle(zip(test, test2)
+
-        heapSort(test[0..len], test2[0..len]);
-        assert(isSorted(test[0..len]));
@@ -769 +788 @@
-
-unittest {
-    gen.seed(unpredictableSeed);
-    uint[] test = new uint[100], stability = new uint[100];
-    foreach(ref e; test) {
-gen, 
+
-    }
-    foreach(i; 0..1_000) {
@@ -779 +797 @@
-            e = j;
-        }
-MultiShuffle(gen, 
-gen, 
+Shuffle(
+
-        // Testing bubble sort distance against bubble sort,
-        // since bubble sort distance computed by bubble sort
@@ -932 +950 @@
-
-unittest {
-    gen.seed(unpredictableSeed);
-    enum n = 1000;
-    uint[] test = new uint[n];
@@ -938 +955 @@
-    uint[] lockstep = new uint[n];
-    foreach(ref e; test) {
-gen, 
+
-    }
-    foreach(i; 0..1_000) {
-        randomShuffle(test, gen);
-        test2[] = test[];
-        lockstep[] = test[];
-gen, 
+
-        qsort!("a > b")(test2[0..len]);
-gen, 
+
-        auto qsRes = partitionK!("a > b")(test[0..len], lockstep[0..len], k);
-        assert(qsRes == test2[k]);
@@ -960 +976 @@
-}
-
-addElement function,
+put function, this output range
- * maintains the invariant that the top N according to compFun will be
- * contained in the data structure.  Uses a heap internally, O(log N) insertion
@@ -977 +993 @@
- * randomShuffle(nums, gen);
- * foreach(n; nums) {
-addElemen
-addElemen
+pu
+pu
- * }
- *  assert(less.getSorted == [0U, 1,2,3,4,5,6,7,8,9]);
@@ -999 +1015 @@
-
-    /** Insert an element into the topN struct.*/
-addElemen
+pu
-        if(nAdded < n) {
-            nodes[nAdded] = elem;
@@ -1041 +1057 @@
-        randomShuffle(nums, gen);
-        foreach(n; nums) {
-addElemen
-addElemen
+pu
+pu
-        }
-        assert(less.getSorted == [0U, 1,2,3,4,5,6,7,8,9]);
@@ -1052 +1068 @@
-        randomShuffle(nums, gen);
-        foreach(n; nums[0..5]) {
-addElemen
-addElemen
+pu
+pu
-        }
-        assert(less.getSorted == qsort!("a < b")(nums[0..5]));

trunk/stats.d

@@ -1 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-
-module dstats.stats;

trunk/summary.d

@@ -1 @@
-
+
+
+
- *
- * Bugs:  This whole module assumes that input will be reals or types implicitly
- *        convertible to real.  No allowances are made for user-defined numeric
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -33 +70 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-
-
-module dstats.summary;
-
-
+
+
-
-import dstats.sort, dstats.base, dstats.alloc;
-
-version(unittest) {
-
-
-
+
-
-    void main() {
@@ -51 +88 @@
-}
-
-in O(N) time on average.  In the case of an even number of
-lements, the mean of the two middle elements is returned.  This is a
-convenience founction designed specifically for numeric types, where the
-averaging of the two middle elements is desired.  A more general selection
-algorithm that can handle any type with a total ordering, as well as
-
+of an input range in O(N) time on average.  In the case of an
+ven number of elements, the mean of the two middle elements is returned.
+This is a convenience founction designed specifically for numeric types,
+where the averaging of the two middle elements is desired.  A more general
+selection algorithm that can handle any type with a total ordering, as well
+as 
- * dstats.sort.quickSelect() and dstats.sort.partitionK().
- * Allocates memory, does not reorder input data.*/
-
-
-
+
+
+
+
+
+
+
-    return medianPartition(dataDup);
-}
@@ -69 +110 @@
- * median, and elements larger than the median will have larger indices than
- * that of the median. Useful both for its partititioning and to avoid
-
-
+
+
+
+
+
+
+
-in {
-    assert(data.length > 0);
@@ -77 +123 @@
-    // with an index larger than the lower median, after the array is
-    // partially sorted.
-    static T min(T[] data) {
-        T min = data[0];
-        foreach(d; data[1..$]) {
-            if(d < min)
-                min = d;
-        }
-        return min;
-    }
-
-    if(data.length == 1) {
-        return data[0];
@@ -92 +129 @@
-    } else {
-        auto lower = partitionK(data, data.length / 2 - 1);
-
+
+
+
+
+
+
+
+
-        return lower * 0.5L + upper * 0.5L;
-    }
@@ -105 +149 @@
-    }
-
-    gen.seed(unpredictableSeed);
-    float[] test = new float[1000];
-    uint upperBound, lowerBound;
-    foreach(testNum; 0..1000) {
-        foreach(ref e; test) {
-gen, 
+
-        }
-        do {
-gen, 
-gen, 
+
+
-        } while(lowerBound == upperBound);
-        if(lowerBound > upperBound) {
@@ -126 +169 @@
-        assert(approxEqual(quickRes, accurateRes));
-    }
+
+    // Make sure everything works with lowest common denominator range type.
+    struct Count {
+        uint num;
+        uint upTo;
+        uint front() {
+            return num;
+        }
+        void popFront() {
+            num++;
+        }
+        bool empty() {
+            return num >= upTo;
+        }
+    }
+
+    Count a;
+    a.upTo = 100;
+    assert(approxEqual(median(a), 49.5));
-    writeln("Passed median unittest.");
-}
-
-
-
+
+
+
+
-    OnlineMean meanCalc;
-    foreach(element; data) {
-addElemen
+pu
-    }
-    return meanCalc.mean;
-}
-
-Struct
+Output range
- * check for N == 0.  This struct uses O(1) space and does *NOT* store the
- * individual elements.
@@ -145 +209 @@
- * ---
- * OnlineMean summ;
-addElemen
-addElemen
-addElemen
-addElemen
-addElemen
+pu
+pu
+pu
+pu
+pu
- * assert(summ.mean == 3);
- * ---*/
@@ -157 +221 @@
-    real k = 0;
-public:
-
-
+
+
+
+
+
-        result += (element - result) / ++k;
-    }
@@ -178 +245 @@
-}
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * on large array summing operations.  However, by default, return type is
- * T (same as input type).*/
-
+
+
-    U sum = 0;
-    foreach(value; data) {
@@ -198 +305 @@
-}
-
-/**User has option of making U a different type than T
- * to prevent overflows on large array summing operations.
- * However, by default, return type is T (same as input type).*/
-U absSum(T, U = Mutable!(T))(const T[] data) {
-    U sum=0;
-    foreach(value; data) {
-        sum+=abs(value);
-    }
-    return sum;
-}
-
-unittest {
-
-
-
-
-
-
-
+
+
+
+
+
-    writefln("Passed sum/mean unittest.");
-}
-
-
-Allows computation of
+Outpu range to compute
- * for stdev, var cost a few floating point ops.  Getter for mean costs
- * a single branch to check for N == 0.  Relatively expensive floating point
@@ -230 +324 @@
- * ---
- * OnlineMeanSD summ;
-addElemen
-addElemen
-addElemen
-addElemen
-addElemen
+pu
+pu
+pu
+pu
+pu
- * assert(summ.mean == 3);
- * assert(summ.stdev == sqrt(2.5));
@@ -246 +340 @@
-public:
-    ///
-addElemen
+pu
-        real kNeg1 = 1.0L / ++_k;
-        _var += (element * element - _var) * kNeg1;
@@ -292 +386 @@
-}
-
-
-
+
+
+
+
-    return meanVariance(data).SD;
-}
-
-
-
+
+
+
+
-    OnlineMeanSD meanSDCalc;
-    foreach(element; data) {
-addElemen
+pu
-    }
-
@@ -307 +405 @@
-}
-
-
-
+
+
+
+
-    auto ret = meanVariance(data);
-    ret.SD = sqrt(ret.SD);
@@ -314 +414 @@
-}
-
-
-
+
+
+
+
-    return meanStdev(data).SD;
-}
-
-unittest {
-
+cast(int[]) 
-    assert(approxEqual(res.SD, 1.7078));
-    assert(approxEqual(res.mean, 3.25));
-
+cast(double[]) 
-    assert(approxEqual(res.SD, 1.5811));
-    assert(approxEqual(res.mean, 3));
@@ -329 +431 @@
-}
-
-
-
-
-
-
-
-
-
+
- * max online. Using this struct is relatively expensive, so if you just need
- * mean and/or stdev, try OnlineMeanSD or OnlineMean. Getter methods for stdev,
@@ -347 +442 @@
- * ---
- * OnlineSummary summ;
-addElemen
-addElemen
-addElemen
-addElemen
-addElemen
+pu
+pu
+pu
+pu
+pu
- * assert(summ.N == 5);
- * assert(summ.mean == 3);
@@ -371 +466 @@
-public:
-    ///
-addElemen
-nvariant
+pu
+mmutable
-        _min = (element < _min) ? element : _min;
-        _max = (element > _max) ? element : _max;
@@ -438 +533 @@
- * the variance is due to infrequent, large deviations from the mean.  Low
- * kurtosis means that the variance is due to frequent, small deviations from
-
-
+
+
+
+
-    OnlineSummary kCalc;
-    foreach(elem; data) {
-addElemen
+pu
-    }
-    return kCalc.kurtosis;
@@ -449 +546 @@
-unittest {
-    // Values from Octave.
-
-
-
+.dup
+.dup
+.dup
-    writefln("Passed kurtosis unittest.");
-}
@@ -458 +555 @@
- * means that the right tail is longer/fatter than the left tail.  Negative
- * skewness means the left tail is longer/fatter than the right tail.  Zero
-
-
+
+
+
+
-    OnlineSummary sCalc;
-    foreach(elem; data) {
-addElemen
+pu
-    }
-    return sCalc.skewness;
@@ -469 +568 @@
-unittest {
-    // Values from Octave.
-
-
-
+
+
+
+
+
+
+
+
-    writeln("Passed skewness test.");
-}
@@ -510 +614 @@
-}
-
-
-
-
+
+
+
+
+
-    OnlineSummary summ;
-    foreach(elem; data) {
-addElemen
+pu
-    }
-    return Summary(data.length, summ.mean, summ.var, summ.stdev, summ.skewness,

trunk/tests.d

@@ -1 @@
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- * All rights reserved.
- *
@@ -29 +62 @@
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-
+
+
-module dstats.tests;
-
-import dstats.base, dstats.distrib, dstats.alloc, dstats.summary, dstats.sort,
-
+, std.range
-
-version(unittest) {
@@ -61 +95 @@
- * Alt.GREATER, meaning mean(data) > mean, and Alt.TWOSIDE, meaning mean(data)
- * != mean.
+ *
- * Returns:  The p-value against the given alternative.*/
-
+
+
-    auto meanSd = meanStdev(data);
-    real t = (meanSd.mean - mean) / (meanSd.SD / sqrt(cast(real) data.length));
@@ -75 +111 @@
-
-unittest {
-
-
-
+.dup
+.dup
+.dup
-    writeln("Passed 1-sample studentsTTest test.");
-}
@@ -86 +122 @@
- * mean(sample2), and Alt.TWOSIDE, meaning mean(sample1) != mean(sample2).
- * Returns:  The p-value against the given alternative.*/
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-                 (n1 + n2 - 2));
-    real t = (s1summ.mean - s2summ.mean) /
@@ -107 +150 @@
-unittest {
-    // Values from R.
-, [1,3,4,5,7,9]
-, [1,3,4,5,7,9]
+.dup, [1,3,4,5,7,9].dup
+.dup, [1,3,4,5,7,9].dup
-           0.1173));
-, [1,3,4,5,7,9]
+.dup, [1,3,4,5,7,9].dup
-           0.8827));
-, [2,2,1,3,4]
-, [2,2,1,3,4]
+.dup, [2,2,1,3,4].dup
+.dup, [2,2,1,3,4].dup
-           0.965));
-, [2,2,1,3,4]
+.dup, [2,2,1,3,4].dup
-           0.03492));
-    writeln("Passed 2-sample studentsTTest test.");
@@ -126 +169 @@
- * != mean(sample2).
- * Returns:  The p-value against the given alternative.*/
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    real t = (s1summ.mean - s2summ.mean) / sx1x2;
-s1summ.SD / n1 + s2summ.SD
+v1 / n1 + v2
-    numerator *= numerator;
-s1summ.SD
+v1
-    denom1 = denom1 * denom1 / (n1 - 1);
-s2summ.SD
+v2
-    denom2 = denom2 * denom2 / (n2 - 1);
-    real df = numerator / (denom1 + denom2);
@@ -151 +202 @@
-unittest {
-        // Values from R.
-, [1,3,4,5,7,9]
-, [1,3,4,5,7,9]
+.dup, [1,3,4,5,7,9].dup
+.dup, [1,3,4,5,7,9].dup
-           0.1079));
-, [1,3,4,5,7,9]
+.dup, [1,3,4,5,7,9].dup
-           0.892));
-, [2,2,1,3,4]
-, [2,2,1,3,4]
+.dup, [2,2,1,3,4].dup
+.dup, [2,2,1,3,4].dup
-           0.967));
-, [2,2,1,3,4]
+.dup, [2,2,1,3,4].dup
-           0.03308));
-    writeln("Passed welchTTest test.");
@@ -170 +221 @@
- * greater than testMean, and Alt.TWOSIDE, meaning the true mean difference is not
- * equal to testMean.
+ *
- * Returns:  The p-value against the given alternative.*/
-
-
-
-
-
+
+
-    OnlineMeanSD msd;
-
-
-
-
-
+
+
+
+
+
+
+
+
+
-
-    if(alt == Alt.LESS) {
-before.length
+len
-    } else if(alt == Alt.GREATER) {
-before.length
+len
-    } else if(t > 0) {
-before.length
+len
-    } else {
-before.length
+len
-    }
-}
@@ -196 +249 @@
-unittest {
-    // Values from R.
-, [2,3,5,5,6]
-, [2,3,5,5,6]
-, [2,3,5,5,6]
-, [2,3,5,5,6]
-, [2,3,5,5,6]
-, [2,3,5,5,6]
+.dup, [2,3,5,5,6].dup
+.dup, [2,3,5,5,6].dup
+.dup, [2,3,5,5,6].dup
+.dup, [2,3,5,5,6].dup
+.dup, [2,3,5,5,6].dup
+.dup, [2,3,5,5,6].dup
-    writeln("Passed pairedTTest unittest.");
-}
@@ -210 +263 @@
- * dereference that can be used in wilcoxonRankSumPval() to adjust for ties in
- * the input data.
+ *
+ * Input ranges for this function must define a length.
+ *
- * Returns:  The Wilcoxon test statistic W.*/
-
-
-
-
-
+
+
+
+
+
+
-
-    float[] ranks = newStack!(float)(N);
-    rankSort(combined, ranks);
-sum!(float, real)(
+reduce!("a + b")(0.0L, 
-    TempAlloc.free;  // Free ranks.
-
@@ -248 +305 @@
-
-unittest {
-, [2, 4, 6, 8, 10]
-, [1, 2, 3, 4, 5]
-, [2, 4, 6, 8, 10]
+.dup, [2, 4, 6, 8, 10].dup
+.dup, [1, 2, 3, 4, 5].dup
+.dup, [2, 4, 6, 8, 10].dup
-    writeln("Passed wilcoxonRankSum test.");
-}
@@ -298 +355 @@
-    assert(approxEqual(wilcoxonRankSumPval(1500, 50, 50, Alt.LESS), .957903));
-    assert(approxEqual(wilcoxonRankSumPval(8500, 100, 200, Alt.LESS), .01704));
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6), 0.9273));
-    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6, Alt.GREATER), 0.4636));
@@ -325 +382 @@
- * In these cases, exactThresh will be ignored.
- *
+ * Input ranges for this function must define a length.
+ *
- * Returns:  The p-value against the given alternative.*/
-
-
+
+
+
-    real tieSum;
-    real W = wilcoxonRankSum(sample1, sample2, &tieSum);
@@ -337 +397 @@
-     // Values from R.  Simple stuff (no ties) first.  Testing approximate
-     // calculation first.
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-           Alt.TWOSIDE, 0), 0.9273));
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-           Alt.LESS, 0), 0.6079));
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-           Alt.GREATER, 0), 0.4636));
-, [3,5,7,8,13,15]
+.dup, [3,5,7,8,13,15].dup
-            Alt.TWOSIDE, 0), 0.4113));
-, [3,5,7,8,13,15]
+.dup, [3,5,7,8,13,15].dup
-            Alt.LESS, 0), 0.2057));
-, [3,5,7,8,13,15]
+.dup, [3,5,7,8,13,15].dup
-            Alt.GREATER, 0), 0.8423));
-, [2,4,6,8,10]
+.dup, [2,4,6,8,10].dup
-            Alt.TWOSIDE, 0), .6745));
-, [2,4,6,8,10]
+.dup, [2,4,6,8,10].dup
-            Alt.LESS, 0), .3372));
-, [2,4,6,8,10]
+.dup, [2,4,6,8,10].dup
-            Alt.GREATER, 0), .7346));
-
-    // Now, lots of ties.
-, [2,3,4,5,6]
+.dup, [2,3,4,5,6].dup
-           Alt.TWOSIDE, 0), 0.3976));
-, [2,3,4,5,6]
+.dup, [2,3,4,5,6].dup
-           Alt.LESS, 0), 0.1988));
-, [2,3,4,5,6]
+.dup, [2,3,4,5,6].dup
-           Alt.GREATER, 0), 0.8548));
-, [1,2,3,1,1]
+.dup, [1,2,3,1,1].dup
-           Alt.TWOSIDE, 0), 0.9049));
-, [1,2,3,1,1]
+.dup, [1,2,3,1,1].dup
-           Alt.LESS, 0), 0.4524));
-, [1,2,3,1,1]
+.dup, [1,2,3,1,1].dup
-           Alt.GREATER, 0), 0.64));
-
-    // Now, testing the exact calculation on the same data.
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-       Alt.TWOSIDE), 0.9307));
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-           Alt.LESS), 0.6039));
-, [1,3,5,7,11,9]
+.dup, [1,3,5,7,11,9].dup
-           Alt.GREATER), 0.4654));
-, [3,5,7,8,13,15]
+.dup, [3,5,7,8,13,15].dup
-            Alt.TWOSIDE), 0.4286));
-, [3,5,7,8,13,15]
+.dup, [3,5,7,8,13,15].dup
-            Alt.LESS), 0.2143));
-, [3,5,7,8,13,15]
+.dup, [3,5,7,8,13,15].dup
-            Alt.GREATER), 0.8355));
-, [2,4,6,8,10]
+.dup, [2,4,6,8,10].dup
-            Alt.TWOSIDE), .6905));
-, [2,4,6,8,10]
+.dup, [2,4,6,8,10].dup
-            Alt.LESS), .3452));
-, [2,4,6,8,10]
+.dup, [2,4,6,8,10].dup
-            Alt.GREATER), .7262));
-    writeln("Passed wilcoxonRankSumPval test.");
@@ -501 +561 @@
- * dereference that can be used in wilcoxonRankSumPval() to adjust for ties in
- * the input data.
+ *
+ * Input ranges for this function must define a length.
+ *
- * Returns:  The Wilcoxon test statistic W.*/
-
-
+
+
+
-    mixin(newFrame);
-    float[] diffRanks = newStack!(float)(before.length);
@@ -517 +581 @@
-    }
-
-
-
-
-
+
+
+
+
+
+
+
+
-    }
-    rankSort(diffs, diffRanks);
@@ -556 +624 @@
-
-unittest {
-, [2,1,4,5,3]
-, [2,7,1,8,2]
-, [0,9,8,6,7]
+.dup, [2,1,4,5,3].dup
+.dup, [2,7,1,8,2].dup
+.dup, [0,9,8,6,7].dup
-    writeln("Passed wilcoxonSignedRank unittest.");
-}
@@ -583 +651 @@
- * or equal pairs.
- *
+ * The input ranges for this function must define a length and must be
+ * forward ranges.
+ *
- * Returns:  The p-value against the given alternative.*/
-const T[] before, const T[]
+T before, T
-      Alt alt = Alt.TWOSIDE, uint exactThresh = 50)
+if(realInput!(T) && dstats.base.hasLength!(T) && isForwardRange!(T))
-in {
-    assert(before.length == after.length);
-} body {
-  
-  
-  
-  foreach(i; 0..before.length
-  if(before[i] == after[i]
+
+
+
+while(!before.empty && !after.empty
+if(before.front == after.front
-            N--;
-
-
+
+
+
+
+
-}
-
@@ -603 +678 @@
-    alias approxEqual ae;
-    // With ties, normal approx.
-, [2,1,4,5,3]
-, [2,7,1,8,2]
-, [0,9,8,6,7]
-, [2,1,4,5,3]
-, [2,7,1,8,2]
-, [0,9,8,6,7]
-, [2,1,4,5,3]
-, [2,7,1,8,2]
-, [0,9,8,6,7]
+.dup, [2,1,4,5,3].dup
+.dup, [2,7,1,8,2].dup
+.dup, [0,9,8,6,7].dup
+.dup, [2,1,4,5,3].dup
+.dup, [2,7,1,8,2].dup
+.dup, [0,9,8,6,7].dup
+.dup, [2,1,4,5,3].dup
+.dup, [2,7,1,8,2].dup
+.dup, [0,9,8,6,7].dup
-
-    // Exact.
-, [2,-4,-8,16,32]
-, [2,-4,-8,16,32]
-, [2,-4,-8,16,32]
-, [2,-4,-8,-16,32]
-, [2,-4,-8,-16,32]
-, [2,-4,-8,-16,32]
+.dup, [2,-4,-8,16,32].dup
+.dup, [2,-4,-8,16,32].dup
+.dup, [2,-4,-8,16,32].dup
+.dup, [2,-4,-8,-16,32].dup
+.dup, [2,-4,-8,-16,32].dup
+.dup, [2,-4,-8,-16,32].dup
-    writeln("Passed wilcoxonSignedRankPval unittest.");
-}
@@ -755 +830 @@
- * elements of after, and Alt.TWOSIDE, meaning that there is a significant
- * difference in either direction.*/
-
-
-
-
+
+
-    uint greater, less;
-foreach(i; 0..before.length
-[i] < after[i]
+while(!before.empty && !after.empty
+.front < after.front
-            less++;
-[i] < before[i]
+.front < before.front
-            greater++;
-        // Ignore equals.
+        before.popFront;
+        after.popFront;
-    }
-    if(alt == Alt.LESS) {
@@ -780 +855 @@
-unittest {
-    alias approxEqual ae;
-, [1,2,4,8,16]
-, [1,2,4,8,16]
-, [1,2,4,8,16]
-, [1,2,3,4,5]
-, [1,2,3,4,5]
-, [1,2,3,4,5]
+.dup, [1,2,4,8,16].dup
+.dup, [1,2,4,8,16].dup
+.dup, [1,2,4,8,16].dup
+.dup, [1,2,3,4,5].dup
+.dup, [1,2,3,4,5].dup
+.dup, [1,2,3,4,5].dup
-}
-
@@ -795 +870 @@
- * elements of after, and Alt.TWOSIDE, meaning that there is a significant
- * difference in either direction.*/
-const T[]
-!isArray
+T
+realInput!(T) && isNumeric
-    uint greater, less;
-i; 0..data.length
-data[i]
+elem; data
+elem
-            less++;
-data[i]
+elem
-            greater++;
-        // Ignore equals.
@@ -817 +892 @@
-
-unittest {
-
+.dup
-    writeln("Passed signTest unittest.");
-}
@@ -828 +903 @@
- * To get the D value used in standard notation,
- * simply take the absolute value of this D value.*/
-
+
+
-    auto TAState = TempAlloc.getState;
-    scope(exit) {
@@ -834 +910 @@
-        TempAlloc.free(TAState);
-    }
-F.tempdup(TAState), Fprime.tempdup(TAStat
+tempdup(F), tempdup(Fprim
-}
-
-/**Same as ksTest, but sorts input data in place instead of duplicating
-
-
+
+
+
+
+
-    qsort(F);
-    qsort(Fprime);
@@ -865 +944 @@
-unittest {
-    // Values from R.
-, [1,2,3,4,5]
-, [1,2,2,3,5]
-, [1,2,2,3,5]
-, [1,2,2,3,5,7,8]
-
-
+.dup, [1,2,3,4,5].dup
+.dup, [1,2,2,3,5].dup
+.dup, [1,2,2,3,5].dup
+.dup, [1,2,2,3,5,7,8].dup
+.dup
+.dup
-    writeln("Passed 2-sample ksTestDestructive.");
-}
@@ -885 +964 @@
- * ---
- */
-const T[]
-is(Func == function) || is(Func == delegate
+T
+realInput!(T) && (is(Func == function) || is(Func == delegate)
-    scope(exit) TempAlloc.free;
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-    qsort(Femp);
-    real D = 0;
@@ -909 +990 @@
-    // Testing against values from R.
-    auto stdNormal = parametrize!(normalCDF)(0.0L, 1.0L);
-
-
-
+.dup
+.dup
+.dup
-    assert(approxEqual(ksTestDestructive(lotsOfTies, stdNormal), -0.8772));
-    writeln("Passed 1-sample ksTestDestructive.");
@@ -940 +1021 @@
-/**KS-test, returns 2-sided P-val instead of D.
- * Bugs:  Exact calculation not implemented.  Uses asymptotic approximation.*/
-
+
+
-    return ksPvalD(F.length, Fprime.length, ksTest(F, Fprime));
-}
-
-unittest {
-, [1, 2, 3, 4, 5, 5, 5]
+.dup, [1, 2, 3, 4, 5, 5, 5].dup
-           .9375));
-    writeln("Passed ksPval 2-sample test.");
@@ -952 +1034 @@
-/**One-sided K-S test, returns P-val instead of D.
- * Bugs:  Exact calculation not implemented.  Uses asymptotic approximation.*/
-const T[]
-is(Func == function) || is(Func == delegate
+T
+realInput!(T) && (is(Func == function) || is(Func == delegate)
-    return ksPvalD(Femp.length, ksTest(Femp, F));
-}
@@ -959 +1041 @@
-unittest {
-    auto stdNormal = parametrize!(normalCDF)(0.0L, 1.0L);
-
+.dup
-
-    auto uniform01 = parametrize!(uniformCDF)(0, 1);
-
+.dup
-
-    writeln("Passed ksPval 1-sample test.");
@@ -1014 +1096 @@
-    alias contingencyTable c;
-
-nvariant
+mmutable
-                   n2 = c[1][0] + c[1][1],
-                   n  = c[0][0] + c[1][0];
-
-nvariant
+mmutable
-        cast(uint) ((cast(real) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
-nvariant
-nvariant
+mmutable
+mmutable
-
-    if(approxEqual(pExact, pMode, 1e-7)) {
-        return 1;
-    } else if(c[0][0] < mode) {
-nvariant
+mmutable
-
-        // Special case to prevent binary search from getting stuck.
@@ -1039 +1121 @@
-                    (cast(ulong) max + cast(ulong) min) / 2UL;
-
-nvariant
+mmutable
-            if(pGuess <= pExact &&
-                hypergeometricPMF(guess - 1, n1, n2, n) > pExact) {
@@ -1053 +1135 @@
-               hypergeometricCDFR(guess, n1, n2, n), 1.0L);
-    } else {
-nvariant
+mmutable
-
-        // Special case to prevent binary search from getting stuck.
@@ -1086 +1168 @@
-    // Simple, naive impl. of two-sided to test against.
-    static real naive(const uint[][] c) {
-nvariant
+mmutable
-                   n2 = c[1][0] + c[1][1],
-                   n  = c[0][0] + c[1][0];
-nvariant
+mmutable
-            cast(uint) ((cast(real) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
-nvariant
-nvariant
+mmutable
+mmutable
-        if(approxEqual(pExact, pMode, 1e-7))
-            return 1;
@@ -1107 +1189 @@
-
-    foreach(i; 0..1000) {
-gen, 
-gen, 
-gen, 
-gen, 
+
+
+
+
-        real naiveAns = naive(c);
-        real fastAns = fisherExact(c);
@@ -1173 +1255 @@
- * Bugs:  No exact calculation of the P-value.  Asymptotic approximation only.
- */
-
-
+
+
+
-    foreach(elem; obs) {
-addElemen
+pu
-    }
-    return r.pVal(alt);
@@ -1185 +1268 @@
-    // hard-coded as the equivalent to positive().  The median of this data
-    // is 0.5, so everything works.
-nvariant
+mmutable
-    assert(approxEqual(runsTest(data), 0.3581));
-    assert(approxEqual(runsTest(data, Alt.LESS), 0.821));
@@ -1206 +1289 @@
-
-    ///
-addElemen
+pu
-        bool curPos = pos(elem);
-        if(nRun == 0) {
@@ -1345 +1428 @@
- * Returns:   An array of Q-values with indices
- * corresponding to the indices of the p-values passed in.*/
-
+
+
-    // Not optimized at all because I can't imagine anyone writing code where
-    // FDR calculations are the main bottleneck.
-pVals.tempdup
+tempdup(pVals)
-    auto perm = newStack!(uint)(pVals.length);  scope(exit) TempAlloc.free;
-    foreach(i, ref elem; perm)
@@ -1375 +1459 @@
-unittest {
-    // Comparing results to R's qvalue package.
-
+.dup
-    auto qVals = falseDiscoveryRate(pVals);
-    assert(approxEqual(qVals[0], .9));