From d21559aacf3b902fdc236c7e9de5d056f622b48f Mon Sep 17 00:00:00 2001 From: MatthewStuber <33096280+MatthewStuber@users.noreply.github.com> Date: Mon, 26 Mar 2018 09:10:08 -0400 Subject: [PATCH] Direct Contractor Update, Website Start --- .gitignore | 3 + docs/make.jl | 3 + docs/mkdocs.yml | 26 ++ docs/src/FDI.jl | 344 ++++++++++++++++++++++ docs/src/contractors.md | 5 + docs/src/index.md | 23 ++ docs/src/tests.md | 0 docs/src/utilities.md | 0 lib/OldParent.jl | 56 ---- src/EAGOParametricInterval.jl | 3 +- src/build/assets/Documenter.css | 18 -- src/build/assets/mathjaxhelper.js | 25 -- src/build/index.md | 13 - src/docs/build/assets/Documenter.css | 18 -- src/docs/build/assets/mathjaxhelper.js | 25 -- src/docs/build/index.md | 13 - src/docs/make.jl | 9 - src/docs/src/index.md | 7 - src/lib/Parametric_Contractor.jl | 381 ++++++++++++++++++++++++- src/lib/Sparse_Conditioner.jl | 19 ++ test/ParametricContractor_Tests.jl | 16 ++ test/SparseCntr_Tests.jl | 17 ++ 22 files changed, 833 insertions(+), 191 deletions(-) create mode 100644 docs/make.jl create mode 100644 docs/mkdocs.yml create mode 100644 docs/src/FDI.jl create mode 100644 docs/src/contractors.md create mode 100644 docs/src/index.md create mode 100644 docs/src/tests.md create mode 100644 docs/src/utilities.md delete mode 100644 lib/OldParent.jl delete mode 100644 src/build/assets/Documenter.css delete mode 100644 src/build/assets/mathjaxhelper.js delete mode 100644 src/build/index.md delete mode 100644 src/docs/build/assets/Documenter.css delete mode 100644 src/docs/build/assets/mathjaxhelper.js delete mode 100644 src/docs/build/index.md delete mode 100644 src/docs/make.jl delete mode 100644 src/docs/src/index.md diff --git a/.gitignore b/.gitignore index 8c960ec..372ef70 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ *.jl.cov *.jl.*.cov *.jl.mem + +docs/build/ +docs/site/ diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..e618b5d --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,3 @@ +using Documenter, EAGOParametricInterval + +makedocs(pages = ["Home"]) diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml new file mode 100644 index 0000000..0745396 --- /dev/null +++ b/docs/mkdocs.yml @@ -0,0 +1,26 @@ +site_name: EAGOParametricInterval.jl +repo_url: https://github.com/MatthewStuber/EAGOParametricInterval.jl +site_description: Parametric Interval Methods for Bounding Functions +site_author: MatthewStuber + +theme: readthedocs + +extra_css: + - assets/Documenter.css + +extra_javascript: + - https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML + - assets/mathjaxhelper.js + +markdown_extensions: + - extra + - tables + - fenced_code + +docs_dir: 'build' + +pages: + - Parametric Interval Methods: index.md + - Contractors: contractors.md + - Utilities: utilities.md + - Tests: tests.md diff --git a/docs/src/FDI.jl b/docs/src/FDI.jl new file mode 100644 index 0000000..85ef19d --- /dev/null +++ b/docs/src/FDI.jl @@ -0,0 +1,344 @@ +workspace() +#using EAGOSemiInfinite +#using EAGOGlobalSolver +using IntervalArithmetic +using EAGOBranchBound +#using EAGOSmoothMcCormickGrad +#using StaticArrays + + +function FDI_invApprox(x::Interval) + Ap = 0.00005 + g = 9.8 + low = -(x^2)/(2*g*Ap^2) + high = (x^2)/(2*g*Ap^2) + if x.lo>0 + return high + elseif x.hi<0 + return low + else + return Interval(0,max(low.hi,high.hi)) + end +end + +function FDI_invApprox(x::Float64) + Ap = 0.00005 + g = 9.8 + if x>=0 + return (x^2)/(2*g*Ap^2) + else x<0 + return -(x^2)/(2*g*Ap^2) + end +end + +function xbnds(u,p) + + # parameters + Ap = 0.00005 + g = 9.8 + z = [] + # bound z via McCormick operators + z1 = (((u[1]/p[3])/Ap)^2)/(2.0*g) + z3 = (1.0/(2.0*g))*((u[1]+u[2])/(p[4]*Ap))^2 + z2 = (-((u[2]-p[4]*Ap*sqrt(2*g*z3))/(p[5])/Ap)^2)/(2.0*g) + z4 = (((u[1]*p[1]/p[3])/Ap)^2)/(2.0*g) + z6 = (1.0/(2.0*g))*((p[1]*u[1]+u[2])/(p[4]*Ap))^2 + z5 = (-((u[2]-p[4]*Ap*sqrt(2*g*z6))/(p[5])/Ap)^2)/(2.0*g) + z7 = (((u[1]/p[3])/Ap)^2)/(2.0*g) + z9 = (1.0/(2.0*g))*((u[1]+u[2])/(p[4]*(Ap+3.14159*p[2]^2)))^2 + z8 = FDI_invApprox((u[2]-p[4]*(Ap+3.14159*p[2]^2)*sqrt(2.0*g*z9))/(p[5])) + push!(z,z1,z2,z3,z4,z5,z6,z7,z8,z9) + + # rotate coordinates back to X & intersects + + x = [] + x1 = min(max(z1-z2+z3,0.0),0.75) + x2 = min(max(z3,0.0),0.75) + x3 = min(max(z3-z2,0.0),0.75) + x4 = min(max(z4-z5+z6,0.0),0.75) + x5 = min(max(z6,0.0),0.75) + x6 = min(max(z6-z5,0.0),0.75) + x7 = min(max(z7-z8+z9,0.0),0.75) + x8 = min(max(z9,0.0),0.75) + x9 = min(max(z9-z8,0.0),0.75) + push!(x,x1,x2,x3,x4,x5,x6,x7,x8,x9) + + return x +end + +# sets up SIP constraint with imbedded x(p) calculation +function SIP_g(un,p) + u = un[1:2] + xbnd = xbnds(u,p) + r_xbnds = [zero(xbnd[1]) for i=1:3,j=1:3] + r_xbnds[:,1] = xbnd[1:3] + r_xbnds[:,2] = xbnd[4:6] + r_xbnds[:,3] = xbnd[7:9] + out = un[3] + for i=1:3 + for j=1:2 + for k=(j+1):3 + out = out - (r_xbnds[i,j]-r_xbnds[i,k])^2 + end + end + end + return out +end + +# sets up SIP objective +function SIP_f(un) + return -un[3] +end +# sets up SIP +ex_SIP_X = [(1E-5)..(1E-4),(1E-5)..(1E-4),(-10.0)..(10.0)] +ex_SIP_P = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)] + +#Intv = SIP_g(ex_SIP_X,ex_SIP_P) +#Trial1 = [1E-5,1E-5,0.0] +#Intv1 = SIP_g(Trial1,ex_SIP_P) +#Trial2 = [1E-5,1E-4,0.0] +#Intv2 = SIP_g(Trial2,ex_SIP_P) +#Trial3 = [1E-4,1E-5,0.0] +#Intv3 = SIP_g(Trial3,ex_SIP_P) +#Trial4 = [1E-4,1E-4,0.0] +#Intv4 = SIP_g(Trial4,ex_SIP_P) +#Trial5 = [(1E-5+1E-4)/2,(1E-5+1E-4)/2,0.0] +#Intv5 = SIP_g(Trial5,ex_SIP_P) +#Trial6 = [(1E-5+1E-4)/3,(1E-5+1E-4)/3,0.0] +#Intv6 = SIP_g(Trial6,ex_SIP_P) +#Trial7 = [2*(1E-5+1E-4)/3,2*(1E-5+1E-4)/3,0.0] +#Intv7 = SIP_g(Trial7,ex_SIP_P) +#Trial8 = [(1E-5+1E-4)/2,(1E-5+1E-4)/3,0.0] +#Intv8 = SIP_g(Trial8,ex_SIP_P) +#Trial9 = [(1E-5+1E-4)/2,2*(1E-5+1E-4)/3,0.0] +#Intv9 = SIP_g(Trial9,ex_SIP_P) +#Trial10 = [(1E-5+1E-4)/3,(1E-5+1E-4)/2,0.0] +#Intv10 = SIP_g(Trial10,ex_SIP_P) +Trial11 = [2*(1E-5+1E-4)/3,(1E-5+1E-4)/2,0.0] +Intv11 = SIP_g(Trial11,ex_SIP_P) + +ex_SIP_P1a = [0.66,(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)] +ex_SIP_P2a = [(0.54)..(0.66),0.005,(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)] +ex_SIP_P3a = [(0.54)..(0.66),(0.0005)..(0.005),1.15,(0.65)..(0.95),(0.85)..(1.15)] +ex_SIP_P4a = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.95,(0.85)..(1.15)] +ex_SIP_P5a = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),1.15] +Intv11a1 = SIP_g(Trial11,ex_SIP_P1a) +Intv11b1 = SIP_g(Trial11,ex_SIP_P2a) +Intv11c1 = SIP_g(Trial11,ex_SIP_P3a) +Intv11d1 = SIP_g(Trial11,ex_SIP_P4a) +Intv11e1 = SIP_g(Trial11,ex_SIP_P5a) + +ex_SIP_P1b = [0.54,(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)] +ex_SIP_P2b = [(0.54)..(0.66),0.0005,(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)] +ex_SIP_P3b = [(0.54)..(0.66),(0.0005)..(0.005),0.85,(0.65)..(0.95),(0.85)..(1.15)] +ex_SIP_P4b = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.65,(0.85)..(1.15)] +ex_SIP_P5b = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),0.85] +Intv11a2 = SIP_g(Trial11,ex_SIP_P1b) +Intv11b2 = SIP_g(Trial11,ex_SIP_P2b) +Intv11c2 = SIP_g(Trial11,ex_SIP_P3b) +Intv11d2 = SIP_g(Trial11,ex_SIP_P4b) +Intv11e2 = SIP_g(Trial11,ex_SIP_P5b) +ex_SIP_P4ca = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.70,(0.85)..(1.15)] +ex_SIP_P4c = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.75,(0.85)..(1.15)] +ex_SIP_P4da = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.80,(0.85)..(1.15)] +ex_SIP_P4d = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.85,(0.85)..(1.15)] +Intv11d2a = SIP_g(Trial11,ex_SIP_P4ca) +Intv11d3 = SIP_g(Trial11,ex_SIP_P4c) +Intv11d3a = SIP_g(Trial11,ex_SIP_P4da) +Intv11d4 = SIP_g(Trial11,ex_SIP_P4d) + + +disc1 = [0.6,0.0005,1.0,0.65,1.00] +ex_SIP_Xt = [(1E-5)..(1E-4),(1E-5)..(1E-4),(0.0)..(0.0)] +Xtest = SIP_g(ex_SIP_Xt,disc1) +m = BnBModel([(1E-5)..(1E-4),(1E-5)..(1E-4),(-10.0)..(10.0)]) +s = BnBSolver() +set_to_default!(s) +s.BnB_atol = 1E-2 + +# solve optimization problem via interval extension +function test_lower(X,k,pos,opt_inner,UBD) + FInt = -X[3] + disc1 = [0.6,0.005,1.0,0.65,1.00] + GInt = SIP_g(X,disc1) + feas = true + for i=1:length(GInt) + if (GInt[i].lo>0.0) + feas = false + break + end + end + val = FInt.lo + pnt = mid.(X) + return pnt, val, feas, X, [FInt] +end + +function test_upper(X,k,pos,opt_inner,temp) + feas = true + val = temp[1].hi + pnt = mid.(X) + return pnt, val, feas, X +end + +s.Lower_Prob = test_lower +s.Upper_Prob = test_upper +solveBnB!(s,m) + +mi = BnBModel([(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)]) +si = BnBSolver() +set_to_default!(si) +si.BnB_atol = 5E-2 + +function inner_lower(X,k,pos,opt_inner,UBD) + disc1 = [5.5E-5,5.5E-5,4.99] + FInt = -SIP_g(disc1,X) + feas = true + val = FInt.lo + pnt = mid.(X) + return pnt, val, feas, X, [FInt] +end + +function inner_upper(X,k,pos,opt_inner,temp) + feas = true + val = temp[1].hi + pnt = mid.(X) + return pnt, val, feas, X +end + +si.Lower_Prob = inner_lower +si.Upper_Prob = inner_upper +#solveBnB!(si,mi) + +m1 = BnBModel([(1E-5)..(1E-4),(1E-5)..(1E-4),(-10.0)..(10.0)]) +s1 = BnBSolver() +set_to_default!(s1) +s1.BnB_atol = 1E-2 + +# solve optimization problem via interval extension +function test_lower(X,k,pos,opt_inner,UBD) + FInt = -X[3] + disc1 = [0.6,0.005,1.0,0.65,1.00] + GInt = SIP_g(X,disc1) +0.5 + feas = true + for i=1:length(GInt) + if (GInt[i].lo>0.0) + feas = false + break + end + end + val = FInt.lo + pnt = mid.(X) + return pnt, val, feas, X, [FInt] +end + +function test_upper(X,k,pos,opt_inner,temp) + feas = true + val = temp[1].hi + pnt = mid.(X) + return pnt, val, feas, X +end + +s1.Lower_Prob = test_lower +s1.Upper_Prob = test_upper +solveBnB!(s1,m1) + +#= +SIPopt = SIP_opts() +set_to_default!(SIPopt) +SIPopt.Verbosity = "Full" #"Normal" + +# solves sample LBP +vals = 1 +count = 1 +p_sto = zeros(vals^5,5) +for i=1:vals + for j=1:vals + for k=1:vals + for p=1:vals + for q=1:vals + t1 = ex_SIP_P[1].lo + (ex_SIP_P[1].hi - ex_SIP_P[1].lo)*(i-1)/(vals-1) + t2 = ex_SIP_P[2].lo + (ex_SIP_P[2].hi - ex_SIP_P[2].lo)*(j-1)/(vals-1) + t3 = ex_SIP_P[3].lo + (ex_SIP_P[3].hi - ex_SIP_P[3].lo)*(k-1)/(vals-1) + t4 = ex_SIP_P[4].lo + (ex_SIP_P[4].hi - ex_SIP_P[4].lo)*(p-1)/(vals-1) + t5 = ex_SIP_P[5].lo + (ex_SIP_P[5].hi - ex_SIP_P[5].lo)*(q-1)/(vals-1) + p_sto[count,:] = [t1,t2,t3,t4,t5] + count += 1 + end + end + end + end +end + +bnd_z2 = ex_SIP_X[2]-ex_SIP_P[4]*0.00005*sqrt((ex_SIP_P[1]*ex_SIP_X[1]+ex_SIP_X[2])/(ex_SIP_P[5]*(ex_SIP_P[4]*0.00005)^2)) +bnd_z5 = ex_SIP_X[2]-ex_SIP_P[4]*0.00005*sqrt((ex_SIP_P[1]*ex_SIP_X[1]+ex_SIP_X[2])/(ex_SIP_P[5]*(ex_SIP_P[4]*0.00005)^2)) +bnd_z8 = (ex_SIP_X[2]-ex_SIP_P[4]*(0.00005+3.14159*ex_SIP_P[2]^2)*sqrt((((ex_SIP_X[1]+ex_SIP_X[2])/(ex_SIP_P[4]*(0.00005+3.14159*ex_SIP_P[2]^2)))^2)))/(ex_SIP_P[5]) + +SIPopt.tol = 1E-1 +SIPopt.eps_g0 = 1.0 +SIPopt.P_LBD = [p_sto[i,:] for i=1:vals^5] +SIPopt.P_UBD = [p_sto[i,:] for i=1:vals^5] +SIPopt.LBP_Opt.DAG_depth = -1 +SIPopt.LBP_Opt.f = SIP_f +SIPopt.LBP_Opt.g = SIP_g +SIPopt.LBP_Opt.X = ex_SIP_X +SIPopt.LBP_Opt.g = x -> BndProb_reform(x,[],SIP_g,[],0.0) +#SIPopt.LBP_Opt.LBD_func_relax = "Diff2-MV-OFF" +#SIPopt.LBP_Opt.LBD_problem_relax = "NLP2" +#SIPopt.LBP_Opt.LBD_problem_solver = "Ipopt" +SIPopt.LBP_Opt.LBD_func_relax = "Interval" +SIPopt.LBP_Opt.LBD_problem_relax = "Interval" +SIPopt.LBP_Opt.LBD_problem_solver = "Interval" +SIPopt.LBP_Opt.UBD_func_relax = "Interval" +SIPopt.LBP_Opt.UBD_problem_relax = "Interval" +SIPopt.LBP_Opt.UBD_problem_solver = "Interval" +SIPopt.LBP_Opt.BnB_options.Verbosity = "Normal" #"Full" +SIPopt.LBP_Opt.BnB_options.BnB_tol = 1E-2 +#LBDg, xbar, feas, tLBP1,tLBP2 = EAGO_Global_Explicit(SIPopt.LBP_Opt) +#time_lower = @time EAGO_Global_Explicit(SIPopt.LBP_Opt) +#time_lower = @time EAGO_Global_Explicit(SIPopt.LBP_Opt) + +xbar = mid.(ex_SIP_X) +SIPopt.LLP_Opt.DAG_depth = -1 +SIPopt.LLP_Opt.f = p -> -SIP_g(xbar,p) +SIPopt.LLP_Opt.X = ex_SIP_P +#SIPopt.LLP_Opt.LBD_func_relax = "Diff2-MV-OFF" +#SIPopt.LLP_Opt.LBD_problem_relax = "NLP2" +#SIPopt.LLP_Opt.LBD_problem_solver = "Ipopt" +SIPopt.LLP_Opt.LBD_func_relax = "Interval" +SIPopt.LLP_Opt.LBD_problem_relax = "Interval" +SIPopt.LLP_Opt.LBD_problem_solver = "Interval" +SIPopt.LLP_Opt.UBD_func_relax = "Interval" +SIPopt.LLP_Opt.UBD_problem_relax = "Interval" +SIPopt.LLP_Opt.UBD_problem_solver = "Interval" +SIPopt.LLP_Opt.BnB_options.Verbosity = "Normal" +SIPopt.LLP_Opt.BnB_options.BnB_tol = 1E-2 +#INNg1, pbar,feas,tLLP1,tLLP2 = EAGO_Global_Explicit(SIPopt.LLP_Opt) +#time_inner = @time EAGO_Global_Explicit(SIPopt.LLP_Opt) +#time_inner = @time EAGO_Global_Explicit(SIPopt.LLP_Opt) + +SIPopt.UBP_Opt.DAG_depth = -1 +SIPopt.UBP_Opt.f = SIP_f +SIPopt.UBP_Opt.g = SIP_g +SIPopt.UBP_Opt.X = ex_SIP_X +SIPopt.UBP_Opt.g = x -> BndProb_reform(x,[],SIP_g,[],1.0) +#SIPopt.UBP_Opt.LBD_func_relax = "Diff2-MV-OFF" +#SIPopt.UBP_Opt.LBD_problem_relax = "NLP2" +#SIPopt.UBP_Opt.LBD_problem_solver = "Ipopt" +SIPopt.UBP_Opt.LBD_func_relax = "Interval" +SIPopt.UBP_Opt.LBD_problem_relax = "Interval" +SIPopt.UBP_Opt.LBD_problem_solver = "Interval" +SIPopt.UBP_Opt.UBD_func_relax = "Interval" +SIPopt.UBP_Opt.UBD_problem_relax = "Interval" +SIPopt.UBP_Opt.UBD_problem_solver = "Interval" +SIPopt.UBP_Opt.BnB_options.Verbosity = "Normal" +SIPopt.UBP_Opt.BnB_options.BnB_tol = 1E-2 +#UBD_temp, xbar, feas,tUBP1,tUBP2 = EAGO_Global_Explicit(SIPopt.UBP_Opt) +#time_upper = @time EAGO_Global_Explicit(SIPopt.UBP_Opt) +#time_upper = @time EAGO_Global_Explicit(SIPopt.UBP_Opt) +#pbar = [0.66, 0.00275, 1.13595, 0.95, 1.15] +#ubar = [5.5e-05, 0.000100001] × [5.5e-05, 0.000100001] × [9.9781, 9.97815] + +xbnd_range1 = xbnds(ex_SIP_X[1:2],ex_SIP_P) +SIP_range = SIP_g(ex_SIP_X,ex_SIP_P) +ex_out = Explicit_SIP_Solve(SIP_f,[],SIP_g,ex_SIP_X,ex_SIP_P,SIPopt) +=# diff --git a/docs/src/contractors.md b/docs/src/contractors.md new file mode 100644 index 0000000..b6e78c4 --- /dev/null +++ b/docs/src/contractors.md @@ -0,0 +1,5 @@ +# Parametric Contractor Methods + +```@docs +PIn_KrawczykCW +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..5ed2af4 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,23 @@ +# Parametric Interval Methods + +```@contents +``` + +```@docs +PI_NewtonGS +``` + +```@docs +PIn_NewtonGS +``` + +```@docs +PI_KrawczykCW +``` + +```@docs +PIn_KrawczykCW +``` + +```@index +``` diff --git a/docs/src/tests.md b/docs/src/tests.md new file mode 100644 index 0000000..e69de29 diff --git a/docs/src/utilities.md b/docs/src/utilities.md new file mode 100644 index 0000000..e69de29 diff --git a/lib/OldParent.jl b/lib/OldParent.jl deleted file mode 100644 index 8d83775..0000000 --- a/lib/OldParent.jl +++ /dev/null @@ -1,56 +0,0 @@ -#__precompile__() -#= -module EAGOParametricInterval - -using IntervalArithmetic -using EAGODomainReduction - -type Param_Bisect_Opts - DAGflag::Bool - LPflag::Bool - kmax_main::Int64 - kmax_cntr::Int64 - style::String - display::String - ptol::Float64 - etol::Float64 - rtol::Float64 - DAGpass::Int64 - p_rel_bisect::Bool - DAGh - DAGg - DAGgL - DAGgU - DAGsym -end -Param_Bisect_Opts() = Param_Bisect_Opts(false, #DAGflag::Bool - false, #LPflag::Bool - 1E2, #kmax_main - 1E2, #kmax_cntr - "KrawczykCW", #style - "Full", #display - 1E-3, #ptol - 1E-3, #etol - 1E-6, #rtol - 5, #DAGpass - false, # p_rel_bisect - [], - [], - [], - [], - []) - -include("src/Parametric_Utility.jl") -include("src/Parametric_Test.jl") -include("src/Parametric_Contractor_STD.jl") -#include("src/Parametric_Contractor_Inplace.jl") -include("src/Parametric_Bisection.jl") -include("src/Parametric_Main.jl") - -#function __init__() -#end - -export Generalized_Param_Bisection, Param_Bisect_Opts, - GenerateH, setprec, PI_NewtonGS, PI_KrawczykCW -end -=# diff --git a/src/EAGOParametricInterval.jl b/src/EAGOParametricInterval.jl index 2165f56..6e068b1 100644 --- a/src/EAGOParametricInterval.jl +++ b/src/EAGOParametricInterval.jl @@ -45,8 +45,9 @@ export Param_Bisect_Opts,setprec, Miranda, MirandaExc, partialIncTop, partialIncBot, Strict_XinY, isEqual, extDivide, extProcess, =# -export PI_NewtonGS, PI_KrawczykCW, +export PI_NewtonGS, PI_KrawczykCW, PIn_NewtonGS, PIn_KrawczykCW, + PId_NewtonGS, PId_KrawczykCW, Sparse_Precondition!, SparseInSto end diff --git a/src/build/assets/Documenter.css b/src/build/assets/Documenter.css deleted file mode 100644 index f7ae84a..0000000 --- a/src/build/assets/Documenter.css +++ /dev/null @@ -1,18 +0,0 @@ -div.wy-menu-vertical ul.current li.toctree-l3 a { - font-weight: bold; -} - -a.documenter-source { - float: right; -} - -.documenter-methodtable pre { - margin-left: 0px; - margin-right: 0px; - margin-top: 0px; - padding: 0px; -} - -.documenter-methodtable pre.documenter-inline { - display: inline; -} diff --git a/src/build/assets/mathjaxhelper.js b/src/build/assets/mathjaxhelper.js deleted file mode 100644 index 3561b10..0000000 --- a/src/build/assets/mathjaxhelper.js +++ /dev/null @@ -1,25 +0,0 @@ -MathJax.Hub.Config({ - "tex2jax": { - inlineMath: [['$','$'], ['\\(','\\)']], - processEscapes: true - } -}); -MathJax.Hub.Config({ - config: ["MMLorHTML.js"], - jax: [ - "input/TeX", - "output/HTML-CSS", - "output/NativeMML" - ], - extensions: [ - "MathMenu.js", - "MathZoom.js", - "TeX/AMSmath.js", - "TeX/AMSsymbols.js", - "TeX/autobold.js", - "TeX/autoload-all.js" - ] -}); -MathJax.Hub.Config({ - TeX: { equationNumbers: { autoNumber: "AMS" } } -}); diff --git a/src/build/index.md b/src/build/index.md deleted file mode 100644 index ad7f6c9..0000000 --- a/src/build/index.md +++ /dev/null @@ -1,13 +0,0 @@ - - - -# EAGOParametricInterval Documentation - - - - -## Functions - - -'''@docs GenerateJacobianX(ExprArr,SymXArr,SymPArr) ''' - diff --git a/src/docs/build/assets/Documenter.css b/src/docs/build/assets/Documenter.css deleted file mode 100644 index f7ae84a..0000000 --- a/src/docs/build/assets/Documenter.css +++ /dev/null @@ -1,18 +0,0 @@ -div.wy-menu-vertical ul.current li.toctree-l3 a { - font-weight: bold; -} - -a.documenter-source { - float: right; -} - -.documenter-methodtable pre { - margin-left: 0px; - margin-right: 0px; - margin-top: 0px; - padding: 0px; -} - -.documenter-methodtable pre.documenter-inline { - display: inline; -} diff --git a/src/docs/build/assets/mathjaxhelper.js b/src/docs/build/assets/mathjaxhelper.js deleted file mode 100644 index 3561b10..0000000 --- a/src/docs/build/assets/mathjaxhelper.js +++ /dev/null @@ -1,25 +0,0 @@ -MathJax.Hub.Config({ - "tex2jax": { - inlineMath: [['$','$'], ['\\(','\\)']], - processEscapes: true - } -}); -MathJax.Hub.Config({ - config: ["MMLorHTML.js"], - jax: [ - "input/TeX", - "output/HTML-CSS", - "output/NativeMML" - ], - extensions: [ - "MathMenu.js", - "MathZoom.js", - "TeX/AMSmath.js", - "TeX/AMSsymbols.js", - "TeX/autobold.js", - "TeX/autoload-all.js" - ] -}); -MathJax.Hub.Config({ - TeX: { equationNumbers: { autoNumber: "AMS" } } -}); diff --git a/src/docs/build/index.md b/src/docs/build/index.md deleted file mode 100644 index ad7f6c9..0000000 --- a/src/docs/build/index.md +++ /dev/null @@ -1,13 +0,0 @@ - - - -# EAGOParametricInterval Documentation - - - - -## Functions - - -'''@docs GenerateJacobianX(ExprArr,SymXArr,SymPArr) ''' - diff --git a/src/docs/make.jl b/src/docs/make.jl deleted file mode 100644 index 2148037..0000000 --- a/src/docs/make.jl +++ /dev/null @@ -1,9 +0,0 @@ -using Documenter, EAGOParametricInterval - -makedocs(modules=[EAGOParametricInterval], - doctest=true) - -deploydocs(deps = Deps.pip("mkdocs", "python-markdown-math"), - repo = "github.com/GITHUBNAME/GITHUBREPO.git", - julia = "0.6.0", - osname = "windows") diff --git a/src/docs/src/index.md b/src/docs/src/index.md deleted file mode 100644 index adf7d90..0000000 --- a/src/docs/src/index.md +++ /dev/null @@ -1,7 +0,0 @@ -# EAGOParametricInterval Documentation - -## Functions - -'''@docs -GenerateJacobianX(ExprArr,SymXArr,SymPArr) -''' diff --git a/src/lib/Parametric_Contractor.jl b/src/lib/Parametric_Contractor.jl index 3ea24aa..afdabf5 100644 --- a/src/lib/Parametric_Contractor.jl +++ b/src/lib/Parametric_Contractor.jl @@ -383,14 +383,12 @@ function PI_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}}, return X,Eflag,Iflag,inclusionLow,inclusionHigh end - """ - PIn_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}},hj!::Function, - h!::Function,opt::Vector{Any},Eflag::Bool,Iflag::Bool,eDflag::Bool) + PIn_NewtonGS Uses the interval arithmetic presented in ValidatedNumerics to bound unique implicit functions x:P->X defined by a system of equations h(z,p)=0 via -Gauss-Siedel Newton. The jacobian of h w.r.t. z is calculated in placeInputs: +Gauss-Siedel Newton. The jacobian of h w.r.t. z is calculated in place. Inputs: * `X0::Vector{Interval{T}}`: Bounds for dependent variables * `P::Vector{Interval{T}}`: Bounds for independent variables * `hj!::Function`: Jacobian of h(x,p) with respect to x (takes input J,x,p) @@ -602,8 +600,7 @@ function PIn_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}}, end """ - PIn_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}},hj!::Function, - h!::Function,opt::Vector{Any},Eflag::Bool,Iflag::Bool) + PIn_KrawczykCW Uses the interval arithmetic presented in ValidatedNumerics to bound unique implicit functions x:P->X defined by a system of equations h(z,p)=0 via @@ -775,3 +772,375 @@ function PIn_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}}, end return X,Eflag,Iflag,inclusionLow,inclusionHigh end + +""" + PId_NewtonGS + +Uses the interval arithmetic presented in ValidatedNumerics to bound unique +implicit functions x:P->X defined by a system of equations h(z,p)=0 via +Gauss-Siedel Newton. The jacobian of h w.r.t. z is calculated in place. Inputs: +* `X0::Vector{Interval{T}}`: Bounds for dependent variables +* `P::Vector{Interval{T}}`: Bounds for independent variables +* `hj!::Function`: Jacobian of h(x,p) with respect to x (takes input J,x,p) +* `h!::Function`: Equations defining potential implicit function (takes inputs H,x,p) +* `opt::Vector{Any}`: `opt[1]::Int64` is the number of iterations, `opt[2]::Float64` + is the equality tolerance for intervals, `opt[3]::Float64` + is the amount added to M[i,i] when processing extended + interval division. +* `Eflag::Bool`: True if no implicit function exists in box +* `Iflag::Bool`: True if unique implicit function exists in box +* `eDflag::Bool`: True if extended division occured in iteration +Returns a tuple `(X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh)` where +* `X::Vector{Interval{T}}`: Output interval box +* `Xtemp::Vector{Interval{T}}`: Extra box formed from extended division +* `Eflag::Bool`: True if no implicit function exists in box +* `Iflag::Bool`: True if unique implicit function exists in box +* `inclusionLow::Vector{Bool}`: Each component is true if corresponding variable + contracted from lower bound +* `inclusionHigh::Vector{Bool}`: Each component is true if corresponding variable + contracted from upper bound +""" +function PId_NewtonGS(X0::Vector{Interval{T}},P::Vector{Interval{T}}, + hj!::Function,h!::Function, + opt::Vector{Any},Eflag::Bool, + Iflag::Bool,eDflag::Bool) where {T<:AbstractFloat} + # unpacks option file + kmax::Int64 = opt[1] + etol::Float64 = opt[2] + rtol::Float64 = opt[3] + + # Initializes variables + nx::Int64 = length(X0) + S1::Interval{T} = Interval(0.0) + S2::Interval{T} = Interval(0.0) + S3::Interval{T} = Interval(0.0) + inclusion::Vector{Bool} = [false for i=1:nx] + inclusionHigh::Vector{Bool} = [false for i=1:nx] + inclusionLow::Vector{Bool} = [false for i=1:nx] + eD::Int64 = 0 + exclusion::Bool = false + Iflag::Bool = false + Eflag::Bool = false + eDflag::Bool = false + k::Int64 = 1 + X::Vector{Interval{T}} = copy(X0) + Xi::Vector{Interval{T}} = copy(X0) + N::Vector{Interval{T}} = copy(X) + + # creates storage for h(xm,P), J_x(X,P), and preconditioner Y + H::Vector{Interval{T}} = [Interval(0.0) for i=1:nx] + J::Array{Interval{T},2} = zeros(Interval{T},nx,nx) + + # calculates extension of h(xm,P) and J_x(X,P) + x_mid::Vector{T} = mid.(X) + h!(H,x_mid,P) + hj!(J,X,P) + Dense_Precondition!(H,J,mid.(J),nx) + + for i=1:nx + S1 = Interval(0.0) + S2 = Interval(0.0) + for j=1:nx + if (ji) + S2 += J[i,j]*(X[j]-x_mid[j]) + end + end + if J[i,i].lo*J[i,i].hi > 0.0 + N[i] = x_mid[i] - (H[i]+S1+S2)/J[i,i] + else + Ntemp = copy(N) + eD,N[i],Ntemp[i] = extProcess(N[i],X[i],J[i,i],S1,S2,H[i],rtol) + if eD == 1 + eDflag = true + Xtemp = copy(X) + Xtemp[i] = Ntemp[i] ∩ X[i] + X[i] = N[i] ∩ X[i] + return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh + end + end + if Strict_XinY(N[i],X[i]) + inclusion[i] = true + inclusionHigh[i] = true + inclusionLow[i] = true + else + inclusion[i] = false + inclusionLow[i] = false + inclusionHigh[i] = false + if (N[i].lo>X[i].lo) + inclusionLow[i] = true + elseif (N[i].hii) + S2 += J[i,j]*(X[j]-x_mid[j]) + end + end + if J[i,i].lo*J[i,i].hi > 0.0 + N[i] = x_mid[i] - (H[i]+S1+S2)/J[i,i] + else + Ntemp = copy(N) + eD,N[i],Ntemp[i] = extProcess(N[i],X[i],J[i,i],S1,S2,H[i],rtol) + if eD == 1 + eDflag = true + Xtemp = copy(X) + Xtemp[i] = Ntemp[i] ∩ X[i] + X[i] = N[i] ∩ X[i] + return X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh + end + end + if Strict_XinY(N[i],X[i]) + inclusion[i] = true + inclusionHigh[i] = true + inclusionLow[i] = true + else + inclusion[i] = false + inclusionLow[i] = false + inclusionHigh[i] = false + if (N[i].lo>X[i].lo) + inclusionLow[i] = true + elseif (N[i].hiX defined by a system of equations h(z,p)=0 via +a component-wise Krawczyk method. Inputs: +* `X0::Vector{Interval{T}}`: Bounds for dependent variables +* `P::Vector{Interval{T}}`: Bounds for independent variables +* `hj!::Function`: Jacobian of h(x,p) with respect to x (takes input J,x,p) +* `h!::Function`: Equations defining potential implicit function (takes inputs H,x,p) +* `opt::Vector{Any}`: `opt[1]::Int64` is the number of iterations, `opt[2]::Float64` + is the equality tolerance for intervals, `opt[3]::Float64` + is the amount added to M[i,i] when processing extended + interval division. +* `Eflag::Bool`: True if no implicit function exists in box +* `Iflag::Bool`: True if unique implicit function exists in box +Returns a tuple `(X,Xtemp,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh)` where +* `X::Vector{Interval{T}}`: Output interval box +* `Xtemp::Vector{Interval{T}}`: Extra box formed from extended division +* `Eflag::Bool`: True if no implicit function exists in box +* `Iflag::Bool`: True if unique implicit function exists in box +* `inclusionLow::Vector{Bool}`: Each component is true if corresponding variable + contracted from lower bound +* `inclusionHigh::Vector{Bool}`: Each component is true if corresponding variable + contracted from upper bound +""" +function PId_KrawczykCW(X0::Vector{Interval{T}},P::Vector{Interval{T}}, + hj!::Function,h!::Function, + opt::Vector{Any},Eflag::Bool,Iflag::Bool) where {T<:AbstractFloat} + + # unpacks option file + kmax::Int64 = opt[1] + etol::Float64 = opt[2] + + # Initializes variables + nx::Int64 = length(X0) + S::Vector{Interval{T}} = [Interval(0.0) for i=1:nx] + inclusion::Vector{Bool} = [false for i=1:nx] + inclusionHigh::Vector{Bool} = [false for i=1:nx] + inclusionLow::Vector{Bool} = [false for i=1:nx] + exclusion::Bool = false + Iflag::Bool = false + Eflag::Bool = false + eDflag::Bool = false + X::Vector{Interval{T}} = copy(X0) + Xi::Vector{Interval{T}} = copy(X) + N::Vector{Interval{T}} = [Interval(0.0) for i=1:nx] + x_mid::Vector{T} = mid.(X) + k::Int64 = 1 + # creates storage for h(xm,P), J_x(X,P), and preconditioner Y + H::Vector{Interval{T}} = [Interval(0.0) for i=1:nx] + J::Array{Interval{T},2} = zeros(Interval{T},nx,nx) + + # calculates extension of h(xm,P) and J_x(X,P) + h!(H,x_mid,P) + hj!(J,X,P) + Dense_Precondition!(H,J,mid.(J),nx) + + for i=1:nx + S1 = Interval(0.0) + S2 = Interval(0.0) + S3 = Interval(0.0) + for j=1:nx + if (ji) + S2 += -J[i,j]*(X[j]-x_mid[j]) + else + S3 += (1.0-J[i,j])*(X[j]-x_mid[j]) + end + end + N[i] = x_mid[i] - H[i] + S1 + S2 +S3 + + if Strict_XinY(N[i],X[i]) + inclusion[i] = true + inclusionHigh[i] = true + inclusionLow[i] = true + else + inclusion[i] = false + inclusionLow[i] = false + inclusionHigh[i] = false + if (N[i].lo>X[i].lo) + inclusionLow[i] = true + elseif (N[i].hii) + S2 += -J[i,j]*(X[j]-x_mid[j]) + else + S3 += (1.0-J[i,j])*(X[j]-x_mid[j]) + end + end + N[i] = x_mid[i] - H[i] + S1 + S2 +S3 + + if Strict_XinY(N[i],X[i]) + inclusion[i] = true + inclusionHigh[i] = true + inclusionLow[i] = true + else + inclusion[i] = false + inclusionLow[i] = false + inclusionHigh[i] = false + if (N[i].lo>X[i].lo) + inclusionLow[i] = true + elseif (N[i].hi