Merge pull request #14 from Vaibhavdixit02/sobolmethod

Sobol Method implementation

src/sobol_sensitivity.jl

@@ -0,0 +1,143 @@
+function give_rand_p(p_range,p_fixed=nothing,indices=nothing)
+    if p_fixed == nothing
+        p = [(p_range[j][2] -p_range[j][1])*rand() + p_range[j][1] for j in 1:length(p_range)]
+    else
+        p =  zeros(length(p_range))
+        j = 1
+        for i in 1:length(p_range)
+            if i in indices
+                p[i] = p_fixed[j]
+                j += 1
+            else
+                p[i] = (p_range[i][2] -p_range[i][1])*rand() + p_range[i][1]
+            end
+        end
+    end
+    p
+end
+
+function calc_mean_var(f,p_range,N)
+    y1 = Array(f(give_rand_p(p_range)))
+    y0 = zeros(length(y1))
+    v = zeros(length(y1))
+    if length(size(y1)) != 1
+        y0 = reshape(y0,size(y1)[1],size(y1)[2])
+        v = reshape(v,size(y1)[1],size(y1)[2])
+    end
+    for i in 1:N
+        y1 = Array(f(give_rand_p(p_range)))
+        @. y0 += y1
+        @. v += y1^2
+    end
+    y0 = @. y0/N
+    y0_sq = [i.^2 for i in y0]
+    v = @. v/N - y0_sq
+    y0,v
+end
+
+function first_order_var(f,p_range,N,y0)
+    ys = []
+    for i in 1:length(p_range)
+        y = zeros(length(y0))
+        if length(size(y0)) != 1
+            y = reshape(y,size(y0)[1],size(y0)[2])
+        end
+        for j in 1:N
+            p2 = give_rand_p(p_range)
+            p1 = give_rand_p(p_range,[p2[i]],[i])
+            yer =  Array(f(p1)) .* Array(f(p2))
+            @. y += yer
+        end
+        y = @. y/N - y0^2
+        push!(ys,copy(y))
+    end
+    ys
+end
+
+function second_order_var(f,p_range,N,y0)
+    ys = []
+    for i in 1:length(p_range)
+        for j in i+1:length(p_range)
+            y = zeros(length(y0))
+            if length(size(y0)) != 1
+                y = reshape(y,size(y0)[1],size(y0)[2])
+            end
+            for k in 1:N
+                p2 = give_rand_p(p_range)
+                p1 = give_rand_p(p_range,[p2[i],p2[j]],[i,j])
+                @. y +=  Array(f(p1)) * Array(f(p2)) 
+            end
+            y = @. y/N - y0^2
+            push!(ys,copy(y))
+        end
+    end
+    ys_frst_order = first_order_var(f,p_range,N,y0)
+    j = 1
+    for i in 1:length(p_range)
+        for k in i+1:length(p_range)
+            ys[j] = @. ys[j] - ( ys_frst_order[i] + ys_frst_order[k] )
+            j += 1
+        end
+    end
+    ys
+end
+
+
+function total_var(f,p_range,N,y0)
+    ys = []
+    for i in 1:length(p_range)
+        y = zeros(length(y0))
+        if length(size(y0)) != 1
+            y = reshape(y,size(y0)[1],size(y0)[2])
+        end
+        for j in 1:N
+            p_fixed_all = []
+            p_fixed_indices = []
+            p2 = give_rand_p(p_range)
+            for j in 1:length(p2)
+                if j != i
+                    push!(p_fixed_all,p2[j])
+                    push!(p_fixed_indices,j)
+                end
+            end
+            p1 = give_rand_p(p_range,p_fixed_all,p_fixed_indices) 
+            yer =  Array(f(p1)) .* Array(f(p2))
+            @. y += yer
+        end
+        y = @. y/N - y0^2
+        push!(ys,copy(y))
+    end
+    ys
+end
+
+function sobol_sensitivity(f,p_range,N,order=2)
+    y0,v = calc_mean_var(f,p_range,N)
+    if order == 1
+        first_order = first_order_var(f,p_range,N,y0)
+        for i in 1:length(first_order)
+            first_order[i] = @. first_order[i] / v
+        end
+        first_order[2:end]
+    elseif order == 2
+        second_order = second_order_var(f,p_range,N,y0)
+        for i in 1:length(second_order)
+            second_order[i] = @. second_order[i] / v
+        end
+        second_order[2:end]
+    else
+        total_indices = total_var(f,p_range,N,y0)
+        for i in 1:length(total_indices)
+            total_indices[i] = @. 1 - (total_indices[i] / v)
+        end
+        total_indices[2:end]
+    end
+end
+
+function sobol_sensitivity(prob::DEProblem,alg,t,p_range,N,order=2)
+    f = function (p)
+        prob1 = remake(prob;p=p)
+        Array(solve(prob1,alg;saveat=t))
+    end
+    @assert length(prob.p) == length(p_range)
+    sobol_sensitivity(f,p_range,N,order)
+end

test/chemical_model.jl

@@ -0,0 +1,37 @@
+using DiffEqSensitivity, OrdinaryDiffEq, ParameterizedFunctions
+
+# f = @ode_def_nohes chem_model begin
+#     dy1 = -k1*y1
+#     dy2 = k1*y1 - k2*(y2^2)
+#     dy3 = k2*(y2^2)
+# end k1 k2
+
+# p = [1.0,1.0]
+# u0 = [1.0,0.0,0.0]
+# prob1 = ODEProblem(f,u0,(0.0,10.0),p)
+# prob1
+
+A = [1,0,2,3]
+function f1(p)
+    x = A[1]*p[1] + A[2]*p[2]
+    y = A[3]*p[1] + A[4]*p[2]
+    [x,y] 
+end
+p_range = [[0.5,1.5],[0.5,1.5]]
+N = 100000
+total = sobol_sensitivity(f1,p_range,N)
+first_order = sobol_sensitivity(f1,p_range,N,1)
+second_order = sobol_sensitivity(f1,p_range,N,2)
+
+@test [total[1][1],total[2][1]] ≈ [first_order[1][1],first_order[2][1]] + second_order[1][1] atol=1e-1
+@test [total[1][2],total[2][2]] ≈ [first_order[1][2],first_order[2][2]] + second_order[2][2] atol=1e-1
+
+# fitz = @ode_def_nohes FitzhughNagumo begin
+#     dV = c*(V - V^3/3 + R)
+#     dR = -(1/c)*(V -  a - b*R)
+# end a b c
+# u0 = [-1.0;1.0]
+# tspan = (0.0,20.0)
+# p = [0.2,0.2,3.0]
+# prob2 = ODEProblem(fitz,u0,tspan,p)
+# prob2