One Hat Cyber Team
Your IP :
216.73.216.14
Server IP :
194.44.31.54
Server :
Linux zen.imath.kiev.ua 4.18.0-553.77.1.el8_10.x86_64 #1 SMP Fri Oct 3 14:30:23 UTC 2025 x86_64
Server Software :
Apache/2.4.37 (Rocky Linux) OpenSSL/1.1.1k
PHP Version :
5.6.40
Buat File
|
Buat Folder
Eksekusi
Dir :
~
/
usr
/
share
/
Macaulay2
/
Edit File:
NumericalImplicitization.m2
newPackage("NumericalImplicitization", Headline => "numerical invariants of images of varieties", AuxiliaryFiles => true, Version => "2.2.0", Date => "November 24, 2020", Authors => { {Name => "Justin Chen", Email => "justin.chen@math.gatech.edu", HomePage => "https://people.math.gatech.edu/~jchen646/"}, {Name => "Joe Kileel", Email => "jkileel@math.princeton.edu", HomePage => "https://web.math.princeton.edu/~jkileel/"} }, Keywords => {"Numerical Algebraic Geometry"}, PackageExports => {"NumericalAlgebraicGeometry"}, Certification => { "journal name" => "The Journal of Software for Algebra and Geometry", "journal URI" => "https://msp.org/jsag/", "article title" => "Numerical implicitization", "acceptance date" => "11 April 2019", "published article URI" => "https://msp.org/jsag/2019/9-1/p07.xhtml", "published article DOI" => "10.2140/jsag.2019.9.55", "published code URI" => "https://msp.org/jsag/2019/9-1/jsag-v9-n1-x07-NumericalImplicitization.m2", "release at publication" => "2f801d123692462f4a65ccb135d411be425c28bd", -- git commit number in hex "version at publication" => "2.1.0", "volume number" => "9", "volume URI" => "https://msp.org/jsag/2019/9-1/" } ) export { "numericalSourceSample", "numericalImageSample", "numericalEval", "numericalNullity", "Precondition", "SVDGap", "numericalImageDim", "numericalHilbertFunction", "ConvertToCone", "NumericalInterpolationTable", "hilbertFunctionArgument", "hilbertFunctionValue", "UseSLP", "imagePoints", "interpolationBasis", "interpolationSVD", "interpolationMatrix", "extractImageEquations", "AttemptZZ", "numericalImageDegree", "pseudoWitnessSet", "DoRefinements", "DoTraceTest", "MaxAttempts", "MaxPoints", "MaxThreads", "Repeats", "TraceThreshold", -- "Endgame", "PseudoWitnessSet", "isCompletePseudoWitnessSet", "sourceEquations", "sourceSlice", "generalCombinations", "imageSlice", "witnessPointPairs", "isOnImage", "realPoint", "optimizeNelderMead", "lineSearch", "Initial" } -- Point sampling code load "NumericalImplicitization/approxPoint.m2" -- software options: default is M2engine throughout NumericalInterpolationTable = new Type of HashTable NumericalInterpolationTable.synonym = "numerical interpolation table" globalAssignment NumericalInterpolationTable net NumericalInterpolationTable := T -> ( (net ofClass class T | ", indicating") || ("the space of degree " | (toString T.hilbertFunctionArgument) | " forms in the ideal of the image has dimension " | (toString T.hilbertFunctionValue)) ) PseudoWitnessSet = new Type of HashTable PseudoWitnessSet.synonym = "pseudo-witness set" globalAssignment PseudoWitnessSet net PseudoWitnessSet := W -> ( (net ofClass class W | ", indicating") || ("the image has degree " | (toString W.degree)) ) checkRings = method(Options => {symbol ConvertToCone => false}) -- checks if the rings of F and I agree and have floating point arithmetic, and converts F, I, pts to the affine cone if ConvertToCone is false checkRings (Matrix, Ideal, List) := Sequence => opts -> (F, I, pts) -> ( k := coefficientRing ring I; if not numrows F == 1 then error "Expected map to be given by a 1-row matrix of polynomials"; if not ring F === ring I then error "Expected same rings for ideal and map"; if not instance(class(1_k), InexactFieldFamily) then error "Expected coefficient field with floating point arithmetic"; if opts.ConvertToCone then ( JJ := getSymbol "JJ"; S := k(monoid[append(gens ring I, JJ)]); toS := map(S, ring I); ((last gens S)*(toS F | matrix{{1_S}}), toS(I), pts/(p -> point{append(p#Coordinates, 1_k)})) ) else (F, I, pts) ) numericalSourceSample = method(Options => {Software => M2engine}) numericalSourceSample (Ideal, Thing, ZZ) := List => opts -> (I, W, sampleSize) -> ( R := ring I; if I == 0 then ( k := coefficientRing R; return (entries random(k^(sampleSize), k^(#gens R)))/(p -> {p})/point; ); -- samplePoints := if instance(W, Point) and not I.cache.?WitnessSet then ( -- d := first numericalDimensions(vars R, I, W); -- squaredUpSource := randomSlice(gens I, #gens R - d, {}); -- startSys := squaredUpSource | randomSlice(vars R, d, {W, "source"}); -- flatten apply(sampleSize, i -> track(startSys, squaredUpSource | randomSlice(vars R, d, {}), {W}, opts))) -- track not reliable, often fails if instance(opts.Software, FunctionClosure) then return apply(sampleSize, i -> (opts.Software) I); if not I.cache.?WitnessSet then I.cache.WitnessSet = if instance(W, WitnessSet) then W else first components(numericalIrreducibleDecomposition(I, opts)); samplePoints := apply(sampleSize, i -> sample I.cache.WitnessSet); if precision R <= precision ring samplePoints#0#Coordinates#0 then samplePoints else refine(polySystem(I_*), samplePoints, Bits => precision R) ) numericalSourceSample (Ideal, WitnessSet) := List => opts -> (I, W) -> numericalSourceSample(I, W, 1, opts) numericalSourceSample (Ideal, Point) := List => opts -> (I, p) -> numericalSourceSample(I, p, 1, opts) numericalSourceSample (Ideal, ZZ) := List => opts -> (I, sampleSize) -> numericalSourceSample(I, null, sampleSize, opts) numericalSourceSample Ideal := List => opts -> I -> numericalSourceSample(I, 1, opts) numericalImageSample = method(Options => options numericalSourceSample) numericalImageSample (Matrix, Ideal, List, ZZ) := List => opts -> (F, I, pts, sampleSize) -> ( samplePoints := if #pts > 0 then numericalSourceSample(I, pts#0, sampleSize-#pts, opts) else numericalSourceSample(I, sampleSize, opts); numericalEval(F, samplePoints, false) /point ) numericalImageSample (Matrix, Ideal, ZZ) := List => opts -> (F, I, sampleSize) -> numericalImageSample(F, I, {}, sampleSize, opts) numericalImageSample (Matrix, Ideal) := List => opts -> (F, I) -> numericalImageSample(F, I, {}, 1, opts) numericalImageSample (List, Ideal, List, ZZ) := List => opts -> (F, I, pts, sampleSize) -> numericalImageSample(matrix{F}, I, pts, sampleSize, opts) numericalImageSample (List, Ideal, ZZ) := List => opts -> (F, I, sampleSize) -> numericalImageSample(matrix{F}, I, {}, sampleSize, opts) numericalImageSample (List, Ideal) := List => opts -> (F, I) -> numericalImageSample(matrix{F}, I, {}, 1, opts) numericalImageSample (RingMap, Ideal, List, ZZ) := List => opts -> (F, I, pts, sampleSize) -> numericalImageSample(F.matrix, I, pts, sampleSize, opts) numericalImageSample (RingMap, Ideal, ZZ) := List => opts -> (F, I, sampleSize) -> numericalImageSample(F.matrix, I, {}, sampleSize, opts) numericalImageSample (RingMap, Ideal) := List => opts -> (F, I) -> numericalImageSample(F.matrix, I, {}, 1, opts) numericalEval = method() numericalEval (Matrix, List, Boolean) := List => (F, upstairsPoints, includeUpstairs) -> ( -- returns a list of either matrices, or pairs of the form (Point, Matrix) evalPts := upstairsPoints/(p -> (p, sub(F, matrix p))); if includeUpstairs then evalPts else evalPts/last ) numericalDimensions = method(Options => options numericalSourceSample) numericalDimensions (Matrix, Ideal, Point) := List => opts -> (F, I, p) -> ( --outputs {dim(V(I)), dim(F(V(I))} (F, I, p) = checkRings(F, I, {p}); p0 := 1/norm(2, matrix p#0)*(matrix p#0); dF := sub(transpose jacobian F, p0); if I == 0 then return {#gens ring I, #gens ring I - numericalNullity(dF, false)}; sourceJacobian := sub(transpose jacobian I, p0); sourceDim := numericalNullity(sourceJacobian, false); {sourceDim, sourceDim - numericalNullity(sourceJacobian || dF, false)} ) numericalDimensions (Matrix, Ideal) := ZZ => opts -> (F, I) -> numericalDimensions(F, I, first numericalSourceSample(I, Software => opts.Software), opts) numericalImageDim = method(Options => options numericalSourceSample) numericalImageDim (Matrix, Ideal, Point) := ZZ => opts -> (F, I, p) -> last numericalDimensions(F, I, p, opts) numericalImageDim (Matrix, Ideal) := ZZ => opts -> (F, I) -> last numericalDimensions(F, I, opts) numericalImageDim (List, Ideal, Point) := ZZ => opts -> (F, I, p) -> last numericalDimensions(matrix{F}, I, p, opts) numericalImageDim (List, Ideal) := ZZ => opts -> (F, I) -> last numericalDimensions(matrix{F}, I, opts) numericalImageDim (RingMap, Ideal, Point) := ZZ => opts -> (F, I, p) -> last numericalDimensions(F.matrix, I, p, opts) numericalImageDim (RingMap, Ideal) := ZZ => opts -> (F, I) -> last numericalDimensions(F.matrix, I, opts) -- converts M to a list of 1-element lists of row matrices (to normalize rows easily) -- listForm satisfies M == matrix listForm M (if numrows M, numcols M > 0), and this conversion is fast listForm Matrix := A -> apply(entries A, r -> {matrix{r}}) rowScale := (L, s) -> matrix flatten apply(L, r -> if r#0 == 0 then {} else {(s/norm(2,r#0))*r}) -- deletes any zero rows -- doubleScale := L -> transpose rowScale((entries transpose rowScale(L,1))/(r -> {matrix{r}}), sqrt(#L/(numcols(L#0#0)))) numericalNullity = method(Options => {symbol SVDGap => 1e5, Verbose => false, symbol Precondition => false}) numericalNullity (List, Boolean) := List => opts -> (M, keepSVD) -> ( if matrix M == 0 then return if keepSVD then {numcols M#0#0, 0} else numcols M#0#0; if opts.Verbose then print "Performing normalization preconditioning ..."; T := timing A := if opts.Precondition then rowScale(M, 1) else matrix M; if opts.Verbose then print(" -- used " | toString(T#0) | " seconds"); if opts.Verbose then print "Computing numerical kernel ..."; T = timing (S, U, Vt) := SVD A; -- do not use DivideConquer => true! if opts.Verbose then print(" -- used " | toString(T#0) | " seconds"); largestGap := (#S, opts.SVDGap); for i from 1 to #S-1 do ( if S#i == 0 then ( largestGap = (i, "infinity"); break; ) else if S#(i-1)/S#i > last largestGap then ( largestGap = (i, S#(i-1)/S#i); break; ); ); if keepSVD then {numcols A - first largestGap, (S, U, Vt)} else numcols A - first largestGap ) numericalNullity (Matrix, Boolean) := ZZ => opts -> (M, keepSVD) -> if numrows M == 0 then numcols M else numericalNullity(listForm M, keepSVD, opts) numericalNullity Matrix := ZZ => opts -> M -> numericalNullity(M, false, opts) debug needsPackage "SLPexpressions" monomialGate = method() monomialGate (RingElement, List, List) := ProductGate => (m, varList, expList) -> ( productGate flatten apply(#gens ring m, i -> apply(expList#i, j -> varList#i)) ) monomialGate (RingElement, List) := ProductGate => (m, varList) -> monomialGate(m, varList, first exponents m) makeInterpolationMatrix = method() makeInterpolationMatrix (Matrix, List) := List => (mons, pts) -> ( X := apply(#gens ring mons, i -> inputGate ("x"|i)); Y := matrix{apply(flatten entries mons, m -> monomialGate(m, X))}; -- E := makeEvaluator(Y, matrix{X}); E := makeSLProgram(matrix{X}, Y); out := mutableMatrix(ring pts#0, numrows Y, numcols Y); apply(pts/mutableMatrix, p -> ( evaluate(E, p, out); {matrix out} )) ) numericalHilbertFunction = method(Options => { symbol ConvertToCone => false, symbol Precondition => true, Software => M2engine, symbol SVDGap => 1e5, symbol UseSLP => false, Verbose => true}) numericalHilbertFunction (Matrix, Ideal, List, ZZ) := NumericalInterpolationTable => opts -> (F, I, sampleImagePoints, d) -> ( --outputs a degree d interpolation table for F(V(I)) (F, I, sampleImagePoints) = checkRings(F, I, sampleImagePoints, ConvertToCone => opts.ConvertToCone); y := getSymbol "y"; allMonomials := basis(d, (coefficientRing ring I)(monoid[y_0..y_(numcols F-1)])); N := numcols allMonomials; if #sampleImagePoints < N then ( if opts.Verbose then print "Sampling image points ..."; T := timing sampleImagePoints = sampleImagePoints | numericalImageSample(F, I, sampleImagePoints, N, Software => opts.Software); if opts.Verbose then print(" -- used " | toString(T#0) | " seconds"); ); sampleImagePoints = apply(sampleImagePoints/matrix, p -> 1/norm(2,p)*p); if opts.Verbose then print "Creating interpolation matrix ..."; T = timing A := if opts.UseSLP then makeInterpolationMatrix(allMonomials, sampleImagePoints) else apply(sampleImagePoints, p -> {sub(allMonomials, p)}); if opts.Verbose then print(" -- used " | toString(T#0) | " seconds"); interpolationData := numericalNullity(A, true, Precondition => opts.Precondition, SVDGap => opts.SVDGap, Verbose => opts.Verbose); new NumericalInterpolationTable from { symbol hilbertFunctionArgument => d, symbol hilbertFunctionValue => first interpolationData, symbol imagePoints => VerticalList sampleImagePoints, symbol interpolationBasis => allMonomials, symbol interpolationSVD => last interpolationData, symbol interpolationMatrix => matrix A, symbol map => F } ) numericalHilbertFunction (Matrix, Ideal, ZZ) := NumericalInterpolationTable => opts -> (F, I, d) -> numericalHilbertFunction(F, I, {}, d, opts) numericalHilbertFunction (List, Ideal, List, ZZ) := NumericalInterpolationTable => opts -> (F, I, sampleImagePoints, d) -> numericalHilbertFunction(matrix{F}, I, sampleImagePoints, d, opts) numericalHilbertFunction (List, Ideal, ZZ) := NumericalInterpolationTable => opts -> (F, I, d) -> numericalHilbertFunction(matrix{F}, I, {}, d, opts) numericalHilbertFunction (RingMap, Ideal, List, ZZ) := NumericalInterpolationTable => opts -> (F, I, sampleImagePoints, d) -> numericalHilbertFunction(F.matrix, I, sampleImagePoints, d, opts) numericalHilbertFunction (RingMap, Ideal, ZZ) := NumericalInterpolationTable => opts -> (F, I, d) -> numericalHilbertFunction(F.matrix, I, {}, d, opts) realPartMatrix := A -> matrix apply(entries A, r -> r/realPart) imPartMatrix := A -> if class ring A === RealField then 0 else matrix apply(entries A, r -> r/imaginaryPart) extractImageEquations = method(Options => {symbol Threshold => 5, symbol AttemptZZ => false}) extractImageEquations NumericalInterpolationTable := Matrix => opts -> T -> ( n := opts.Threshold; (V, mons) := (last T.interpolationSVD, T.interpolationBasis); A := clean(10.0^(-n), conjugate transpose V^{numrows V-T.hilbertFunctionValue..numrows V-1}); if not opts.AttemptZZ === false then ( if opts.AttemptZZ === 2 then ( B := if class ring A === ComplexField then matrix table(numrows A, numcols A, (i,j) -> matrix{{realPart A_(i,j),imaginaryPart A_(i,j)}}) else A; C := matrix apply(entries B, r -> r/(e -> lift(round(10^(1+n)*round(n, e)), ZZ))); D := submatrix(LLL(C), numcols A..numcols C-1); E := mons*colReduce(sub(D, ring mons), 10.0^(-n)); ) else ( A = T.interpolationMatrix; B = random(RR)*realPartMatrix A + random(RR)*imPartMatrix A; C = matrix apply(entries B, r -> r/(e -> lift(round(10^(1+n)*round(n, e)), ZZ))); D = submatrix(LLL(id_(ZZ^(numcols C)) || C), toList (0..<numcols mons), toList(0..<T.hilbertFunctionValue)); E = mons*sub(D, ring mons); ); val := sub(E, T.imagePoints#0); if clean(10.0^(-n), val) != 0 then ( << "Warning: some of the integer equations may be inexact. Their values at a sample image point are " << val << endl; ); E ) else mons*sub(A, ring mons) ) extractImageEquations (Matrix, Ideal, ZZ) := Matrix => opts -> (F, I, d) -> extractImageEquations(numericalHilbertFunction(F, I, d), opts) extractImageEquations (List, Ideal, ZZ) := Matrix => opts -> (F, I, d) -> extractImageEquations(numericalHilbertFunction(matrix{F}, I, d), opts) extractImageEquations (RingMap, Ideal, ZZ) := Matrix => opts -> (F, I, d) -> extractImageEquations(numericalHilbertFunction(F.matrix, I, d), opts) round (ZZ, ZZ) := ZZ => (n, x) -> x round (ZZ, CC) := CC => (n, x) -> round(n, realPart x) + ii*round(n, imaginaryPart x) round (ZZ, BasicList) := BasicList => (n, L) -> L/round_n round (ZZ, Matrix) := Matrix => (n, M) -> matrix(entries M/round_n) round (ZZ, RingElement) := RingElement => (n, f) -> ( C := coefficients f; ((C#0)*round(n, lift(C#1, coefficientRing ring f)))_(0,0) ) pseudoWitnessSet = method(Options => { symbol DoRefinements => false, symbol DoTraceTest => true, symbol MaxAttempts => 5, symbol MaxPoints => infinity, symbol MaxThreads => 1, Software => M2engine, symbol Repeats => 3, symbol TraceThreshold => 1e-5, symbol Threshold => 5, -- symbol Endgame => false, Verbose => true}) pseudoWitnessSet (Matrix, Ideal, List, Thing) := PseudoWitnessSet => opts -> (F, I, pointPairs, sliceMatrix) -> ( --outputs a pseudo-witness set for F(V(I)) local imagePointString, local pairTable, local startSystem; y := getSymbol "y"; k := coefficientRing ring I; targetRing := k(monoid[y_1..y_(numcols F)]); if #pointPairs == 0 then error "Expected source point"; sourcePoint := pointPairs#0#0; dims := numericalDimensions(F, I, sourcePoint); numAttempts := 0; traceResult := opts.TraceThreshold + 1; (fiberSlice, fiberdim) := ({}, first dims - last dims); while not traceResult < opts.TraceThreshold and numAttempts < opts.MaxAttempts do ( if numAttempts > 0 then sourcePoint = first numericalSourceSample(I, sourcePoint, Software => opts.Software); pullbackSlice := if sliceMatrix === null then randomSlice(F, last dims, {sourcePoint, "source"}) else ( if numAttempts == 0 and not all(pointPairs, pair -> clean((10.0)^(-opts.Threshold), sub(sliceMatrix, matrix pair#0)) == 0) then error "Expected input points to lie on input slice"; flatten entries sliceMatrix ); squaredUpSource := if I == 0 then {} else randomSlice(gens I, #gens ring I - first dims, {}); if fiberdim > 0 then ( fiberSlice = randomSlice(vars ring I, fiberdim, {sourcePoint, "source"}); if numAttempts == 0 and #pointPairs > 1 then pointPairs = numericalEval(F, {sourcePoint} | flatten apply(toList(1..#pointPairs-1), i -> ( codimSlice := randomSlice(F - sub(matrix pointPairs#i#1, ring F), first dims - fiberdim, {}); localFiberSlice := codimSlice | squaredUpSource | randomSlice(vars ring I, fiberdim, {pointPairs#i#0, "source"}); globalFiberSlice := codimSlice | squaredUpSource | fiberSlice; myTrack(localFiberSlice, globalFiberSlice, {pointPairs#i#0}) )), true); ); newStartSystem := squaredUpSource | fiberSlice | pullbackSlice; newPairs := if numAttempts > 0 then numericalEval(F, myTrack(startSystem, newStartSystem, (values pairTable)/first, opts), true) else pointPairs/(pair -> (pair#0, matrix pair#1)); if #newPairs == 0 then ( if opts.Verbose then print "Failed to track old points to new slice. Retrying..."; numAttempts = numAttempts + 1; continue; ); pairTable = new MutableHashTable; for pair in newPairs do ( imagePointString = toString round(opts.Threshold, last pair); if not pairTable#?imagePointString then pairTable#imagePointString = pair; ); startSystem = newStartSystem; pointPairs = monodromyLoop(F, last dims, startSystem, pairTable, opts); if not opts.DoTraceTest then break; if opts.DoRefinements then ( if opts.Verbose then print "Refining solutions..."; pointPairs = numericalEval(F, refine(startSystem, pointPairs/first, Bits => precision ring I), true); ); if opts.Verbose then print("Running trace test ..."); traceResult = traceTest(F, last dims, pointPairs, startSystem, opts); if not traceResult < opts.TraceThreshold and opts.Verbose then print("Failed trace test! Trace: " | toString traceResult); numAttempts = numAttempts + 1; ); if opts.Verbose then ( if traceResult > opts.TraceThreshold then ( print("Degree of image should be at least " | #pointPairs); print("Consider changing parameters (Repeats, MaxAttempts, Threshold) or reparametrizing for a better result."); -- Alternatively, consider increasing precision (e.g. changing ground field to CC_100). ); ); new PseudoWitnessSet from { symbol isCompletePseudoWitnessSet => traceResult < opts.TraceThreshold, symbol degree => #pointPairs, symbol map => F, symbol sourceEquations => I, symbol generalCombinations => matrix{squaredUpSource}, symbol sourceSlice => matrix{fiberSlice}, symbol imageSlice => matrix{pullbackSlice}, symbol witnessPointPairs => VerticalList apply(pointPairs, pair -> (pair#0, point pair#1)), symbol trace => traceResult } ) pseudoWitnessSet(Matrix, Ideal, Point) := PseudoWitnessSet => opts -> (F, I, p) -> ( (F, I, p) = checkRings(F, I, {p}); pseudoWitnessSet(F, I, numericalEval(F, p, true), null, opts) ) pseudoWitnessSet (Matrix, Ideal) := PseudoWitnessSet => opts -> (F, I) -> ( if opts.Verbose then print "Sampling point in source ..."; pseudoWitnessSet(F, I, first numericalSourceSample I, opts) ) pseudoWitnessSet(List, Ideal, List, Thing) := PseudoWitnessSet => opts -> (F, I, pointPairs, L) -> pseudoWitnessSet(matrix{F}, I, pointPairs, L, opts) pseudoWitnessSet(List, Ideal, Point) := PseudoWitnessSet => opts -> (F, I, p) -> pseudoWitnessSet(matrix{F}, I, p, opts) pseudoWitnessSet (List, Ideal) := PseudoWitnessSet => opts -> (F, I) -> pseudoWitnessSet(matrix{F}, I, opts) pseudoWitnessSet(RingMap, Ideal, List, Thing) := PseudoWitnessSet => opts -> (F, I, pointPairs, L) -> pseudoWitnessSet(F.matrix, I, pointPairs, L, opts) pseudoWitnessSet(RingMap, Ideal, Point) := PseudoWitnessSet => opts -> (F, I, p) -> pseudoWitnessSet(F.matrix, I, p, opts) pseudoWitnessSet (RingMap, Ideal) := PseudoWitnessSet => opts -> (F, I) -> pseudoWitnessSet(F.matrix, I, opts) numericalImageDegree = method(Options => options pseudoWitnessSet) numericalImageDegree PseudoWitnessSet := ZZ => opts -> W -> W.degree numericalImageDegree (Matrix, Ideal) := ZZ => opts -> (F, I) -> (pseudoWitnessSet(F, I, opts)).degree numericalImageDegree (List, Ideal) := ZZ => opts -> (F, I) -> (pseudoWitnessSet(matrix{F}, I, opts)).degree numericalImageDegree (RingMap, Ideal) := ZZ => opts -> (F, I) -> (pseudoWitnessSet(F.matrix, I, opts)).degree myTrack = method(Options => options pseudoWitnessSet) myTrack (List, List, List) := List => opts -> (startSystem, targetSystem, startSolutions) -> ( k := coefficientRing ring startSystem#0; randomGamma := random k; if #startSolutions > max(10, 2*opts.MaxThreads) and opts.MaxThreads > 1 then ( -- prints many errors, but continues to run --setIOExclusive(); -- buggy: causes isReady to indefinitely hang startSolutionsList := pack(ceiling(#startSolutions/opts.MaxThreads), startSolutions); threadList := {}; for paths in startSolutionsList do ( threadList = append(threadList, schedule(x -> timing track x, (startSystem, targetSystem, paths, gamma => randomGamma, Software => opts.Software))); ); while not all(threadList, isReady) do sleep 1; results := delete(null, threadList/taskResult); targetSolutions := flatten(results/last); if opts.Verbose then print("Finished tracking " | #targetSolutions | " paths in parallel, in " | toString sum(results/first) | " seconds"); ) else ( -- if startSolutions is empty then error is thrown! T := timing targetSolutions = track(startSystem, targetSystem, startSolutions, gamma => randomGamma, Software => opts.Software); if opts.Verbose and T#0 > 1 then print (" -- used " | toString(T#0) | " seconds"); ); goodSols := select(targetSolutions, p -> p.cache#?SolutionStatus and status p == Regular); if opts.Verbose and #goodSols < #startSolutions then print("Paths going to infinity: " | #startSolutions - #goodSols | " out of " | #startSolutions); if opts.DoRefinements then goodSols = apply(refine(polySystem targetSystem, goodSols, Bits => precision k), p -> point sub(matrix p, k)); goodSols ) randomSlice = method() -- returns a list of c random linear combinations of polys (row matrix) passing through (optional source or target) point, via translation randomSlice (Matrix, ZZ, List) := List => (polys, c, pointData) -> ( R := ring polys; coeffs := random(R^(numcols polys), R^c); G := polys*coeffs; flatten entries(G - if #pointData == 0 then 0 else sub(if pointData#1 == "source" then sub(G, matrix pointData#0) else (matrix pointData#0)*coeffs, R)) ) monodromyLoop = method(Options => options pseudoWitnessSet) monodromyLoop (Matrix, ZZ, List, MutableHashTable) := List => opts -> (F, imageDim, startSystem, pairTable) -> ( numRepetitiveMonodromyLoops := 0; numPts := {#values pairTable}; if opts.Verbose then print "Tracking monodromy loops ..."; while numRepetitiveMonodromyLoops < opts.Repeats do ( intermediateSystem1 := drop(startSystem, -imageDim) | randomSlice(F | matrix{{10_(ring F)}}, imageDim, {}); startSols := (values pairTable)/first; -- increment := if opts.Endgame and #startSols > 100 and #numPts > 5 and all((firstDifference firstDifference numPts)_{-3..-1}, d -> d < 0) then ( -- startSols = startSols_(randomInts(#startSols, max(100, #startSols//10))); -- 1/4 -- ) else 1; intermediateSolutions1 := myTrack(startSystem, intermediateSystem1, startSols, opts); if #intermediateSolutions1 > 0 then ( endSolutions := myTrack(intermediateSystem1, startSystem, intermediateSolutions1, opts); if #endSolutions > 0 then ( candidatePairs := numericalEval(F, endSolutions, true); for pair in candidatePairs do ( imagePointString := toString round(opts.Threshold, last pair); if not pairTable#?imagePointString then pairTable#imagePointString = pair; ); ); ); if numPts#-1 < #values pairTable then numRepetitiveMonodromyLoops = 0 else numRepetitiveMonodromyLoops = numRepetitiveMonodromyLoops + 1; numPts = append(numPts, #values pairTable); if opts.Verbose then print ("Points found: " | numPts#-1); if numPts#-1 >= opts.MaxPoints then break; ); values pairTable ) traceTest = method(Options => options pseudoWitnessSet) traceTest (Matrix, ZZ, List, List) := RR => opts -> (F, imageDim, intersectionPointPairs, startSystem) -> ( C := coefficientRing ring F; startUpstairsPoints := intersectionPointPairs/first; startDownstairsPoints := intersectionPointPairs/last; for translationMagnitude in {0,1,3,2,-1,5,-2,6} do ( randomTranslation := 10^(translationMagnitude)*flatten entries(map(C^1, C^(#startSystem - imageDim), 0) | random(C^1, C^imageDim)); gammas := {random C, random C}; firstStepSystem := startSystem + (first gammas)*randomTranslation; secondStepSystem := startSystem + (last gammas)*randomTranslation; firstStepUpstairsPoints := myTrack(startSystem, firstStepSystem, startUpstairsPoints, opts); if #firstStepUpstairsPoints == #startUpstairsPoints then ( secondStepUpstairsPoints := myTrack(startSystem, secondStepSystem, startUpstairsPoints, opts); if #secondStepUpstairsPoints == #startUpstairsPoints then ( firstStepDownstairsPoints := numericalEval(F, firstStepUpstairsPoints, false); secondStepDownstairsPoints := numericalEval(F, secondStepUpstairsPoints, false); traceList := (1/first gammas)*(firstStepDownstairsPoints - startDownstairsPoints) - (1/last gammas)*(secondStepDownstairsPoints - startDownstairsPoints); return norm(2,sum traceList); ); ); ); infinity ) isOnImage = method(Options => { MaxThreads => 1, Software => M2engine, Threshold => 5, Verbose => true}) isOnImage (PseudoWitnessSet, Point) := Boolean => opts -> (W, q) -> ( q = matrix q; if not W.isCompletePseudoWitnessSet then print "Warning: not a complete pseudo-witness set! May return false negative."; F := W.map; I := W.sourceEquations; if not ring q === coefficientRing ring I then error "Point must have coordinates in the coefficient ring of the ideal."; fiberSlice := flatten entries W.sourceSlice; pullbackSlice := flatten entries W.imageSlice; squaredUpSource := flatten entries W.generalCombinations; startUpstairsPoints := W.witnessPointPairs /first; newPullbackSlice := randomSlice(F, #pullbackSlice, {q, "target"}); targetUpstairsPoints := myTrack(squaredUpSource | fiberSlice | pullbackSlice, squaredUpSource | fiberSlice | newPullbackSlice, startUpstairsPoints, opts); imagePointTable := hashTable apply(numericalEval(F, targetUpstairsPoints, false), p -> round(opts.Threshold, p) => 0); imagePointTable#?(round(opts.Threshold, q)) ) isOnImage (Matrix, Ideal, Point) := Boolean => opts -> (F, I, q) -> isOnImage(pseudoWitnessSet(F, I, opts), q, opts) isOnImage (List, Ideal, Point) := Boolean => opts -> (F, I, q) -> isOnImage(matrix{F}, I, q, opts) isOnImage (RingMap, Ideal, Point) := Boolean => opts -> (F, I, q) -> isOnImage(F.matrix, I, q, opts) isWellDefined NumericalInterpolationTable := Boolean => T -> ( -- CHECK DATA STRUCTURE -- CHECK KEYS K := keys T; expectedKeys := set { symbol hilbertFunctionArgument, symbol hilbertFunctionValue, symbol imagePoints, symbol interpolationBasis, symbol interpolationSVD, symbol interpolationMatrix, symbol map }; if set K =!= expectedKeys then ( if debugLevel > 0 then ( added := toList(K - expectedKeys); missing := toList(expectedKeys - K); if #added > 0 then << "-- unexpected key(s): " << toString added << endl; if #missing > 0 then << "-- missing keys(s): " << toString missing << endl; ); return false ); -- CHECK TYPES if not instance(T.hilbertFunctionArgument, ZZ) then ( if debugLevel > 0 then << "-- expected `hilbertFunctionArgument' to be an integer" << endl; return false ); if not instance(T.hilbertFunctionValue, ZZ) then ( if debugLevel > 0 then << "-- expected `hilbertFunctionValue' to be an integer" << endl; return false ); if not instance(T.map, Matrix) then ( if debugLevel > 0 then << "-- expected `map' to be a matrix" << endl; return false ); if not instance(T.interpolationBasis, Matrix) then ( if debugLevel > 0 then << "-- expected `interpolationBasis' to be a matrix" << endl; return false ); if not instance(T.interpolationMatrix, Matrix) then ( if debugLevel > 0 then << "-- expected `interpolationMatrix' to be a matrix" << endl; return false ); if not instance(T.interpolationSVD, Sequence) then ( if debugLevel > 0 then << "-- expected `interpolationSVD' to be a sequence" << endl; return false ); if not instance(first T.interpolationSVD, List) then ( if debugLevel > 0 then << "-- expected first element of `interpolationSVD' to be a list" << endl; return false ); if not all(first T.interpolationSVD, s -> instance(s, RR)) then ( if debugLevel > 0 then << "-- expected first element of `interpolationSVD' to be a list of singular values" << endl; return false ); if not all(drop(T.interpolationSVD, 1), M -> instance(M, Matrix)) then ( if debugLevel > 0 then << "-- expected second and third elements of `interpolationSVD' to be matrices" << endl; return false ); -- CHECK MATHEMATICAL STRUCTURE if not unique flatten last degrees T.interpolationBasis === {T.hilbertFunctionArgument} then ( if debugLevel > 0 then << ("-- expected `interpolationBasis' to consist of monomials of degree " | T.hilbertFunctionArgument) << endl; return false ); if not all({coefficientRing ring T.interpolationBasis, ring(T.interpolationSVD#2)}/class, C -> C === ComplexField) then ( if debugLevel > 0 then << "-- expected ground field to be complex numbers" << endl; return false ); numMonomials := binomial(numcols T.map + T.hilbertFunctionArgument - 1, T.hilbertFunctionArgument); if not #gens ring T.interpolationBasis === numcols T.map or not numcols T.interpolationBasis === numMonomials then ( if debugLevel > 0 then << ("-- expected `interpolationBasis' to have " | numMonomials | " monomials in " | #T.map | " variables") << endl; return false ); true ) isWellDefined PseudoWitnessSet := Boolean => W -> ( -- CHECK DATA STRUCTURE -- CHECK KEYS K := keys W; expectedKeys := set { symbol isCompletePseudoWitnessSet, symbol degree, symbol map, symbol sourceEquations, symbol sourceSlice, symbol generalCombinations, symbol imageSlice, symbol witnessPointPairs, symbol trace }; if set K =!= expectedKeys then ( if debugLevel > 0 then ( added := toList(K - expectedKeys); missing := toList(expectedKeys - K); if #added > 0 then << "-- unexpected key(s): " << toString added << endl; if #missing > 0 then << "-- missing keys(s): " << toString missing << endl; ); return false ); -- CHECK TYPES if not instance(W.isCompletePseudoWitnessSet, Boolean) then ( if debugLevel > 0 then << "-- expected `isCompletePseudoWitnessSet' to be a Boolean" << endl; return false ); if not instance(W.degree, ZZ) then ( if debugLevel > 0 then << "-- expected `degree' to be an integer" << endl; return false ); if not instance(W.map, Matrix) then ( if debugLevel > 0 then << "-- expected `map' to be a matrix" << endl; return false ); if not instance(W.sourceEquations, Ideal) then ( if debugLevel > 0 then << "-- expected `sourceEquations' to be an ideal" << endl; return false ); if not instance(W.sourceSlice, Matrix) then ( if debugLevel > 0 then << "-- expected `sourceSlice' to be a matrix" << endl; return false ); if not instance(W.generalCombinations, Matrix) then ( if debugLevel > 0 then << "-- expected `generalCombinations' to be a matrix" << endl; return false ); if not instance(W.imageSlice, Matrix) then ( if debugLevel > 0 then << "-- expected `imageSlice' to be a matrix" << endl; return false ); if not instance(W.witnessPointPairs, List) then ( if debugLevel > 0 then << "-- expected `witnessPointPairs' to be a list" << endl; return false ); if not all(W.witnessPointPairs, pair -> instance(pair, Sequence)) then ( if debugLevel > 0 then << "-- expected `witnessPointPairs' to be a list of sequences" << endl; return false ); if not all(W.witnessPointPairs, pair -> all(pair, p -> instance(p, Point))) then ( if debugLevel > 0 then << "-- expected `witnessPointPairs' to be a list of sequences of points" << endl; return false ); if not instance(W.trace, RR) then ( if debugLevel > 0 then << "-- expected `trace' to be a real number" << endl; return false ); -- CHECK MATHEMATICAL STRUCTURE R := ring W.sourceEquations; if not R === ring W.map then ( if debugLevel > 0 then << "-- expected `map' and `sourceEquations' to have the same ring" << endl; return false ); if not instance(class 1_(coefficientRing R), InexactFieldFamily) then ( if debugLevel > 0 then << "-- expected ground field to have floating point arithmetic" << endl; return false ); if not all(W.witnessPointPairs, pair -> #(pair#0#Coordinates) === #gens R and #(pair#1#Coordinates) === numcols W.map) then ( if debugLevel > 0 then << "-- number of coordinates in `witnessPointPairs' do not match" << endl; return false ); if not all(W.witnessPointPairs/first, p -> clean(1e-5, sub(gens W.sourceEquations, matrix p)) == 0) then ( if debugLevel > 0 then << " -- expected first components of `witnessPointPairs' to satisfy `sourceEquations'" << endl; return false ); if not all(W.witnessPointPairs, pair -> clean(1e-5, matrix last pair - sub(W.map, matrix first pair)) == 0) then ( if debugLevel > 0 then << " -- expected components `witnessPointPairs' to correspond under `map'" << endl; return false ); if not all(W.witnessPointPairs/first/matrix, p -> clean(1e-5, sub(W.imageSlice, p)) == 0) then ( if debugLevel > 0 then << " -- expected second components of `witnessPointPairs' to lie on `imageSlice'" << endl; return false ); true ) beginDocumentation() load "NumericalImplicitization/doc.m2" load "NumericalImplicitization/tests.m2" end-- restart needsPackage "NumericalImplicitization" loadPackage("NumericalImplicitization", Reload => true) uninstallPackage "NumericalImplicitization" installPackage "NumericalImplicitization" installPackage("NumericalImplicitization", RemakeAllDocumentation => true) viewHelp "NumericalImplicitization" check "NumericalImplicitization" -- Defaults: 10 in monodromyLoop for affine term, 20 in myTrack (for parallelization), {0,1,3,2,-1,5,-2,6} for translationMagnitude in traceTest -- Future: Alpha-certify option, improvements to interpolation (using non-monomial bases, Jacobi SVD) -- Guarantee correct number of sample points, starting from initial point p goodPts = {} sampleSize = 100 eps = 1e-10 elapsedTime while #goodPts < sampleSize do ( q = first numericalSourceSample(I, p); if q#?SolutionStatus and q#SolutionStatus === Regular and clean(eps, sub(gens I, matrix q)) == 0 then goodPts = goodPts | {q}; ) -- high degree rational normal curve R = CC[s,t],; F = basis(40,R); I = ideal 0_R; numericalImageDim(F, I) time tests = toList(1..100)/(i -> pseudoWitnessSet(F,I,Repeats=>2,Verbose=>false)); -- Generic Pinched Veronese R = CC[x_0..x_3] F = toList(1..5)/(i -> random(10,R)); allowableThreads = maxAllowableThreads pseudoWitnessSet(F,ideal 0_R,Repeats=>2) -- Trifocal variety R=CC[a00,a01,a02,a03,a10,a11,a12,a13,a20,a21,a22,a23,b10,b11,b12,b13,b20,b21,b22,b23],;A = transpose genericMatrix(R,a00,4,3),;B = matrix{{0,0,0,1},{b10,b11,b12,b13},{b20,b21,b22,b23}},;C = matrix{{1_R,0,0,0},{0,1,0,0},{0,0,1,0}},;M = A||B||C,;F = flatten flatten apply(3, i-> apply(3, j-> apply(reverse subsets(3,2), k->det submatrix(M,{i}|{j+3}|(k+{6,6}) , ) ) )); allowableThreads = 4 elapsedTime pseudoWitnessSet(F,ideal 0_R,Repeats=>2,MaxThreads=>allowableThreads) -- Tensor product surface (a,b) = (3,1) R=CC[s,t,u,v, Degrees=>{{1,0},{1,0},{0,1},{0,1}}] Ix=intersect(ideal(s,u),ideal(t,v)) B=super basis({a,b},Ix) C=matrix{{1_R,1,0,0,0,0},{0,1,1,0,0,0},{0,0,1,1,0,0},{0,0,0,1,1,1}} F = C*transpose(B) I = ideal 0_R numericalImageDim(F,I) W = pseudoWitnessSet(F,I) T = numericalHilbertFunction(F,I,W.degree) extractImageEquations(T, AttemptZZ => true) -- Undirected graphical model on 4 variables setRandomSeed 0 loadPackage "GraphicalModels" G = graph({1,2,3,4},{{1,2},{1,3},{2,3},{3,4}}) R = CC[x_(1,1),x_(1,2),x_(1,4),x_(2,2),x_(2,4),x_(3,3),x_(3,4),x_(4,4)] M = matrix{{x_(1,1),x_(1,2),0,x_(1,4)},{x_(1,2),x_(2,2),0,x_(2,4)},{0,0,x_(3,3),x_(3,4)},{x_(1,4),x_(2,4),x_(3,4),x_(4,4)}} F = flatten(for i from 1 to 4 list ( for j from i to 4 list ( det(submatrix'(M, {i-1}, {j-1})) ) )) I = ideal 0_R numericalImageDim(F, I) pseudoWitnessSet(F, I, Repeats => 2) T = numericalHilbertFunction(F, I, 2) extractImageEquations(T, AttemptZZ => true) -- Check approximate equations: T = numericalHilbertFunction(F, ideal 0_R, 2); E = extractImageEquations(T, AttemptZZ => true); all((toList T.imagePoints)/(p -> clean(1e-11, sub(E, toList(0..<#(p#Coordinates))/(i -> (gens ring E)#i => (p#Coordinates)#i)))), v -> v == 0) --------------- Implicitization Challenge + variants -- (line) secant of (P^1)^n in P^31, n = 5: degree 3256 n = 5 R = CC[a_1..a_n,b_1..b_n,s,t]; F = s*(terms product apply(toList(1..n), i-> 1 + a_i)) + t*(terms product apply(toList(1..n), i-> 1 + b_i)); allowableThreads = maxAllowableThreads time W = pseudoWitnessSet(F, ideal 0_R, Repeats => 1) elapsedTime W = pseudoWitnessSet(F, ideal 0_R, Repeats => 1, MaxThreads => allowableThreads) -- Challenge: Hadamard square of line secant of (P^1)^4 in P^15, degree 110, passed in 188.084 seconds t = symbol t; n = 4 R = CC[a_1..a_n,b_1..b_n, c_1..c_n, d_1..d_n,t_0..t_3]; F1 = t_0*(terms product apply(toList(1..n), i->(1 + a_i))) + t_1*(terms product apply(toList(1..n), i->(1 + b_i))); F2 = t_2*(terms product apply(toList(1..n), i->(1 + c_i))) + t_3*(terms product apply(toList(1..n), i->(1 + d_i))); F = apply(toList(0..15), i -> F1#i * F2#i); allowableThreads = maxAllowableThreads time W = pseudoWitnessSet(F, ideal 0_R, Repeats => 1, MaxThreads => allowableThreads) -- precision tests R = CC_54[s,t]; I = ideal 0_R; W = pseudoWitnessSet(basis(3, R), I) toList W.witnessPointPairs /first/(p -> p#Coordinates )/first/ring prec = 500 setDefault(Precision => prec) R = CC_prec[s,t]; I = ideal 0_R; F = basis(3, R); W = pseudoWitnessSet(F, I) R = CC[s,t]; F = basis(4, R); I = ideal 0_R T = numericalHilbertFunction(F, I, 2) A = matrix T.interpolationMatrix prec = 5 printingPrecision = 16 setDefault(Precision => prec) R = CC_prec[x_0..x_3] R = CC[x_0..x_2] I = ideal random(R^1,R^{-2,-3}) I = ideal(random(2,R), random(3,R)) F = random(R^1,R^{3:-3}) F = matrix{toList(1..3)/(i -> random(3,R))} d = 18 -- ConvertToCone tests R = CC[x_1..x_5] F = vars R I = ideal(x_1*x_2^2 - x_3^2, x_1*x_4^2 - x_5^2) elapsedTime I.cache.WitnessSet = first components numericalIrreducibleDecomposition I time pseudoWitnessSet(F,I) F = (gens R)_(toList(0,1,2)) I = ideal(x_1-x_3*x_4,x_2-x_3*x_5,x_4*x_5-1) -- SO(5) n = 5 R = CC[x_0..x_(n^2-1)] A = genericMatrix(R,n,n) I = ideal(A*transpose A - id_(R^n)); F = vars R p = point id_((coefficientRing R)^n) q = first numericalSourceSample(I, p) allowableThreads = 4 elapsedTime PW = pseudoWitnessSet(F,I,p,Threshold=>3,MaxThreads=>allowableThreads, MaxPoints=>384) -- passed in 153.496 seconds -- Dim challenge example jacI=(d,l,n)->(S=CC[x_(0,1)..x_(n,l),c_0..c_(binomial(l,n)-1)];R= S[X_0..X_n]; M=for i from 1 to l list matrix{toList (x_(0,i)..x_(n,i))}; H=for b in(for i from 0 to#subsets(M,n)-1 list for a in(subsets(M,n))_i list{a}) list matrix b; P=for t from 0 to #H-1 list for j from 0 to n list(-1)^(j)*(minors(n,H_t))_(n-j); F=sum for i from 0 to #P-1 list c_(i)*(sum for j from 0 to n list P_i_j*X_j)^d; I=transpose substitute((coefficients F)#1,S)) t=13 time F = jacI(t,t+1,2); time numericalImageDim(F, ideal 0_S)
Simpan