prep for v0.12.0

also: - bug fix in RVASInv - phyLiNC: no bang, un-export
JuliaPhylo · Jul 14, 2020 · 54a1d72 · 54a1d72 · cecileane · Jul 14, 2020
1 parent 3dc4e7f
commit 54a1d72
Show file tree

Hide file tree

Showing 7 changed files with 55 additions and 17 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "PhyloNetworks"
 uuid = "33ad39ac-ed31-50eb-9b15-43d0656eaa72"
 license = "MIT"
-version = "0.11.0"
+version = "0.12.0"
 
 [deps]
 BioSequences = "7e6ae17a-c86d-528c-b3b9-7f778a29fe59"

diff --git a/docs/src/lib/public.md b/docs/src/lib/public.md
@@ -139,5 +139,4 @@ nstates
 stationary
 empiricalDNAfrequencies
 mapindividuals
-phyLiNC!
 ```
diff --git a/src/PhyloNetworks.jl b/src/PhyloNetworks.jl
@@ -150,7 +150,7 @@ module PhyloNetworks
         empiricalDNAfrequencies,
         nni!,
         mapindividuals,
-        phyLiNC!,
+        # phyLiNC,
         # neighbor joining
         nj
 

diff --git a/src/phyLiNCoptimization.jl b/src/phyLiNCoptimization.jl
@@ -70,7 +70,7 @@ struct CacheLengthLiNC
 end
 
 """
-    phyLiNC!(net::HybridNetwork, fastafile::String, substitutionModel::Symbol)
+    phyLiNC(net::HybridNetwork, fastafile::String, substitutionModel::Symbol)
 
 Estimate a phylogenetic network from concatenated DNA data using
 maximum likelihood, ignoring incomplete lineage sorting
@@ -173,7 +173,7 @@ Regardless of these arguments, once a final topology is chosen, branch lenghts
 are optimized using stricter tolerances (1e-10, 1e-12, 1e-10, 1e-10) for better
 estimates.
 """
-function phyLiNC!(net::HybridNetwork, fastafile::String, modSymbol::Symbol,
+function phyLiNC(net::HybridNetwork, fastafile::String, modSymbol::Symbol,
                   rvsymbol=:noRV::Symbol, rateCategories=4::Int;
                   maxhybrid=1::Int, no3cycle=true::Bool,
                   nohybridladder=true::Bool,
@@ -184,6 +184,7 @@ function phyLiNC!(net::HybridNetwork, fastafile::String, modSymbol::Symbol,
     if !isempty(speciesfile)
         net, constraints = mapindividuals(net, speciesfile)
     else
+        net = deepcopy(net)
         constraints = TopologyConstraint[]
     end
     if !isempty(cladefile)
@@ -218,6 +219,15 @@ function phyLiNC!(net::HybridNetwork, fastafile::String, modSymbol::Symbol,
     phyLiNC!(obj; maxhybrid=maxhybrid, no3cycle=no3cycle, nohybridladder=nohybridladder,
             constraints=constraints, verbose=verbose, kwargs...)
 end
+
+"""
+    phyLiNC!(obj::SSM; kwargs...)
+
+Called by [`phyLiNC`](@ref) after `obj` is created (containing both the data
+and the model) and after checks are made to start from a network that satisfies
+all the constraints. Different runs are distributed to different processors,
+if more than one are available.
+"""
 function phyLiNC!(obj::SSM;
                   maxhybrid=1::Int, no3cycle=true::Bool,
                   nohybridladder=true::Bool, maxmoves=100::Int, nreject=75::Int,
@@ -431,7 +441,7 @@ If the number of processors is > 1, this will be false because workers can't
 write on streams opened by master. `logfile` will be stdout if `writelog_1proc`
 is false. Otherwise, it will be the log file created by `phyLiNC!`.
 
-See [`phyLiNC!`](@ref) for other arguments.
+See [`phyLiNC`](@ref) for other arguments.
 """
 function phyLiNCone!(obj::SSM, maxhybrid::Int, no3cycle::Bool,
                     nohybridladder::Bool, maxmoves::Int, nrejectmax::Int,
@@ -511,7 +521,7 @@ end
         constraints=TopologyConstraint[]::Vector{TopologyConstraint},
         verbose::Bool=false)
 
-Check that `net` is an adequate starting network before phyLiNC:
+Check that `net` is an adequate starting network before `phyLiNC!`:
 remove nodes of degree 2 (possibly including the root);
 check that `net` meets the topological `constraints`,
 has no polytomies (except at species constraints),
@@ -604,7 +614,7 @@ to be performed is in `PhyloNetworks.moveweights_LiNC`.
 - `nrejectmax`: the search stops when there has been this number of moves that
   have been rejected in a row (ignoring root changes)
 
-For a description of other arguments, see [`phyLiNC!`](@ref).
+For a description of other arguments, see [`phyLiNC`](@ref).
 
 Assumptions:
 - `checknetworkbeforeLiNC` and `discrete_corelikelihood!` have been called on
@@ -691,7 +701,7 @@ is accepted.
 Return true if move accepted, false if move rejected. Return nothing if there
 are no nni moves possible in the network.
 
-For arguments, see [`phyLiNC!`](@ref).
+For arguments, see [`phyLiNC`](@ref).
 
 Called by [`optimizestructure!`](@ref), which is called by [`phyLiNC!`](@ref).
 
@@ -760,7 +770,7 @@ optimizes branch lengths and gammas locally as part of PhyLiNC optimization.
 Return true if accepted add hybrid move. If move not accepted, return false.
 If cannot add a hybrid, return nothing.
 
-For arguments, see [`phyLiNC!`](@ref).
+For arguments, see [`phyLiNC`](@ref).
 Called by [`optimizestructure!`](@ref).
 """
 function addhybridedgeLiNC!(obj::SSM, currLik::Float64,
@@ -856,7 +866,7 @@ This creates a problem if the user asked for `nohybridladder`:
 this request may not be met.
 fixit: In future, we could check for this case and prevent it.
 
-For a description of arguments, see [`phyLiNC!`](@ref).
+For a description of arguments, see [`phyLiNC`](@ref).
 Called by [`optimizestructure!`](@ref), which does some checks.
 """
 function deletehybridedgeLiNC!(obj::SSM, currLik::Float64,

diff --git a/src/substitutionModels.jl b/src/substitutionModels.jl
@@ -1130,8 +1130,10 @@ struct RVASInv <: RateVariationAcrossSites
     lograteweight::StaticArrays.MVector{2,Float64}
 end
 function RVASInv(pinv=0.05::Float64)
+    r = StaticArrays.MVector{2,Float64}(undef) # rates
+    r[1] = 0.0 # invariable category
     obj = RVASInv(StaticArrays.MVector{1,Float64}(pinv),
-            StaticArrays.MVector{2,Float64}(undef), # rates
+            r,
             StaticArrays.MVector{2,Float64}(undef)) # log weights
     setpinv!(obj, pinv) # checks for 0 <= pinv < 1
     return obj
@@ -1148,10 +1150,12 @@ end
 function RVASGammaInv(pinv::Float64, alpha::Float64, ncat::Int)
     @assert ncat > 1 "ncat must be 2 or more for the Gamma+I model"
     s = 1+ncat
+    r = StaticArrays.MVector{s,Float64}(undef) # rates
+    r[1] = 0.0 # invariable category
     obj = RVASGammaInv{s}(
             StaticArrays.MVector{1,Float64}(pinv),
             StaticArrays.MVector{1,Float64}(alpha), ncat,
-            StaticArrays.MVector{s,Float64}(undef), # rates
+            r,
             StaticArrays.MVector{s,Float64}(undef)) # log weights
     setpinvalpha!(obj, pinv, alpha) # checks for α >= 0 and 0 <= pinv < 1
     return obj

diff --git a/src/traitsLikDiscrete.jl b/src/traitsLikDiscrete.jl
@@ -109,6 +109,30 @@ const SSM = StatisticalSubstitutionModel
 # fasta constructor: from net, fasta filename, modsymbol, and maxhybrid
 # Works for DNA in fasta format. Probably need different versions for
 # different kinds of data (snp, amino acids). Similar to fitdiscrete()
+"""
+    StatisticalSubstitutionModel(model::SubstitutionModel,
+            ratemodel::RateVariationAcrossSites,
+            net::HybridNetwork, trait::AbstractVector,
+            siteweight=nothing::Union{Nothing, Vector{Float64}},
+            maxhybrid=length(net.hybrid)::Int)
+
+Inner constructor. Makes a deep copy of the input model, rate model.
+Warning: does *not* make a deep copy of the network:
+modification of the `object.net` would modify the input `net`.
+Assumes that the network has valid gamma values (to extract displayed trees).
+
+    StatisticalSubstitutionModel(net::HybridNetwork, fastafile::String,
+            modsymbol::Symbol, rvsymbol=:noRV::Symbol,
+            ratecategories=4::Int;
+            maxhybrid=length(net.hybrid)::Int)
+
+Constructor from a network and a fasta file.
+The model symbol should be one of `:JC69`, `:HKY85`, `:ERSM` or `:BTSM`.
+The `rvsymbol` should be as required by [`RateVariationAcrossSites`](@ref).
+
+The network's gamma values are modified if they are missing. After that,
+a deep copy of the network is passed to the inner constructor.
+"""
 function StatisticalSubstitutionModel(net::HybridNetwork, fastafile::String,
         modsymbol::Symbol, rvsymbol=:noRV::Symbol, ratecategories=4::Int;
         maxhybrid=length(net.hybrid)::Int)
@@ -124,6 +148,7 @@ function StatisticalSubstitutionModel(net::HybridNetwork, fastafile::String,
     model = defaultsubstitutionmodel(net, modsymbol, data, siteweights)
     ratemodel = RateVariationAcrossSites(rvsymbol, ratecategories)
     dat2 = traitlabels2indices(view(data, :, 2:size(data,2)), model)
+    # check_matchtaxonnames makes a deep copy of the network
     o, net = check_matchtaxonnames!(data[:,1], dat2, net) # calls resetNodeNumbers, which calls preorder!
     trait = dat2[o]
     obj = StatisticalSubstitutionModel(model, ratemodel, net, trait, siteweights,

diff --git a/test/test_phyLiNCoptimization.jl b/test/test_phyLiNCoptimization.jl
@@ -239,12 +239,12 @@ end
 
 @testset "phyLiNC no constraints: HKY, rate variation" begin
 net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");
-obj = @test_nowarn PhyloNetworks.phyLiNC!(net, fastasimple, :JC69, :G, 2; maxhybrid=2, # no missing BLs, so they're not re-estimated
+obj = @test_nowarn PhyloNetworks.phyLiNC(net, fastasimple, :JC69, :G, 2; maxhybrid=2, # no missing BLs, so they're not re-estimated
                     no3cycle=true, nohybridladder=true, maxmoves=2,
                     nreject=1, nruns=1, filename="", verbose=false, seed=105)
 @test obj.loglik > -27.27
 net = readTopology("(((A:2.0,(B:1.0)#H1:0.1::0.9):1.5,(C:0.6,#H1:1.0::0.1):1.0):0.5,D:2.0);");
-obj = @test_nowarn PhyloNetworks.phyLiNC!(net, fastasimple, :HKY85; maxhybrid=2,
+obj = @test_nowarn PhyloNetworks.phyLiNC(net, fastasimple, :HKY85; maxhybrid=2,
                     no3cycle=true, nohybridladder=true, maxmoves=2, probST=1.0, # not enough moves to get back to a good topology
                     nreject=1, nruns=1, filename="phyLiNC2", verbose=false, seed=0)
 @test obj.loglik > -24.21
@@ -259,7 +259,7 @@ addprocs(1) # multiple cores
 #using Distributed; @everywhere begin; using Pkg; Pkg.activate("."); using PhyloNetworks; end
 originalstdout = stdout  # verbose=true below
 redirect_stdout(open("/dev/null", "w")) # not portable to Windows
-obj = PhyloNetworks.phyLiNC!(net, fastasimple, :JC69; maxhybrid=2, no3cycle=true,
+obj = PhyloNetworks.phyLiNC(net, fastasimple, :JC69; maxhybrid=2, no3cycle=true,
                         nohybridladder=true, maxmoves=2, nreject=1, nruns=2,
                         filename="phyLiNCmult", verbose=true, seed=106)
 redirect_stdout(originalstdout)
@@ -304,7 +304,7 @@ lcache = PhyloNetworks.CacheLengthLiNC(obj, 1e-2,1e-2,1e-2,1e-2, 5)
         1e-2, 1e-2, 0.0,50.0, 0.01,.9, γcache, lcache)
 @test obj.loglik > -65.0
 
-obj = phyLiNC!(net_level1_s, # missing BLs, so BLs are re-estimated before starting
+obj = PhyloNetworks.phyLiNC(net_level1_s, # missing BLs, so BLs are re-estimated before starting
             fastaindiv, :JC69, :Inv; maxhybrid=2, no3cycle=true, nohybridladder=true,
             verbose=false, filename="", speciesfile=mappingfile, seed=106, nruns=1,
             maxmoves=10, nreject=2)