Merge pull request #18 from BioJulia/constructors

Expansion of constructors methods for Nucleic Acid Substitution Models

src/SubstitutionModels.jl

@@ -5,36 +5,15 @@ module SubstitutionModels
     StaticArrays,
     LinearAlgebra
 
-  import
-    Base.getindex,
-    Base.show,
-    Base.setindex!,
-    Base.checkbounds
-
   include("core.jl")
   include("indexing.jl")
   include("nucleic_acid/nucleic_acid.jl")
-  include("nucleic_acid/jc69/abstract.jl")
-  include("nucleic_acid/jc69/absolute.jl")
-  include("nucleic_acid/jc69/relative.jl")
-  include("nucleic_acid/k80/abstract.jl")
-  include("nucleic_acid/k80/absolute.jl")
-  include("nucleic_acid/k80/relative.jl")
-  include("nucleic_acid/f81/abstract.jl")
-  include("nucleic_acid/f81/absolute.jl")
-  include("nucleic_acid/f81/relative.jl")
-  include("nucleic_acid/f84/abstract.jl")
-  include("nucleic_acid/f84/absolute.jl")
-  include("nucleic_acid/f84/relative.jl")
-  include("nucleic_acid/hky85/abstract.jl")
-  include("nucleic_acid/hky85/absolute.jl")
-  include("nucleic_acid/hky85/relative.jl")
-  include("nucleic_acid/tn93/abstract.jl")
-  include("nucleic_acid/tn93/absolute.jl")
-  include("nucleic_acid/tn93/relative.jl")
-  include("nucleic_acid/gtr/abstract.jl")
-  include("nucleic_acid/gtr/absolute.jl")
-  include("nucleic_acid/gtr/relative.jl")
+
+  for m in ["jc69"; "k80"; "f81"; "f84"; "hky85"; "tn93"; "gtr"]
+    for f in ["abstract"; "absolute"; "relative"; "outer_constructors"]
+      include("nucleic_acid/$(m)/$(f).jl")
+    end
+  end
 
   export
     SubstitutionModel, SM,

src/core.jl

@@ -19,3 +19,18 @@ const Qmatrix = SMatrix{4, 4, Float64}
 
 
 const Pmatrix = SMatrix{4, 4, Float64}
+
+
+function Base.convert(::Type{T}, θ::F...; safe::Bool=true) where {T <: NASM, F <: Float64}
+  return T(θ..., safe)
+end
+
+
+function Base.convert(::Type{T}, θ_vec::A; safe::Bool=true) where {T <: NASM, A <: AbstractArray}
+  return T(θ_vec, safe)
+end
+
+
+function Base.convert(::Type{T}, θ_vec::A, π_vec::A; safe::Bool=true) where {T <: NASM, A <: AbstractArray}
+  return T(θ_vec, π_vec, safe)
+end

src/indexing.jl

@@ -1,38 +1,38 @@
-@inline function checkbounds(a::AbstractArray, i::T, j::T) where T <: NucleicAcid
+@inline function Base.checkbounds(a::AbstractArray, i::T, j::T) where T <: NucleicAcid
   if iscertain(i) & iscertain(j)
     return true
   end
   checkbounds(a, trailing_zeros(i) + 1, trailing_zeros(j) + 1)
 end
 
 
-@inline function checkbounds(a::AbstractArray, i::NucleicAcid)
+@inline function Base.checkbounds(a::AbstractArray, i::NucleicAcid)
   if iscertain(i)
       return true
   end
   checkbounds(a, trailing_zeros(i) + 1)
 end
 
 
-@inline function getindex(a::AbstractArray, i::T, j::T) where T <: NucleicAcid
+@inline function Base.getindex(a::AbstractArray, i::T, j::T) where T <: NucleicAcid
   @boundscheck checkbounds(a, i, j)
   @inbounds return a[trailing_zeros(i) + 1, trailing_zeros(j) + 1]
 end
 
 
-@inline function getindex(a::AbstractArray, i::NucleicAcid)
+@inline function Base.getindex(a::AbstractArray, i::NucleicAcid)
   @boundscheck checkbounds(a, i)
   @inbounds return a[trailing_zeros(i) + 1]
 end
 
 
-@inline function setindex!(a::AbstractArray, x, i::T, j::T) where T <: NucleicAcid
+@inline function Base.setindex!(a::AbstractArray, x, i::T, j::T) where T <: NucleicAcid
   @boundscheck checkbounds(a, i, j)
   @inbounds return setindex!(a, x, trailing_zeros(i) + 1, trailing_zeros(j) + 1)
 end
 
 
-@inline function setindex!(a::AbstractArray, x, i::NucleicAcid)
+@inline function Base.setindex!(a::AbstractArray, x, i::NucleicAcid)
   @boundscheck checkbounds(a, i)
   @inbounds return setindex!(a, x, trailing_zeros(i) + 1)
 end

src/nucleic_acid/f81/absolute.jl

