startingBL: no unzipping option

src/phyLiNCoptimization.jl

@@ -156,10 +156,9 @@ function phyLiNC!(net::HybridNetwork, fastafile::String, substitutionModel::Symb
        Warning: need to call updateSSM after using checknetwork_LiNC =#
     updateSSM!(obj, true; constraints=constraints)
     if any([e.length < 0.0 for e in obj.net.edge]) # check for missing branch lengths
-        startingBL!(obj.net, true, obj.trait, obj.siteweight) # true: to unzip
-    else # keep existing branch lengths in input network, only unzip reticulations
-        unzip_canonical!(obj.net)
+        startingBL!(obj.net, obj.trait, obj.siteweight)
     end
+    unzip_canonical!(obj.net)
     phyLiNC!(obj; maxhybrid=maxhybrid, no3cycle=no3cycle, nohybridladder=nohybridladder,
             constraints=constraints, verbose=verbose, kwargs...)
 end
@@ -932,15 +931,15 @@ end
 
 ## Optimize Branch Lengths and Gammas ##
 """
-    startingBL!(net::HybridNetwork, unzip::Bool,
+    startingBL!(net::HybridNetwork,
                 trait::AbstractVector{Vector{Union{Missings.Missing,Int}}},
                 siteweight=ones(length(trait[1]))::AbstractVector{Float64})
 
 Calibrate branch lengths in `net` by minimizing the mean squared error
 between the JC-adjusted pairwise distance between taxa, and network-predicted
 pairwise distances, using [`calibrateFromPairwiseDistances!`](@ref).
 `siteweight[k]` gives the weight of site (or site pattern) `k` (default: all 1s).
-`unzip` = true sets all edges below a hybrid node to length zero.
+Note: the network is not "unzipped", as required by PhyLiNC.
 
 Assumptions:
 
@@ -953,7 +952,7 @@ Assumptions:
   If not, all pairwise hamming distances are scaled by `.75/(m*1.01)` where `m`
   is the maximum observed hamming distance, to make them all < 0.75.
 """
-function startingBL!(net::HybridNetwork, unzip::Bool,
+function startingBL!(net::HybridNetwork,
         trait::AbstractVector{Vector{Union{Missings.Missing,Int}}},
         siteweight=ones(length(trait[1]))::AbstractVector{Float64})
     nspecies = net.numTaxa
@@ -997,9 +996,6 @@ function startingBL!(net::HybridNetwork, unzip::Bool,
             e.length = 0.0001
         end
     end
-    if unzip
-        unzip_canonical!(net)
-    end
     return net
 end
 

test/test_phyLiNCoptimization.jl

@@ -49,7 +49,8 @@ PhyloNetworks.checknetwork_LiNC!(obj.net, 3, true, true, emptyconstraint)
 # checknetwork removes degree-2 nodes (including root) and 2- and 3-cycles
 # and requires that the network is preordered.
 PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
-PhyloNetworks.startingBL!(obj.net, true, obj.trait, obj.siteweight)
+PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
+PhyloNetworks.unzip_canonical!(obj.net)
 ## Local BL
 lengthe = obj.net.edge[27].length
 @test_nowarn PhyloNetworks.optimizelocalBL_LiNC!(obj, obj.net.edge[27], 1e-6,1e-6,1e-2,1e-3)
@@ -127,7 +128,8 @@ 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):
 obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69, 1)
 PhyloNetworks.checknetwork_LiNC!(obj.net, 1, true, true)
 PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
-PhyloNetworks.startingBL!(obj.net, true, obj.trait, obj.siteweight)
+PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
+PhyloNetworks.unzip_canonical!(obj.net)
 PhyloNetworks.discrete_corelikelihood!(obj)
 @test obj.loglik ≈ -29.7762035
 maxmoves = 2
@@ -153,7 +155,8 @@ for nohybridladder in [true, false]
     obj = PhyloNetworks.StatisticalSubstitutionModel(net, fastasimple, :JC69, 1)
     PhyloNetworks.checknetwork_LiNC!(obj.net, 1, no3cycle, nohybridladder)
     PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
-    PhyloNetworks.startingBL!(obj.net, true, obj.trait, obj.siteweight)
+    PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
+    PhyloNetworks.unzip_canonical!(obj.net)
     obj.loglik = -Inf # missing otherwise, which would cause an error below
     nullio = open("/dev/null", "w")
     γcache = PhyloNetworks.CacheGammaLiNC(obj)
@@ -218,7 +221,8 @@ c_species[1] = PhyloNetworks.TopologyConstraint(0x01, c_species[1].taxonnames, o
 PhyloNetworks.updateSSM!(obj, true; constraints=emptyconstraint)
 
 for e in obj.net.edge e.length = 0.1; end # was -1.0 for missing
-PhyloNetworks.startingBL!(obj.net, true, obj.trait, obj.siteweight) # true: to unzip
+PhyloNetworks.startingBL!(obj.net, obj.trait, obj.siteweight)
+PhyloNetworks.unzip_canonical!(obj.net)
 obj.loglik = -Inf # actual likelihood -56.3068141288164. Need something non-missing
 seed = 103
 nullio = open("/dev/null", "w")

test/test_traitLikDiscrete.jl

@@ -617,18 +617,18 @@ dna_net = readTopology("((((((((((((((Ae_caudata_Tr275,Ae_caudata_Tr276),Ae_caud
 dat2 = PhyloNetworks.traitlabels2indices(dna_dat[!,2:end], JC69([0.5]))
 o, dna_net = @test_logs (:warn, "the network contains taxa with no data: those will be pruned") match_mode=:any PhyloNetworks.check_matchtaxonnames!(copy(dna_dat[1]), dat2, dna_net)
 trait = view(dat2, o)
-PhyloNetworks.startingBL!(dna_net, false, trait, dna_weights)
+PhyloNetworks.startingBL!(dna_net, trait, dna_weights)
 @test maximum(e.length for e in dna_net.edge) > 0.03
-@test_logs PhyloNetworks.startingBL!(dna_net, false, trait) # no dna_weights
+@test_logs PhyloNetworks.startingBL!(dna_net, trait) # no dna_weights
 
 dna_dat, dna_weights = readfastatodna(fastafile, true);
 dna_net = readTopology("((((((((((((((Ae_caudata_Tr275,Ae_caudata_Tr276),Ae_caudata_Tr139))#H1,#H2),(((Ae_umbellulata_Tr266,Ae_umbellulata_Tr257),Ae_umbellulata_Tr268),#H1)),((Ae_comosa_Tr271,Ae_comosa_Tr272),(((Ae_uniaristata_Tr403,Ae_uniaristata_Tr357),Ae_uniaristata_Tr402),Ae_uniaristata_Tr404))),(((Ae_tauschii_Tr352,Ae_tauschii_Tr351),(Ae_tauschii_Tr180,Ae_tauschii_Tr125)),(((((((Ae_longissima_Tr241,Ae_longissima_Tr242),Ae_longissima_Tr355),(Ae_sharonensis_Tr265,Ae_sharonensis_Tr264)),((Ae_bicornis_Tr408,Ae_bicornis_Tr407),Ae_bicornis_Tr406)),((Ae_searsii_Tr164,Ae_searsii_Tr165),Ae_searsii_Tr161)))#H2,#H4))),(((T_boeoticum_TS8,(T_boeoticum_TS10,T_boeoticum_TS3)),T_boeoticum_TS4),((T_urartu_Tr315,T_urartu_Tr232),(T_urartu_Tr317,T_urartu_Tr309)))),(((((Ae_speltoides_Tr320,Ae_speltoides_Tr323),Ae_speltoides_Tr223),Ae_speltoides_Tr251))H3,((((Ae_mutica_Tr237,Ae_mutica_Tr329),Ae_mutica_Tr244),Ae_mutica_Tr332))#H4))),Ta_caputMedusae_TB2),S_vavilovii_Tr279),Er_bonaepartis_TB1),H_vulgare_HVens23);");
 dat2 = PhyloNetworks.traitlabels2indices(dna_dat[!,2:end], HKY85([0.5], [0.25, 0.25, 0.25, 0.25], true))
 o, dna_net = @test_logs (:warn, "the network contains taxa with no data: those will be pruned") match_mode=:any PhyloNetworks.check_matchtaxonnames!(copy(dna_dat[1]), dat2, dna_net)
 trait = view(dat2, o)
-PhyloNetworks.startingBL!(dna_net, false, trait, dna_weights)
+PhyloNetworks.startingBL!(dna_net, trait, dna_weights)
 @test maximum(e.length for e in dna_net.edge) > 0.03
-@test_logs PhyloNetworks.startingBL!(dna_net, false, trait) # no dna_weights
+@test_logs PhyloNetworks.startingBL!(dna_net, trait) # no dna_weights
 end # of startingBL!
 
 @testset "testing prep and wrapper functions" begin