From 442cd53075dedbf484a1b29f91546cae07b64281 Mon Sep 17 00:00:00 2001 From: Thomas Pellet Date: Sat, 30 Mar 2019 13:00:15 -0400 Subject: [PATCH 1/5] TP 1 --- setup | 16 +++++ src/HWconstrained.jl | 140 ++++++++++++++++++++++++------------------- 2 files changed, 94 insertions(+), 62 deletions(-) create mode 100644 setup diff --git a/setup b/setup new file mode 100644 index 0000000..3bd5b1b --- /dev/null +++ b/setup @@ -0,0 +1,16 @@ +# Set up + +cd(joinpath(DEPOT_PATH[1],"dev","HWconstrained")) # go to package location +] activate . +] instantiate + +] activate # no args goes back to main +using HWconstrained # precompiles +data() # errors: function is incomplete! go and complete it in your editor! + +# if in the package directory: +] activate . +] test + +# from any other location +] test HWconstrained diff --git a/src/HWconstrained.jl b/src/HWconstrained.jl index 3b6faa4..b65b268 100644 --- a/src/HWconstrained.jl +++ b/src/HWconstrained.jl @@ -5,24 +5,51 @@ greet() = print("Hello World!") using JuMP, NLopt, DataFrames, Ipopt using LinearAlgebra - export data, table_NLopt, table_JuMP + export data, table_NLopt, table_JuMP, max_JuMP, max_NLopt function data(a=0.5) - - - - - - - - - + na = 3 + p = ones(3) + e = [2.0, 0.0, 0.0] + z = zeros(16,3) + z_2 = [0.72, 0.92, 1.12, 1.32] + z_3 = [0.86, 0.96, 1.06, 1.16] + z[:,1] .= 1.0 + for i=1:4 + z[1+(4*(i-1)):4+(4*(i-1)), 2] .= z_2[i] + end + for j = 1:4 + z[j,3] = z_3[j] + z[j+4,3] = z_3[j] + z[j+4*2,3] = z_3[j] + z[j+4*3,3] = z_3[j] + end + ns = length(z_2) + nss = length(z_2) * length(z_3) + pi = 1/nss + nc = 1 + + # na number of assets + # ns number of states per asset + # nss total number of assets + # nc number of constraints + + return Dict("a"=>a,"na"=>na,"nc"=>nc,"ns"=>ns,"nss"=>nss,"e"=>e,"p"=>p,"z"=>z,"pi"=>pi) end +dat = data() function max_JuMP(a=0.5) + m = Model(with_optimizer(Ipopt.Optimizer)) + @variable(m, omega[1:3]) + @variable(m, c >= 0.0) + @NLconstraint(m, c + sum(data()["p"][i] * (omega[i] - data()["e"][i]) for i in 1:3) == 0.0) + @NLobjective(m, Max, + -exp(-a*c) + data()["pi"] * sum(-exp(-a * + sum(data()["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data()["nss"])) + JuMP.optimize!(m) return Dict("obj"=>objective_value(m),"c"=>value(c),"omegas"=>[value(omega[i]) for i in 1:length(omega)]) end @@ -39,66 +66,55 @@ greet() = print("Hello World!") return d end - - - - - - + ## Does not work function obj(x::Vector,grad::Vector,data::Dict) - - - - - - - - - - - - - - - - - - - - - - - - + if length(grad) > 0 + # grad wrt x1 + a = data["a"] + # grad wrt c + grad[4] = (-a) * exp(-a * c) + # grad with respect to omega 1 + grad[1] = (-a) * data["pi"] * sum(data["z"][i,1] * exp(-a * + sum(data["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data["nss"]) + + # grad wrt to omega 2 + grad[2] = (-a) * data["pi"] * sum(data["z"][i,2] * exp(-a * + sum(data["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data["nss"]) + + # grad wrt to omega 3 + grad[3] = (-a) * data["pi"] * sum(data["z"][i,3] * exp(-a * + sum(data["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data["nss"]) + + end + return (-1) * (-exp(-a*x[4]) + data["pi"] * sum(-exp(-a * + sum(data["z"][j,i] * x[i] for i in 1:3)) for j in 1:data["nss"])) end + + function constr(x::Vector,grad::Vector,data::Dict) - - - - - - - - - - - - - - + if length(grad) > 0 + grad[4] = 1 + grad[1] = data["p"][1] + grad[2] = data["p"][2] + grad[3] = data["p"][3] + + end + return x[4] + sum(data["p"][i] * (x[i] - data["e"][i]) for i in 1:3) + end + function max_NLopt(a=0.5) - - - - - - - - + opt = Opt(:LD_MMA, 4) + lower_bounds!(opt, [-Inf, -Inf, -Inf, 0.]) + xtol_rel!(opt, 1e-4) + min_objective!(opt,(x,grad) -> obj(x,grad,data())) + inequality_constraint!(opt, (x,grad) -> constr(x,grad,data())) + ftol_rel!(opt,1e-9) + vector_init = vcat(data()["e"], 1.0) + NLopt.optimize(opt, vector_init) end function table_NLopt() From 7eb47d40a759393319bfc05522ca03de870c4b8a Mon Sep 17 00:00:00 2001 From: Thomas Pellet Date: Sat, 30 Mar 2019 16:52:36 -0400 Subject: [PATCH 2/5] TP 2 --- setup | 7 ++++++- src/HWconstrained.jl | 40 ++++++++++++++++++++++++++-------------- 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/setup b/setup index 3bd5b1b..8e6d548 100644 --- a/setup +++ b/setup @@ -1,7 +1,7 @@ # Set up cd(joinpath(DEPOT_PATH[1],"dev","HWconstrained")) # go to package location -] activate . +] activate . ] instantiate ] activate # no args goes back to main @@ -14,3 +14,8 @@ data() # errors: function is incomplete! go and complete it in your editor! # from any other location ] test HWconstrained + +# commit changes +git add . # adds everything. edit if that's not what you want +git commit -m 'my homework' +git push diff --git a/src/HWconstrained.jl b/src/HWconstrained.jl index b65b268..a9ec3be 100644 --- a/src/HWconstrained.jl +++ b/src/HWconstrained.jl @@ -69,26 +69,33 @@ dat = data() ## Does not work function obj(x::Vector,grad::Vector,data::Dict) + a = data["a"] + println("$a") + if length(grad) > 0 - # grad wrt x1 - a = data["a"] # grad wrt c - grad[4] = (-a) * exp(-a * c) + grad[4] = (-a) * exp(-a * x[4]) + # grad with respect to omega 1 grad[1] = (-a) * data["pi"] * sum(data["z"][i,1] * exp(-a * - sum(data["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data["nss"]) + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"]) # grad wrt to omega 2 grad[2] = (-a) * data["pi"] * sum(data["z"][i,2] * exp(-a * - sum(data["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data["nss"]) + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"]) # grad wrt to omega 3 grad[3] = (-a) * data["pi"] * sum(data["z"][i,3] * exp(-a * - sum(data["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data["nss"]) + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"]) end + # test = (-1) * (-exp(-a*x[4]) + data["pi"] * sum(-exp(-a * + # (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"])) + # println(grad) + + return (-1) * (-exp(-a*x[4]) + data["pi"] * sum(-exp(-a * - sum(data["z"][j,i] * x[i] for i in 1:3)) for j in 1:data["nss"])) + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"])) end @@ -107,13 +114,17 @@ dat = data() function max_NLopt(a=0.5) + count = 0 # keep track of # function evaluations opt = Opt(:LD_MMA, 4) lower_bounds!(opt, [-Inf, -Inf, -Inf, 0.]) - xtol_rel!(opt, 1e-4) - min_objective!(opt,(x,grad) -> obj(x,grad,data())) - inequality_constraint!(opt, (x,grad) -> constr(x,grad,data())) + xtol_rel!(opt, 1e-7) + grad = zeros(4) + println(length(grad)) + d = data(a) + min_objective!(opt,(x,grad) -> obj(x,grad,d)) + inequality_constraint!(opt, (x,grad) -> constr(x,grad,d)) ftol_rel!(opt,1e-9) - vector_init = vcat(data()["e"], 1.0) + vector_init = vcat(d["e"], 0) NLopt.optimize(opt, vector_init) end @@ -121,10 +132,11 @@ dat = data() d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3)) for i in 1:nrow(d) xx = max_NLopt(d[i,:a]) - for j in 2:ncol(d)-1 - d[i,j] = xx[2][j-1] + for j in 1:ncol(d)-3 + d[i,j+2] = xx[2][j] end - d[i,end] = xx[1] + d[i,2] = xx[2][end] + d[i,end] = -xx[1] end return d end From c4f44a7b33d2de0f7dffe53a3d070c484566b68b Mon Sep 17 00:00:00 2001 From: Thomas Pellet Date: Sat, 30 Mar 2019 19:03:02 -0400 Subject: [PATCH 3/5] TP 3 --- src/HWconstrained.jl | 2 - src/HWconstrained_working.jl | 145 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 49 ++++++++---- 3 files changed, 181 insertions(+), 15 deletions(-) create mode 100644 src/HWconstrained_working.jl diff --git a/src/HWconstrained.jl b/src/HWconstrained.jl index a9ec3be..a1b5013 100644 --- a/src/HWconstrained.jl +++ b/src/HWconstrained.jl @@ -70,7 +70,6 @@ dat = data() ## Does not work function obj(x::Vector,grad::Vector,data::Dict) a = data["a"] - println("$a") if length(grad) > 0 # grad wrt c @@ -119,7 +118,6 @@ dat = data() lower_bounds!(opt, [-Inf, -Inf, -Inf, 0.]) xtol_rel!(opt, 1e-7) grad = zeros(4) - println(length(grad)) d = data(a) min_objective!(opt,(x,grad) -> obj(x,grad,d)) inequality_constraint!(opt, (x,grad) -> constr(x,grad,d)) diff --git a/src/HWconstrained_working.jl b/src/HWconstrained_working.jl new file mode 100644 index 0000000..0319fd0 --- /dev/null +++ b/src/HWconstrained_working.jl @@ -0,0 +1,145 @@ +module HWconstrained + +greet() = print("Hello World!") + + using JuMP, NLopt, DataFrames, Ipopt + using LinearAlgebra + + export data, table_NLopt, table_JuMP, max_JuMP, max_NLopt + + function data(a=0.5) + na = 3 + p = ones(3) + e = [2.0, 0.0, 0.0] + z = zeros(16,3) + z_2 = [0.72, 0.92, 1.12, 1.32] + z_3 = [0.86, 0.96, 1.06, 1.16] + z[:,1] .= 1.0 + for i=1:4 + z[1+(4*(i-1)):4+(4*(i-1)), 2] .= z_2[i] + end + for j = 1:4 + z[j,3] = z_3[j] + z[j+4,3] = z_3[j] + z[j+4*2,3] = z_3[j] + z[j+4*3,3] = z_3[j] + end + ns = length(z_2) + nss = length(z_2) * length(z_3) + pi = 1/nss + nc = 1 + + # na number of assets + # ns number of states per asset + # nss total number of assets + # nc number of constraints + + + return Dict("a"=>a,"na"=>na,"nc"=>nc,"ns"=>ns,"nss"=>nss,"e"=>e,"p"=>p,"z"=>z,"pi"=>pi) + end + +dat = data() + + function max_JuMP(a=0.5) + + m = Model(with_optimizer(Ipopt.Optimizer)) + @variable(m, omega[1:3]) + @variable(m, c >= 0.0) + @NLconstraint(m, c + sum(data()["p"][i] * (omega[i] - data()["e"][i]) for i in 1:3) == 0.0) + @NLobjective(m, Max, + -exp(-a*c) + data()["pi"] * sum(-exp(-a * + sum(data()["z"][j,i] * omega[i] for i in 1:3)) for j in 1:data()["nss"])) + JuMP.optimize!(m) + return Dict("obj"=>objective_value(m),"c"=>value(c),"omegas"=>[value(omega[i]) for i in 1:length(omega)]) + end + + function table_JuMP() + d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3)) + for i in 1:nrow(d) + xx = max_JuMP(d[i,:a]) + d[i,:c] = xx["c"] + d[i,:omega1] = xx["omegas"][1] + d[i,:omega2] = xx["omegas"][2] + d[i,:omega3] = xx["omegas"][3] + d[i,:fval] = xx["obj"] + end + return d + end + + + ## Does not work + function obj(x::Vector,grad::Vector,data::Dict) + a = data["a"] + + if length(grad) > 0 + # grad wrt c + grad[4] = (-a) * exp(-a * x[4]) + + # grad with respect to omega 1 + grad[1] = (-a) * data["pi"] * sum(data["z"][i,1] * exp(-a * + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"]) + + # grad wrt to omega 2 + grad[2] = (-a) * data["pi"] * sum(data["z"][i,2] * exp(-a * + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"]) + + # grad wrt to omega 3 + grad[3] = (-a) * data["pi"] * sum(data["z"][i,3] * exp(-a * + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"]) + + end + # test = (-1) * (-exp(-a*x[4]) + data["pi"] * sum(-exp(-a * + # (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"])) + # println(grad) + + + return (-1) * (-exp(-a*x[4]) + data["pi"] * sum(-exp(-a * + (data["z"][i,1] * x[1] + data["z"][i,2] * x[2] + data["z"][i,3] * x[3])) for i in 1:data["nss"])) + end + + + + function constr(x::Vector,grad::Vector,data::Dict) + if length(grad) > 0 + grad[4] = 1 + grad[1] = data["p"][1] + grad[2] = data["p"][2] + grad[3] = data["p"][3] + + end + return x[4] + sum(data["p"][i] * (x[i] - data["e"][i]) for i in 1:3) + + end + + + function max_NLopt(a=0.5) + count = 0 # keep track of # function evaluations + opt = Opt(:LD_MMA, 4) + lower_bounds!(opt, [-Inf, -Inf, -Inf, 0.]) + xtol_rel!(opt, 1e-7) + grad = zeros(4) + println(length(grad)) + d = data(a) + min_objective!(opt,(x,grad) -> obj(x,grad,d)) + inequality_constraint!(opt, (x,grad) -> constr(x,grad,d)) + ftol_rel!(opt,1e-9) + vector_init = vcat(d["e"], 0) + NLopt.optimize(opt, vector_init) + end + + function table_NLopt() + d = DataFrame(a=[0.5;1.0;5.0],c = zeros(3),omega1=zeros(3),omega2=zeros(3),omega3=zeros(3),fval=zeros(3)) + for i in 1:nrow(d) + xx = max_NLopt(d[i,:a]) + for j in 1:ncol(d)-3 + d[i,j+2] = xx[2][j] + end + d[i,2] = xx[2][end] + d[i,end] = -xx[1] + end + return d + end + + + +end # module diff --git a/test/runtests.jl b/test/runtests.jl index 1793ffc..e54ffc7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,22 +5,44 @@ using Test @testset "testing components" begin @testset "tests gradient of objective function" begin - - - - - - + + truth = HWconstrained.DataFrame(a=[0.5;1.0;5.0], + c = [1.008;1.004;1.0008], + omega1=[-1.41237;-0.20618;0.758763], + omega2=[0.801455;0.400728;0.0801455], + omega3=[1.60291;0.801455;0.160291], + fval=[-1.20821;-0.732819;-0.013422]) + x = zeros(3,4) + x[:,1] = truth[:omega1] + x[:,2] = truth[:omega2] + x[:,3] = truth[:omega3] + x[:,4] = truth[:c] + grad = zeros(4) + d = data() + # Extract gradient + f(x) = grad -> HWconstrained.obj(x[1,:],grad,d) + for i in 1:HWconstrained.nrow(truth) + d = data(truth[:a][i]) + @test isa(HWconstrained.obj(x[i,:],grad,d) , Real) == true + + end + + + + + + + end @testset "tests gradient of constraint function" begin - - - - - - + + + + + + end end @@ -34,7 +56,7 @@ using Test fval=[-1.20821;-0.732819;-0.013422]) @testset "checking result of NLopt maximization" begin - + tol2 = 1e-1 t1 = table_NLopt() for c in names(truth) @test all(maximum.(abs.(t1[c].-truth[c])) .< tol2) @@ -43,6 +65,7 @@ using Test @testset "checking result of NLopt maximization" begin + tol2 = 1e-1 t1 = table_JuMP() for c in names(truth) @test all(maximum.(abs.(t1[c].-truth[c])) .< tol2) From 907a8ea4d46c4415988156de2e53e9ba9b1b6317 Mon Sep 17 00:00:00 2001 From: Thomas Pellet Date: Sun, 31 Mar 2019 08:46:59 -0400 Subject: [PATCH 4/5] TP/JS --- test/runtests.jl | 37 +++++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e54ffc7..6c1e0bd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,27 +17,40 @@ using Test x[:,2] = truth[:omega2] x[:,3] = truth[:omega3] x[:,4] = truth[:c] - grad = zeros(4) + grad = rand(size(x,2)) + testvec = deepcopy(grad) d = data() # Extract gradient - f(x) = grad -> HWconstrained.obj(x[1,:],grad,d) - for i in 1:HWconstrained.nrow(truth) - d = data(truth[:a][i]) - @test isa(HWconstrained.obj(x[i,:],grad,d) , Real) == true - - end - - - - - + grad! = (x,grad) -> HWconstrained.obj(x,grad,d) + r = grad!(x[1,:],grad) + # Cannot test gradient is zero since constrained regression + @test grad != testvec end @testset "tests gradient of constraint function" begin + truth = HWconstrained.DataFrame(a=[0.5;1.0;5.0], + c = [1.008;1.004;1.0008], + omega1=[-1.41237;-0.20618;0.758763], + omega2=[0.801455;0.400728;0.0801455], + omega3=[1.60291;0.801455;0.160291], + fval=[-1.20821;-0.732819;-0.013422]) + x = zeros(3,4) + x[:,1] = truth[:omega1] + x[:,2] = truth[:omega2] + x[:,3] = truth[:omega3] + x[:,4] = truth[:c] + grad = rand(size(x,2)) + testvec = deepcopy(grad) + d = data() + # Extract gradient + gradconstr! = (x,grad) -> HWconstrained.constr(x,grad,d) + r = gradconstr!(x[1,:],grad) + # Cannot test whether the gradient is zero since it is a constraint regression + @test grad != testvec From 32c2325fca67ec750074c315c3c490d06963cb4a Mon Sep 17 00:00:00 2001 From: Thomas Pellet Date: Sun, 31 Mar 2019 08:49:05 -0400 Subject: [PATCH 5/5] TP/JS 2 --- setup | 21 --------------------- 1 file changed, 21 deletions(-) delete mode 100644 setup diff --git a/setup b/setup deleted file mode 100644 index 8e6d548..0000000 --- a/setup +++ /dev/null @@ -1,21 +0,0 @@ -# Set up - -cd(joinpath(DEPOT_PATH[1],"dev","HWconstrained")) # go to package location -] activate . -] instantiate - -] activate # no args goes back to main -using HWconstrained # precompiles -data() # errors: function is incomplete! go and complete it in your editor! - -# if in the package directory: -] activate . -] test - -# from any other location -] test HWconstrained - -# commit changes -git add . # adds everything. edit if that's not what you want -git commit -m 'my homework' -git push