CAD via Projection and Lifting (v7)Dr. Matthew England (University of Bath)28th Nov 2012restart:with(RegularChains):This file is an introduction to the module "ProjectionCAD".read("ProjectionCAD.mm"):
with(ProjectionCAD):The code in this module is an implementation of algorithms for constructing CAD via projection and lifting.The two stages can be done individually via the commands CADProjection and CADLifting, or in one step via CADFull.Each of these commands has two mandatory inputs, first a set of polynomials and second a list of variables (in decreasing order from the main variable). Optional arguments can be used to specify the algorithms used and output formats, as discussed further below.
<Text-field style="Heading 1" layout="Heading 1">The main commands</Text-field>ProjectionCADProjection outputs a complete set of projection polynomials (which includes the input polynomials themselves).Two projection operators are possible, Collins and McCallum. If not specified, then McCallum's will be used.For example:f:=y^3+y^2+y*x^2-1;psetC:=CADProjection( {f}, [y,x], method=Collins ); nops(%);psetM:=CADProjection( {f}, [y,x], method=McCallum ); nops(%);CADProjection( {f}, [y,x] ); LiftingCADLifting takes a set of projection polynomials and outputs a CAD.There are currently four different output formats, discussed further below.Sample points are defined using a regular chain and an isolating cube. To work with them you should load the RegularChains package. Note that the lifting algorithms again depend on whether the method is Collins or McCallum, (with McCallum the default). If you were to lift a set of projections polynomials from one algorithm with the other then the output is no longer guaranteed.CADLifting( psetC, [y,x], method=Collins): nops(%);cad:=CADLifting( psetM, [y,x], method=McCallum): nops(%);CADThe command FullCAD does both projection and lifting and is the recommended choice for most users.CADFull( {f}, [y,x], method=McCallum ): nops(%);
CADFull( {f}, [y,x], method=Collins ): nops(%);
<Text-field style="Heading 1" layout="Heading 1">CAD output formats</Text-field>Output formatsThere are currently four style of output: list, listwithrep, rootof and piecewise. The default is list, which gives a list of cells, each containing an index and a samplepoint. Note that this is the same format as RegularChains:-SemiAlgebraicSetTools:-CylindricalAlgebraicDecompose with output=list.cadl:=CADFull( {f}, [y,x], method=McCallum, output=list ); nops(%);The choice output=rootof will give a list of cells, each represented by conditions on the variables.Note that this is the same format as RegularChains:-SemiAlgebraicSetTools:-CylindricalAlgebraicDecompose with output=rootof.cadr:=CADFull( {f}, [y,x], method=McCallum, output=rootof ); nops(%);The choice output=listwithrep will give a list of cells which contain all the information from the previous two formats.Note that there is no comparable format within RegularChains:-SemiAlgebraicSetTools:-CylindricalAlgebraicDecompose.cadlwr:=CADFull( {f}, [y,x], method=McCallum, output=listwithrep ); nops(%);The choice output=piecewise will show the cad arranged cylindrically using Maple's piecewise construction.This format is comparable with RegularChains:-SemiAlgebraicSetTools:-CylindricalAlgebraicDecompose with output=piecewise, however, here radicals are used instead of RootOf constructions for easier examples.cadp:=CADFull( {f}, [y,x], method=McCallum, output=piecewise ); Note that we cannot measure the number of cells in a piecewise cad using the nops functions. Use the function CADNumCellsInPiecewise insteadnops(cadp); CADNumCellsInPiecewise(cadp);
<Text-field style="Heading 1" layout="Heading 1">Possible failure and userinfo level</Text-field>FailureCollins algorithm will always succeed (given sufficient time and memory).McCallum's algorithm can fail when the input polynomials are not well oriented. (Note: Unlike qepcad, they do not need to be very well oriented.)McCallum's CADW algorithm has been adapted to maximise success by guaranteeing only a sign-invariant CAD (rather than order-invariant) as a final output.(This is because if failure to be well-oriented occurs at the final lift then the output is still a sign-invariant CAD).Hence to give an example of failure we need at least 5-dimensions (This is because 2d is always OI, 3d can only fail on dim zero cell for which we have McCallum's delineating polynomial method and 4d can still be guaranteed sign-invariant.)A trivial example of failure:f:=a*e+b*d+c*e+d+e;CADFull( {f}, [a,b,c,d,e], method=McCallum ): nops(%);In this case the algorithm continued, but the warning implies that the output may be meaningless. We can change what happens in this situation with an optional argument. Using failure=err will generate an error message as opposed to the default warning. Using failure=giveFAIL will return the boolean value FAIL.CADFull( {f}, [a,b,c,d,e], method=McCallum, failure=err );CADFull( {f}, [a,b,c,d,e], method=McCallum, failure=giveFAIL );Note that we can get a guaranteed correct output by using Collins:CADFull( {f}, [a,b,c,d,e], method=Collins ): nops(%);RegularChains:-SemiAlgebraicSetTools:-CylindricalAlgebraicDecompose( [f], RegularChains:-PolynomialRing([a,b,c,d,e]) ): CADNumCellsInPiecewise(%);Note the number of cells in Collins guaranteed SI CAD was the same as in the "failed" McCallum cad. Hints that McCallum's algorithm does not really fail? The MapleCAD is far smaller however.UserInfoWe can get more information of the progress of the algorithm by changing the userinfo level from its default of 1. Level 3 will describe what is happening while level 4 will actually print some of the polynomials. infolevel[CADFull]:=4:CADFull( {f}, [a,b,c,d,e], method=McCallum ): nops(%);
<Text-field style="Heading 1" layout="Heading 1">Order Invariance</Text-field>Order-invarianceAs discussed above, McCallum's CADW algorithm has been modified to maximise success - this is the default.To run the original CADW algorithm, guaranteeing order-invariance at the final step, (or failure), then the optional argument finalCAD=OI must be given.Consider for example:f:=x^2-z*y; infolevel[CADFull]:=3:cad0:=CADFull( {f}, [z,y,x], method=Collins, output=listwithrep ): nops(%);cad1:=CADFull( {f}, [z,y,x], method=McCallum, output=listwithrep ): nops(%);cad2:=CADFull( {f}, [z,y,x], method=McCallum, finalCAD=OI, output=listwithrep ): nops(%);The first cad only has 21 cells, the second 23. Look at the difference by considering the indicesconvert( map(X->op(1,X), cad2),set) minus convert(map(X->op(1,X),cad1),set);select(X->op(1,X)=[2,2,1], cad2);
select(X->op(1,X)=[2,2,2], cad2);
select(X->op(1,X)=[2,2,3], cad2);select(X->op(1,X)=[2,2,1], cad1);The difference occurs when x=0 and y=0. The first cad has just one cell here, with z free. The second splits into three cells, with z<0,z=0 and z>0.This is clearly needed for OI, since the order of f at zero is 2 while the order elsewhere is 1 or 0.By considering the user info statements we see why the outputs differ. In both cases there is a nullification at the final lift. The first time this was ignored since it was the final lift and it was assumed only a sign-invariant cad is required.The second time this nullification was examined. Since the cell in question was zero-dimensional we could use McCallum's delineating polynomial method to guarantee order invariance.
<Text-field style="Heading 1" layout="Heading 1">Generate Stack</Text-field>GenerateStackPart of the lifting algorithm involves building a stack over a cell such that a set of polynomials are delineable. To do this we use the internal RegularChains algorithms, RegularChains:-TRD_generate_stack. The polynomials passed to this algorithm should separate above the cell, hence to ensure correctness, we must first process our polynomials so that this is true. The user algorithm, CADGenerateStack first error checks, then does this pre-processing before finally passing to the TRD algorithm. Consider for example:f:=z^2+y^2+x^2-1;
fullcad:=CADFull([f], [z,y,x], method=McCallum): nops(%);We can force CADFull to break halfway and return the cad of over (x,y)-space using the optional argument retcad=ixycad:=CADFull([f], [z,y,x], method=McCallum, retcad=2): nops(%);Consider a cell of this cad.cell:=xycad[3];If we generate the stack over this cell wrt to f, we get a list of new cells.CADGenerateStack( cell, [f], [z,y,x] );Note that each cell has an index which starts with the index of the previous cell. If we were to do this to all cells and put them together we get the full cad:cad:=[]:
for i from 1 to 13 do
cad:=[op(cad), op(CADGenerateStack(xycad[i],[f],[z,y,x]))]:
od: nops(cad);The two cads can be shown to be identical under inspection