infogap_results/opportuneness_vs_robustness.png added

examples/learn_julia/julia_tutorial.jl

@@ -82,8 +82,8 @@ PyPlot.loglog(S2,linewidth=3,"go--")
 PyPlot.close("all")
 
 PyPlot.figure()
-PyPlot.plot(S2, label = "sv's of A",color="green", linestyle="dashed", marker="o",markersize=12)
-PyPlot.plot(S, label = "sv's of B",color="magenta", linestyle="dashed", marker="o",markersize=12)
+PyPlot.plot(S2, label = "sv's of A",color="green", linestyle="dashed", marker="o", markersize=12)
+PyPlot.plot(S, label = "sv's of B",color="magenta", linestyle="dashed", marker="o", markersize=12)
 PyPlot.legend()
 
 PyPlot.figure()

examples/model_analysis/infogap-jumplin.jl

@@ -1,3 +1,3 @@
 include(Pkg.dir("Mads") * "/src/MadsInfoGap.jl")
 
-infogap_jumplin(horizons=[0.05, 0.1, 0.2, 0.5], retries=1, maxiter=1000, verbosity=0, seed=2015)
+min, max = infogap_jump_polinomial(model=1, horizons=[0.1, 0.2, 0.5], retries=1, maxiter=1000, verbosity=0, seed=2015)

examples/model_analysis/infogap.jl

@@ -1,9 +1,24 @@
-import Mads
+# import Mads
 import Gadfly
 include(Pkg.dir("Mads") * "/src/MadsInfoGap.jl")
 
-md = Mads.loadmadsfile("models/internal-polynomial.mads")
+# md = Mads.loadmadsfile("models/internal-polynomial.mads")
+
+if !isdir("infogap_results")
+	mkdir("infogap_results")
+end
 
 info("Information Gap analysis")
 
-infogap_jump(md, retries=2, seed=2016)
+h = [0.001, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1]
+lmin = Array(Any, 3)
+lmax = Array(Any, 3)
+colors = ["blue", "red", "green"]
+for i = 1:3
+	min, max = infogap_jump_polinomial(model=i, horizons=h, retries=10, maxiter=1000, verbosity=0, seed=2015)
+	lmin[i] = Gadfly.layer(x=min, y=h, Gadfly.Geom.line, Gadfly.Theme(default_color=parse(Colors.Colorant, colors[i])))
+	lmax[i] = Gadfly.layer(x=max, y=h, Gadfly.Geom.line, Gadfly.Theme(default_color=parse(Colors.Colorant, colors[i])))
+end
+f = Gadfly.plot(lmin..., lmax..., Gadfly.Guide.xlabel("o5"), Gadfly.Guide.ylabel("Horizon of uncertainty"), Gadfly.Guide.title("Opportuneness vs. Robustness"), Gadfly.Guide.manual_color_key("Models", ["y = a * t + c", "y = a * t^(1.1) + b * t + c", "y = a * t^n + b * t + c"], colors))
+Gadfly.draw(Gadfly.PNG("infogap_results/opportuneness_vs_robustness.png", 6Gadfly.inch, 4Gadfly.inch), f)
+display(f)

src/MadsInfoGap.jl

@@ -138,23 +138,27 @@ function infogap_jump_polinomial(madsdata::Associative=Dict(); horizons::Vector=
 		pmin = [-10,-10,-5,-3]
 		pmax = [10,10,5,3]
 	end
+	pi = similar(pinit)
 	par_best = Array(Float64, 0)
 	obs_best = Array(Float64, 0)
 	hmin = Array(Float64, 0)
 	hmax = Array(Float64, 0)
 	for h in horizons
 		for mm = ("Min", "Max")
-			phi_best = (mm == "Min") ? 0 : Inf
+			phi_best = (mm == "Max") ? 0 : Inf
 			for r = 1:retries
 				m = JuMP.Model(solver=Ipopt.IpoptSolver(max_iter=maxiter, print_level=verbosity))
-				if r > 1 || random
+				if r == 1 && !random
 					for i = 1:np
-						p[i] = rand() * (pmax[i] - pmin[i]) + pmin[i]
+						pi[i] = pinit[i]
 					end
 				else
-					p = pinit
+					for i = 1:np
+						pi[i] = rand() * (pmax[i] - pmin[i]) + pmin[i]
+					end
 				end
-				@JuMP.variable(m, p[i = 1:np], start = p[i])
+				# @show pi
+				@JuMP.variable(m, p[i = 1:np], start = pi[i])
 				@JuMP.variable(m, o[1:no])
 				@JuMP.constraint(m, p[i = 1:np] .>= pmin[i = 1:np])
 				@JuMP.constraint(m, p[i = 1:np] .<= pmax[i = 1:np])
@@ -181,7 +185,8 @@ function infogap_jump_polinomial(madsdata::Associative=Dict(); horizons::Vector=
 				end
 				JuMP.solve(m)
 				phi = JuMP.getobjectivevalue(m)
-				if mm == "Min"
+				# !quiet && println("OF = $(phi)")
+				if mm == "Max"
 					if phi_best < phi
 						phi_best = phi
 						par_best = JuMP.getvalue(p)