diff --git a/M2/Macaulay2/packages/=distributed-packages b/M2/Macaulay2/packages/=distributed-packages index 25c732be3b0..ea5fc4cc301 100644 --- a/M2/Macaulay2/packages/=distributed-packages +++ b/M2/Macaulay2/packages/=distributed-packages @@ -295,3 +295,4 @@ MRDI EliminationTemplates WittVectors Padic +EuclideanDistanceDegree \ No newline at end of file diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree.m2 b/M2/Macaulay2/packages/EuclideanDistanceDegree.m2 new file mode 100644 index 00000000000..fc274504319 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree.m2 @@ -0,0 +1,701 @@ +newPackage( + "EuclideanDistanceDegree", + Version => "1.1", + Date => "April 2026", + Authors => { + {Name => "Jose Israel Rodriguez", + Email => "Jose@Math.wisc.edu", + HomePage => "http://www.math.wisc.edu/~jose/"}, + {Name => "Will Huang", + Email => "whuang259@wisc.edu", + HomePage => ""} + }, + Headline => "produce critical equations and compute ED degrees", + DebuggingMode => false, + AuxiliaryFiles => true, + PackageImports => {"Elimination","MonodromySolver"}, + PackageExports => {"Bertini","NumericalAlgebraicGeometry"}, + Configuration => { "Continuation"=>"Bertini" }, + OptionalComponentsPresent => (readPackage "Bertini").OptionalComponentsPresent, + CacheExampleOutput => true +) + +randomCC=()->random CC +randCC=()->random CC +randomRR=()->((-1)^(random(1,2)) *random RR) +randomZZ=()->random(1,30103) +randomValue=(kk)-> if kk===CC then randomCC() else if kk===RR then randomRR() else randomZZ() +randomVector=method(Options=>{ }) +randomVector(ZZ,Thing):= o->(n,R) ->apply(n,i->randomValue(R))--list of length n of randomValue + +load "./EuclideanDistanceDegree/Determinantal.m2" +load "./EuclideanDistanceDegree/LeftKernel.m2" +load "./EuclideanDistanceDegree/Numerical.m2" +load "./EuclideanDistanceDegree/Parameterization.m2" + +export { + "TempDirectory", + "ReturnCriticalIdeal", + "UseMonodromy", + "Submodel", + "SampleGenerator", + "symbolicWeightEDDegree", + "determinantalUnitEDDegree", + "determinantalGenericEDDegree", + "leftKernelWeightEDDegree", + "leftKernelUnitEDDegree", + "leftKernelGenericEDDegree", + "newNumericalComputationOptions", + "NumericalComputationOptions", + "homotopyEDDegree", + "numericWeightEDDegree", + "numericGenericEDDegree", + "numericUnitEDDegree", + "parameterizedWeightEDDegree", + "parameterizedGenericEDDegree", + "parameterizedUnitEDDegree", + "averageNumericEDDegree" +} + +--##########################################################################-- +-- INTERNAL METHODS +--##########################################################################-- +parString = (aString) -> ("("|toString(aString)|")"); +addSlash = (aString) -> ( + if #aString =!= 0 then ( + if aString_-1 === " " then error(aString | " cannot end with whitespace."); + if aString_-1 =!= "/" then aString = aString | "/"; + ); + aString +); +makeJac = (system,unknowns) -> ( + -- it is a list of lists of partial derivatives of a polynomial + for i in system list for j in unknowns list diff(j,i) +) +checkZero = (aSol, eps) -> if aSol/abs//min < eps then false else true +sortPointFunction = (aSol) -> (if not (apply(aSol,i->{realPart i,imaginaryPart i}/abs//max)//min<1e-8) then true else false); + +--##########################################################################-- +-- DOCUMENTATION +--##########################################################################-- + +beginDocumentation() + +doc /// -- Package + Key + EuclideanDistanceDegree + Headline + produce critical equations and compute ED degrees + Description + Text + The Euclidean distance (ED) degree of a varieties arises from the following problem: given a (generic) data point, find the point on a + given variety which minimizes the (possibly weighted) Euclidean distance function. The number of nonsingular critical points of the + distance function is the ED degree of the variety. + Text + This package provides several routines for determining the (generic or unit) Euclidean distance degree of an algebraic variety. These + routines include symbolic methods and numerical methods for determining these degrees; they can be grouped into four main types: + Code + UL { + TO symbolicWeightEDDegree, + TO leftKernelWeightEDDegree, + TO numericWeightEDDegree, + TO parameterizedWeightEDDegree + } + Text + As an example, this code computes the (unit) ED degree of a circle both symbolically and numerically + Example + R = QQ[x,y]; + F = {x^2 + y^2 - 1}; + determinantalUnitEDDegree(F) + leftKernelUnitEDDegree(F) + Text + When the variety is an affine cone, one is able to compute ED degrees using ED degree homotopies, i.e., a structured parameter + homotopy. The easiest case is when the variety is a hypersurface (or more generally, a complete intersection) + Example + R = QQ[x1,x2,x3,x4]; + F = G = {det genericMatrix(R,2,2)}; + numericGenericEDDegree(F, G) + Text + When a $V(F)$ is not a complete intersection we incorporate a membership test to filter out residual critical points. Here $V(F)$ + is an irreducible component of $V(G)$ (a reducible variety) and `#G===codim ideal F`. These methods employ an equation by equation + method called regeneration. + Example + R = QQ[x1,x2,x3,x4,x5,x6]; + F = (minors(2, genericMatrix(R,3,2)))_*; + G = drop(F, -1); + #G == codim ideal F; + numericGenericEDDegree(F, G) + Text + One may also determine ED degrees using a parameter homotopy called a Weight-ED Degree Homotopy. This code computes a unit ED degree: + Example + R = QQ[x1,x2,x3,x4,x5,x6]; + F = (minors(2, genericMatrix(R,3,2)))_*; + G = drop(F, -1); + NCO = newNumericalComputationOptions(F, G); + NCO#"TargetWeight" = apply(#gens R, i->1); + homotopyEDDegree(NCO, "Weight", true, true) +/// + +doc /// -- symbolic + Key + ReturnCriticalIdeal + [symbolicWeightEDDegree, ReturnCriticalIdeal] + [determinantalUnitEDDegree, ReturnCriticalIdeal] + [determinantalGenericEDDegree, ReturnCriticalIdeal] + Headline + whether to return the ideal of critical points + Usage + ICP = symbolicWeightEDDegree(F, U, W, ReturnCriticalIdeal => true) + Description + Text + Symbolic methods compute the Euclidean Distance degree of a variety by computing the degree of its critical ideal, which defines the + critical points of the distance function on the variety. Setting this option will cause symbolic methods to return this ideal rather + than its degree. + Example + R = QQ[x,y]; + F = {x^2 + y^2 - 1}; + (U,W) = ({12, 23}, {15, 331}); + ICP = symbolicWeightEDDegree(F, U, W, ReturnCriticalIdeal => true) +/// + +doc /// + Key + symbolicWeightEDDegree + (symbolicWeightEDDegree, List, List, List) + determinantalUnitEDDegree + (determinantalUnitEDDegree, List) + determinantalGenericEDDegree + (determinantalGenericEDDegree, List) + Headline + symbolically compute ED degrees of affine varieties + Usage + UED = determinantalUnitEDDegree(F) + GED = determinantalGenericEDDegree(F) + GED = symbolicWeightEDDegree(F,U,W) + Inputs + F:List + a system of polynomials (system need not be square) + U:List + a (generic) data vector + W:List + a (generic) weight vector + Outputs + ED:ZZ + a generic/unit Euclidean distance degree + ICP:Ideal + ideal of critical points (if ReturnCriticalIdeal is set) + Description + Text + This method computes Euclidean distance (ED) degrees for the variety defined by the system $F$ via symbolic computation. The ideal of + critical points is formed by saturating the defining ideal of the variety with minors of the Jacobian and Augmented Jacobian. The + degree of this ideal is the ED degree of the variety. The unit variant of this method computes an ED degree using random (integer) + data and unit weights, whereas the generic variant will use random data and random weights. + Example + R = QQ[x,y]; + F = {x^2 + y^2 - 1}; + (U,W) = ({12, 23}, {15, 331}); + UED = determinantalUnitEDDegree F + GED = determinantalGenericEDDegree F + ICP = symbolicWeightEDDegree(F, U, W, ReturnCriticalIdeal => true) +/// + +doc /// -- parameterized + Key + UseMonodromy + [parameterizedWeightEDDegree, UseMonodromy] + [parameterizedUnitEDDegree, UseMonodromy] + [parameterizedGenericEDDegree, UseMonodromy] + Headline + whether to use monodromy methods to compute degrees + Usage + GED = parameterizedWeightEDDegree(F, U, W, UseMonodromy => true) + Description + Text + If this option is set, the parameterized ED degree methods will use the MonodromySolver to numerically solve the critical equations + and report the number of regular solutions as the ED degree of the parameterized variety. + Example + R = QQ[x,y]; + F = {x^2 + 1, x * y, y - 1}; + GED = parameterizedGenericEDDegree(F, UseMonodromy => true) + Caveat + Using monodromy may result in an ED degree that is less than expected. In particular, if the parameterization map is $n$-to-one, + using monodromy to solve for the degree may result in a degree than is off by a factor of $n$. +/// + +doc /// + Key + parameterizedWeightEDDegree + (parameterizedWeightEDDegree, List, List, List) + parameterizedUnitEDDegree + (parameterizedUnitEDDegree, List) + parameterizedGenericEDDegree + (parameterizedGenericEDDegree, List) + Headline + compute ED degrees of parameterized varieties + Usage + UED = parameterizedUnitEDDegree(F) + GED = parameterizedGenericEDDegree(F) + GED = parameterizedWeightEDDegree(F,U,W) + Inputs + F:List + a system of polynomials in $d$ variables parameterizing a $d$-dimensional variety + U:List + a (generic) data vector + W:List + a (generic) weight vector + Outputs + ED:ZZ + a generic/unit Euclidean distance degree + Description + Text + This method computes Euclidean distance (ED) degrees for the variety parameterized by the a set of polynomials $F$ in $d$ variables. If + the resulting variety is $d$-dimensional, then by finding a global description for the kernel of the transpose of the Jacobian map, the + critical equations of the variety can be computed. The unit variant of this method computes an ED degree using random (integer) data + and unit weights, whereas the generic variant will use random data and random weights. + Example + R = QQ[x,y]; + F = {x^2 + 1, x * y, y - 1}; + (U,W) = ({12, 23, 25}, {15, 331, 1}); + UED = parameterizedUnitEDDegree F + GED = parameterizedGenericEDDegree F + GED = parameterizedWeightEDDegree(F, U, W) +/// + +doc /// -- TempDirectory + Key + TempDirectory + [leftKernelWeightEDDegree, TempDirectory] + [leftKernelUnitEDDegree, TempDirectory] + [leftKernelGenericEDDegree, TempDirectory] + [numericWeightEDDegree, TempDirectory] + [numericUnitEDDegree, TempDirectory] + [numericGenericEDDegree, TempDirectory] + [averageNumericEDDegree, TempDirectory] + Headline + set directory for Bertini runs + Usage + GED = leftKernelGenericEDDegree(F, U, W, TempDirectory => dir) + Description + Text + This option is used to change the directory from which Bertini files are written to and run from for the numerical methods. If unset, + a temp directory will be created and used. + Example + R = QQ[x,y]; + F = {x^2 + y^2 - 1}; + dir = temporaryFileName(); + GED = leftKernelGenericEDDegree(F, TempDirectory => dir) +/// + +doc /// -- leftKernel + Key + leftKernelWeightEDDegree + (leftKernelWeightEDDegree, List, List, List) + leftKernelUnitEDDegree + (leftKernelUnitEDDegree, List) + leftKernelGenericEDDegree + (leftKernelGenericEDDegree, List) + Headline + numerically compute Euclidean distance degrees of complete intersections + Usage + UED = leftKernelUnitEDDegree(F) + GED = leftKernelGenericEDDegree(F) + GED = leftKernelWeightEDDegree(F,U,W) + Inputs + F:List + a system of polynomials (need not be square) defining an affine variety that is a complete intersection + U:List + a (generic) data vector + W:List + a (generic) weight vector + Outputs + ED:ZZ + a generic/unit Euclidean distance degree + Description + Text + This method computes Euclidean distance (ED) degrees for affine varieties defined by the system $F$ numerically using Lagrange + multipliers. If the resulting variety is a complete intersection, the left kernel of the augmented Jacobian is used to derive a set of + critical erquations which are passed into Bertini. The resulting number of critical points is returned as the ED degree. The unit + variant of this method computes an ED degree using random (complex) data and unit weights, whereas `leftKernelGenericEDDegree` will + use random data and random weights. + Example + R = QQ[x,y]; + F = {x^2 + y^2 - 1}; + (U,W) = ({.12, .23}, {.15, .331}); + UED = leftKernelUnitEDDegree F + GED = leftKernelGenericEDDegree F + GED = leftKernelWeightEDDegree(F, U, W) + Caveat + The computed ED degree may be lower than expected due to path tracking. +/// + +doc /// -- NumericalComputationOptions + Key + NumericalComputationOptions + Headline + define homotopy options and configurations + Description + Text + This object stores all the options needed compute an ED degree using homotopy continuation with Bertini. It is created using the + @TO newNumericalComputationOptions@ method. + Text + Keys available to customize the homotopy include: + "Directory" to change the directory Bertini is run from, the directory need not exist initially. + "StartData" and "TargetData" if executing a Data homotopy, defaults to random complex data. + "StartWeight" and "TargetWeight" if executing a Weight homotopy, defaults to random complex and unit weights respectively. + "Infinity" to add a hyperplane at infinity for homogenization, one is created by default in @TO homotopyEDDegree@ if not present. + "PrimalCoordinates" to add additional variables to the ambient space. By default this field will already contain the variables of + `ring F`, thus it is safer to append extra variables rather than overwriting the less, otherwise subsequent methods may fail. + Text + Keys are also available to customize the Bertini run: + "FinerRestriction" a list of polynomials to filter down critical points. Critical points will only be kept if they vanish + for every polynomial in this list. + "BertiniStartFiberSolveConfiguration" a list of options (e.g. `FinalTol -> 1e-8`) that will be passed into the Berni inputs file. By + default, the options `{"TrackType"=>0, "PrintPathProgress"=>1000}` are included. +/// + +doc /// + Key + newNumericalComputationOptions + (newNumericalComputationOptions, List, List) + [newNumericalComputationOptions, TempDirectory] + [newNumericalComputationOptions, Submodel] + Headline + define homotopy options and configurations + Usage + NCO = newNumericalComputationOptions(F, G) + Inputs + F:List + a system of polynomials (need not be square) defining the variety + G:List + a system of polynomials (complete intersection) such that V(F) is an irreducible component of V(G), i.e. a witness model + Submodel => List + a system of linear polynomials to slice the variety + Outputs + NCO:NumericalComputationOptions + a MutableHashTable to keep track of options and configurations for the homotopy methods + Description + Text + Creates a @TO NumericalComputationOptions@ object. At mimumum it requires the model $F$ for which the ED Degree will be computed and + a witness model $G$. A submodel $L$ may be passed in to "slice" the variety, by default it will be the empty set. A String indicating + the temporary directory from which Bertini will read/write files during the run may also be passed in as an optional argument, by + default it will be a random temporary file name. The temporary directory will be created if it does not exist. + Example + R = QQ[x,y]; + F = G = {x^2 + y^2 - 1}; + dir = temporaryFileName(); + NCO = newNumericalComputationOptions(F, G, TempDirectory => dir); + NCO#"TargetWeight" = apply(#gens R, i -> 1_R); + UED = homotopyEDDegree(NCO, "Weight", true, true) +/// + +doc /// -- homotopy + Key + homotopyEDDegree + (homotopyEDDegree, NumericalComputationOptions, String, Boolean, Boolean) + Headline + numerically compute ED degrees of affine cones using homotopy continuation + Usage + GED = homotopyEDDegree(NCO, "Weight", true, false) + Inputs + NCO:NumericalComputationOptions + a MutableHashTable to keeps track of the options and configurations for the homotopy methods + ht:String + one of "Weight" or "Data" which defines the type of homotopy to execute + isStageOne:Boolean + if set, runs stage 1 to solve the start system + isStageTwo:Boolean + if set, runs stage 2 to solve the target system via path tracking + Outputs + ED:ZZ + a generic/unit Euclidean distance degree + Description + Text + Executes the homotopy defined by the passed in @TO NumericalComputationOptions@ object. A weight homotopy will vary the weights of the + distance function whereas a data homotopy will vary the data point. A stage1 homotopy solves the start system as defined in NCO whereas + stage2 will perform the weight/data homotopy to solve the defined target system. It is possible to perform stage2 multiple times for + different target systems without rerunning stage1, but at least one stage1 homotopy must be executed in the directory specified in NCO + for stage2 to executed properly. + Example + R = QQ[x,y]; + F = G = {x^2+y^2-1}; + NCO = newNumericalComputationOptions(F, G); + NCO#"TargetWeight" = apply(#gens R, i -> random RR); + GED = homotopyEDDegree(NCO, "Weight", true, true) + + NCO#"TargetWeight" = apply(#gens R, i -> 1_R); + UED = homotopyEDDegree(NCO, "Weight", false, true) + Caveat + Inaccurate results may be returned if $V(F)$ is contained in $V(L)$. The computed ED degree may be lower than expected due to path tracking. +/// + +doc /// -- numeric + Key + numericWeightEDDegree + (numericWeightEDDegree, List, List, List, List) + numericUnitEDDegree + (numericUnitEDDegree, List, List) + numericGenericEDDegree + (numericGenericEDDegree, List, List) + [numericWeightEDDegree, Submodel] + [numericUnitEDDegree, Submodel] + [numericGenericEDDegree, Submodel] + Headline + numerically compute ED degrees of affine cones using homotopy continuation + Usage + UED = numericUnitEDDegree(F,G) + GED = numericGenericEDDegree(F,G) + GED = homotopyEDDegree(F,G,U,V) + Inputs + F:List + a system of polynomials (need not be square) defining the variety + G:List + a system of polynomials (complete intersection) such that $V(F)$ is an irreducible component of $V(G)$, i.e. a witness model + U:List + a (generic) data vector + W:List + a (generic) weight vector + Submodel => List + list of linear polynomials to slice the variety + Outputs + ED:ZZ + a generic/unit Euclidean distance degree + Description + Text + Special instances of the @TO homotopyEDDegree@ method that uses a homotopy on weights. The unit variant of this method computes an ED + degree using random (complex) data and unit weights, whereas `numericGenericEDDegree` will use random data and random weights. + Example + R = QQ[x,y]; + F = G = {x^2+y^2-1}; + (U,W) = ({.12, .23}, {.15, .331}); + UED = numericUnitEDDegree(F, G) + GED = numericGenericEDDegree(F, G) + GED = numericWeightEDDegree(F, G, U, W) + Caveat + Inaccurate results may be returned if $V(F)$ is contained in $V(L)$. The computed ED degree may be lower than expected due to path tracking. +/// + +doc /// -- average + Key + averageNumericEDDegree + (averageNumericEDDegree, List, List, ZZ) + [averageNumericEDDegree, Submodel] + [averageNumericEDDegree, SampleGenerator] + [averageNumericEDDegree, Tolerance] + Headline + compute average ED degrees using sampled data + Usage + aED = averageNumericEDDegree(F, G, 10) + Inputs + F: List + a system of polynomials (need not be square) defining the variety + G: List + list of polynomials (complete intersection) such that $V(F)$ is an irreducible component of $V(G)$, i.e. a witness model + n: ZZ + number of samples to take + Submodel => List + list of linear polynomials to slice the variety + SampleGenerator => FunctionClosure + a function which generates data samples + Tolerance => RR + tolerance to use when checking if a critical point is real + Outputs + aED: RR + average ED degree after n trials + Description + Text + Generate data samples using the given function and uses homotopy continuation to find critical points of the distance function. The + average number of real critical points after $n$ trials is returned. This method creates a @TO NumericalComputationOptions@ object and + computes critical points using the @TO homotopyEDDegree@ method. By default, `random(RR)` is used to generate data samples. Points are + tested using the @TO realPoints@ function from the @TO NumericalAlgebraicGeometry@ package, a tolerance can be passed along using the + `Tolerance` option, by default it is 1e-6. + Text + + Example + R = QQ[x,y]; + F = G = {x^2 + y^2 - 1}; + sampleGen = () -> apply(#gens R, i -> random(RR)); + aED = averageNumericEDDegree(F, G, 10, SampleGenerator => sampleGen) +/// + +--##########################################################################-- +-- END DOCUMENTATION +--##########################################################################-- + +TEST /// -- basic examples + setRandomSeed(123); + R = QQ[x0, x1, x2]; + F = {x0^2*x2 - x1^2*(x1 + x2)}; + assert(determinantalGenericEDDegree(F) === 7); + assert(determinantalUnitEDDegree(F) === 7); + + R = QQ[jj]/ideal(jj^2+1)[x0,x1,x2]; + F = {x0^2*x1 -(x1 - jj*x2)^2*x2}; + assert(determinantalGenericEDDegree(F) === 2*7); + assert(determinantalUnitEDDegree(F) === 2*7); + + R = QQ[x0,x1,x2,x3]; + F = {x0^2*x1-x2*x3^2}; + assert(determinantalGenericEDDegree(F) === 10); + assert(determinantalUnitEDDegree(F) === 10); + + R = QQ[x,y]; + F = {x^2 + y^2 - 1}; + (U, W) = ({12, 23}, {15, 331}); + assert(symbolicWeightEDDegree(F, U, W) === 4); + assert(determinantalGenericEDDegree F === 4); + assert(determinantalUnitEDDegree F === 2); + + R = QQ[x,y]; + F = G = {x^2+y^2-1}; + (U, W) = ({.12, .23}, {.15, .331}) + assert(leftKernelWeightEDDegree(F, U, W) === 4); + assert(leftKernelGenericEDDegree F === 4); + assert(leftKernelUnitEDDegree F === 2); + + assert(numericWeightEDDegree(F, G, U, W) === 4); + assert(numericGenericEDDegree(F, G) === 4); + assert(numericUnitEDDegree(F, G) === 2); +/// + +TEST /// -- Selected surfaces from Herwig Hauser's algebraic surfaces gallery + setRandomSeed(123); + R = QQ[x,y,z]; + surfaces = { + 3*x^2 + 3*y^2 + z^2 - 1, + x^2 - y^2*z, + x^2 + y^2 + z^2 - 1, + x^2 + z^2 - y^3*(y-1)^3, + x^2*y*z + x*y^2 + y^3 + y^3*z - x^2*z^2 + }; + + for surface in surfaces do ( + F = {surface}; + + GED_symb = determinantalGenericEDDegree F; + GED_left = leftKernelGenericEDDegree F; + GED_numeric = numericGenericEDDegree(F, F); + + assert(GED_symb === GED_left); + assert(GED_left === GED_numeric); + ) +/// + +TEST /// -- PNN function space + d = (3,1,1); + r = 2; + + R = QQ[W_0..W_(d_1 * d_0), V_0..V_(d_2 * d_1)]; + W = genericMatrix(R, W_0, d_1, d_0); + V = genericMatrix(R, V_0, d_2, d_1); + + T = R[x_0..x_(d_0 - 1)]; + X = transpose matrix{apply(d_0, i -> x_i)}; + Z = W * X; + A = matrix table(d_1, 1, (i, j) -> (Z_(i,j))^r); + Phi = V * A; + + mons = flatten entries basis(r, T); + S = QQ[c_0..c_(d_2 * #mons - 1)]; + im = flatten apply(d_2, i -> ( + f = Phi_(i,0); + apply(mons, m -> coefficient(m, f)) + )); + + paramMap = map(R, S, im); + I = kernel paramMap; + F = flatten entries gens I; + c = codim I; + G = apply(c, i -> sum apply(#F, j -> (random(QQ) * F_j))); + assert(leftKernelUnitEDDegree(G) === numericUnitEDDegree(G, G)); + assert(leftKernelGenericEDDegree(G) === numericGenericEDDegree(G, G)); +/// + +TEST /// -- parameterization: basic test + setRandomSeed(123456); + R = QQ[x,y]; + F = {x^2 + 1, x * y, y - 1}; + GED_param = parameterizedGenericEDDegree F; + + n = #F; + numX = #gens R; + S = QQ[gens R] ** QQ[y_1..y_(n-numX), u_1..u_n]; + M = sub(matrix{F}, S); + imageModel = eliminate(support M, ideal(M - matrix{for i from 1 to n list u_i})); + GED_implicit = determinantalUnitEDDegree((sub(imageModel, QQ[support imageModel]))_*); + assert(GED_param === GED_implicit); +/// + +TEST /// -- parameterization: spurious critical points + setRandomSeed(123456); + R = QQ[x]; + F = {x^2, x^3}; + GED_param = parameterizedGenericEDDegree F; + + n = #F; + numX = #gens R; + S = QQ[gens R] ** QQ[y_1..y_(n-numX), u_1..u_n]; + M = sub(matrix{F}, S); + imageModel = eliminate(support M, ideal(M - matrix{for i from 1 to n list u_i})); + GED_implicit = determinantalUnitEDDegree((sub(imageModel, QQ[support imageModel]))_*); + assert(GED_param === GED_implicit); +/// + +TEST /// -- parameterization: monodromy vs symbolic + setRandomSeed(123); + R = QQ[x1,x2,x3,x4,x5]; + F = {x1^2+x4^2, x2^2+x5^2, x3^2+1, x1*x2+x4*x5, x1*x3+x4, x2*x3+x5}; + U = {1,2,3,4,5,6}; + W = {1,1,1,1,1,1}; + GED1 = parameterizedWeightEDDegree(F,U,W); + GED2 = parameterizedWeightEDDegree(F,U,W, UseMonodromy => true); + assert(GED1 === GED2); +/// + +end + + + +uninstallPackage "EuclideanDistanceDegree" +restart +loadPackage "EuclideanDistanceDegree" +installPackage "EuclideanDistanceDegree" +check "EuclideanDistanceDegree" + +viewHelp + +-- Debugging surface test seed issue +loadPackage "EuclideanDistanceDegree" +setRandomSeed(123); +R = QQ[x,y,z]; + surfaces = { + 3*x^2 + 3*y^2 + z^2 - 1, + x^2 + y^2 + z^3 - z^2, + x^2 + y^3 + z^5, + x^2 - y^2*z, + x^2 + y^2*z^3, + x^2 + y^2 + z^2 - 1, + x^2 + z^2 - y^3*(y-1)^3, + x^2 - y^3*z^3, + x^2 + y^2 + z, + x^2*y*z + x*y^2 + y^3 + y^3*z - x^2*z^2 + }; + +GED1 = {}; +GED2 = {}; +GED3 = {}; + +for surface in surfaces do ( + F = {surface}; + print(toString surface); + + GED1 = GED1 | {determinantalGenericEDDegree F}; + GED2 = GED2 | {leftKernelGenericEDDegree F}; + GED3 = GED3 | {numericGenericEDDegree(F, F)}; +) + +print("checking degrees"); +scan(#GED1, i -> ( + if GED1_i =!= GED2_i then print("symb-left F: " | toString surfaces_i | "; " | GED1_i | ", " | GED2_i); + if GED2_i =!= GED3_i then print("left-num F: " | toString surfaces_i | "; " | GED2_i | ", " | GED3_i); + if GED1_i =!= GED3_i then print("sym-num F: " | toString surfaces_i | "; " | GED1_i | ", " | GED3_i); +)) \ No newline at end of file diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/Determinantal.m2 b/M2/Macaulay2/packages/EuclideanDistanceDegree/Determinantal.m2 new file mode 100644 index 00000000000..5357cd050d7 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/Determinantal.m2 @@ -0,0 +1,49 @@ +rand := randomZZ --This is the function used to compute random numbers. +--random(QQ) is not a good choice in practice. + +--This function computes weighted ED degrees +symbolicOptions = { ReturnCriticalIdeal => false } +symbolicWeightEDDegree = method(Options => symbolicOptions) +symbolicWeightEDDegree(List, List, List) := o-> (F, data, weight) -> ( + I := ideal F; + xList := gens ring I; + jac := matrix makeJac(F, xList); + topRow := apply(#weight, i -> 2 * weight_i * (xList_i-data_i)); + M := matrix{topRow} || jac; -- augmented Jacobian + c := codim I; + win := I + minors(c+1, M); -- critical ideal + sl := ideal singularLocus I; + win = saturate(win, sl); + if o.ReturnCriticalIdeal then return win else return degree win +) + +determinantalUnitEDDegree = method(Options => symbolicOptions) +determinantalUnitEDDegree(List) := o-> (F) -> symbolicWeightEDDegree( + F, + apply(#gens ring first F, i->rand()), + apply(#gens ring first F, i->1_(ring first F)), + o +) + +determinantalGenericEDDegree = method(Options => symbolicOptions) +determinantalGenericEDDegree(List) := o-> (F) -> symbolicWeightEDDegree( + F, + apply(#gens ring first F, i->rand()), + apply(#gens ring first F, i->rand()), + o +) + +end + +-* +restart +load "./EDD_Determinantal.m2" +makeJac = (system, unknowns) -> (for i in system list for j in unknowns list diff(j, i)) + +R = QQ[x, y]; +F = {x^2 + y^2 - 1} +(U,W) = ({12, 23}, {15, 331}) +UED = determinantalUnitEDDegree F +GED = determinantalGenericEDDegree F +ICP = symbolicWeightEDDegree(F, U, W, ReturnCriticalIdeal => true) +*- \ No newline at end of file diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/LeftKernel.m2 b/M2/Macaulay2/packages/EuclideanDistanceDegree/LeftKernel.m2 new file mode 100644 index 00000000000..304eabffcae --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/LeftKernel.m2 @@ -0,0 +1,72 @@ +leftKernelOptions = { TempDirectory => null } +leftKernelWeightEDDegree = method(Options => leftKernelOptions) +leftKernelWeightEDDegree(List, List, List) := o-> (F, data, weight) -> ( + dir := ""; + if o.TempDirectory === null then dir = temporaryFileName() else dir = o.TempDirectory; + if not fileExists dir then mkdir dir; + + R := ring first F; + numX := #gens R; + coef := coefficientRing R; + S := R ** coef[apply(#F+1, i->"L"|i)] ** coef[apply(numX, i->"u"|i)] ** coef[apply(numX, i->"w"|i)]; + + xList := flatten entries basis({1,0,0,0}, S); + lamList := flatten entries basis({0,1,0,0}, S); + uList := flatten entries basis({0,0,1,0}, S); + wList := flatten entries basis({0,0,0,1}, S); + + c := #lamList - 1; + jac := sub(matrix makeJac(apply(F, i->sub(i,S)), xList), S); + topRow := apply(#weight, i->2*wList_i*(xList_i-uList_i)); + M := matrix{topRow} || jac; -- augmented Jacobian + critEq := flatten entries((matrix{lamList} * sub(M, S))); + restrictLam := apply(#lamList - 1 - #F, i->makeB'Section(drop(lamList, 1))); + + win := restrictLam | F | critEq; + consts := (transpose{uList, data}) | (transpose{wList, weight}); + --sl = ideal singularLocus I; + --win = saturate(win, sl); + + makeB'InputFile( + dir, + AffVariableGroup => {xList}, + HomVariableGroup => {lamList}, + BertiniInputConfiguration => {"UseRegeneration" => 1, "TrackType" => 0, "PrintPathProgress" => 1000}, + B'Polynomials => win, + B'Constants => consts + ); + + dir = addSlash(dir); + runBertini(dir); + --readFile(dir); + numSols := null; + scanLines(ell -> (numSols = value ell; break), dir | "nonsingular_solutions"); + return numSols +) + +leftKernelUnitEDDegree = method(Options => leftKernelOptions) +leftKernelUnitEDDegree(List) := o -> (F) -> leftKernelWeightEDDegree( + F, + apply(#gens ring first F, i->randCC()), + apply(#gens ring first F, i->1_(ring first F)), + o +) + +leftKernelGenericEDDegree = method(Options => leftKernelOptions) +leftKernelGenericEDDegree(List) := o-> (F) -> leftKernelWeightEDDegree( + F, + apply(#gens ring first F, i->randCC()), + apply(#gens ring first F, i->randCC()), + o +) + +end + +-* +R = QQ[x,y]; +F = {x^2+y^2-1}; +(U,W) = ({.12, .23}, {.15, .331}) +GED = leftKernelWeightEDDegree(F, U, W); +GED = leftKernelGenericEDDegree F; +UED = leftKernelUnitEDDegree F; +*- \ No newline at end of file diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/Numerical.m2 b/M2/Macaulay2/packages/EuclideanDistanceDegree/Numerical.m2 new file mode 100644 index 00000000000..7c1a0b296a3 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/Numerical.m2 @@ -0,0 +1,553 @@ +--Projective formulation for intersections with linear spaces +rand := randomValue +(stageOne,stageTwo) = (1,2); + +--Assume ring is a complex inexact field +--G is a subset of F. + +NumericalComputationOptions = new Type of MutableHashTable + +parameterKeys = {"StartWeight", "TargetWeight", "StartData", "TargetData", "GammaVector"} +jacKeys= {"JacobianWitnessModel", "JacobianStartSubmodel", "JacobianTargetSubmodel"} +modelKeys = {"Model","WitnessModel", "StartSubmodel", "TargetSubmodel"} +degreeKeys = {"DegreeWitnessModel", "DegreeSubmodel"} + +bertiniKeys = {"BertiniStartFiberSolveConfiguration", "BertiniMembershipTestConfiguration", "BertiniSubstitute", "BertiniConstants"} +coordinateKeys = {"PrimalCoordinates", "HomogeneousVariableGroups", "AffineVariableGroups"} +planeKeys = {"Infinity", "PairGeneralHyperplaneList"} +variableKeys = {"JacobianVars", "GradientVars", "ScaleVars", "DataVars", "WeightVars"} + +directoryKeys = {"Directory"} +solutionKeys = {"TrackSolutions"} +outputKeys = {"OutputType", "FinerRestriction"} +fixValues = {"FixedData", "FixedWeight", "FixedSubmodel", "FixedJacobianSubmodel"} + +nocKeys = parameterKeys|jacKeys|modelKeys|degreeKeys|bertiniKeys|coordinateKeys|planeKeys|variableKeys +nocKeys = nocKeys|directoryKeys|solutionKeys|outputKeys|fixValues + +newNumericalComputationOptions = method(Options => { TempDirectory => null, Submodel => {}}) +newNumericalComputationOptions(List, List) := o -> (F, G) -> ( + NCO := new NumericalComputationOptions from apply(nocKeys, i -> i=>null); + + -- Temp directory for Bertini input file + if o.TempDirectory =!= null then dir := o.TempDirectory else dir = temporaryFileName(); + NCO#"Directory" = dir; + + -- Model keys + NCO#"Model"=F; + NCO#"WitnessModel"=G; + L := o.Submodel; + NCO#"StartSubmodel"=L; + NCO#"TargetSubmodel"=L; + + -- Degrees of models + NCO#"DegreeSubmodel"=L/degree/first; + NCO#"DegreeWitnessModel"=G/degree/first; + + -- Jacobians of models + NCO#"JacobianWitnessModel"=diff(matrix{gens ring first G}, transpose matrix{G}); + NCO#"JacobianTargetSubmodel"=NCO#"JacobianStartSubmodel"=diff(matrix{gens ring first G}, transpose matrix {L}); + + -- Data keys (used for other homotopy types) + numX:=#gens ring first G; + NCO#"TargetData"=NCO#"StartData"=apply(numX,i->random CC); + NCO#"TargetWeight"=apply(numX,i->1); + NCO#"StartWeight"=apply(numX,i->random CC); + --NCO#"GammaVector"=apply(numX-1,i->random CC); + + -- Variable group keys for Bertini + scan(bertiniKeys,i->NCO#i={}); + NCO#"HomogeneousVariableGroups"={gens ring first F}; + NCO#"AffineVariableGroups"={}; + NCO#"PrimalCoordinates"=gens ring first F; -- This is different when working with a parameterization + + -- Fixed keys (which model to use for which run) + NCO#"FixedData"="StartData"; + NCO#"FixedWeight"="StartWeight"; + NCO#"FixedSubmodel"="StartSubmodel"; + NCO#"FixedJacobianSubmodel"="JacobianStartSubmodel"; + + -- Keys defining variable block names + NCO#"JacobianVars"="jv"; + NCO#"GradientVars"="gv"; + NCO#"ScaleVars"="sv"; + NCO#"DataVars"="dv"; + NCO#"WeightVars"="wv"; + + -- Hyperplanes + NCO#"Infinity"=null; + NCO#"PairGeneralHyperplaneList"=null; + + -- Flags + NCO#"IsProjective"=false; -- currently not used + NCO#"OutputType"="Standard"; + NCO#"FinerRestriction"={}; + + return NCO +) + +possibleHT = {"Weight", "Data", "Submodel"} +stageOne=1 +stageTwo=2 +ht="Weight" +isStageOne=true +isStageTwo=true + +homotopyEDDegree = method() +homotopyEDDegree(NumericalComputationOptions, String, Boolean, Boolean) := (NCO, ht, isStageOne, isStageTwo) -> ( + if not member(ht,possibleHT) then error("Argument 1 is not in "|toString possibleHT); + weightHomotopy := ht === possibleHT_0; + dataHomotopy := ht === possibleHT_1; + submodelHomotopy := ht === possibleHT_2; + + if not fileExists NCO#"Directory" then mkdir NCO#"Directory"; + + jacL := NCO#(NCO#"FixedJacobianSubmodel"); + L := NCO#(NCO#"FixedSubmodel"); + startWeight := NCO#"StartWeight"; + targetWeight := NCO#"TargetWeight"; + startData := NCO#"StartData"; + targetData := NCO#"TargetData"; + + -- Symbol names for special variables used in the homotopy machinery/inputs + (lagMult,numerHB,denomQ,primal,tWeight) := ("lagMult","numerHB","denomQ","primal","tWeight"); + (jv,gv,sv) := (NCO#"JacobianVars",NCO#"GradientVars",NCO#"ScaleVars"); + (dv,wv) := (NCO#"DataVars",NCO#"WeightVars"); + + -- F is the "model" (target variety). V(G) ∩ V(L) is a CI inside V(F) ∩ V(L). + (F,G,startL,targetL,jacG) := (NCO#"Model",NCO#"WitnessModel",NCO#"StartSubmodel",NCO#"TargetSubmodel",NCO#"JacobianWitnessModel"); + + -- Optional gamma vector for randomization (not used here; keep for extensibility) + -- randomGamma:=NCO#"GammaVector"; + + -- List of primal coordinates (ambient variables) + xList := NCO#"PrimalCoordinates"; + + -- Ensure a hyperplane at infinity section is set; create it if missing. + if NCO#"Infinity" === null then NCO#"Infinity" = makeB'Section(xList, NameB'Section=>"HX"); + + vRing := (numVars,s,kk) -> kk[apply(numVars,i->s|i)]; + + -- Build extensic symbolic ring + nc := #xList; + kk0 := QQ; + extrinsicRing := kk0[flatten transpose apply(#G+#L,i->apply(nc,j->jv|i|"v"|j))]; + scan({sv,gv,dv,wv}, {#G+#L+1,nc,nc,nc}, (s,numVars) -> extrinsicRing = extrinsicRing ** vRing(numVars,s,kk0)); + + symbolicScaleMatrix := basis({0,1,0,0,0},extrinsicRing); -- the sv's + symbolicGradient := basis({0,0,1,0,0},extrinsicRing); -- the gv's + symbolicJac := genericMatrix(extrinsicRing,#G+#L,nc); + + -- symbolicSystem is the model system: (Scale) * (Gradient || Jacobian) that we want to solve + -- to be specialized by pairing functions below. + symbolicSystem := symbolicScaleMatrix*(symbolicGradient||symbolicJac); + + jacZero:={}; + pairJac:={}; + pairRowFunction := (M1, M2, hom) -> ( + arg := flatten entries M1; + val := flatten entries M2; + scan(#arg, i -> if val_i==0 then ( + jacZero = jacZero|{arg_i=>0}; + symbolicSystem = sub(symbolicSystem, {arg_i=>0}) + ) + else pairJac = pairJac | {makeB'Section( + {val_i}, + B'NumberCoefficients=>{1}, + B'Homogenization=>hom, + NameB'Section=>arg_i + )} + ) + ); + + pairParameterFunction := (p0,p1,r1,r2,sym,bool) -> ( + if bool then apply(#p0, + --- for each i we encode an equation like sym_i = p_0_i* r1+p1_i*r2, + --- x = 2*r1+ 7*r2; Typically, r1 = 1-T and r2 = T. + i->makeB'Section( + {p0_i,p1_i}, + B'NumberCoefficients=>{r1,r2}, + NameB'Section=>sym_i + )) + else apply(#p0, + --- for each i we encode an equation like sym_i = p0_i * 1, + --- x = 2*1; + i->makeB'Section( + {p0_i}, + B'NumberCoefficients=>{1}, + NameB'Section=>sym_i + )) + ); + + -- Pair symbolic vars with their start/target values and set up the homotopy between them + weight := symbolicWeight := gens vRing(nc,wv,kk0); + data := symbolicData := gens vRing(nc,dv,kk0); + pairData := pairParameterFunction(startData,targetData, "(1-TData)","TData", symbolicData, dataHomotopy); + pairWeight := pairParameterFunction(startWeight,targetWeight, "(1-TWeight)","TWeight", symbolicWeight, weightHomotopy); + + kk2 := ring first startWeight; + topS := kk2[numerHB,denomQ,tWeight]; + (topNumerHB,topDenomQ,topTWeight) := toSequence flatten entries basis({1},topS); + + pairGradient := apply(#xList, i -> makeB'Section( + {xList_i}, + ContainsPoint=>{data_i}, + B'NumberCoefficients=>{weight_i}, + NameB'Section=>symbolicGradient_(0,i), + B'Homogenization=>"HX" + )); + + -- Homogenize and pair the Jacobian + jacLG := jacL||jacG; + kk3 := coefficientRing ring first F; + jacRing := kk3[gens ring first F|{"HX"}]; + HX := last gens jacRing; + + homogLG := homogenize(sub(matrix{L|G},jacRing), HX) // entries // flatten; + homogJac := matrix apply(numrows jacLG,i->apply(numcols jacLG,j->diff((gens jacRing)_j,homogLG_i))); + pairRowFunction(symbolicJac, homogJac, "HX"); + + -- Pair lagrange multipliers for column homogenization + (degSubmodel,degWitnessModel) := (NCO#"DegreeSubmodel",NCO#"DegreeWitnessModel"); + degAugJac := {1} | apply(degSubmodel|degWitnessModel, i -> i-1); + maxDegree := degAugJac//max; + degRescale := degAugJac/(i -> maxDegree-i); -- how much to rescale each column + bLagrangeVars := lagList := apply(#degRescale, i -> "L"|i); + rescaling := new MutableList from apply(#degRescale, i -> lagList_i); + + -- Build generic linear forms H(i,v,j) used to homogenize per-column degrees + --Homogenize cols by multiplying by a diagonal matrix of linear products on the left. + generalHyperplaneList:={}; + scan(#degRescale, i -> scan( + degRescale_i, + j -> ( + hg:="H"|i|"v"|j; --wants to be both + rescaling#i = (rescaling#i)|"*"|hg; + generalHyperplaneList = generalHyperplaneList|{hg} + ) + )); + + -- Reuse previously paired hyperplanes if provided; otherwise pair new ones now. + if NCO#"PairGeneralHyperplaneList"=!=null then pairGeneralHyperplanes:=NCO#"PairGeneralHyperplaneList" + else ( + pairGeneralHyperplanes=apply(#generalHyperplaneList, i-> + makeB'Section(xList|{"HX"}, + NameB'Section=>generalHyperplaneList_i) + ); + NCO#"PairGeneralHyperplaneList"=pairGeneralHyperplanes + ); + + -- Connect each scale slot with the corresponding homogenizing product. + pairScale := apply(flatten entries symbolicScaleMatrix, rescaling, (i,j) -> i=>j); + + bModelVars := gens ring first F|{"HX"}; + bPoly := homogLG|flatten entries symbolicSystem; + bConfiguration := {"TrackType"=>0, "PrintPathProgress"=>1000} | (NCO#"BertiniStartFiberSolveConfiguration"); + BF := pairData|pairWeight|pairJac|pairGradient|pairGeneralHyperplanes|pairScale; + + writeSolveInputFunction := (stage,nif) -> ( + if stage===stageOne then ( + PG:={"adfadfdf"}; + BC:={"TData"=>0,"TWeight"=>0} + ) + else if stage===stageTwo then ( + BC={}; + if dataHomotopy then PG={"TData"} + else if weightHomotopy then PG={"TWeight"} + ); + + if stage===stageOne then bConfiguration = bConfiguration | {"UseRegeneration"=>1}; + if stage===stageTwo then bConfiguration = bConfiguration | {"UseRegeneration"=>0}; + + makeB'InputFile( + NCO#"Directory", + NameB'InputFile=>nif, + HomVariableGroup=>{bLagrangeVars,bModelVars}, + BertiniInputConfiguration=>bConfiguration|{"ParameterHomotopy"=>stage}, + B'Polynomials=>bPoly, + B'Constants=>jacZero|BC, + ParameterGroup=>PG, + B'Functions=>BF + ) + ); + + --our solution file will always be member points. + criticalPointName := "criticalPointFile"; + runSolveInputFunction := (stage,nif) -> ( + writeSolveInputFunction(stage,nif); + + if stage==stageTwo then( + writeParameterFile(NCO#"Directory",{0},NameParameterFile=>"start_parameters"); + writeParameterFile(NCO#"Directory",{1},NameParameterFile=>"final_parameters") + ); + + runBertini(NCO#"Directory",NameB'InputFile=>nif); + readFile(NCO#"Directory"); + + if stage==stageOne then( + moveB'File(NCO#"Directory","bertini_session.log","stageOne_log",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions","stageOne_solutions",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions","start",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions","member_points",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions",criticalPointName,CopyB'File=>true) + ); + + if stage==stageTwo then( + moveB'File(NCO#"Directory","bertini_session.log","stageTwo_log",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions","stageTwo_solutions",CopyB'File=>true); + moveB'File(NCO#"Directory","main_data","stageTwo_main_data",CopyB'File=>true); + --moveB'File(NCO#"Directory","nonsingular_solutions","start",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions","member_points",CopyB'File=>true); + moveB'File(NCO#"Directory","nonsingular_solutions",criticalPointName,CopyB'File=>true) + ); + ); + + ------------------------------------------------------------------------ + -- (FUNCTIONS 3) Membership tests and incidence matrices. + -- We create and run two inputs per hypersurface: + -- - TrackType=1 (positive-dimensional solve) + -- - TrackType=3 (membership test to decide if points lie on the hypersurface) + -- Output: an incidence matrix (list of lists) from importIncidenceMatrix. + ------------------------------------------------------------------------ + ttOne:=1; + ttThree:=3; + nameFileFunction:=(stage,case,indexCase,hypersurface,theTrackType)->("input_first_MT_"|case|"_"|indexCase|"_"|theTrackType); + + writeIsMembershipFunction := (stage,case,indexCase,hypersurface,theTrackType) -> ( + nif := nameFileFunction(stage,case,indexCase,hypersurface,theTrackType); + if stage===stageOne then BC:={"TData"=>0,"TWeight"=>0}; + if stage===stageTwo then BC={"TData"=>1,"TWeight"=>1}; + if not member(stage,{1,2}) then error"stage is in {1,2}"; + + makeB'InputFile( + NCO#"Directory", + NameB'InputFile=>nif, + AffVariableGroup=>flatten flatten {bLagrangeVars,bModelVars}, + BertiniInputConfiguration=>bConfiguration|{"TrackType"=>theTrackType}, + B'Polynomials=>{hypersurface}, + B'Constants=>jacZero|BC, + --ParameterGroup=>PG, + B'Functions=>BF + ) + ); + + -- Example test: isMembershipFunction(stageOne, "TT", 0, "x1*x2-x3*x4") + isMembershipFunction:=(stage,case,indexCase,hypersurface)->( + --Pos dim solve TrackType=>1 + writeIsMembershipFunction(stage,case,indexCase,hypersurface,ttOne); + nif:=nameFileFunction(stage,case,indexCase,hypersurface,ttOne); + runBertini(NCO#"Directory",NameB'InputFile=>nif); + moveB'File(NCO#"Directory","bertini_session.log","bertini_session_"|nif|".log",CopyB'File => false); + + --MT TrackType=>3 + writeIsMembershipFunction(stage,case,indexCase,hypersurface,ttThree); + nif=nameFileFunction(stage,case,indexCase,hypersurface,ttThree); + runBertini(NCO#"Directory",NameB'InputFile=>nif); + moveB'File(NCO#"Directory","bertini_session.log","bertini_session_"|nif|".log",CopyB'File => false); + + outIM := importIncidenceMatrix(NCO#"Directory"); + return outIM + ); + + -- Filter points using incidence matrix + filterSolutionFunction:=(nsf,kp,ns,numCoords)->( + -- member_points structure assumption: first line header, then groups + -- of (1 + numCoords) lines per solution. Only selected solutions are kept. + + firstLine := true; + countSol := 0; + countLine := 0; + groupSize := 1+numCoords; + isSelected := null; + + sf := openOut(NCO#"Directory"|"/"|nsf); + scanLineSolutionFunction := (ell) -> ( + if firstLine then ( + firstLine=false; + sf << toString(#kp) << endl + ) + else if countSol < ns then ( + if countLine==0 then isSelected = member(countSol,kp); + countLine = countLine+1; + if isSelected then sf << ell << endl; + if countLine == groupSize then ( + countLine = 0; + countSol = countSol+1; + ) + ) + ); + scanLines(scanLineSolutionFunction,(NCO#"Directory")|"/"|"member_points"); + close sf; + return (nsf) + ); + + -- positionFilterFunction: + -- stage : stage1 or stage2 (affects BC and file naming) + -- case : label (e.g. "SaturateH", "IntersectF") + -- indexCase : index within the polyList + -- hypersurface : polynomial to test + -- bin : "typeA" (drop if ON hypersurface) or "typeB" (keep only if ON) + positionFilterFunction := (stage,case,indexCase,hypersurface,bin) -> ( + if bin==="typeA" then isOffHypersurface := (m->(m==={})) + else if bin==="typeB" then isOffHypersurface = (m->(m=!={})) + else error"last argument is typeA or typeB"; + + imMT := isMembershipFunction(stage,case,indexCase,hypersurface); + kp := {}; + scan(#imMT, i -> if isOffHypersurface(imMT_i) then kp=kp|{i}); + + ns:= #imMT; + (nsf,nc) := ("filterFile", #flatten {bLagrangeVars,bModelVars}); + + filterSolutionFunction("filterFile",kp,ns,nc); + moveB'File(NCO#"Directory","filterFile","member_points",CopyB'File=>true); + return #kp + ); + + -- stageEDDegBound holds a label and the counts after stage1 and stage2 filtering. + stageEDDegBound:=new MutableList from {"GEDvUED",null,null}; + + -- Drop solutions that lie on offPolyList + runSaturateUnionFunction:=(polyList,stage)->( + (case,bin) := ("SaturateH","typeA"); + scan(#polyList,i->( + stageEDDegBound#stage = positionFilterFunction(stageOne,case,i,polyList_i,bin); + )); + ); + + -- Keep only solutions that lie on onPolyList + runRestrictIntersectionFunction:=(polyList,stage)->( + (case, bin) := ("IntersectF", "typeB"); + scan(#polyList, i -> ( + stageEDDegBound#stage = positionFilterFunction(stageOne, case, i, polyList_i, bin); + )); + ); + + runComputationStage:=(stage,offPolyList,onPolyList)->( + if stage==stageOne then runSolveInputFunction(stageOne,"input_first_solve") + else runSolveInputFunction(stageTwo,"input_second_solve"); + + if stage===stageOne or NCO#"OutputType"=!="TestHomotopyConjectureGEDvUED" then runSaturateUnionFunction(offPolyList,stage); + runRestrictIntersectionFunction(onPolyList,stage); + if stage===stageTwo and NCO#"FinerRestriction"=!={} then runRestrictIntersectionFunction(NCO#"FinerRestriction",stage); + ); + + -- Polynomials that should vanish (onPolyList) or should not vanish (offPolyList). + -- offPolyList includes: + -- - HX (hyperplane at infinity) + -- - "L0" (first Lagrange multiplier) to remove trivial solutions at infinity + -- - all generic hyperplanes used in rescaling + offPolyList := {HX,"L0"}|((pairGeneralHyperplanes/(i->i#NameB'Section))); + + -- onPolyList are the homogenized F equations (in the Jacobian ring with HX) + onPolyList := F/(i->homogenize(sub(i,jacRing),HX)); + + -- Execute the requested stages, then return the count for the last executed stage. + if isStageOne then runComputationStage(stageOne,offPolyList,onPolyList); + if isStageTwo then runComputationStage(stageTwo,offPolyList,onPolyList); + if isStageTwo then return stageEDDegBound#2 else if isStageOne then return stageEDDegBound#1 +) + +numericalOptions = { TempDirectory => null, Submodel => {} } +numericWeightEDDegree = method(Options => numericalOptions) +numericWeightEDDegree(List, List, List, List) := o -> (F, G, data, weight) -> ( + NCO := newNumericalComputationOptions(F, G, o); + NCO#"StartWeight" = weight; + NCO#"TargetData" = NCO#"StartData" = data; + homotopyEDDegree(NCO, "Weight", true, false) +) + +numericUnitEDDegree = method(Options => numericalOptions) +numericUnitEDDegree(List, List) := o -> (F, G) -> numericWeightEDDegree( + F, G, + apply(#gens ring first F, i->randCC()), + apply(#gens ring first F, i->1_(ring first F)), + o +) + +numericGenericEDDegree = method(Options => numericalOptions) +numericGenericEDDegree(List, List) := o -> (F, G) -> numericWeightEDDegree( + F, G, + apply(#gens ring first F, i->randCC()), + apply(#gens ring first F, i->randCC()), + o +) + +aEDOptions = { TempDirectory => null, Tolerance => 1e-6, Submodel => {}, SampleGenerator => null } +averageNumericEDDegree = method(Options => aEDOptions) +averageNumericEDDegree(List, List, ZZ) := o -> (F, G, n) -> ( + if o.Submodel =!= null then L := o.Submodel else L = {}; + if o.TempDirectory =!= null then dir := o.TempDirectory else dir = temporaryFileName(); + if o.SampleGenerator =!= null then + SampleGenerator := o.SampleGenerator + else + SampleGenerator = () -> apply(#gens ring first F, i -> randomRR()); + + -- Initial run + R := ring first F; + NCO := newNumericalComputationOptions(F, G, Submodel => o.Submodel, TempDirectory => o.TempDirectory); + NCO#"StartData" = SampleGenerator(); + NCO#"TargetData" = apply(#gens R, i -> randomRR()); + NCO#"StartWeight" = NCO#"TargetWeight" = apply(#gens R, i -> 1_R); + NCO#"BertiniStartFiberSolveConfiguration" = {"FinalTol" => 1e-12}; + homotopyEDDegree(NCO, "Data", true, false); + + -- path to sampled data and average real critical points + avgEDDeg := 0; + for i to n-1 do ( + NCO#"StartData" = SampleGenerator(); + homotopyEDDegree(NCO, "Data", false, true); + critPoints := importMainDataFile(NCO#"Directory", NameMainDataFile => "stageTwo_main_data"); + avgEDDeg = avgEDDeg + #realPoints(critPoints, Tolerance => o.Tolerance); + ); + numeric(avgEDDeg/n) +) + +end + +restart +loadPackage "EuclideanDistanceDegree" +loadPackage "Probability" +n = 100; +R = QQ[x,y]; +F = G = {x^2 + 4*y^2 - 4}; +Z = normalDistribution(); +sampleGen = () -> apply(#gens R, i -> random Z); +averageNumericEDDegree({F,G}, sampleGen, n) + +restart +loadPackage "EuclideanDistanceDegree" +R = QQ[x,y]; +F = G = {x^2 + 4*y^2 - 4}; +dir = temporaryFileName(); +NCO = newNumericalComputationOptions(dir, F, G); +NCO#"BertiniStartFiberSolveConfiguration" = {"FinalTol" => 1e-12}; + +NCO#"StartData" = NCO#"TargetData" = apply(#gens R, i -> random(RR)); +homotopyEDDegree(NCO, "Data", true, false); + +NCO#"StartData" = apply(#gens R, i -> random(RR)); +homotopyEDDegree(NCO, "Data", false, true); + +limitPoints := importMainDataFile(NCO#"Directory", NameMainDataFile => "stageTwo_main_data"); + +cnt1 = number(limitPoints, P -> ( + coords := drop(drop(P#Coordinates, #G+1), -1); + HX := last P#Coordinates; + affineCoords := apply(coords, x -> x / HX); + all(affineCoords, x -> abs(imaginaryPart x) <= 1e-4) +)); +cnt2 = realPoints(limitPoints); + +affinePoints = apply(limitPoints, P -> ( + coords = drop(drop(P#Coordinates, #G+1), -1); + HX = last P#Coordinates; + apply(coords, x -> x / HX) +)); +netList limitPoints + +netList readDirectory dir +scanLines(l -> print(l), dir | "/stageTwo_solutions"); + +numericGenericEDDegree(F, G, TempDirectory => dir) +averageNumericEDDegree({F,G}, sampleGen, 100); \ No newline at end of file diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/Parameterization.m2 b/M2/Macaulay2/packages/EuclideanDistanceDegree/Parameterization.m2 new file mode 100644 index 00000000000..f9d65d40564 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/Parameterization.m2 @@ -0,0 +1,93 @@ +rand := () -> random QQ; + +parameterOptions = { UseMonodromy => false } +parameterizedWeightEDDegree = method(Options => parameterOptions) +parameterizedWeightEDDegree(List, List, List) := o -> (F, U, W) -> ( + -- print("U: " | toString U | " W: " | toString W); + R := ring first F; + coef := coefficientRing R; + numX := #gens R; + n := #F; + + S := R ** coef[apply(n - numX, i->"L"|i)]; + xList := flatten entries basis({1,0}, S); + lamList := flatten entries basis({0,1}, S); + M := sub(matrix{F}, S); + + -- Find a spanning set for ker(jacM) + jacM := transpose sub(matrix makeJac(apply(F, i->sub(i,S)), xList), S); + assert(rank jacM == numX); -- ensure dim X = d + columnVectors := gens kernel jacM; -- bottleneck + evalColumnVectors := sub(columnVectors, apply(xList, x -> x => random(1, 100))); + + A := matrix for i to n-1 list {}; + count := 0; + selectColumns := {}; + while #selectColumns < n-numX and count < #columnVectors do ( + B := A | matrix(evalColumnVectors_count); + if rank B > rank A then ( + A = B; + selectColumns = selectColumns | {count} + ); + count = count + 1 + ); + assert(#selectColumns == n - numX); + unscaledMatrix := columnVectors_selectColumns; + outMatrix := unscaledMatrix * transpose matrix{lamList}; + + gradObjective := 2 * diagonalMatrix(W) * transpose(M - matrix{U}); -- 2w(f - u) + criteqs := outMatrix - gradObjective; + + if o.UseMonodromy then ( + sys := polySystem flatten entries criteqs; + sols := sparseMonodromySolve sys; + number(points sols, p -> status p === Regular) + ) + else ( + I := ideal(criteqs); + degree I + ) +) + +parameterizedUnitEDDegree = method(Options => parameterOptions) +parameterizedUnitEDDegree(List) := o -> (F) -> parameterizedWeightEDDegree( + F, + apply(#F, i -> rand()), + apply(#F, i -> 1_(ring first F)), + o +) + +parameterizedGenericEDDegree = method(Options => parameterOptions) +parameterizedGenericEDDegree(List) := o -> (F) -> parameterizedWeightEDDegree( + F, + apply(#F, i -> rand()), + apply(#F, i -> rand()), + o +) + +end + +loadPackage "MonodromySolver" +loadPackage "EuclideanDistanceDegree" +setRandomSeed(123456); +makeJac = (system,unknowns) -> for i in system list for j in unknowns list diff(j,i) +R = QQ[x1,x2,x3,x4,x5]; +F = {x1^2+x4^2, x2^2+x5^2, x3^2+1, x1*x2+x4*x5, x1*x3+x4, x2*x3+x5}; +U = {1,2,3,4,5,6}; +W = {1,1,1,1,1,1}; +parameterizedWeightEDDegree(F,U,W) +parameterizedWeightEDDegree(F,U,W, UseMonodromy => true) + +restart +loadPackage "EuclideanDistanceDegree" +makeJac = (system,unknowns) -> for i in system list for j in unknowns list diff(j,i) +R = QQ[x]; +F = {x^2 - 1, x^3 - x}; +parameterizedGenericEDDegree F + +n = #F; +numX = #gens R; +S = QQ[gens R] ** QQ[y_1..y_(n-numX), u_1..u_n]; +M = sub(matrix{F}, S); +imageModel = eliminate(support M, ideal(M - matrix{for i from 1 to n list u_i})); +determinantalUnitEDDegree((sub(imageModel, QQ[support imageModel]))_*) \ No newline at end of file diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Euclidean__Distance__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Euclidean__Distance__Degree.out new file mode 100644 index 00000000000..df5d101fdf3 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Euclidean__Distance__Degree.out @@ -0,0 +1,52 @@ +-- -*- M2-comint -*- hash: 399917589567742851 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + y^2 - 1}; + +i3 : determinantalUnitEDDegree(F) + +o3 = 2 + +i4 : leftKernelUnitEDDegree(F) +Warning: The HomVariableGroup is written first and then the AffVariableGroup is written second. + +o4 = 2 + +i5 : R = QQ[x1,x2,x3,x4]; + +i6 : F = G = {det genericMatrix(R,2,2)}; + +i7 : numericGenericEDDegree(F, G) +-- warning: experimental computation over inexact field begun +-- results not reliable (one warning given per session) + +o7 = 6 + +i8 : R = QQ[x1,x2,x3,x4,x5,x6]; + +i9 : F = (minors(2, genericMatrix(R,3,2)))_*; + +i10 : G = drop(F, -1); + +i11 : #G == codim ideal F; + +i12 : numericGenericEDDegree(F, G) + +o12 = 10 + +i13 : R = QQ[x1,x2,x3,x4,x5,x6]; + +i14 : F = (minors(2, genericMatrix(R,3,2)))_*; + +i15 : G = drop(F, -1); + +i16 : NCO = newNumericalComputationOptions(F, G); + +i17 : NCO#"TargetWeight" = apply(#gens R, i->1); + +i18 : homotopyEDDegree(NCO, "Weight", true, true) + +o18 = 2 + +i19 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Return__Critical__Ideal.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Return__Critical__Ideal.out new file mode 100644 index 00000000000..057afcf5b52 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Return__Critical__Ideal.out @@ -0,0 +1,16 @@ +-- -*- M2-comint -*- hash: 17266523904631417738 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + y^2 - 1}; + +i3 : (U,W) = ({12, 23}, {15, 331}); + +i4 : ICP = symbolicWeightEDDegree(F, U, W, ReturnCriticalIdeal => true) + + 2 2 +o4 = ideal (316x*y - 7613x + 180y, x + y - 1) + +o4 : Ideal of R + +i5 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Temp__Directory.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Temp__Directory.out new file mode 100644 index 00000000000..199f68e88a4 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Temp__Directory.out @@ -0,0 +1,14 @@ +-- -*- M2-comint -*- hash: 17265651732443607604 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + y^2 - 1}; + +i3 : dir = temporaryFileName(); + +i4 : GED = leftKernelGenericEDDegree(F, TempDirectory => dir) +Warning: The HomVariableGroup is written first and then the AffVariableGroup is written second. + +o4 = 4 + +i5 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Use__Monodromy.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Use__Monodromy.out new file mode 100644 index 00000000000..7fb65a77989 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/___Use__Monodromy.out @@ -0,0 +1,11 @@ +-- -*- M2-comint -*- hash: 13689930562824927558 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + 1, x * y, y - 1}; + +i3 : GED = parameterizedGenericEDDegree(F, UseMonodromy => true) + +o3 = 7 + +i4 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_average__Numeric__E__D__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_average__Numeric__E__D__Degree.out new file mode 100644 index 00000000000..40b5ca3caf5 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_average__Numeric__E__D__Degree.out @@ -0,0 +1,15 @@ +-- -*- M2-comint -*- hash: 15222729530691853413 + +i1 : R = QQ[x,y]; + +i2 : F = G = {x^2 + y^2 - 1}; + +i3 : sampleGen = () -> apply(#gens R, i -> random(RR)); + +i4 : aED = averageNumericEDDegree(F, G, 10, SampleGenerator => sampleGen) + +o4 = 1.9 + +o4 : RR (of precision 53) + +i5 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_homotopy__E__D__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_homotopy__E__D__Degree.out new file mode 100644 index 00000000000..36b4534d4aa --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_homotopy__E__D__Degree.out @@ -0,0 +1,23 @@ +-- -*- M2-comint -*- hash: 15398408180617996408 + +i1 : R = QQ[x,y]; + +i2 : F = G = {x^2+y^2-1}; + +i3 : NCO = newNumericalComputationOptions(F, G); + +i4 : NCO#"TargetWeight" = apply(#gens R, i -> random RR); + +i5 : GED = homotopyEDDegree(NCO, "Weight", true, true) +-- warning: experimental computation over inexact field begun +-- results not reliable (one warning given per session) + +o5 = 4 + +i6 : NCO#"TargetWeight" = apply(#gens R, i -> 1_R); + +i7 : UED = homotopyEDDegree(NCO, "Weight", false, true) + +o7 = 2 + +i8 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_left__Kernel__Weight__E__D__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_left__Kernel__Weight__E__D__Degree.out new file mode 100644 index 00000000000..9d4168c4a1f --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_left__Kernel__Weight__E__D__Degree.out @@ -0,0 +1,24 @@ +-- -*- M2-comint -*- hash: 10893047696075686640 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + y^2 - 1}; + +i3 : (U,W) = ({.12, .23}, {.15, .331}); + +i4 : UED = leftKernelUnitEDDegree F +Warning: The HomVariableGroup is written first and then the AffVariableGroup is written second. + +o4 = 2 + +i5 : GED = leftKernelGenericEDDegree F +Warning: The HomVariableGroup is written first and then the AffVariableGroup is written second. + +o5 = 4 + +i6 : GED = leftKernelWeightEDDegree(F, U, W) +Warning: The HomVariableGroup is written first and then the AffVariableGroup is written second. + +o6 = 4 + +i7 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_new__Numerical__Computation__Options.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_new__Numerical__Computation__Options.out new file mode 100644 index 00000000000..7ab75ed942b --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_new__Numerical__Computation__Options.out @@ -0,0 +1,19 @@ +-- -*- M2-comint -*- hash: 12902145464090925420 + +i1 : R = QQ[x,y]; + +i2 : F = G = {x^2 + y^2 - 1}; + +i3 : dir = temporaryFileName(); + +i4 : NCO = newNumericalComputationOptions(F, G, TempDirectory => dir); + +i5 : NCO#"TargetWeight" = apply(#gens R, i -> 1_R); + +i6 : UED = homotopyEDDegree(NCO, "Weight", true, true) +-- warning: experimental computation over inexact field begun +-- results not reliable (one warning given per session) + +o6 = 2 + +i7 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_numeric__Weight__E__D__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_numeric__Weight__E__D__Degree.out new file mode 100644 index 00000000000..dd339bcf08d --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_numeric__Weight__E__D__Degree.out @@ -0,0 +1,23 @@ +-- -*- M2-comint -*- hash: 6775569773143953516 + +i1 : R = QQ[x,y]; + +i2 : F = G = {x^2+y^2-1}; + +i3 : (U,W) = ({.12, .23}, {.15, .331}); + +i4 : UED = numericUnitEDDegree(F, G) + +o4 = 2 + +i5 : GED = numericGenericEDDegree(F, G) +-- warning: experimental computation over inexact field begun +-- results not reliable (one warning given per session) + +o5 = 4 + +i6 : GED = numericWeightEDDegree(F, G, U, W) + +o6 = 4 + +i7 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_parameterized__Weight__E__D__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_parameterized__Weight__E__D__Degree.out new file mode 100644 index 00000000000..00f0fd2b769 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_parameterized__Weight__E__D__Degree.out @@ -0,0 +1,21 @@ +-- -*- M2-comint -*- hash: 4180291041213606841 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + 1, x * y, y - 1}; + +i3 : (U,W) = ({12, 23, 25}, {15, 331, 1}); + +i4 : UED = parameterizedUnitEDDegree F + +o4 = 7 + +i5 : GED = parameterizedGenericEDDegree F + +o5 = 7 + +i6 : GED = parameterizedWeightEDDegree(F, U, W) + +o6 = 7 + +i7 : diff --git a/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_symbolic__Weight__E__D__Degree.out b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_symbolic__Weight__E__D__Degree.out new file mode 100644 index 00000000000..496e55a6d70 --- /dev/null +++ b/M2/Macaulay2/packages/EuclideanDistanceDegree/examples/_symbolic__Weight__E__D__Degree.out @@ -0,0 +1,24 @@ +-- -*- M2-comint -*- hash: 8919791488280722637 + +i1 : R = QQ[x,y]; + +i2 : F = {x^2 + y^2 - 1}; + +i3 : (U,W) = ({12, 23}, {15, 331}); + +i4 : UED = determinantalUnitEDDegree F + +o4 = 2 + +i5 : GED = determinantalGenericEDDegree F + +o5 = 4 + +i6 : ICP = symbolicWeightEDDegree(F, U, W, ReturnCriticalIdeal => true) + + 2 2 +o6 = ideal (316x*y - 7613x + 180y, x + y - 1) + +o6 : Ideal of R + +i7 :