Merge pull request #33 from himoto/develop

Release v0.4.1

Author: Hiroaki Imoto

.travis.yml

@@ -1,7 +1,7 @@
 name = "BioMASS"
 uuid = "324734c7-f323-4536-9335-775d9be9d101"
 authors = ["Hiroaki Imoto <himoto@protein.osaka-u.ac.jp>"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

Project.toml

@@ -15,8 +15,9 @@ This module provides a Julia interface to the [BioMASS](https://github.com/okada
 
 BioMASS.jl supports:
 
-- Parameter Estimation of ODE/DDE models
-- Bifurcation Analysis
+- parameter estimation of ODE/DDE models
+- visualization of simulation results
+- bifurcation analysis
 
 ## Usage
 
@@ -32,7 +33,7 @@ optimize(model, 1, max_generation=20000, allowable_error=0.5)
 visualize(model, viz_type="best", show_all=true)
 ```
 
-## Convert optimized parameters into BioMASS format
+## Convert optimized parameters to BioMASS format
 
 ![](docs/src/assets/conversion.png)
 

README.md

@@ -9,6 +9,7 @@ This module provides a Julia interface to the [BioMASS](https://github.com/okada
 ## Features
 
 - Parameter estimation
+- Visualization of simulation results
 - Bifurcation analysis
 
 ## Installation

docs/src/index.md

@@ -53,6 +53,10 @@ Find a parameter set that reproduces experimental observations.
 
     - The maximum number of generations over which the entire population is evolved.
 
+  - `initial_threshold`::Float64 (default: 1e12)
+
+    - Allowable error used to generate initial population. Default is 1e12 (numerically solvable).
+
   - `allowable_error`::Float64 (default: 0.0)
 
     - Optimization stops when Best Fitness <= allowable_error.
@@ -92,9 +96,13 @@ Save simulation results with optimized parameter values.
     - Whether to show all simulation results.
 
   - `stdev`::Bool (default: `false`)
+
     - If True, the standard deviation of simulated values will be shown
       (only available for `"average"` visualization type).
 
+  - `save_format`::String (default: `"pdf"`)
+    - Either "png" or "pdf", indicating whether to save figures as png or pdf format.
+
 ---
 
 **param2biomass**(`path_to_model`::String)

docs/src/usage/parameter_estimation.md

@@ -95,18 +95,18 @@ function get_search_region()::Matrix{Float64}
     search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()
     search_param::Vector{Float64} = init_search_param(search_idx, p, u0)
 
-    search_rgn::Matrix{Float64} = zeros(2, length(p)+length(u0))
+    search_rgn::Matrix{Float64} = zeros(2, length(p) + length(u0))
 
     # Default: 0.1 ~ 10x
-    for (i,j) in enumerate(search_idx[1])
-        search_rgn[1,j] = search_param[i]*0.1  # lower bound
-        search_rgn[2,j] = search_param[i]*10.0  # upper bound
+    for (i, j) in enumerate(search_idx[1])
+        search_rgn[1,j] = search_param[i] * 0.1  # lower bound
+        search_rgn[2,j] = search_param[i] * 10.0  # upper bound
     end
 
     # Default: 0.5 ~ 2x
-    for (i,j) in enumerate(search_idx[2])
-        search_rgn[1,j+length(p)] = search_param[i+length(search_idx[1])]*0.5  # lower bound
-        search_rgn[2,j+length(p)] = search_param[i+length(search_idx[1])]*2.0  # upper bound
+    for (i, j) in enumerate(search_idx[2])
+        search_rgn[1,j + length(p)] = search_param[i + length(search_idx[1])] * 0.5  # lower bound
+        search_rgn[2,j + length(p)] = search_param[i + length(search_idx[1])] * 2.0  # upper bound
     end
 
     # search_rgn[:, C.param_name] = [lower_bound, upper_bound]
@@ -200,11 +200,11 @@ function update_param(indiv::Vector{Float64})::Tuple{Array{Float64,1},Array{Floa
 
     search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index()
 
-    for (i,j) in enumerate(search_idx[1])
+    for (i, j) in enumerate(search_idx[1])
         @inbounds p[j] = indiv[i]
     end
-    for (i,j) in enumerate(search_idx[2])
-        @inbounds u0[j] = indiv[i+length(search_idx[1])]
+    for (i, j) in enumerate(search_idx[2])
+        @inbounds u0[j] = indiv[i + length(search_idx[1])]
     end
 
     # constraints --------------------------------------------------------------
@@ -240,7 +240,7 @@ function decode_gene2val(indiv_gene::Vector{Float64})::Vector{Float64}
         )
     end
 
-    return round.(indiv,sigdigits=7)
+    return round.(indiv, sigdigits=7)
 end
 
 
@@ -261,18 +261,18 @@ end
 
 
 function encode_bestIndivVal2randGene(
-        idx::Int64,
+        gene_idx::Int64,
         best_indiv::Vector{Float64},
         p0_bounds::Vector{Float64})::Float64
     search_rgn::Matrix{Float64} = get_search_region()
     rand_gene::Float64 = (
         log10(
-            best_indiv[idx]*10^(
-                rand() * log10(p0_bounds[2]/p0_bounds[1]) + log10(p0_bounds[1])
+            best_indiv[gene_idx] * 10^(
+                rand() * log10(p0_bounds[2] / p0_bounds[1]) + log10(p0_bounds[1])
             )
-        ) - search_rgn[1,idx]
+        ) - search_rgn[1,gene_idx]
     ) / (
-        search_rgn[2,idx] - search_rgn[1,idx]
+        search_rgn[2,gene_idx] - search_rgn[1,gene_idx]
     )
     return rand_gene
 end
@@ -284,15 +284,15 @@ function init_search_param(
         u0::Vector{Float64})::Vector{Float64}
     duplicate::Vector{String} = []
     if length(search_idx[1]) != length(unique(search_idx[1]))
-        for idx in findall([count(x->x==i,search_idx[1]) 
+        for idx in findall([count(x -> x == i, search_idx[1]) 
                             for i in unique(search_idx[1])] .!= 1)
             push!(duplicate, C.NAMES[search_idx[1][idx]])
         end
         error(
             "Duplicate parameters (C.): $duplicate"
         )
     elseif length(search_idx[2]) != length(unique(search_idx[2]))
-        for idx in findall([count(x->x==i,search_idx[2]) 
+        for idx in findall([count(x -> x == i, search_idx[2]) 
                             for i in unique(search_idx[2])] .!= 1)
             push!(duplicate, V.NAMES[search_idx[2][idx]])
         end
@@ -303,11 +303,11 @@ function init_search_param(
     search_param = zeros(
         length(search_idx[1]) + length(search_idx[2])
     )
-    for (i,j) in enumerate(search_idx[1])
+    for (i, j) in enumerate(search_idx[1])
         @inbounds search_param[i] = p[j]
     end
-    for (i,j) in enumerate(search_idx[2])
-        @inbounds search_param[i+length(search_idx[1])] = u0[j]
+    for (i, j) in enumerate(search_idx[2])
+        @inbounds search_param[i + length(search_idx[1])] = u0[j]
     end
 
     if any(x -> x == 0.0, search_param)
@@ -339,35 +339,35 @@ end
 function conv_lin2log!(
         search_rgn::Matrix{Float64},
         search_idx::Tuple{Array{Int64,1},Array{Int64,1}})::Matrix{Float64}
-    for i=1:size(search_rgn,2)
+    for i = 1:size(search_rgn, 2)
         if minimum(search_rgn[:,i]) < 0.0
             msg = "search_rgn[lower_bound,upper_bound] must be positive.\n"
             if i <= C.NUM
                 error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
             else
-                error(@sprintf("`V.%s` ", V.NAMES[i-C.NUM]) * msg)
+                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
             end
         elseif minimum(search_rgn[:,i]) == 0.0 && maximum(search_rgn[:,i]) != 0.0
             msg = "lower_bound must be larger than 0.\n"
             if i <= C.NUM
                 error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
             else
-                error(@sprintf("`V.%s` ", V.NAMES[i-C.NUM]) * msg)
+                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
             end
         elseif search_rgn[2,i] - search_rgn[1,i] < 0.0
             msg = "lower_bound must be smaller than upper_bound.\n"
             if i <= C.NUM
                 error(@sprintf("`C.%s` ", C.NAMES[i]) * msg)
             else
-                error(@sprintf("`V.%s` ", V.NAMES[i-C.NUM]) * msg)
+                error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg)
             end
         end
     end
 
     nonzero_idx::Vector{Int} = []
-    for i=1:size(search_rgn,2)
+    for i = 1:size(search_rgn, 2)
         if search_rgn[:,i] != [0.0,0.0]
-            push!(nonzero_idx,i)
+            push!(nonzero_idx, i)
         end
     end
     difference::Vector{Int} = collect(
@@ -381,7 +381,7 @@ function conv_lin2log!(
             if idx <= C.NUM
                 println(@sprintf("`C.%s`", C.NAMES[Int(idx)]))
             else
-                println(@sprintf("`V.%s`", V.NAMES[Int(idx)-C.NUM]))
+                println(@sprintf("`V.%s`", V.NAMES[Int(idx) - C.NUM]))
             end
         end
         error(

examples/fos_model/set_search_param.jl

@@ -23,7 +23,7 @@ for observable in observables
 end
 
 const dt = 1.0
-t = collect(0.0:dt:5400.0)  # 0, 1, 2, ..., 5400 [sec.]
+const t = collect(0.0:dt:5400.0)  # 0, 1, 2, ..., 5400 [sec.]
 
 const conditions = ["EGF", "HRG"]
 
@@ -42,7 +42,11 @@ function solveode(
         prob = ODEProblem(f, u0, (t[1], t[end]), p)
         sol = solve(
             prob,CVODE_BDF(),
-            abstol=ABSTOL,reltol=RELTOL,saveat=dt,verbose=false
+            abstol=ABSTOL,
+            reltol=RELTOL,
+            saveat=dt,
+            dtmin=eps(),
+            verbose=false
         )
         is_successful = ifelse(sol.retcode === :Success, true, false)
     catch
@@ -67,9 +71,13 @@ function get_steady_state(
         sol = solve(
             prob,
             DynamicSS(
-                CVODE_BDF();abstol=ABSTOL,reltol=RELTOL
+                CVODE_BDF();
+                abstol=ABSTOL,
+                reltol=RELTOL
             ),
-            dt=dt,verbose=false
+            dt=dt,
+            dtmin=eps(),
+            verbose=false
         )
         is_successful = ifelse(sol.retcode === :Success, true, false)
     catch