Merge pull request #3 from JuliaMusic/weights

allow weights in sampling

CHANGELOG.md

@@ -7,5 +7,8 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 # master
 work in progress changes are contained in this section.
 
+# v0.2.0
+Keyword argument `weights` allows you to choose weights for the motifs that are used in the initial sequence.
+
 # v0.1.0 - Initial Release
 Changelog is kept with respect to this release.

REQUIRE

@@ -1,2 +1,3 @@
 julia 0.7.0-beta2
 Combinatorics 0.6.0
+StatsBase 0.23.0

src/MotifSequenceGenerator.jl

@@ -6,7 +6,7 @@ All main functionality is given by the function [`random_sequence`](@ref).
 """
 module MotifSequenceGenerator
 
-using Combinatorics, Random
+using Combinatorics, Random, StatsBase
 
 export random_sequence, all_possible_sums
 
@@ -23,11 +23,12 @@ Base.showerror(io::IO, e::DeadEndMotifs) = print(io,
 The algorithm works as follows: First a random sequence of motifs is created,
 so that it has length of `q - δq ≤ ℓ ≤ q - δq + maximum(motiflengths)`.
 The possible tries of random sequences is set by the `tries` keyword (default `5`).
+The sequence is optionally sampled given a probability vector.
 
-For each random try, it is first check whether the sequence is already correct.
+For each random try, it is first checked whether the sequence is already correct.
 If not, the last entry of the sequence is dropped. Then, since the sequence is now
 already smaller than `q`, all possible sums of `summands` out of the motif pool
-are checked. If some combination of `summands` sums to exactly the difference,
+are checked. If some combination of `summands` sums to the difference,
 they are added to the sequence.
 For multiple satisfactory combinations, a random one is picked.
 
@@ -50,7 +51,9 @@ a proper one, an error is thrown.
 Create a random sequence of motifs of type `M`, under the constraint that the
 sequence has "length" `ℓ` **exactly** within `q - δq ≤ ℓ ≤ q + δq`.
 Return the sequence itself as well as the
-sequence of indices of `motifs` used to create it.
+sequence of indices of `motifs` used to create it. A vector of probabilities `weights`
+can be given as a keyword argument, which then dictates the sampling probability
+for each entry of `motifs` for the initial sequence created.
 
 "length" here means an abstracted length defined by the struct `M`,
 based on the `limits` and `translate` functions.
@@ -64,7 +67,7 @@ It does **not** refer to the amount of elements!
   motif which is translated by `t` (either negative or positive), with
   respect to the same units as `q`.
 
-## Keywords
+## Other Keywords
 Please see the source code (use `@which`) for a full description of the algorithm.
 
 * `tries = 5` : Up to how many initial random sequences are accepted.
@@ -74,7 +77,10 @@ Please see the source code (use `@which`) for a full description of the algorith
 """
 function random_sequence(motifs::Vector{M}, q,
     limits, translate, δq = 0;
-    tries = 5, summands = 3, tailcut = 2) where {M}
+    tries = 5, summands = 3, tailcut = 2,
+    weights = ones(length(motifs))) where {M}
+
+    ws = _toweight(weights)
 
     idxs = 1:length(motifs)
     motifs0, motiflens = _motifs_at_origin(motifs, limits, translate)
@@ -94,7 +100,7 @@ function random_sequence(motifs::Vector{M}, q,
     while worked == false
         count > tries && throw(DeadEndMotifs(tries, summands, tailcut))
 
-        seq, seq_length = _random_sequence_try(motiflens, q, δq)
+        seq, seq_length = _random_sequence_try(motiflens, q, δq, ws)
 
         worked = _complete_sequence!(seq, motiflens, q, δq, summands, tailcut)
 
@@ -104,6 +110,8 @@ function random_sequence(motifs::Vector{M}, q,
     return _instantiate_sequence(motifs0, motiflens, seq, translate), seq
 end
 
+_toweight(a) = (s = sum(a); ProbabilityWeights(a./s, 1))
+
 """
     _motifs_at_origin(motifs, limits, translate) -> (motifs0, motiflens)
 Bring all motifs to the origin and compute the motif lengths.
@@ -121,15 +129,15 @@ function _motifs_at_origin(motifs::Vector{M}, limits, translate) where M
 end
 
 """
-    _random_sequence_try(motiflens, q) -> seq, seq_length
+    _random_sequence_try(motiflens, q, δq [, ws]) -> seq, seq_length
 Return a random sequence of motif indices
 so that the total sequence is *guaranteed* to have total length of
 `q - δq ≤ ℓ ≤ q - δq + maximum(motiflens)`.
 """
-function _random_sequence_try(motiflens, q, δq)
+function _random_sequence_try(motiflens, q, δq, ws = defaultweights(motiflens))
     seq = Int[]; seq_length = 0; idxs = 1:length(motiflens)
     while seq_length < q - δq
-        i = rand(idxs)
+        i = sample(idxs, ws)
         push!(seq, i)
         seq_length += motiflens[i]
     end
@@ -174,13 +182,6 @@ function _complete_sequence_remainder!(seq, motiflens, q, δq, summands, tailcut
         pop!(seq)
         isempty(seq) && return false
 
-        # I Think the following if is unecessary?...
-        # if q - δq - sum(motiflens[k] for k in seq) < 0
-        #     ok = _complete_sequence_extra!(seq, motiflens, q, δq)
-        #     isempty(seq) && return false
-        #     ok && return true
-        # end
-
         # At this point ℓ is guaranteed less than q - δq
         remainder = q - δq - sum(motiflens[k] for k in seq)
         @assert remainder > 0

test/float_length_tests.jl

@@ -35,6 +35,16 @@ end
     end
 end
 
+@testset "Float Length, Weights, δq=$(δq)" for δq in [1.0, 2.0]
+    weights = rand(1:5, N)
+    for j in 1:N
+        r, s = random_sequence(shouts, q, shoutlimits, shouttranslate, δq;
+        weights = weights)
+        ℓ = shoutlens(r)
+        @test q - δq ≤ ℓ ≤ q + δq
+    end
+end
+
 using MotifSequenceGenerator: DeadEndMotifs
 @test_throws ArgumentError random_sequence(shouts, q, shoutlimits, shouttranslate, 0.0)
 @test_throws DeadEndMotifs random_sequence(shouts, q, shoutlimits, shouttranslate, 0.000001)

test/integer_length_tests.jl

@@ -34,6 +34,16 @@ end
     end
 end
 
+@testset "Integer Length, Weights, δq=$(δq)" for δq in [0, 2]
+    weights = rand(1:5, N)
+    for j in 1:N
+        r, s = random_sequence(shouts, q, shoutlimits, shouttranslate, δq;
+        tries = 10, weights = weights)
+        ℓ = shoutlens(r)
+        @test q - δq ≤ ℓ ≤ q + δq
+    end
+end
+
 using MotifSequenceGenerator: DeadEndMotifs
 @test_throws DeadEndMotifs random_sequence(shouts, 7, shoutlimits, shouttranslate, 0)