Merge pull request YangLab-um#39 from biphy/bugfix/dependency

Fix bugs related to newer versions of dependencies

Project.toml

@@ -4,48 +4,30 @@ authors = ["Franco Tavella"]
 version = "0.0.2"
 
 [deps]
-CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
-JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
-LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LatinHypercubeSampling = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Luxor = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-CSV = "0.10"
 Catalyst = "13"
 DSP = "0.7"
 DataFrames = "1"
 DifferentialEquations = "7"
-Documenter = "0.27"
-DocumenterTools = "0.1"
-JLD2 = "0.4"
-LaTeXStrings = "1"
 Latexify = "0.15"
 LatinHypercubeSampling = "1"
 ModelingToolkit = "8"
 Peaks = "0.4"
-Plots = "1"
-PyPlot = "2"
-Revise = "3"
 julia = "1.8"
 
 [extras]

src/feature_calculation.jl

@@ -74,11 +74,20 @@ function calculate_amplitude(ode_solution::ODESolution)
         # Peaks
         peak_positions, ~ = findmaxima(ode_solution[i,:])
         peak_positions, ~ = peakproms(peak_positions, ode_solution[i,:], minprom=rough_amp/3.0)
-        peak_values = ode_solution[i, peak_positions]
+        if isempty(peak_positions)
+            # This case must be explicitly handled as ode_solution[i, Int64[]] raises an ArgumentError
+            peak_values = Float64[]
+        else
+            peak_values = ode_solution[i, peak_positions]
+        end
         # Troughs
         trough_positions, ~ = findminima(ode_solution[i,:])
         trough_positions, ~ = peakproms(trough_positions, ode_solution[i,:], minprom=rough_amp/3.0)
-        trough_values = ode_solution[i, trough_positions]
+        if isempty(trough_positions)
+            trough_values = Float64[]
+        else
+            trough_values = ode_solution[i, trough_positions]
+        end
         if length(peak_values) > 0 && length(trough_values) > 0
             amplitude[i] = mean(peak_values) - mean(trough_values)
             peak_variation[i] = std(peak_values)/amplitude[i]

src/gene_regulatory_network.jl

@@ -17,8 +17,8 @@ function gene_regulatory_network(connectivity::AbstractMatrix)
     # Number of edges
     E = count(!iszero, connectivity)
     @variables t
-    @species (m(t))[collect(1:N)]
-    @species (X(t))[collect(1:N)]
+    @species (m(t))[1:N]
+    @species (X(t))[1:N]
     @parameters α[1:N], β[1:N], δ[1:N], γ[1:E], κ[1:E], η[1:E]
     rxs = Reaction[]
     # Node-specific reactions for mRNA and Protein
@@ -352,4 +352,4 @@ function grn_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target
         end
         return hit_rate
     end
-end
+end

src/protein_interaction_network.jl

@@ -17,7 +17,7 @@ function protein_interaction_network(connectivity::AbstractMatrix)
     # Number of edges
     E = count(!iszero, connectivity)
     @variables t
-    @species (X(t))[collect(1:N)]
+    @species (X(t))[1:N]
     @parameters α[1:N], β[1:N], γ[1:E], κ[1:E], η[1:E]
     rxs = Reaction[]
     # Node-specific reactions
@@ -344,4 +344,4 @@ function pin_hit_rate(connectivity::AbstractMatrix, initial_samples::Int; target
         end
         return hit_rate
     end
-end
+end

test/gene_regulatory_network/grn_tests.jl

@@ -11,8 +11,8 @@ using BiologicalOscillations, Catalyst, ModelingToolkit, DifferentialEquations
 # Test that gene_regulatory_network creates the correct ReactionSystems
 # Repressilator
 @variables t
-@species (m(t))[collect(1:3)]
-@species (X(t))[collect(1:3)]
+@species (m(t))[1:3]
+@species (X(t))[1:3]
 @parameters α[1:3], β[1:3], δ[1:3], γ[1:3], κ[1:3], η[1:3]
 
 true_reactions = [Reaction(1 - α[1]*m[1], nothing, [m[1]]), 
@@ -127,4 +127,4 @@ for (idx, row) in enumerate(eachrow(all_params))
     params_1 = parameter_sets[idx, :]
     params_2 = collect(row)
     @test params_1 == params_2
-end
+end

test/protein_interaction_network/pin_tests.jl

@@ -11,7 +11,7 @@ using DataFrames, Peaks
 @test_throws DomainError protein_interaction_network([0 1; 3 0])
 ## Correct reactions - Repressilator
 @variables t
-@species (X(t))[collect(1:3)]
+@species (X(t))[1:3]
 @parameters α[1:3], β[1:3], γ[1:3], κ[1:3], η[1:3]
 
 true_reactions = [Reaction(α[1]*(1.0 - X[1]), nothing, [X[1]]), 
@@ -323,4 +323,4 @@ hyperparameters["simulation_output"] = sim_output_config
 result = find_pin_oscillations(connectivity_T0, samples; 
                                hyperparameters=hyperparameters)
 @test result["simulation_result"] isa DataFrame
-@test size(result["simulation_result"], 2) == 21
+@test size(result["simulation_result"], 2) == 21

test/simulation_tests.jl

@@ -6,6 +6,11 @@ model = protein_interaction_network(connectivity)
 parameter_values = [1.0, 44.3, 0.13, 75.3, 2.3, 0.10, 2472.4, 0.27, 3.0, 442.2, 0.91, 4.8, 4410.0, 0.29, 3.3]
 parameter_sets = transpose(repeat(parameter_values, 1, 3))
 initial_conditions = [0.5*ones(3) for i=1:3]
+
+α = parameter_values[1:2:5]
+β = parameter_values[2:2:6]
+γ = parameter_values[7:3:13]
+
 equilibration_times = [pin_timescale(α, β, γ) for i=1:3]
 
 equilibration_data = equilibrate_ODEs(model, parameter_sets, initial_conditions, equilibration_times)