(* ::Package:: *) (* :Title: InequalityGraphics *) (* :Author: Roger Germundsson, 1998-2000 *) (* :Summary: Provides InequalityPlot and InequalityPlot3D for visualizing logical combinations of polynomial and algebraic inequalities over the reals. And ComplexInequalityPlot for visualizing regions over the complexes. *) (* :Context: Graphics`InequalityGraphics` *) (* :Package Version: 1.1 *) (* :History: V0.9 Roger Germundsson, April 1998 Initial package and examples. V1.0 Roger Germundsson, December 1998 Second implementation and cleanup V1.1 Roger Germundsson, March 2000 Added support for algebraic numbers and functions. Added support for numeric inequalities. Packaged in single file for StandardPackages format. Added support for complex regions. Added message trapping and better filtering of messages. Added new gridding strategy for 3D surface patches. *) (* :Keywords: Real Regions, Complex Regions, Analytic Geometry, Plane Geometry, Solid Geometry, Constructive Plane Geometry, Constructive Solid Geometry, Implicit Regions, Boundary Representation, Polynomial Inequalities, Semi-Algebraic Sets, Graphics, Visualization *) (* :Sources: Roger Germundsson *) (* :Warning: *) (* :Mathematica Version: 4.1 *) (* :Copyright: Copyright 1998-2007, Wolfram Research, Inc.*) (* :Limitation: InequalityPlot, InequalityPlot3D and ComplexInequalityPlot relies on CylindricalDecomposition and as such can deal with polynomial and algebraic inequalities, or semi-algebraic sets. *) Message[General::obspkg, "Graphics`InequalityGraphics`"] Quiet[ BeginPackage[ "Graphics`InequalityGraphics`", {"Utilities`FilterOptions`"} ] , {General::obspkg, General::newpkg}] InequalityPlot::usage = "InequalityPlot[ineqs, {x,xmin,xmax}, {y,ymin,ymax}] plots \ region defined by ineqs within the box bounded by {xmin, ymin} and {xmax, ymax}."; InequalityPlot3D::usage = "InequalityPlot3D[ineqs, {x,xmin,xmax}, {y,ymin,ymax}, \ {z,zmin,zmax}] plots the set described by the inequalities ineqs within the box bounded \ by {xmin,ymin,zmin} and {xmax,ymax,zmax}."; ComplexInequalityPlot::usage = "ComplexInequalityPlot[ineqs, {z, zmin, zmax}] plots \ the the region defined by ineqs within the box bounded by {Re[zmin], Im[zmin]} and \ {Re[zmax], Im[zmax]}. The functions that occur within the inequality need to be real valued \ functions of a complex argument, e.g. Abs, Re and Im."; InequalityPlot::pllim = "The range specification `1` is not of the form {x, xmin, xmax} or {x}."; InequalityPlot3D::pllim = "The range specification `1` is not of the form {x, xmin, xmax} or {x}."; ComplexInequalityPlot::pllim = "The range specification `1` is not of the form {z, zmin, zmax} or {z}."; InequalityPlot::region = "The region defined by `1` could not be broken down into cylinders."; InequalityPlot3D::region = "The region defined by `1` could not be broken down into cylinders."; ComplexInequalityPlot::region = "The region defined by `1` could not be broken down into cylinders."; InequalityPlot::unbounded = "The region defined by `1` is unbounded. Try adding bounds for your variables."; InequalityPlot3D::unbounded = "The region defined by `1` is unbounded. Try adding bounds for your variables."; ComplexInequalityPlot::unbounded = "The region defined by `1` is unbounded. Try adding bounds for your variable."; ComplexInequalityPlot::notreal = "The region defined by `1` does not appear to be a logical combination of real valued functions."; InequalityPlot3D::method = "The value of Method -> `1`, should be one of Automatic or \ \"FixedSampleCountGridding\"."; Begin[ "`Private`" ]; issueObsoleteFunMessage[fun_, context_] := Message[General::obspkgfn, fun, context]; Unprotect[ InequalityPlot, InequalityPlot3D, ComplexInequalityPlot ]; (* InequalityPlot[ ineqs, xrange, yrange, opts ] * Input: * ineqs -- Logic combination of equalities and inequalities. * xrange -- {x, xmin, xmax} range for x-axis * yrange -- {y, ymin, ymax} range for y-axis * Output: * Graphics representing the full dimensional components of ineqs. *) Options[ InequalityPlot ] = Join[ Options[Plot], {BoundaryStyle -> Black} ]; SetOptions[ InequalityPlot, Axes -> True, AspectRatio -> Automatic, Filling -> {1 -> {2}}, FillingStyle -> Cyan, PlotRange -> All, PlotStyle -> Black, Method -> {"AxesInFront" -> True}]; Fix2DRange[ {x_} ] := {x, -Infinity, Infinity}; Fix2DRange[ {x_, xmin_, xmax_} ] := {x, xmin, xmax}; Fix2DRange[ xr_ ] := ( Message[ InequalityPlot::pllim, xr ]; xr ); InequalityPlot[ ineqs_List, rest___ ] := InequalityPlot[ Apply[ And, ineqs ], rest ]; InequalityPlot[ ineqs_, xr_List, yr_List, opts___?OptionQ ] := (issueObsoleteFunMessage[InequalityPlot,"Graphics`InequalityGraphics`"]; plot /; FreeQ[ plot = inequalityPlot[ ineqs, {xr, yr}, opts ], $Failed ]); inequalityPlot[ ineqs_, plimits_List, opts___?OptionQ ] := Module[{plims, nineqs, vars, nests, plots, cad, bstyle}, bstyle = BoundaryStyle /. {opts} /. Options[ InequalityPlot ]; If[!ListQ[bstyle], bstyle = {bstyle}]; Check[ plims = Map[ Fix2DRange, plimits ], Return[ $Failed ] ]; vars = Map[ First, plims ]; nineqs = ineqs && FromNestedLists[ {plims} ]; Check[ Block[ {$CADMsgTrap=True}, cad = InequalityDecompose[ nineqs, vars ]], ( Message[ InequalityPlot::region, nineqs ]; Return[ $Failed ] ), CylindricalDecomposition::nrpti, GenericCylindricalDecomposition::nrtpi ]; If[ cad === False, (* empty region *) Return[ Show[ Graphics[ {}, FilterOptions[ Graphics, opts, Sequence @@ Options[InequalityPlot]] ] ] ] ]; nests = ToNestedLists[ cad, vars ]; If[ !FreeQ[ nests, _DirectedInfinity | Infinity, Infinity ], ( Message[ InequalityPlot::unbounded, nineqs ]; Return[ $Failed ] ) ]; plots = Map[ NestedPlot[#, DisplayFunction->Identity, PlotStyle->bstyle, FilterOptions[ Plot, opts, Sequence @@ Options[InequalityPlot] ] ]&, nests ]; Show[ plots, FilterOptions[ Graphics, opts, Sequence @@ Options[InequalityPlot] ] ] ]; NestedPlot[ {{x1_, l1_, u1_}, {x2_, l2_, u2_}}, opts___ ] := Plot[ Evaluate[{l2, u2}], {x1, l1, u1}, Evaluate[FilterOptions[ Plot, opts ]], Filling -> {1 -> {2}}, FillingStyle -> Cyan, PlotRange -> All, Method -> {"AxesInFront" -> True}]; (* ComplexInequalityPlot[ ineqs, zrange, opts ] * Input: * ineqs -- Logic combination of equalities and inequalities. * zrange -- {x, zmin, zmax} defines the bounding box: * {Re[zmin],Im[zmin]} through {Re[zmax], Im[zmax]} * Output: * Graphics representing the full dimensional components of ineqs. * Notes: * In order to generate messages associated with ComplexInequalityPlot * rather then InequalityPlot there is some code duplication. This also * gives us the freedom to diagnose slightly differently since the likely * errors that come up are somewhat different. *) Options[ ComplexInequalityPlot ] = Join[ Options[Plot], {BoundaryStyle -> {}} ]; SetOptions[ ComplexInequalityPlot, Axes -> True, AspectRatio -> Automatic, Filling -> {1 -> {2}}, FillingStyle -> Cyan, PlotRange -> All, PlotStyle -> Black, Method -> {"AxesInFront" -> True}]; FixComplexRange[ {x_} ] := {x, {Re[x], -Infinity, Infinity}, {Im[x], -Infinity, Infinity}}; FixComplexRange[ {x_, xmin_, xmax_} ] := {x, {Re[x], Re[xmin], Re[xmax]}, {Im[x], Im[xmin], Im[xmax]}}; FixComplexRange[ xr_ ] := ( Message[ ComplexInequalityPlot::pllim, xr ]; xr ); ComplexInequalityPlot[ ineqs_List, rest___ ] := (issueObsoleteFunMessage[ComplexInequalityPlot,"Graphics`InequalityGraphics`"]; InequalityPlot[ Apply[ And, ineqs ], rest ]); ComplexInequalityPlot[ ineqs_, zr_List, opts___?OptionQ ] := (issueObsoleteFunMessage[ComplexInequalityPlot,"Graphics`InequalityGraphics`"]; plot /; FreeQ[ plot = complexInequalityPlot[ ineqs, zr, opts ], $Failed ]); complexInequalityPlot[ ineqs_, zr_List, opts___?OptionQ ] := Module[{plims, nineqs, vars, nests, plots, cad, x, y, xr, yr, z, zmin, zmax, cineqs, bstyle}, bstyle = BoundaryStyle /. {opts} /. Options[ ComplexInequalityPlot ]; If[!ListQ[bstyle], bstyle = {bstyle}]; Check[ plims = FixComplexRange[ zr ], Return[ $Failed ] ]; {z, xr, yr} = plims; vars = {x, y}; cineqs = ineqs && FromNestedLists[ { {xr, yr} } ]; nineqs = ComplexToRealRegion[ cineqs, z, {x, y} ]; If[ !FreeQ[ nineqs, Complex[0,1]| _Re| _Im| _Conjugate, Infinity ], (* I'm guessing this is a typical mistake *) Message[ ComplexInequalityPlot::notreal, cineqs ]; Return[ $Failed ] ]; Check[ Block[ {$CADMsgTrap=True}, cad = InequalityDecompose[ nineqs, vars ]], ( Message[ ComplexInequalityPlot::region, cineqs ]; (* Complex inequalities used here *) Return[ $Failed ] ), CylindricalDecomposition::nrpti, GenericCylindricalDecomposition::nrtpi ]; If[ cad === False, (* empty region *) Return[ Show[ Graphics[ {}, FilterOptions[ Graphics, opts, Sequence @@ Options[InequalityPlot]] ] ] ] ]; nests = ToNestedLists[ SelectCylinders[ cad, vars, 2 ], vars ]; If[ !FreeQ[ nests, _DirectedInfinity|Infinity, Infinity ], (* unbounded region *) ( Message[ ComplexInequalityPlot::unbounded, cineqs ]; Return[ $Failed ] ) ]; plots = Show[ Map[ NestedPlot[#, DisplayFunction->Identity, FilterOptions[ Plot, opts, Sequence @@ Options[InequalityPlot] ] ]&, nests ], DisplayFunction -> Identity ]; Show[ plots /. l_Line :> Append[ bstyle, l ], FilterOptions[ Graphics, opts, Sequence @@ Options[InequalityPlot] ] ] ]; (* ComplexToRealRegion[ ineqs, z, {x,y} ] * Input: * ineqs -- Logical combination of inequalities involving polynomials, * rational and algebraic functions as well as Abs, Re and Im. * z -- Complex variable. * {x,y} -- Real variables to use: z = x + I y * Output: * A logical combination of inequalities involving only x and y. * Notes: This works somewhat more reliably then say doing only ComplexExpand, e.g. * * ComplexExpand[Inequality[Abs[z], LessEqual, Re[z^2], Less, 2] /. z -> x + I*y, TargetFunctions -> {Re, Im}] * ==> Inequality[Sqrt[x^2 + y^2], LessEqual, Re[(x + I*y)^2], Less, 2] * * whereas * * ComplexToRealRegion[Inequality[Abs[z], LessEqual, Re[z^2], Less, 2], z, {x, y}] * ==> Inequality[Sqrt[x^2 + y^2], LessEqual, x^2 - y^2, Less, 2] *) ComplexToRealRegion[ ineq_, z_, {x_, y_} ] := Module[ {ce, nineq, i, f, g, less, lessEqual, greater, greaterEqual, equal, unequal}, ce = ComplexExpand[#, TargetFunctions -> {Re, Im}] &; nineq = ReplaceRepeated[ ineq /. z -> x + \[ImaginaryI] y, { i_Inequality :> Apply[ Inequality, Table[ If[ OddQ[j], ce[i[[j]]], i[[j]]], {j, Length[i]}]], i_Less :> less @@ ce /@ i, i_LessEqual :> lessEqual @@ ce /@ i, i_Greater :> greater @@ ce /@ i, i_GreaterEqual :> greaterEqual @@ ce /@ i, i_Equal :> equal @@ ce /@ i, i_Unequal :> unequal @@ ce /@ i } ]; mineq = nineq /. {less -> Less, lessEqual -> LessEqual, greater -> Greater, greaterEqual -> GreaterEqual, equal -> Equal, unequal -> Unequal}; mineq /. { Re[x]->x, Im[x]->x, Re[y]->y, Im[y]->y, Abs[x]->x, Abs[y]->y, Conjugate[x]->x, Conjugate[y]->y } (* Postprocessing for some bugs in ComplexExpand *) ]; (* InequalityPlot3D[ ineqs, xrange, yrange, zrange, opts ] * Input: * ineqs -- Logic combination of equalities and inequalities. * xrange -- {x, xmin, xmax} range for x-axis * yrange -- {y, ymin, ymax} range for y-axis * zrange -- {z, zmin, zmax} range for z-axis * Output: * Graphics3D representing the full dimensional components of ineqs. * Notes: * o InequalityPlot3D currently implement a cylindrical algebraic decomposition * *) Options[ InequalityPlot3D ] = Join[ Developer`Graphics3DOptions[], {PlotPoints -> 8, BoundaryStyle -> {}} ]; SetOptions[ InequalityPlot3D, Axes -> True, Method -> Automatic ]; Fix3DRange[ {x_} ] := {x, -Infinity, Infinity}; Fix3DRange[ {x_, xmin_, xmax_} ] := {x, xmin, xmax}; Fix3DRange[ xr_ ] := ( Message[ InequalityPlot3D::pllim, xr ]; xr ) InequalityPlot3D[ ineqs_List, rest___ ] := InequalityPlot3D[ Apply[ And, ineqs ], rest ]; InequalityPlot3D[ ineqs_, xr_List, yr_List, zr_List, opts___?OptionQ ] := (issueObsoleteFunMessage[InequalityPlot3D,"Graphics`InequalityGraphics`"]; plot /; FreeQ[ plot = inequalityPlot3D[ ineqs, {xr, yr, zr}, opts ], $Failed ]); inequalityPlot3D[ ineqs_, plimits_List, opts___?OptionQ ] := Module[{plims, vars, nineqs, nests, plots, cad, bstyle, method, nxy, nx, ny, dx, dy, xmin, xmax, ymin, ymax}, bstyle = BoundaryStyle /. {opts} /. Options[ InequalityPlot3D ]; If[!ListQ[bstyle], bstyle = {bstyle}]; method = Method /. {opts} /. Options[ InequalityPlot3D ]; Check[ plims = Map[ Fix3DRange, plimits ], Return[ $Failed ] ]; vars = Map[ First, plims ]; nineqs = ineqs && FromNestedLists[ {plims} ]; Check[ Block[ {$CADMsgTrap=True}, cad = InequalityDecompose[ nineqs, vars ]], ( Message[ InequalityPlot3D::region, nineqs ]; (* non semi-algebraic region *) Return[ $Failed ] ), CylindricalDecomposition::nrpti, GenericCylindricalDecomposition::nrtpi ]; If[ cad === False, (* empty region *) Return[ Show[ Graphics3D[ {}, FilterOptions[ Graphics3D, opts, Sequence @@ Options[InequalityPlot3D]] ] ] ] ]; nests = ToNestedLists[ SelectCylinders[ cad, vars, 3 ], vars ]; If[ !FreeQ[ nests, _DirectedInfinity|Infinity, Infinity ], (* unbounded region *) ( Message[ InequalityPlot3D::unbounded, nineqs ]; Return[ $Failed ] ) ]; nxy = PlotPoints /. {opts} /. Options[ InequalityPlot3D ]; If[ ListQ[nxy] && Length[nxy] == 2, {nx, ny} = nxy, nx = ny = nxy ]; Switch[ method, "FixedSampleCountGridding" | Automatic, (* fastest and most robust gridding *) ( plots = Show[ Map[ NestedPlot3D[#, {nx, ny}, DisplayFunction->Identity, FilterOptions[ NestedPlot3D, opts, Sequence @@ Options[ InequalityPlot3D ] ] ]&, nests ], DisplayFunction -> Identity ]; ), "FixedResolutionGridding", (* fixed resolotion, global grid, sub voxel sampling - undocumented *) ( {{xmin, xmax}, {ymin, ymax}} = BoundingBox[nineqs, vars]; {dx, dy} = {(xmax-xmin)/nx, (ymax-ymin)/ny}; plots = Show[ Map[ NestedPlot3D3[#, {dx, dy}, {xmin, ymin}, DisplayFunction->Identity, FilterOptions[ NestedPlot3D, opts, Sequence @@ Options[ InequalityPlot3D ] ] ]&, nests ], DisplayFunction -> Identity ]; ), "FixedResolutionGridding1", (* fixed resolution, no global grid - undocumented *) ( {{xmin, xmax}, {ymin, ymax}} = BoundingBox[nineqs, vars]; {dx, dy} = {(xmax-xmin)/nx, (ymax-ymin)/ny}; plots = Show[ Map[ NestedPlot3D1[#, {dx, dy}, DisplayFunction->Identity, FilterOptions[ NestedPlot3D, opts, Sequence @@ Options[ InequalityPlot3D ] ] ]&, nests ], DisplayFunction -> Identity ]; ), "FixedResolutionGridding2", (* approximately fixed resolution, no global grid - undocumented *) ( {{xmin, xmax}, {ymin, ymax}} = BoundingBox[nineqs, vars]; {dx, dy} = {(xmax-xmin)/nx, (ymax-ymin)/ny}; plots = Show[ Map[ NestedPlot3D2[#, {dx, dy}, DisplayFunction->Identity, FilterOptions[ NestedPlot3D, opts, Sequence @@ Options[ InequalityPlot3D ] ] ]&, nests ], DisplayFunction -> Identity ]; ), _, Message[ InequalityPlot3D::method, method ]; Return[ $Failed ] ]; Show[ plots /. {p__Polygon} :> Join[ bstyle, {p} ], FilterOptions[ Graphics3D, opts, Sequence @@ Options[ InequalityPlot3D ] ] ] ]; (*----------------------- Gridding --------------------------------*) (* Supremum[ obj_, ineqs_, vars_ ] := -Experimental`Infimum[-obj, ineqs, vars ]; *) (* BoundingBox[ ineqs, {x,y,z}] * Input: * ineqs -- Set of inequalities. * {x,y,z} -- Variables * Output: * {{xmin, xmax}, {ymin, ymax}} -- Bounding box for set. * Notes: * o This may be significantly less expensive with generic QE. * o Alternatively we might make it compulsive in the input to * specify bounding box in which case the following is redundant. * o This is used for fixed spatial resolution sampling. *) BoundingBox[ ineqs_, {x_, y_, z_} ] := Module[ {xineqs, yineqs, xmin, xmax, ymin, ymax}, xineqs = Resolve[ Exists[ {y, z}, Element[{x, y, z}, Reals ], ineqs ] ] /. Element[_,Reals] :> True; xmin = Experimental`Infimum[ {x, xineqs}, {x}]; xmax = Experimental`Supremum[ {x, xineqs}, {x}]; yineqs = Resolve[ Exists[ {x, z}, Element[{x, y, z}, Reals ], ineqs ] ] /. Element[_,Reals] :> True; ymin = Experimental`Infimum[ {y, yineqs}, {y}]; ymax = Experimental`Supremum[ {y, yineqs}, {y}]; {{xmin, xmax}, {ymin, ymax}} ]; (* NestedPlot3D[ nest, opts ] * Input: * nest -- Nested representation of volume see pattern below. * Output: * Graphics3D object with graphics primitives to visualize volume. * Notes: * o NestedPlot3D implements a gridding strategy with a fixed number * of sample points nx and ny in the x and y direction respectively. * This means that the spatial resolution will vary, but not the number * of samples. * o The surface patches are guaranteed to fit together without gaps, * but the resulting meshes usually do not. * o There is no global information available to NestedPlot3D. *) Options[ NestedPlot3D ] = {PlotPoints->8, PolygonFunction -> Polygon}; NestedPlot3D[ {{x1_, l1_, u1_}, {x2_, l2_, u2_}, {x3_, l3_, u3_}}, {nx_, ny_}, opts___?OptionQ ] := Module[ {hi, ghi, lo, glo, pfun}, pfun= PolygonFunction /. {opts} /. Options[ NestedPlot3D ]; hi = uniformSampleTable[ {x1, x2, N[u3]}, {x1, N[l1], N[u1], nx}, {x2, N[l2], N[u2], ny}]; ghi = generatePolygons[ hi, pfun, -1 ]; lo = uniformSampleTable[ {x1, x2, N[l3]}, {x1, N[l1], N[u1], nx}, {x2, N[l2], N[u2], ny}]; glo = generatePolygons[ lo, pfun, 1 ]; Show[ {ghi, glo} , FilterOptions[ Graphics3D, opts ] ] ]; (* uniformSampleTable[ expr, xrange, yrange ] * Input: * expr -- Expression to sample * xrange -- {x, lx, ux, xp} sampling description in x direction * yrange -- {y, ly[x], uy[x], yp[x]} sampling description in y direction * Output: * A matrix of sample values. The way this is called the sample values * will be of the form: {x, y, expr[x,y]}, i.e. a matrix of triplets. * * Notes: * -- The input to this function is guaranteed to satisfy: * lx < ux and ly[lx] <= uy[lx] and ly[ux] <= uy[ux] * -- The input specification reads * {x, lx, ux, xp} * x variable * lx lower value for x * ux upper value for x * xp number of sample points * -- The input is guaranteed to be real for all sample points. However using * machine precision evaluation there are cases (typically boundary points) * where this may be violated. We take the real part in those cases. *) uniformSampleTable[ expr_, {x_, lx_, ux_, xp_}, {y_, ly_, uy_, yp_} ] := If[ lx == ux, {{}}, Re[ Table[ If[ Re[ly] == Re[uy], (* First case is when we have end points meet -- padding *) Table[ Evaluate[ReplaceAll[expr, y -> Re[ly]]], {yp + 1} ], Table[ expr, {y, Re[ly], Re[uy], Re[(uy - ly)/yp]} ] ], {x, Re[lx], Re[ux], Re[(ux - lx)/xp]} ] ] ]; (* NestedPlot3D1[ nest, {dx, dy}, opts ] * Input: * nest -- Nested representation of volume see pattern below. * {dx,dy} -- Spatial sampling resolution in the x and y direction. * Output: * Graphics3D object with graphics primitives to visualize volume. * Notes: * o NestedPlot3D1 implements a gridding strategy that samples * with a fixed spatial resolution of dx and dy in the x and y * directions respectively. * o It isn't optimal in that there may be gaps between surface * patches. * o There is no global information available to NestedPlot3D1. *) NestedPlot3D1[ {{x1_, l1_, u1_}, {x2_, l2_, u2_}, {x3_, l3_, u3_}}, {dx_, dy_}, opts___?OptionQ ] := Module[ {nxy, nx, ny, hi, ghi, lo, glo, pfun}, pfun= PolygonFunction /. {opts} /. Options[ NestedPlot3D ]; hi = uniformSampleTable1[ {x1, x2, N[u3]}, {x1, N[l1], N[u1]}, {x2, N[l2], N[u2]}, {N[dx], N[dy]}]; ghi = generatePolygons[ hi, pfun, -1 ]; lo = uniformSampleTable1[ {x1, x2, N[l3]}, {x1, N[l1], N[u1]}, {x2, N[l2], N[u2]}, {N[dx], N[dy]}]; glo = generatePolygons[ lo, pfun, 1 ]; Show[ {ghi, glo} , FilterOptions[ Graphics3D, opts ] ] ]; uniformSampleTable1[ expr_, {x_, lx_, ux_}, {y_, ly_, uy_}, {dx_, dy_} ] := RightPadding[ (* definitely fixed resolution, but leaves gaps *) Table[ Re[expr], {x, lx, ux, dx}, {y, ly, uy, dy} ] ]; RightPadding[ matrix_ ] := Module[ {ny = Max[ Map[ Length, matrix ] ]}, Map[ PadRight[#, ny, {Last[#]}]&, matrix ] ]; (* NestedPlot3D2[ nest, {dx, dy}, opts ] * Input: * nest -- Nested representation of volume see pattern below. * {dx,dy} -- Spatial sampling resolution in the x and y direction. * Output: * Graphics3D object with graphics primitives to visualize volume. * Notes: * o NestedPlot3D2 implements a gridding strategy that samples * with an approximate fixed spatial resolution of dx and dy in the x and y * directions respectively. * o It isn't optimal in that the mesh lines for a multi cylinder region * will not align (this is the approximation). * o There is no global information available to NestedPlot3D1. *) NestedPlot3D2[ {{x1_, l1_, u1_}, {x2_, l2_, u2_}, {x3_, l3_, u3_}}, {dx_, dy_}, opts___?OptionQ ] := Module[ {nxy, nx, ny, hi, ghi, lo, glo, pfun}, pfun= PolygonFunction /. {opts} /. Options[ NestedPlot3D ]; hi = uniformSampleTable2[ {x1, x2, N[u3]}, {x1, N[l1], N[u1]}, {x2, N[l2], N[u2]}, {N[dx], N[dy]}]; ghi = generatePolygons[ hi, pfun, -1 ]; lo = uniformSampleTable2[ {x1, x2, N[l3]}, {x1, N[l1], N[u1]}, {x2, N[l2], N[u2]}, {N[dx], N[dy]}]; glo = generatePolygons[ lo, pfun, 1 ]; Show[ {ghi, glo} , FilterOptions[ Graphics3D, opts ] ] ]; uniformSampleTable2[ expr_, {x_, lx_, ux_}, {y_, ly_, uy_}, {dx_, dy_} ] := Module[ {nx, ny}, nx = Ceiling[ (ux-lx)/dx ]; (* approximately the right x resolution *) ny = Ceiling[ Max[ Table[ Re[uy - ly], {x, lx, ux, (ux-lx)/Ceiling[nx]} ]]/dy ]; (* approximately the right minimal y resolution *) Re[ Table[ If[ ly == uy, (* First case is when we have end points meet -- padding *) Table[ Evaluate[ReplaceAll[expr, y -> ly]], {ny+1} ], Table[ expr, {y, ly, uy, (uy - ly)/ny} ] ], {x, lx, ux, (ux - lx)/nx} ] ] ]; (* NestedPlot3D3[ nest, {dx, dy}, {xmin, ymin}, opts ] * Input: * nest -- Nested representation of volume see pattern below. * {dx,dy} -- Spatial sampling resolution in the x and y direction. * {xmin, ymin} -- Lower left corner of global bounding box. * Output: * Graphics3D object with graphics primitives to visualize volume. * * Notes: * o NestedPlot3D3 implements a gridding strategy that samples * with a fixed spatial resolution of dx and dy in the x and y * directions respectively. * o At or near boundaries sampling below the fixed spatial * resolution is done. This sampling guarantees that surface * patches fit together. * o The global information {xmin, ymin} allows NestedPlot3D3 * to conform to a global sampling grid (except near the boundary * where additional samples may be taken. *) NestedPlot3D3[ {{x1_, l1_, u1_}, {x2_, l2_, u2_}, {x3_, l3_, u3_}}, {dx_, dy_}, {glx_, gly_}, opts___?OptionQ ] := Module[ {nxy, nx, ny, hi, ghi, lo, glo, pfun}, pfun= PolygonFunction /. {opts} /. Options[ NestedPlot3D ]; hi = UniformSpatialSampling[ {x1, x2, Re[N[u3]]}, {x1, N[l1], N[u1], N[dx], N[glx]}, {x2, Re[N[l2]], Re[N[u2]], N[dy], N[gly]}]; ghi = generatePolygons[ hi, pfun, -1 ]; lo = UniformSpatialSampling[ {x1, x2, Re[N[l3]]}, {x1, N[l1], N[u1], N[dx], N[glx]}, {x2, Re[N[l2]], Re[N[u2]], N[dy], N[gly]}]; glo = generatePolygons[ lo, pfun, 1 ]; Show[ {ghi, glo} , FilterOptions[ Graphics3D, opts ] ] ]; PadCenter[ list_, n_, {left_, right_} ] := PadRight[ PadLeft[ list, Floor[(n - Length[list])/2] + Length[list], left], n, right ]; CenterPadding[ matrix_ ] := Module[ {ny = Max[ Map[ Length, matrix ] ]}, Map[ PadCenter[#,ny,{{First[#]}, {Last[#]}}]&, matrix ] ]; uniformSampleTable3[ expr_, {x_, lx_, ux_, dx_, glx_}, {y_, ly_, uy_, dy_, gly_} ] := CenterPadding[ boundaryTable[ (* fixed resolution following global grid, and sub resolution sampling at boundary *) boundaryTable[ Re[expr], {y, Re[ly], Re[uy], dy, gly} ], {x, lx, ux, dx, glx} ] ]; SetAttributes[ boundaryTable, HoldAll ]; boundaryTable[ expr_, {x_, lx_, ux_, dx_, glx_} ] := Module[ {nlx, nux}, (* The minimum and maximum value fitting withing the global grid *) {nlx, nux} = {glx + (Quotient[ lx - glx, dx ]+1)dx, glx + Quotient[ ux - glx, dx ] dx }; Join[ { ReleaseHold[ ReplaceAll[ Hold[ expr ], x->lx ] ] }, Table[ Re[expr], {Evaluate[x], Evaluate[nlx], Evaluate[nux], dx} ], { ReleaseHold[ ReplaceAll[ Hold[ expr ], x->ux ] ] } ] ]; $GridEpsilon = 0.05; (* Keep this settable for now *) $GridDebug = False; UniformSpatialSampling[ expr_, {x_, lx_, ux_, dx_, xc_}, {y_, ly_, uy_, dy_, yc_} ] := Module[ {in, lxb, uxb, lyb, uyb, int, bound, fbound, ymins, ymin, lpads, ymaxs, ymax, rpads, pos, invariant, sample}, in = (* interior points *) Table[ expr, {x, xc + (Quotient[lx - xc, dx] + 1)dx, xc + Quotient[ux - xc, dx] dx, dx}, {y, yc + (Quotient[ly - yc, dy] + 1)dy, yc + Quotient[uy - yc, dy]dy, dy} ]; lxb = (* lower x boundary *) Table[ Evaluate[ expr /. x -> lx], {y, yc + (Quotient[ly - yc, dy] + 1)dy /. x -> lx, yc + Quotient[uy - yc, dy]dy /. x -> lx , dy}]; uxb = (* upper x boundary *) Table[ Evaluate[ expr /. x -> ux], {y, yc + (Quotient[ly - yc, dy] + 1)dy /. x -> ux, yc + Quotient[uy - yc, dy]dy /. x -> ux , dy}]; set = Join[ {lxb}, in, {uxb} ] ; lyb = Join[ (* lower y boundary *) { expr /. y -> ly /. x -> lx }, Table[ Evaluate[ expr /. y -> ly ], {x, xc + (Quotient[lx - xc, dx] + 1)dx, xc + Quotient[ux - xc, dx] dx, dx} ], { expr /. y -> ly /. x -> ux } ]; uyb = Join[ (* upper y boundary *) { expr /. y -> uy /. x -> lx }, Table[ Evaluate[ expr /. y -> uy ], {x, xc + (Quotient[lx - xc, dx] + 1)dx, xc + Quotient[ux - xc, dx] dx, dx} ], { expr /. y -> uy /. x -> ux } ]; bound = MapThread[ Join, { Transpose[{lyb}], set, Transpose[{uyb}] }, 1 ]; (* Pruning points whose y-values are too close in the dy relative measure *) pruneStart[ list_ ] := If[ Abs[ list[[1, 2]] - list[[2, 2]] ] < 2 $GridEpsilon dy, Drop[ list, 1 ], list ]; pruneEnd[ list_ ] := If[ Abs[ list[[-1, 2]] - list[[-2, 2]] ] < 2 $GridEpsilon dy, Drop[ list, -1 ], list ]; fbound = Map[ If[ Length[ # ] == 2, pruneStart[#], pruneEnd[ pruneStart[ # ] ] ] &, bound ]; (* Directional fuzzy rounding, note the critical dependence on $GridEpsilon. To get the actual grid positions in a the global {xmin,ymin} based grid with resolution {dx,dy} *) ymins = Map[ First, Map[ Part[#,2]&, fbound, {2} ]]; ymin = Floor[ Min[ (ymins - yc)/dy ] + $GridEpsilon ]; (* global grid values *) lpads = Map[ Floor[ (# - yc)/dy + $GridEpsilon ] &, ymins ] - ymin; ymaxs = Map[ Last, Map[ Part[#,2]&, fbound, {2} ]]; ymax = Ceiling[ Max[ (ymaxs - yc)/dy ] - $GridEpsilon ]; (* global grid values *) rpads = ymax - Map[ Ceiling[ (# - yc)/dy - $GridEpsilon ] &, ymaxs ]; (* Special case of only one y-value in a row. The directional rounding will not work here. *) Do[ If[ Length[ fbound[[i]] ] == 1, pos = Round[ (fbound[[i,1,2]] - yc)/dy ]; lpads[[i]] = pos - ymin; rpads[[i]] = ymax - pos; ], {i, Length[fbound]} ]; If[ $GridDebug === True, (* An assert check ... off by default. *) invariant = Apply[ Plus, Transpose[{lpads, Map[Length, fbound], rpads }], {1}]; If[ Length[ Union[ invariant ] ] != 1, Print[ "UniformSpatialSampling::Assert\n", "Input: ", HoldForm[ UniformSpatialSampling[ expr, {x, lx, ux, dx, xc}, {y, ly, uy, dy, yc}]], "\n Invariant: ", invariant ]; ]; ]; sample = Table[ Join[ Table[ First[ fbound[[i]] ], {lpads[[i]]}], (* left padding *) fbound[[i]], Table[ Last[ fbound[[i]] ], {rpads[[i]]}] (* right padding *) ], {i, Length[fbound]} ]; (* This will succeed, even if it may not be as optimally aestheatically pleasing *) sample = Map[ Take[ #, Min[ Map[ Length, sample ] ] ]&, sample ]; If[ $GridDebug === True, Polygonize[ sample, LinePolygon ] ]; sample ]; Polygonize[ surf3d_ , pfun_] := Show[ Graphics[ Flatten[ Map[ pfun[ Part[Flatten[#, 1], {1, 2, 4, 3}]] &, Partition[ Take[surf3d, All, All, 2], {2, 2}, {1, 1}], {2} ] ], AspectRatio -> Automatic, PlotRange -> All ] ]; LinePolygon[ {x1_, x2_, x3_, x4_} ] := Line[ {x1, x2, x3, x4, x1}]; (*----------------------- Polygonization ----------------------------*) (* generatePolygons[ surf3d, pfun ] * Input: * surf3d -- A matrix of sample points {xi, yi, zi} as generated by uniformsSampleTable * pfun -- A function to apply to quadrangles {p1, p2, p3, p4} to generate polygons * Output: * Graphics3D object with primitives making up the visualization of surf3d *) generatePolygons[ surf3d_, pfun_ ] := generatePolygons[ surf3d, pfun, 1 ]; (* Orientation: Standard *) generatePolygons[ surf3d_, pfun_, 1 ] := Graphics3D[ Flatten[ Map[ pfun[ Part[ Flatten[#,1], {1,2,4,3}]]&, Partition[ surf3d, {2,2}, {1, 1} ], {2} ] ] ]; (* Orientation: Reversed *) generatePolygons[ surf3d_, pfun_, -1 ] := Graphics3D[ Flatten[ Map[ pfun[ Part[ Flatten[#,1], {1,3,4,2}]]&, Partition[ surf3d, {2,2}, {1, 1} ], {2} ] ] ]; (*-------------------- Region Decomposition -------------------------*) (* InequalityDecompose[ ineq, vars ] * * Will decompose a set of inequalities. It will return the * full dimensional (of dimension == Length[vars]) only. * This indirection will be useful in the future when different * functionality will be used for the decomposition step. *) InequalityDecompose[ ineqs_, vars_ ] := Module[ {oineqs}, (* Optimization1: Use strict inequalities *) oineqs = ineqs /. {LessEqual -> Less, GreaterEqual -> Greater}; (* Optimization2: Numericalize and Rationalize, also improves scope, e.g. x^2 + y^2 < Pi *) oineqs = Rationalize[ N[ oineqs, $MachinePrecision + 1 ], 0 ]; If[ AlgebraicQ[ oineqs, vars ], SelectCylinders[ (* Should use WorkingPrecision -> nnn in 4.1 *) CylindricalDecomposition[ oineqs, vars ], vars, Length[vars] ], First[ GenericCylindricalDecomposition[ oineqs, vars ] ] ] ]; (* This is a very rough first classification *) AlgebraicQ[ expr_, vars_ ] := !FreeQ[ expr, _Sqrt|Power[_,_Rational]|_Root ]; (* Message trapping * Notes: This will trap messages issued issued from either * * CylindricalDecomposition * GenericCylindricalDecomposition * * in such a way that I can invoke my own handler of them. * That is the message will still be issued and picked up by Check. *) Unprotect[ Message ]; $CADMsgTrap = False; msg : Message[ CylindricalDecomposition::nrtpi, args___ ] /; $CADMsgTrap := Block[ {$CADMsgTrap=False, $Messages={}}, msg ]; msg : Message[ GenericCylindricalDecomposition::nrtpi, args___ ] /; $CADMsgTrap := Block[ {$CADMsgTrap=False, $Messages={}}, msg ]; Protect[ Message ]; (* * The following functions implement a set of utilities for dealing with nests * and inequalities in cylindric (or nested) form. * * Main functions: * * ToNestedLists[ ineqs, vars ] ==> nests * FromNestedLists[ nests ] ==> ineqs * * Inequalities in nested or cylindric form: * * a1 < x1 < b1 && a2[x1] < x2 < b2[x1] && ... && an[x1,..,xn-1] < xn < bn[x1,..,xn-1] || * c1 < x1 < d1 && c2[x1] < x2 < d2[x1] && ... && cn[x1,..,xn-1] < xn < dn[x1,..,xn-1] || * ... * * where the inequality < can be either of < or <= and similarly > either of > and >=. * * An equivalent representation in terms of nestings (as used in Integrate, NIntegrate etc) is * of the form: * { * { {x1, a1, b1}, {x2, a2[x1], b2[x1]}, ..., {xn, an[x1,..,xn-1], bn[x1,..,xn-1]} }, * { {x1, c1, d1}, {x2, c2[x1], d2[x1]}, ..., {xn, cn[x1,..,xn-1], dn[x1,..,xn-1]} }, * ... * } * * ToNestedLists and FromNestedLists convert to and from these representations. * * Additional notes: * * o The nesting representation which is essentially a generalization of an Interval, * cannot distinguish between open and closed intervals or equivalently between * strict and non-strict inequalities. By adding a function such as DirectedZero[+1] etc * the nesting representation can distinguish between open and closed intervals. * * o Equalities show up in nestings with lower and upper bound equal. Unequal cannot * be represented. * * o CylindricalDecomposition generates a set of inequalities of appropriate * cylindric or nested form. * * Additional functions: * * SelectCylinders[ ineqs, vars, dim ] ==> ineqs * Selects the inequalities defining a set of dimensions dim. * * The following define certain normal forms for logical combinations: * * ToClauseList[ ineq1 || ... || ineqn ] ==> {ineq1, .., ineqn} * FromClauseList[ {ineq1, .., ineqn} ] ==> ineq1 || ... || ineqn * * ToLiteralLists[ a1 && ... && an || b1 && ... && bm || ... ] ==> * {{a1,...,an}, {b1,...,bm}, ... } * FromLiteralLists[ {{a1,...,an}, {b1,...,bm}, ... } ] ==> * a1 && ... && an || b1 && ... && bm || ... * * *) (* * ToClauseList[ logic ] * FromClauseList[ list ] *) ToClauseList[ ineqs_ ] := toClauseList[ LogicalExpand[ ineqs ] ]; ToClauseList[ ineqs_List ] := Flatten[ Map[ ToClauseList, ineqs ] ]; toClauseList[ ineq_Or ] := Apply[ List, ineq ]; toClauseList[ ineq_ ] := {ineq}; FromClauseList[ list_List ] := Apply[ Or, list ] (* * ToLiteralLists[ logic ] * FromLiteralLists[ list ] *) ToLiteralLists[ ineqs_ ] := Map[ toLiteralList, ToClauseList[ ineqs ] ]; ToLiteralLists[ ineq_List ] := ToLiteralLists[ FromLiteralLists[ ineq ] ]; toLiteralList[ ineq_And ] := Apply[ List, ineq ]; toLiteralList[ ineq_ ] := {ineq}; FromLiteralLists[ list_List ] := Apply[ Or, Apply[ And, list, {1} ] ]; (* ToNestedLists[ ineqs, vars ] * Note: This works only for full dimensional components. * Basic parsing of inequalities generated by CAD. * * 1. [ToLiteralLists] * Compute a normal form which is a list of list corresponding to * Or:s of And:s. To separate the various cylinders. * 2. [toNestedComponentList] * Pick out the various components for each variable, i.e. of the form: * { {x1, ineq[x1]}, {x2, ineq[x1, x2]}, ..., {xn, ineq[x1,..,xn-1]}} * 3. [toBounds] * Get the actual bounds for the main variable in the list of ineqs. * 3.1 [findComponentBounds] * Parse the list of inequalities using findBound and combine the result. * 3.1.1 [findBound] * Parse elementary inequalities. * 3.2 [toStandardBounds] * Compute standard representation, i.e. using Infinity when no bound is present. *) ToNestedLists[ fcad_, vars_ ] := Map[ toNestedList[#, vars]&, ToLiteralLists[ LogicalExpand[ fcad ] ] ]; (* The full R^n *) ToNestedLists[ True, vars_ ] := List[ Map[ {#, -Infinity, Infinity}&, vars ] ]; (* The empty set in R^n -- perhaps it should return $Failed *) ToNestedLists[ False, vars_ ] := List[ Map[ {#, 0, 0}&, vars ] ]; toNestedList[ cyl_, vars_ ] := Map[ toBounds, toNestedComponentList[ cyl, vars ] ]; toNestedComponentList[ cyl_List, vars_List ] := toNestedComponentList[ cyl, vars, {}, {} ]; (* Iterative function that finds bounds for one variable at the time *) toNestedComponentList[ rcomps_List,{xf_, xr___}, dcomps_List, dvars_List ] := Module[ {newc}, newc = Select[ rcomps, FreeQ[#, Alternatives[xr]]& ]; toNestedComponentList[ Complement[ rcomps, newc ], {xr}, Append[ dcomps, { xf, newc } ], Append[ dvars, xf ] ] ]; toNestedComponentList[ rcomps_List, {}, dcomps_List, dvars_List ] := dcomps; toBounds[ {x_, {}} ] := toStandardBounds[ {x, {}, {}} ]; toBounds[ {x_, elineqs_List} ] := toStandardBounds[ findComponentBounds[ {x, elineqs} ] ]; (* toStandardBounds[ {x, {lower}, {upper}} ] --> {x, lower, upper} *) toStandardBounds[ {x_, {}, {}} ] := {x, -Infinity, Infinity}; toStandardBounds[ {x_, {}, hi_} ] := toStandardBounds[{x, {-Infinity}, hi}]; toStandardBounds[ {x_, lo_, {}} ] := toStandardBounds[{x, lo, {Infinity}}]; toStandardBounds[ {x_, {lo_},{ hi_ }} ] := {x, lo, hi}; toStandardBounds[ b:{x_, lo_, hi_} ] := (Message[ cylinder::bounds, x, b ];b); (* findComponentBounds[ {x, {ineq1, ..., ineqn}} ] --> {x, {lower}, {upper}} *) findComponentBounds[ {x_, elineqs_List} ] := Prepend[ Apply[ Join, Transpose[ Map[ findBound[#, x]&, elineqs ] ], {1} ], x ]; (* Notes: * $AlgebraicEpsilon is a way for to simulate open intervals: * * a < x < b ==> {x, a + $AlgebraicEpsilon, b - $AlgebraicEpsilon} * * However for this to really work well the following things need * to be taken into account: * o Safe evaluation near the boundary for a given precision * means: * $AlgebraicEpsilon > apmin * o Correct covering by patches as produced in these nests * demands: * $AlgebraicEpsilon < apmax * * So a full analysis would first compute the aggregate level of * approximation through apmax. Then by increasing precision we * we can make apmin small enough that there is some overlap * in these two inequalities. * * Also the reason for wanting to use open intervals is to avoid * various problems on the boundaries, such as complex numbers and * discontinuities. These result in perturbed graphics. And is most * noticeable for higher degree regions. * * This option is currently not used. *) $AlgebraicEpsilon = 0; (* findBound[ x < expr, x ] --> {{}, {expr}} *) findBound[ Less[x_, expr_], x_ ] := {{}, {expr - $AlgebraicEpsilon}}; findBound[ Less[expr_, x_], x_ ] := {{expr + $AlgebraicEpsilon}, {}}; findBound[ Greater[x_, expr_], x_ ] := {{expr + $AlgebraicEpsilon}, {}}; findBound[ Greater[expr_, x_], x_ ] := {{}, {expr - $AlgebraicEpsilon}}; findBound[ LessEqual[x_, expr_], x_ ] := {{}, {expr}}; findBound[ LessEqual[expr_, x_], x_ ] := {{expr}, {}}; findBound[ GreaterEqual[x_, expr_], x_ ] := {{expr}, {}}; findBound[ GreaterEqual[expr_, x_], x_ ] := {{}, {expr}}; (* FromNestedLists[ ineqs, vars ] * Notes: See above for description. *) FromNestedLists[ nests_List ] := Apply[ Or, Map[ Apply[ And, Map[ boundToInequality, # ] ]&, nests ] ]; boundToInequality[ {x_} ] := True; boundToInequality[ {x_, -Infinity, Infinity} ] := True; boundToInequality[ {x_, -Infinity, ux_} ] := x<= ux; boundToInequality[ {x_, lx_, Infinity} ] := x>=lx; boundToInequality[ {x_, lx_, ux_} ] := lx <= x <= ux; (* SetDimension[ cad ] * *) SetDimension[ cad_, vars_ ] := Max[ Map[ cylinderDimension[#,vars]&, ToClauseList[ LogicalExpand[ cad ] ] ] ]; cylinderDimension[ cyl_, vars_ ] := Length[ vars ] - Count[ cyl, _Equal, {0, Infinity} ]; SelectCylinders[ cad_, vars_, dim_ ] := Apply[ Or, Select[ ToClauseList[LogicalExpand[cad]], (cylinderDimension[#,vars] == dim)& ] ]; SelectCylinders[ True, vars_, dim_ ] := True; SelectCylinders[ False, vars_, dim_ ] := False; Protect[ InequalityPlot, InequalityPlot3D, ComplexInequalityPlot ]; End[ (* "`Private`" *) ]; EndPackage[ (* "Graphics`InequalityGraphics`" *) ]