(* :Title: Interval Plots and Optimization *) (* :Author: Roman E. Maeder *) (* :Summary: 2D and 3D interval plots and global minima *) (* :Context: Intervals` *) (* :Package Version: 1.0 *) (* :Copyright: Copyright 1997, Roman E. Maeder. This package may be used for personal and instructional, noncommercial purposes only. Commercial licensing is available upon request. *) (* :History: Version 1.0 for the Mathematica Journal, July 1997. *) (* :Keywords: interval arithmetic, function plot, optimization, minimization *) (* :Source: Maeder, Roman E. 1997. Interval Plotting and Global Optimisation. The Mathematica Journal, 7(2). *) (* :Mathematica Version:3.0 *) BeginPackage["Intervals`"] IntervalPlot::usage = "IntervalPlot[expr,{x,x0,x1},opts...] plots expr as a function of x." IntervalPlot3D::usage = "IntervalPlot[expr,{x,x0,x1},{y,y0,y1},opts...] plots expr as a function of x and y." IntervalMinimum::usage = "IntervalMinimum[expr,{x,x0,x1},...,opts..] finds the global minimum of expr as a function of the given variables." IntervalMinimumConstrained::usage = "IntervalMinimum[expr,{x,x0,x1},...,constraint, opts..] finds the global minimum of expr as a function of the given variables. Search is restricted to those regions satisfying constraint." RectangleGraph::usage = "RectangleGraph[l] gives a Graphics or Graphics3D of the list of pairs or triplets of intervals." RectangleGraphColor::usage = "RectangleGraphColor[l] gives a Graphics or Graphics3D of the list of pairs or triplets of intervals. Components are colored randomly." EnclosingRectangle::usage = "EnclosingRectangle[l] gives a list of intervals that enclose the given list of lists of intervals." FeasiblePoint::usage = "FeasiblePoint->p is an option of IntervalMinimumConstrained giving an initial feasible point." Options[IntervalPlot] = { Tolerance -> 0.01 } Options[IntervalPlot3D] = { Tolerance -> 0.05 } Options[IntervalMinimum] = { Tolerance -> 0.001 } Options[IntervalMinimumConstrained] = { Tolerance -> 0.001, FeasiblePoint -> Automatic } Begin["`Private`"] Needs["Utilities`FilterOptions`"] Needs["PriorityQueue`"] (* works for intervals or lists of intervals *) length[i_] := Max[i]-Min[i] (* split an interval into two *) split[i_Interval]:= With[{m=(Min[i]+Max[i])/2}, {Interval[{Min[i],m}],Interval[{m,Max[i]}]}] (* split Cartesian products *) split[l_List] := Distribute[split/@l,List,List] (* normalize *) SetAttributes[interval, Listable] interval[i_Interval] := i interval[i_?NumericQ] := Interval[i] interval[_] := Interval[{-Infinity,Infinity}] (* junk *) (* midpoint *) SetAttributes[mid,Listable] mid[i_Interval]:= (Min[i]+Max[i])/2 (* avoid invisible rectangles/cuboids *) rect[l:{xi_,yi_}]/;length[yi]==0 := Line[{Min/@l, Max/@l}] rect[{xi_,yi_,zi_}]/;length[zi]==0 := With[{x0=Min[xi],x1=Max[xi],y0=Min[yi],y1=Max[yi],z=Min[zi]}, Polygon[{{x0,y0,z},{x1,y0,z},{x1,y1,z},{x0,y1,z}}] ] (* infinite boundaries *) infcolor=RGBColor[1,0,0]; huge=10.0^6; rect[{xi__,yi_}]/;Min[yi]==-Infinity||Max[yi]==Infinity := {infcolor, rect[{xi, yi /. DirectedInfinity[d:(-1|1)] :> d huge}]} (* 1D *) rect[l:{x_}] := Line[{{Min[x],0},{Max[x],0}}] (* 2D *) rect[l:{_,_}] := Rectangle[Min/@l, Max/@l] (* 3D *) rect[l:{_,_,_}] := Cuboid[Min/@l, Max/@l] RectangleGraph[l:{{_}..},opts___?OptionQ] := Graphics[{Thickness[0.01],rect/@l}, {opts}] RectangleGraphColor[l:{{_}..},opts___?OptionQ] := Graphics[{Thickness[0.01],{Hue[Random[]],rect[#]}&/@l}, {opts}] RectangleGraph[l:{{_,_}..},opts___?OptionQ] := Graphics[{Thickness[0],rect/@l}, {opts}] RectangleGraphColor[l:{{_,_}..},opts___?OptionQ] := Graphics[{Thickness[0],{Hue[Random[]],rect[#]}&/@l}, {opts}] RectangleGraph[l:{{_,_,_}..,opts___?OptionQ}] := Graphics3D[{EdgeForm[Thickness[0]],rect/@l}, {opts}] EnclosingRectangle[l_List] := Apply[IntervalUnion, Transpose[l], {1}] IntervalPlot[expr_, {x_Symbol,x0_,x1_}, opts___?OptionQ]:= With[{tol = (x1-x0)Tolerance /. {opts} /. Options[IntervalPlot]}, Module[{finals}, finals = refine[Function[x,expr], Interval/@N[{{x0,x1}}], tol]; Show[RectangleGraph[finals], Evaluate[FilterOptions[Graphics,opts]], PlotRange->{{x0,x1},Automatic}, Axes->True] ] ] IntervalPlot3D[expr_,{x_Symbol,x0_,x1_},{y_Symbol,y0_,y1_},opts___?OptionQ]:= With[{tol = (x1-x0)Tolerance /. {opts} /. Options[IntervalPlot3D]}, Module[{finals}, finals = refine[Function[{x,y},expr], Interval/@N[{{x0,x1},{y0,y1}}], tol]; Show[RectangleGraph[finals], Evaluate[FilterOptions[Graphics3D,opts]], PlotRange->{{x0,x1},{y0,y1},Automatic}, Axes->True, BoxRatios->{1,1,0.5}] ] ] refine[f_, init_, tol_] := Module[{new = {init, {}}, finals = h[], i, j}, While[new =!= {}, {i,new} = new; (* dequeue one *) j = interval[f@@i]; If[ length[j] fx, Continue[] ]; (* no longer in the game *) (* split it *) new = split[new[[1]]]; (* update fx, the value of f at a feasible point *) fx = Min[Apply[f, mid/@new, {1}], fx]; (* range intervals *) fs = interval[Apply[f, new, {1}]]; new = Transpose[{new, fs}]; (* stopping criteria *) Scan[ If[ Min[#[[2]]] <= fx, (* requeue *) If[ length[#[[2]]] < tol || length[#[[1,1]]] < tol, fin = h[#, fin], (* no further division *) EnQueue[div, #] (* divide further *) ]]&, new ]; ]; DeleteQueue[div]; fin = List@@Flatten[fin, Infinity, h]; fin = Select[fin, Min[#[[2]]] <= fx&]; (* throw away non-minima *) {new, fs} = Transpose[fin]; {Interval[{Min[fs], fx}], new} ]]] IntervalMinimumConstrained[expr_, ranges:{_Symbol,_,_}.., cond_, opts___?OptionQ] := With[{vars = First/@{ranges},n = Length[{ranges}], tol = Tolerance /. {opts} /. Options[IntervalMinimumConstrained], dom = N[Interval /@ Rest /@ {ranges}]}, With[{f = Function@@Hold[vars, expr], t=Function@@Hold[vars, cond]}, Module[{div = MakeQueue[Min[#1[[2]]] < Min[#2[[2]]]&], fs, fx = Infinity, fin = h[], new, in, out, ok, fp}, fp = FeasiblePoint /. {opts} /. Options[IntervalMinimumConstrained]; If[ t@@fp, fx = f@@fp ]; (* feasible point *) EnQueue[div, {dom, f@@dom}]; While[!EmptyQueue[div], new = DeQueue[div]; (* get smallest one *) If[ Min[new[[2]]] > fx, Continue[] ]; (* no longer in the game *) (* split it *) new = split[new[[1]]]; (* apply constraint *) states = Apply[t, new, {1}]; in = Flatten[Position[states, True, {1}, Heads->False]]; out = Flatten[Position[states, False, {1}, Heads->False]]; ok = Complement[Range[Length[states]], out]; (* not discarded *) (* update fx, the value of f at a feasible point *) fx = Min[Apply[f, mid/@Part[new, in], {1}], fx]; new = Part[new, ok]; (* range intervals *) fs = interval[Apply[f, new, {1}]]; new = Transpose[{new, fs}]; (* stopping criteria *) Scan[ If[ Min[#[[2]]] <= fx, (* requeue *) If[ length[#[[2]]] < tol || length[#[[1,1]]] < tol, fin = h[#, fin], (* no further division *) EnQueue[div, #] (* divide further *) ]]&, new ]; ]; DeleteQueue[div]; fin = List@@Flatten[fin, Infinity, h]; fin = Select[fin, Min[#[[2]]] <= fx&]; (* throw away non-minima *) fin = Select[fin, t@@#[[1]]=!=False&]; (* not outside feasible region *) {new, fs} = Transpose[fin]; {Interval[{Min[fs], fx}], new} ]]] End[] Protect[IntervalPlot, IntervalPlot3D, IntervalMinimum, IntervalMinimumConstrained, RectangleGraph, EnclosingRectangle] EndPackage[]