@@ -5,28 +5,28 @@ struct F81abs <: F81
   πG::Float64
   πT::Float64
   function F81abs(β::Float64,
-                  πA::Float64, πC::Float64, πG::Float64, πT::Float64)
-    if β <= 0.
-      error("F81 parameter β must be positive")
-    elseif sum([πA,πC,πG,πT]) != 1.0
-      error("F81 frequencies must sum to 1.0")
-    elseif any([πA,πC,πG,πT] .<=0.0)
-      error("F81 frequencies must be positive")
+                  πA::Float64, πC::Float64, πG::Float64, πT::Float64,
+                  safe::Bool=true)
+    if safe
+      if β <= 0.
+        error("F81 parameter β must be positive")
+      elseif sum([πA,πC,πG,πT]) != 1.0
+        error("F81 frequencies must sum to 1.0")
+      elseif any([πA,πC,πG,πT] .<=0.0)
+        error("F81 frequencies must be positive")
+      end
     end
     new(β, πA, πC, πG, πT)
   end
 end
 
 
-function show(io::IO, object::F81abs)
+function Base.show(io::IO, object::F81abs)
   print(io, "\r\e[0m\e[1mF\e[0melsenstein 19\e[1m81\e[0m model (absolute rate form)
 β = $(object.β), π = [$(object.πA), $(object.πC), $(object.πG), $(object.πT)]")
 end
 
 
-F81(β, πA, πC, πG, πT) = F81abs(β, πA, πC, πG, πT)
-
-
 @inline function Q(mod::F81abs)
   β = mod.β
   πA = _πA(mod); πC = _πC(mod); πG = _πG(mod); πT = _πT(mod)

src/nucleic_acid/f81/outer_constructors.jl

@@ -0,0 +1,57 @@
+F81(πA::F, πC::F, πG::F, πT::F,
+    safe::Bool=true) where F <: Float64 =
+F81rel(πA, πC, πG, πT)
+
+
+F81(β::F,
+    πA::F, πC::F, πG::F, πT::F,
+    safe::Bool=true) where F <: Float64 =
+F81abs(β, πA, πC, πG, πT, safe)
+
+
+function F81(θ_vec::A,
+             π_vec::A,
+             safe::Bool=true) where A <: AbstractArray
+  if safe && length(π_vec) != 4
+    error("Incorrect base frequency vector length")
+  end
+  if length(θ_vec) == 0
+    return F81rel(π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
+                  safe)
+  elseif length(θ_vec) == 1
+    return F81abs(θ_vec[1],
+                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
+                  safe)
+  else
+    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
+  end
+end
+
+
+function F81rel(θ_vec::A,
+                π_vec::A,
+                safe::Bool=true) where A <: AbstractArray
+  if safe
+    if length(θ_vec) != 0
+      error("Incorrect parameter vector length")
+    elseif length(π_vec) != 4
+      error("Incorrect base frequency vector length")
+    end
+  end
+  return F81rel(π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
+end
+
+
+
+function F81abs(θ_vec::A,
+                π_vec::A,
+                safe::Bool=true) where A <: AbstractArray
+  if safe
+    if length(θ_vec) != 1
+      error("Incorrect parameter vector length")
+    elseif length(π_vec) != 4
+      error("Incorrect base frequency vector length")
+    end
+  end
+  return F81abs(θ_vec[1], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
+end

src/nucleic_acid/f81/relative.jl

@@ -3,26 +3,26 @@ struct F81rel <: F81
   πC::Float64
   πG::Float64
   πT::Float64
-  function F81rel(πA::Float64, πC::Float64, πG::Float64, πT::Float64)
-    if sum([πA,πC,πG,πT]) != 1.0
-      error("F81 frequencies must sum to 1.0")
-    elseif any([πA,πC,πG,πT] .<=0.0)
-      error("F81 frequencies must be positive")
+  function F81rel(πA::Float64, πC::Float64, πG::Float64, πT::Float64,
+                  safe::Bool=true)
+    if safe
+      if sum([πA,πC,πG,πT]) != 1.0
+        error("F81 frequencies must sum to 1.0")
+      elseif any([πA,πC,πG,πT] .<=0.0)
+        error("F81 frequencies must be positive")
+      end
     end
     new(πA, πC, πG, πT)
   end
 end
 
 
-function show(io::IO, object::F81rel)
+function Base.show(io::IO, object::F81rel)
   print(io, "\r\e[0m\e[1mF\e[0melsenstein 19\e[1m81\e[0m model (relative rate form)
 π = [$(object.πA), $(object.πC), $(object.πG), $(object.πT)]")
 end
 
 
-F81(πA, πC, πG, πT) = F81rel(πA, πC, πG, πT)
-
-
 @inline function Q(mod::F81rel)
   πA = _πA(mod); πC = _πC(mod); πG = _πG(mod); πT = _πT(mod)
 

src/nucleic_acid/f84/absolute.jl

@@ -6,30 +6,30 @@ struct F84abs <: F84
   πG::Float64
   πT::Float64
   function F84abs(κ::Float64, β::Float64,
-                  πA::Float64, πC::Float64, πG::Float64, πT::Float64)
-    if κ <= 0.
-      error("F84 parameter α must be positive")
-    elseif β <= 0.
-      error("F84 parameter β must be positive")
-    elseif sum([πA,πC,πG,πT]) != 1.0
-      error("F84 frequencies must sum to 1.0")
-    elseif any([πA,πC,πG,πT] .<= 0.0)
-      error("F84 frequencies must be positive")
+                  πA::Float64, πC::Float64, πG::Float64, πT::Float64,
+                  safe::Bool=true)
+    if safe
+      if κ <= 0.
+        error("F84 parameter α must be positive")
+      elseif β <= 0.
+        error("F84 parameter β must be positive")
+      elseif sum([πA,πC,πG,πT]) != 1.0
+        error("F84 frequencies must sum to 1.0")
+      elseif any([πA,πC,πG,πT] .<= 0.0)
+        error("F84 frequencies must be positive")
+      end
     end
     new(κ, β, πA, πC, πG, πT)
   end
 end
 
 
-function show(io::IO, object::F84abs)
+function Base.show(io::IO, object::F84abs)
   print(io, "\r\e[0m\e[1mF\e[0melsenstein 19\e[1m84\e[0m substitution model (absolute rate form)
 κ = $(object.κ), β = $(object.β), π = [$(object.πA), $(object.πC), $(object.πG), $(object.πT)]")
 end
 
 
-F84(κ, β, πA, πC, πG, πT) = F84abs(κ, β, πA, πC, πG, πT)
-
-
 @inline function Q(mod::F84abs)
   β = mod.β; κ = mod.κ
   πA = _πA(mod); πC = _πC(mod); πG = _πG(mod); πT = _πT(mod)

src/nucleic_acid/f84/outer_constructors.jl

@@ -0,0 +1,58 @@
+F84(κ::F,
+    πA::F, πC::F, πG::F, πT::F,
+    safe::Bool=true) where F <: Float64 =
+  F84rel(κ, πA, πC, πG, πT, safe)
+
+
+F84(κ::F, β::F,
+    πA::F, πC::F, πG::F, πT::F,
+    safe::Bool=true) where F <: Float64 =
+  F84abs(κ, β, πA, πC, πG, πT, safe)
+
+
+function F84(θ_vec::A,
+             π_vec::A,
+             safe::Bool=true) where A <: AbstractArray
+  if safe && length(π_vec) != 4
+    error("Incorrect base frequency vector length")
+  end
+  if length(θ_vec) == 1
+    return F84rel(θ_vec[1],
+                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
+                  safe)
+  elseif length(θ_vec) == 2
+    return F84abs(θ_vec[1], θ_vec[2],
+                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
+                  safe)
+  else
+    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
+  end
+end
+
+
+function F84rel(θ_vec::A,
+                π_vec::A,
+                safe::Bool=true) where A <: AbstractArray
+  if safe
+    if length(θ_vec) != 1
+      error("Incorrect parameter vector length")
+    elseif length(π_vec) != 4
+      error("Incorrect base frequency vector length")
+    end
+  end
+  return F84rel(θ_vec[1], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
+end
+
+
+function F84abs(θ_vec::A,
+                π_vec::A,
+                safe::Bool=true) where A <: AbstractArray
+  if safe
+    if length(θ_vec) != 2
+      error("Incorrect parameter vector length")
+    elseif length(π_vec) != 4
+      error("Incorrect base frequency vector length")
+    end
+  end
+  return F84abs(θ_vec[1], θ_vec[2], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
+end

src/nucleic_acid/f84/relative.jl

@@ -5,28 +5,28 @@ struct F84rel <: F84
   πG::Float64
   πT::Float64
   function F84rel(κ::Float64,
-                  πA::Float64, πC::Float64, πG::Float64, πT::Float64)
-    if κ <= 0.
-      error("F84 parameter κ must be positive")
-    elseif sum([πA,πC,πG,πT]) != 1.0
-      error("F84 frequencies must sum to 1.0")
-    elseif any([πA,πC,πG,πT] .<=0.0)
-      error("F84 frequencies must be positive")
+                  πA::Float64, πC::Float64, πG::Float64, πT::Float64,
+                  safe::Bool=true)
+    if safe
+      if κ <= 0.
+        error("F84 parameter κ must be positive")
+      elseif sum([πA,πC,πG,πT]) != 1.0
+        error("F84 frequencies must sum to 1.0")
+      elseif any([πA,πC,πG,πT] .<=0.0)
+        error("F84 frequencies must be positive")
+      end
     end
     new(κ, πA, πC, πG, πT)
   end
 end
 
 
-function show(io::IO, object::F84rel)
+function Base.show(io::IO, object::F84rel)
   print(io, "\r\e[0m\e[1mF\e[0melsenstein 19\e[1m84\e[0m substitution model (relative rate form)
 κ = $(object.κ), π = [$(object.πA), $(object.πC), $(object.πG), $(object.πT)]")
 end
 
 
-F84(κ, πA, πC, πG, πT) = F84rel(κ, πA, πC, πG, πT)
-
-
 @inline function Q(mod::F84rel)
   κ = mod.κ
   πA = _πA(mod); πC = _πC(mod); πG = _πG(mod); πT = _πT(mod)