Probabilities WIP

src/Information/Information.jl

@@ -5,9 +5,10 @@ module Information
 
   export BLOSUM62_Pa, BLOSUM62_Pab,
 
-  Pseudocount, AdditiveSmoothing, ResidueCount,
+  Pseudocount, AdditiveSmoothing, ResidueCount, ResidueProbability,
 	update!, apply_pseudocount!, count!
 
+
   #Fixed, Pseudofrequencies, ResidueProbabilities, ResiduePairProbabilities
 
   include("BLOSUM62.jl")

src/Information/Probabilities.jl

@@ -246,6 +246,100 @@ end
 
 count!{T, N, UseGap}(n::ResidueCount{T, N, UseGap}, res::AbstractVector{Residue}...) = count!(n, one(T), res...)
 
+## Probabilities
+
+type ResidueProbability{N, UseGap} <: AbstractArray{Float64, N}
+  probabilities::Array{Float64, N}
+  marginals::Array{Float64, 2}
+end
+
+call(::Type{ResidueProbability{1, true}}) = ResidueProbability{1, true}(Array(Float64, 21), Array(Float64, (21,1)))
+call(::Type{ResidueProbability{2, true}}) = ResidueProbability{2, true}(Array(Float64, (21, 21)), Array(Float64, (21,2)))
+
+call(::Type{ResidueProbability{1, false}}) = ResidueProbability{1, false}(Array(Float64, 20), Array(Float64, (20,1)))
+call(::Type{ResidueProbability{2, false}}) = ResidueProbability{2, false}(Array(Float64, (20, 20)), Array(Float64, (20,2)))
+
+function call{N, UseGap}(::Type{ResidueProbability{N, UseGap}})
+  nres = UseGap ? 21 : 20
+  ResidueProbability{N, UseGap}(Array(Float64, (Int[ nres for d in 1:N]...)), Array(Float64, (nres, N)))
+end
+
+zeros(::Type{ResidueProbability{1, true}}) = ResidueProbability{1, true}(zeros(Float64, 21), zeros(Float64, (21, 1)))
+zeros(::Type{ResidueProbability{2, true}}) = ResidueProbability{2, true}(zeros(Float64, (21, 21)), zeros(Float64, (21, 2)))
+
+zeros(::Type{ResidueProbability{1, false}}) = ResidueProbability{1, false}(zeros(Float64, 20), zeros(Float64, (20, 1)))
+zeros(::Type{ResidueProbability{2, false}}) = ResidueProbability{2, false}(zeros(Float64, (20, 20)), zeros(Float64, (20, 2)))
+
+function zeros{N, UseGap}(::Type{ResidueProbability{N, UseGap}})
+  nres = UseGap ? 21 : 20
+  ResidueProbability{N, UseGap}(zeros(Float64, (Int[ nres for d in 1:N]...)), zeros(Float64, (nres, N)))
+end
+
+### Abstract Array Interface
+
+Base.linearindexing(::ResidueProbability)	= Base.LinearFast()
+
+"""`getindex(p::ResidueProbability, i::Int)` gives you access to the probabilities, use `Int()` for indexing with `Residues`"""
+getindex(p::ResidueProbability, i::Int)	= getindex(p.probabilities, i)
+
+"""`setindex!(p::ResidueProbability, v, i::Int)` set a value into probabilities field, but doesn't update the marginals.
+Use `Int()` for indexing with `Residues`. Use `update!` in order to calculate again the marginals."""
+setindex!(p::ResidueProbability, v, i::Int) = setindex!(p.probabilities, v, i)
+
+#### Length & Size
+
+length(p::ResidueProbability{1, true}) = 21
+length(p::ResidueProbability{1,false}) = 20
+
+length(p::ResidueProbability{2, true}) = 441
+length(p::ResidueProbability{2,false}) = 400
+
+length{N,UseGap}(p::ResidueProbability{N, UseGap}) = length(p.probabilities)
+
+size(p::ResidueProbability{1, true}) = (21,)
+size(p::ResidueProbability{1,false}) = (20,)
+
+size(p::ResidueProbability{2, true}) = (21, 21)
+size(p::ResidueProbability{2,false}) = (20, 20)
+
+size{N,UseGap}(p::ResidueProbability{N, UseGap}) = size(p.probabilities)
+
+#### Iteration Interface
+
+start(p::ResidueProbability) = 1
+
+@inbounds next(p::ResidueProbability, state::Int) = (p.probabilities[state], state + 1)
+
+done{N,UseGap}(p::ResidueProbability{N, UseGap}, state) = state > length(p)
+
+eltype{N,UseGap}(::Type{ResidueProbability{N, UseGap}}) = Float64
+
+#### Similar
+
+similar{N,UseGap}(p::ResidueProbability{N,UseGap}) = ResidueProbability{N,UseGap}()
+similar{N,UseGap}(p::ResidueProbability{N,UseGap}, D::Int) = ResidueProbability{D,UseGap}()
+
+### Update!
+
+function update!{UseGap}(p::ResidueProbability{1,UseGap})
+	p.marginals[:] = p.probabilities
+	p
+end
+
+function update!{UseGap}(p::ResidueProbability{2,UseGap})
+	p.marginals[:,1] = sum(p.probabilities, 2) # this is faster than sum(p,2)
+	p.marginals[:,2] = sum(p.probabilities, 1)
+	p
+end
+
+function update!{N, UseGap}(p::ResidueProbability{N, UseGap})
+	for i in 1:N
+		p.marginals[:,i] = sum(p.probabilities, _tuple_without_index(i,N))
+	end
+	p
+end
+
+
 # """Fill the Pab matrix with λ, but doesn't update the marginal probabilities"""
 # function __initialize!(pab::ResiduePairProbabilities)
 #   fill!(pab.Pab, zero(Float64))

test/information.jl

@@ -3,6 +3,8 @@ using MIToS.Information
 using MIToS.Clustering
 using MIToS.MSA
 
+## Test Counts
+
 @test size(ResidueCount{Int, 1, false}().counts) == (20,)
 @test size(ResidueCount{Int, 1, true}().counts) == (21,)
 
@@ -106,6 +108,68 @@ seq = Residue(MSA._to_char)
 @test count!(zeros(ResidueCount{Float64, 2, false}), false_clusters, seq, seq).marginals[:,1] == getweight(false_clusters)
 @test count!(zeros(ResidueCount{Float64, 3, false}), false_clusters, seq, seq, seq).marginals[:,1] == getweight(false_clusters)
 
+## Test Probabilities
+
+
+@test size(ResidueProbability{1, false}().probabilities) == (20,)
+@test size(ResidueProbability{1, true}().probabilities) == (21,)
+
+@test size(ResidueProbability{2, false}().probabilities) == (20,20)
+@test size(ResidueProbability{2, true}().probabilities) == (21,21)
+
+@test size(ResidueProbability{3, false}().probabilities) == (20,20,20)
+@test size(ResidueProbability{3, true}().probabilities) == (21,21,21)
+
+for ndim = 1:4
+	@test size(ResidueProbability{ndim, true}().marginals) == (21, ndim)
+	@test size(ResidueProbability{ndim, false}().marginals) == (20, ndim)
+end
+
+for ndim = 1:4, usegap = Bool[true, false]
+	println("### N: $(ndim) UseGap: $(usegap) ###")
+  nres = usegap ? 21 : 20
+
+  println("Test zeros for ResidueProbability{$(ndim), $(usegap)}")
+  N = zeros(ResidueProbability{ndim, usegap})
+  @test sum(N.probabilities) == zero(Float64)
+
+  println("Test iteration interface for ResidueProbability{$(ndim), $(usegap)}")
+  @test collect(N) == zeros(Float64, nres^ndim)
+
+  println("Test size for ResidueProbability{$(ndim), $(usegap)}")
+  @test size(N) == size(N.probabilities)
+
+  println("Test indexing with i for ResidueProbability{$(ndim), $(usegap)}")
+  @test N[1] == zero(Float64)
+  @test N[Int(Residue('R'))] == zero(Float64)
+
+  println("Test setindex! with i for ResidueProbability{$(ndim), $(usegap)}")
+  N[3] = 1.0
+  @test N[3] == one(Float64)
+end
+
+for usegap = Bool[true, false]
+	println("### UseGap: $(usegap) ###")
+  nres = usegap ? 21 : 20
+
+	N = zeros(ResidueProbability{2, usegap})
+
+  println("Test indexing with ij for ResidueProbability{2, $(usegap)}")
+  @test N[1,1] == zero(Float64)
+  @test N[Int(Residue('R')), Int(Residue('R'))] == zero(Float64)
+
+  println("Test setindex! with ij for ResidueProbability{2, $(usegap)}")
+  N[3,3] = 1.0
+  @test N[3,3] == one(Float64)
+end
+
+for ndim = 1:4
+	@test size(similar(ResidueProbability{ndim, false}()).marginals) == (20, ndim)
+	@test size(similar(ResidueProbability{ndim, true}()).marginals) == (21, ndim)
+end
+
+@test size(similar(ResidueProbability{1, false}(),4).marginals) == (20, 4)
+
 
 # # Copy & deepcopy
 # cnone = copy(none)