diff --git a/python/drafts/nid/Markdown.md b/python/drafts/nid/Markdown.md new file mode 100644 index 000000000..f967bb2e0 --- /dev/null +++ b/python/drafts/nid/Markdown.md @@ -0,0 +1,5 @@ +### Drafts of python versions of code that will make it's way into core + +`SliceMover.ipynb` contains a class that can create homotopies and calculate the slice move on a system using the bertini.nag_algorithm.moving_homotopy(). It contains example driver code as well as a sample graph, other then the driver code, and it only requires access to the bertini library. + +`RegenCascade.ipynb` contains the beginning of the Regeneration Cascade to generate witness points for the system. It requires junk removal to be added still, but it'll start generating the points needed for the next step. Currently it generates start points from codimension 1, and then I'm working on codim 2 and beyond. diff --git a/python/drafts/nid/RegenCascade.ipynb b/python/drafts/nid/RegenCascade.ipynb new file mode 100644 index 000000000..30109684c --- /dev/null +++ b/python/drafts/nid/RegenCascade.ipynb @@ -0,0 +1,593 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "934ab4b8", + "metadata": {}, + "source": [ + "## 2. `RegenCascade.ipynb`\n", + " \n", + "### Purpose\n", + "Drives the regeneration algorithm across an entire system `f1, f2, …, fn`\n", + "one function (\"codimension\") at a time. At each stage it needs a set of\n", + "**start points** for the homotopy that adds the next equation; the\n", + "`BasePointRegen` class builds those start points for the very first stage,\n", + "and `SliceMover` (from the other notebook) is reused to advance the slice\n", + "between stages." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "%run SliceMover.ipynb\n", + "import bertini as b2\n", + "import numpy as np\n", + "from bertini import nag_algorithm\n", + "from bertini import Slice" + ] + }, + { + "cell_type": "markdown", + "id": "031dc474", + "metadata": {}, + "source": [ + "### Class `BasePointRegen`\n", + " \n", + "Stores the system, its first variable group (`self.vars = system.variable_groups()[0]` — assuming a single variable group), and\n", + "`self.dim = system.num_functions()`. Several result attributes\n", + "(`base_points_numeric`, `start_moving_system`, `end_moving_system`,\n", + "`degree`, `start_points`) start as `None` and are filled in by\n", + "`generate_start_points`.\n", + " \n", + "#### `generate_start_points(bottom_slice, *, rng=None)`\n", + " \n", + "Builds the base-level (codimension-1) start points for regeneration by\n", + "using the product of linear forms intersected with a\n", + "generic linear slice, which can be solved directly, without a homotopy.\n", + " \n", + "Step by step:\n", + " \n", + "1. **Degree lookup.** `degree = list(self.system.degrees())[0]` — the\n", + " degree of the system's first function, `f1`. \n", + "2. **Randomize to codimension 1.**\n", + " `self.randomized_sys = self.system.randomize(codimension=1)` collapses\n", + " the system to a single representative equation.\n", + "3. **Build two parallel container systems**, `start_moving` and\n", + " `end_moving`, both starting from an empty system carrying only the\n", + " variable group:\n", + " - `start_moving` gets the product of `degree`\n", + " random linear forms (`add_products_of_linears`), each factor's\n", + " coefficients drawn via `b2.random_complex()` and converted to\n", + " Bertini's symbolic complex type.\n", + " - `end_moving` gets the real (randomized) target function\n", + " `target_fn`.\n", + " - `bottom_slice.add_to(...)` appends the same generic linear slice\n", + " equations to both systems, cutting each down to isolated points.\n", + "4. **Solve the start system in closed form.** Since `start_moving`'s\n", + " defining equation is a product of `degree` linear factors, its\n", + " intersection with the slice is just `degree` separate linear systems. \n", + " For each factor's coefficient row, the code stacks\n", + " that row with the slice's coefficient rows into a square matrix `A`\n", + " and solves `A·x = b` with `np.linalg.solve` — giving the `degree` base\n", + " points directly, with no path tracking needed.\n", + "5. **Convert to Bertini points** via `complex_mp`, and store:\n", + " `self.start_points` (Bertini form), `self.base_points_numeric` (raw\n", + " NumPy), `self.start_moving_system`, `self.end_moving_system`,\n", + " `self.degree`." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "class BasePointRegen:\n", + " def __init__(self, system):\n", + " self.system = system\n", + " self.randomized_sys = None\n", + " self.vars = system.variable_groups()[0]\n", + " self.dim = system.num_functions()\n", + "\n", + " self.start_points = None # base-level start points, set by generate_start_points\n", + "\n", + " self.base_points_numeric = None\n", + " self.start_moving_system = None\n", + " self.end_moving_system = None\n", + " self.degree = None\n", + "\n", + " def generate_start_points(self, bottom_slice, *, rng=None):\n", + " #rng setup\n", + " if rng is None:\n", + " rng = np.random.default_rng()\n", + "\n", + " #var setup\n", + " degree = list(self.system.degrees())[0]\n", + " \n", + " #randomizes R\n", + " self.randomized_sys = self.system.randomize(codimension = 1)\n", + "\n", + " target_fn = self.randomized_sys.function(0)\n", + " bottom_slice_coeffs = bottom_slice.coefficients()\n", + "\n", + " #systems setup\n", + " fixed_sys = b2.System()\n", + " fixed_sys.add_variable_group(self.vars)\n", + "\n", + " start_moving = b2.system.clone(fixed_sys)\n", + " end_moving = b2.system.clone(fixed_sys)\n", + "\n", + " #random coeff for the degree\n", + " factor_coeffs = [[b2.random_complex() for ii in range(self.dim + 1)] for jj in range(degree)]\n", + "\n", + " #converts into symbolics complex numbers\n", + " def convert_to_b2_symbolic_complex(z):\n", + " return b2.symbolics.Complex(str(z.real), str(z.imag))\n", + "\n", + " #add functions to the systems\n", + " start_moving.add_products_of_linears([[[convert_to_b2_symbolic_complex(c) for c in row] for row in factor_coeffs]])\n", + " end_moving.add_function(target_fn)\n", + "\n", + " bottom_slice.add_to(start_moving)\n", + " bottom_slice.add_to(end_moving)\n", + "\n", + " #matrix setup\n", + " base_points_numeric = []\n", + " for row in factor_coeffs:\n", + " row = np.array(row)\n", + " A = np.vstack([row[:-1][None, :], bottom_slice_coeffs[:, :-1]]) # N x N\n", + " b = -np.concatenate([[row[-1]], bottom_slice_coeffs[:, -1]]) # N\n", + "\n", + " base_points_numeric.append(np.linalg.solve(A.astype(complex), b.astype(complex)))\n", + "\n", + " #converts points into a numpy array of complex_mp points\n", + " def to_bertini_point(pt):\n", + " return np.array([b2.complex_mp(str(v.real), str(v.imag)) for v in pt])\n", + "\n", + " base_points_bertini_form = [to_bertini_point(p) for p in base_points_numeric]\n", + "\n", + " #results\n", + " self.start_points = base_points_bertini_form\n", + " self.base_points_numeric = base_points_numeric\n", + " self.start_moving_system = start_moving\n", + " self.end_moving_system = end_moving\n", + " self.degree = degree" + ] + }, + { + "cell_type": "markdown", + "id": "d2f59c45", + "metadata": {}, + "source": [ + "### Driver cell\n", + " \n", + "Builds a concrete 3-variable example and runs the codimension-by-codimension loop:\n", + " \n", + "```python\n", + "x, y, z = b2.Variable('x'), b2.Variable('y'), b2.Variable('z')\n", + "f1 = (y - x**2) * (x**2 + y**2 + z**2 - 1) * (x - 2)\n", + "f2 = (z - x**3) * (x**2 + y**2 + z**2 - 1) * (y - 2)\n", + "f3 = (z - x**3) * (y - x**2) * (x**2 + y**2 + z**2 - 1) * (z - 2)\n", + "```\n", + "`sys` gets the variable group `[x, y, z]` and all three functions;\n", + "`bpr = BasePointRegen(sys)`, and a maximal generic slice\n", + "`bottom_slice = b2.Slice.random_complex(sys.variable_groups(), sys.num_variables() - 1)`\n", + "is built once, up front (2 linear equations, since there are 3 variables).\n", + " \n", + "The main loop, `for codim in range(1, bpr.dim + 1)` (i.e. codim 1, 2, 3):\n", + " \n", + "- **`codim == 1`:**\n", + " Calls `bpr.generate_start_points(bottom_slice)`, pulls out\n", + " `start_points`, `start_system`, `target_system`, and — for\n", + " diagnostics — evaluates `start_system` at every start point and prints\n", + " the residual. An inline comment flags an important numerical caveat:\n", + " **raw residual magnitude isn't scale-invariant**, so a fixed tolerance\n", + " like `1e-13` can't reliably tell a true root from a non-root once the\n", + " system has been rescaled (e.g. by randomization).\n", + "- **`codim > 1`:**\n", + " Takes the slices from the previous stage's target system\n", + " (`target_system.slices()`), drops the first one\n", + " (`target_slices = slices_prev_step[1:]`), and builds a fresh random\n", + " 2-dimensional slice (`Slice.random_complex(vars, 2)`). A new\n", + " `SliceMover(bpr.randomized_sys, slices_prev_step[0], rand_slice, start_points)` is created, then `sm.move_slice()` is called in a loop that runs `bpr.degree` times, each time collecting `results.solutions` into\n", + " `new_start_points` and advancing to a fresh random end slice via `sm.set_end_slice(...)`. \n", + "- **Every iteration (both branches):**\n", + " Builds the actual regeneration homotopy that adds the next equation:\n", + "```python\n", + " hom = b2.nag_algorithm.blend_homotopy(target=target_system, start=start_system)\n", + " homSolver = b2.HomotopySolver(target=target_system, homotopy=hom, start_points=start_points)\n", + " solve_result = homSolver.solve()\n", + "```\n", + " `solve_result` isn't consumed further in the code shown (its print\n", + " statement, and a following call to a not-yet-implemented\n", + " `bpr.set_up_partially_regen_system(codim)`, are both commented out).\n", + " \n", + "---\n", + " \n", + "#### Known issues and loose ends\n", + " \n", + "1. **\"Junk removal\"** after the repeated `move_slice()` calls in the\n", + " `codim > 1` branch is flagged as not yet implemented.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "codim 1\n", + "start_sys:\n", + " [(-1.0421075930853945000020e+00, 1.3575679903465262999980e+00)\n", + " (1.1196502975482482999980e+00, 1.2146505795945861000000e-01)\n", + " (-4.5905432682353849999980e-01, -5.8915688836732570000050e-01)] [(1.2066938925426353000010e-01, 4.6983945322518605000070e-01)\n", + " (-2.5258262307363083000070e-01, 1.3030530841713234999970e+00)\n", + " (5.3674869763359160000100e-01, -6.8625177450125760000060e-01)] [(-2.2410193566864312000040e-01, 1.5548508551481523999990e+00)\n", + " (9.4485329646279010000100e-02, -6.2783987296640339999920e-02)\n", + " (-5.4613936738801480000000e-02, -1.7973703369834137999970e-01)] [(-1.9574953655444255999970e+00, 3.4344881486158657000020e-01)\n", + " (2.3245977432390080000050e+00, 1.3079889924555659999980e+00)\n", + " (-6.2579092806664770000040e-01, -1.5085079986100937999980e+00)] [(-3.9466412111166010000020e-02, 4.7342747385981593000080e-01)\n", + " (-5.4964425796656279999910e-02, 1.2869650186154873000010e+00)\n", + " (4.4236570472771129999990e-01, -7.4186445690234450000060e-01)]\n", + "codim 2\n", + "rand_sys: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "1 function:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + "\n", + "slices_prev_step: linear slice on 3 variables:\n", + "x y z \n", + "\n", + "augmented coefficient matrix (last column is the constant term):\n", + "\n", + "(-0.499635,-0.289307) (-0.269638,0.143842) (0.594747,0.468559) (-0.597089,0.813991)\n", + " (-0.112562,0.566271) (-0.403476,0.482767) (-0.39012,0.344406) (0.779846,0.179663)\n", + "\n", + "slice is not homogeneous\n", + "\n", + "rand_slice: linear slice on 3 variables:\n", + "x y z \n", + "\n", + "augmented coefficient matrix (last column is the constant term):\n", + "\n", + "(-0.346547,-0.461778) (0.393155,0.605749) (0.380316,-0.0228714) (0.15641,0.452016)\n", + " (-0.513664,0.263596) (-0.661661,0.139355) (0.292895,-0.351658) (-0.470703,-0.911298)\n", + "\n", + "slice is not homogeneous\n", + "\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((5479225958647425/8376666887651221,838789883576953/770639950880571)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 15dd7ea74fbc9743, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((5479225958647425/8376666887651221,838789883576953/770639950880571)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 307e7d9f5fbe93f7, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((5479225958647425/8376666887651221,838789883576953/770639950880571)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 3434b77b8aa65fb5, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((5479225958647425/8376666887651221,838789883576953/770639950880571)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 31b25fa3e818e370, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((5479225958647425/8376666887651221,838789883576953/770639950880571)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 18816cd55de9e978, 0 recalled, records at bertini_output)\n", + "new start points: [] [] [] [] []\n", + "start_sys:\n", + " [(-1.0421075930853945000020e+00, 1.3575679903465262999980e+00)\n", + " (1.1196502975482482999980e+00, 1.2146505795945861000000e-01)\n", + " (-4.5905432682353849999980e-01, -5.8915688836732570000050e-01)] [(1.2066938925426353000010e-01, 4.6983945322518605000070e-01)\n", + " (-2.5258262307363083000070e-01, 1.3030530841713234999970e+00)\n", + " (5.3674869763359160000100e-01, -6.8625177450125760000060e-01)] [(-2.2410193566864312000040e-01, 1.5548508551481523999990e+00)\n", + " (9.4485329646279010000100e-02, -6.2783987296640339999920e-02)\n", + " (-5.4613936738801480000000e-02, -1.7973703369834137999970e-01)] [(-1.9574953655444255999970e+00, 3.4344881486158657000020e-01)\n", + " (2.3245977432390080000050e+00, 1.3079889924555659999980e+00)\n", + " (-6.2579092806664770000040e-01, -1.5085079986100937999980e+00)] [(-3.9466412111166010000020e-02, 4.7342747385981593000080e-01)\n", + " (-5.4964425796656279999910e-02, 1.2869650186154873000010e+00)\n", + " (4.4236570472771129999990e-01, -7.4186445690234450000060e-01)]\n", + "codim 3\n", + "rand_sys: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "1 function:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + "\n", + "slices_prev_step: linear slice on 3 variables:\n", + "x y z \n", + "\n", + "augmented coefficient matrix (last column is the constant term):\n", + "\n", + "(-0.499635,-0.289307) (-0.269638,0.143842) (0.594747,0.468559) (-0.597089,0.813991)\n", + " (-0.112562,0.566271) (-0.403476,0.482767) (-0.39012,0.344406) (0.779846,0.179663)\n", + "\n", + "slice is not homogeneous\n", + "\n", + "rand_slice: linear slice on 3 variables:\n", + "x y z \n", + "\n", + "augmented coefficient matrix (last column is the constant term):\n", + "\n", + "(-0.354744,-0.455511) (-0.364772,-0.026187) (0.00124255,0.730014) (-0.236567,-0.621582)\n", + " (-0.298298,0.49432) (-0.118188,0.453867) (-0.65716,0.121832) (0.709366,0.868087)\n", + "\n", + "slice is not homogeneous\n", + "\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((-5570167659392597/424075598455468,4884532633245/2431342296532)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run ee13c3711174fdd1, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((-5570167659392597/424075598455468,4884532633245/2431342296532)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 2a45151bbbe50476, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((-5570167659392597/424075598455468,4884532633245/2431342296532)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 0d77f5d13d414e80, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((-5570167659392597/424075598455468,4884532633245/2431342296532)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 4e02e4bc4d5aebf0, 0 recalled, records at bertini_output)\n", + "hom: 1 variable group:\n", + " group 0: x y z \n", + "\n", + "3 functions:\n", + " f_0 = R . g (R: 1x3 randomization matrix)\n", + " g_0 = (x^2+y^2+z^2-1)*(z-x^3)*(y-x^2)*(z-2)\n", + " g_1 = (x^2+y^2+z^2-1)*(z-x^3)*(y-2)\n", + " g_2 = (x^2+y^2+z^2-1)*(y-x^2)*(x-2)\n", + " R =\n", + " [ 1, (-0.5467+-0.4484*i), (0.2598+-0.6576*i) ]\n", + " f_1..f_2 = (1-t)*A + ((-5570167659392597/424075598455468,4884532633245/2431342296532)*t)*B (blend of 2 systems)\n", + "\n", + "path variable: t\n", + "\n", + "results = SolveResult(0 solutions, run 0754d307a6d9668f, 0 recalled, records at bertini_output)\n", + "new start points: [] [] [] [] []\n", + "start_sys:\n", + " [(-1.0421075930853945000020e+00, 1.3575679903465262999980e+00)\n", + " (1.1196502975482482999980e+00, 1.2146505795945861000000e-01)\n", + " (-4.5905432682353849999980e-01, -5.8915688836732570000050e-01)] [(1.2066938925426353000010e-01, 4.6983945322518605000070e-01)\n", + " (-2.5258262307363083000070e-01, 1.3030530841713234999970e+00)\n", + " (5.3674869763359160000100e-01, -6.8625177450125760000060e-01)] [(-2.2410193566864312000040e-01, 1.5548508551481523999990e+00)\n", + " (9.4485329646279010000100e-02, -6.2783987296640339999920e-02)\n", + " (-5.4613936738801480000000e-02, -1.7973703369834137999970e-01)] [(-1.9574953655444255999970e+00, 3.4344881486158657000020e-01)\n", + " (2.3245977432390080000050e+00, 1.3079889924555659999980e+00)\n", + " (-6.2579092806664770000040e-01, -1.5085079986100937999980e+00)] [(-3.9466412111166010000020e-02, 4.7342747385981593000080e-01)\n", + " (-5.4964425796656279999910e-02, 1.2869650186154873000010e+00)\n", + " (4.4236570472771129999990e-01, -7.4186445690234450000060e-01)]\n" + ] + } + ], + "source": [ + "#driver\n", + "\n", + "#function setup\n", + "x, y, z = b2.Variable('x'), b2.Variable('y'), b2.Variable('z')\n", + "#assuming these are in non-increasing degree order (expose sort so that users can have the system automaticaly sorted)\n", + "f1 = (y - x**2) * (x**2 + y**2 + z**2 - 1) * (x - 2)\n", + "f2 = (z - x**3) * (x**2 + y**2 + z**2 - 1) * (y - 2)\n", + "f3 = (z - x**3) * (y - x**2) * (x**2 + y**2 + z**2 - 1) * (z - 2)\n", + "\n", + "#systems setup\n", + "sys = b2.System()\n", + "vars = b2.VariableGroup([x, y, z])\n", + "sys.add_variable_group(vars)\n", + "clone_sys = sys.clone() #used to copy the base system without the added functions if needed later\n", + "sys.add_functions([f1, f2, f3])\n", + "#sys.sort()\n", + "\n", + "bpr = BasePointRegen(sys)\n", + "\n", + "bottom_slice = b2.Slice.random_complex(sys.variable_groups(), sys.num_variables() - 1)\n", + "\n", + "for codim in range(1, bpr.dim + 1):\n", + " print(f'codim {codim}')\n", + " if codim==1:\n", + " result = bpr.generate_start_points(bottom_slice)\n", + " # print(f\"degree(f{codim + 1}) = {bpr.degree} -> {len(bpr.start_points)} base-level start points\")\n", + " # for p in bpr.base_points_numeric:\n", + " # print(\" \", p)\n", + " \n", + " # print(bpr.start_moving_system)\n", + " # print(bpr.end_moving_system)\n", + "\n", + " start_points = bpr.start_points\n", + " start_system = bpr.start_moving_system\n", + " target_system = bpr.end_moving_system\n", + "\n", + " # residual cannot be used to determine if a point is a solution or not based on scale\n", + " # if x is a solution where x < 10^-11 or similarly close enough to zero where we would assume it is\n", + " # then A * f(x), where A is 10^11 or larger makes x no longer a solution that is close enough to 0,\n", + " # but instead 1. \n", + " # for p in bpr.start_points:\n", + " # print(\"eval:\\n\", *start_system.eval(p)) \n", + "\n", + " # print(\"\\n\\n\")\n", + "\n", + " else:\n", + " # use the endpoints from the previous codim's solve\n", + " # start_system = clone_sys.clone()\n", + " # f1 and product of linears from slicemove\n", + " # code here that takes the linear product of the slices we moved to in previous line\n", + " slices_prev_step = target_system.slices()\n", + " #print(\"target_sys\",target_system)\n", + "\n", + " target_slices = slices_prev_step[1:] # take the top slice off and replace with the real system\n", + "\n", + " #print(f\"len(slices_prev_step): {len(slices_prev_step)}\\nlen(target_slices):{len(target_slices)}\")\n", + "\n", + " new_start_points = []\n", + " print(\"rand_sys:\", bpr.randomized_sys)\n", + " print(\"slices_prev_step:\",slices_prev_step[0])\n", + " rand_slice = Slice.random_complex(vars, 2)\n", + " print(\"rand_slice:\",rand_slice)\n", + "\n", + " sm = SliceMover(bpr.randomized_sys, slices_prev_step[0], rand_slice, start_points)\n", + " \n", + " for ii in range(bpr.degree):\n", + " print(\"hom: \", sm.hom)\n", + " results = sm.move_slice()\n", + " print(\"results = \", results)\n", + " #junk removal here\n", + " new_start_points.append(results.solutions)\n", + " sm.set_end_slice(Slice.random_complex(vars, 2))\n", + "\n", + " print(\"new start points: \", *new_start_points)\n", + "\n", + "\n", + " # do the homotopy solve to solve the first equation in the system\n", + " print(\"start_sys:\\n\", *bpr.start_points)\n", + "\n", + " hom = b2.nag_algorithm.blend_homotopy(target = target_system, start= start_system)\n", + " homSolver = b2.HomotopySolver(target=target_system, homotopy=hom, start_points=start_points)\n", + " solve_result = homSolver.solve()\n", + "\n", + " #print(solve_result)\n", + "\n", + " #when checking if == 0, use abs(x) > 1 * 10^-13\n", + " #bpr.set_up_partially_regen_system(codim)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/python/drafts/nid/SliceMover.ipynb b/python/drafts/nid/SliceMover.ipynb new file mode 100644 index 000000000..46cc76d3a --- /dev/null +++ b/python/drafts/nid/SliceMover.ipynb @@ -0,0 +1,376 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c7a92eee", + "metadata": {}, + "source": [ + "## 1. `SliceMover.ipynb`\n", + " \n", + "### Purpose\n", + "Given a polynomial system and a set of points lying on\n", + "`system ∩ start_slice`, `SliceMover` deforms the linear slice from\n", + "`start_slice` to `end_slice` via a continuation homotopy, tracking each\n", + "witness point to its corresponding point on `system ∩ end_slice`.\n", + " \n", + "### Dependencies\n", + "```python\n", + "import bertini as b2\n", + "from bertini import nag_algorithm\n", + "from bertini import multiprec\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.patches as mpatches\n", + "import sympy\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "8a307586-28f2-4187-98e7-685986d17036", + "metadata": {}, + "outputs": [], + "source": [ + "import bertini as b2\n", + "from bertini import nag_algorithm\n", + "from bertini import multiprec\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.patches as mpatches\n", + "import sympy\n", + "\n", + "#import eigenpy #to find mpfr precision" + ] + }, + { + "cell_type": "markdown", + "id": "f723dd52", + "metadata": {}, + "source": [ + "### Class `SliceMover`\n", + " \n", + "```python\n", + "SliceMover(new_sys, start_slice, end_slice, start_points,\n", + " new_t_start=1, new_t_end=0, new_gamma=None)\n", + "```\n", + " \n", + "| Parameter | Meaning |\n", + "|---|---|\n", + "| `new_sys` | The polynomial system being intersected with the moving slice. |\n", + "| `start_slice`, `end_slice` | Either a `Slice` object (must expose `.as_system()`) or an already-built `System`. Handled polymorphically by `slice_to_sys`. |\n", + "| `start_points` | Points on `new_sys ∩ start_slice`, to be tracked as the slice moves. |\n", + "| `new_t_start`, `new_t_end` | Homotopy path-parameter endpoints (default full range `1 → 0`). |\n", + "| `new_gamma` | Optional fixed gamma constant for the gamma-trick (used to avoid path crossings/singularities during tracking). If omitted, a random rational gamma is drawn via `b2.symbolics.Rational.rand()`. |\n", + " \n", + "**Constructor behavior:** stores the system and points, converts both\n", + "slices to system form, initializes a path-tracking variable\n", + "`self.tracking_var = b2.Variable('t')`, and immediately builds the\n", + "homotopy by calling `set_up_homotopy()`.\n", + " \n", + "#### Methods\n", + " \n", + "- **`slice_to_sys(start_slice, end_slice)`**\n", + " Tries `slice.as_system()` on each argument; if that fails (e.g. the\n", + " argument is already a `System`, not a `Slice`), falls back to using it\n", + " as-is. Sets `self.start_sys` / `self.end_sys`.\n", + "- **`set_sys(new_sys)`**\n", + " Replaces `self.system` and rebuilds the homotopy.\n", + "- **`set_start_slice(start_slice)`** / **`set_end_slice(end_slice)`**\n", + " Update one side of the slice pair, re-derive `start_sys`/`end_sys`, and\n", + " rebuild the homotopy. `set_end_slice` is the one reused in\n", + " `RegenCascade.ipynb` to chain multiple slice moves in sequence.\n", + "- **`set_up_homotopy()`**\n", + " Builds `self.hom = nag_algorithm.moving_homotopy(self.system, self.start_sys, self.end_sys, gamma=self.gamma)` and attaches the tracking\n", + " variable via `self.hom.add_path_variable(self.tracking_var)`. This is\n", + " the homotopy `H(x, t)` that literally moves the slice from start to end.\n", + "- **`move_slice()`**\n", + " Builds the true target system as\n", + " `b2.system.concatenate(self.system, self.end_sys)`, constructs a\n", + " `b2.HomotopySolver(self.hom, self.start_points, target)`, solves it, and\n", + " stores/returns the result (`self.results`).\n", + "- **`__repr__`**\n", + " Prints the system and both slice-systems for debugging." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab4dabf2-66e1-43fd-a277-f827fbdb134a", + "metadata": {}, + "outputs": [], + "source": [ + "class SliceMover():\n", + " def __init__(self, new_sys, start_slice, end_slice, start_points, new_t_start = 1, new_t_end = 0, new_gamma = None):\n", + " self.system = new_sys.clone() #copies so changes outside of SliceMover affect it\n", + " self.start_sys = None\n", + " self.end_sys = None\n", + " \n", + " self.slice_to_sys(start_slice, end_slice)\n", + "\n", + " self.hom = None\n", + " \n", + " self.start_points = start_points\n", + " self.results = None\n", + " \n", + " if new_gamma == None:\n", + " self.gamma = b2.symbolics.Rational.rand()\n", + " else:\n", + " self.gamma = new_gamma\n", + " \n", + " self.t_start = new_t_start\n", + " self.t_end = new_t_end\n", + " self.tracking_var = b2.Variable('t')\n", + " self.set_up_homotopy()\n", + "\n", + " def slice_to_sys(self, start_slice, end_slice):\n", + " # print(\"slice to sys\")\n", + " try:\n", + " self.start_sys = start_slice.as_system()\n", + " except:\n", + " self.start_sys = start_slice\n", + "\n", + " try:\n", + " self.end_sys = end_slice.as_system()\n", + " except:\n", + " self.end_sys = end_slice\n", + "\n", + " def set_sys(self, new_sys):\n", + " self.system = new_sys\n", + " self.set_up_homotopy()\n", + "\n", + " def set_start_slice(self, start_slice):\n", + " self.start_slice = start_slice\n", + " self.slice_to_sys(start_slice, self.end_sys)\n", + " self.set_up_homotopy()\n", + "\n", + " def set_end_slice(self, end_slice):\n", + " self.end_slice = end_slice\n", + " self.slice_to_sys(self.start_sys, end_slice)\n", + " self.set_up_homotopy()\n", + " \n", + " def __repr__(self):\n", + " s = \"system:\\n\" + str(self.system)\n", + " s += \"\\nstart slice:\\n\" + str(self.start_sys)\n", + " s += \"\\nend slice:\\n\" + str(self.end_sys)\n", + " return s\n", + "\n", + " \"\"\"\n", + " Creates a moving_homotopy, stored in self.hom\n", + " \"\"\"\n", + " def set_up_homotopy(self):\n", + " # print(\"setting up hom\")\n", + " self.hom = nag_algorithm.moving_homotopy(self.system, self.start_sys, self.end_sys, gamma=self.gamma)\n", + " self.hom.add_path_variable(self.tracking_var)\n", + "\n", + " \"\"\"\n", + " Solves the slice move using user_homotopy. Results are returned and also stored in self.results\n", + " \"\"\"\n", + " def move_slice(self):\n", + " # print(\"setting up target\")\n", + " target = b2.system.concatenate(self.system, self.end_sys)\n", + " # print(\"creating mover\")\n", + " mover = b2.HomotopySolver(self.hom, self.start_points, target)\n", + " # print(\"solving hom\", self.hom)\n", + " self.results = mover.solve()\n", + " # print(\"returning\")\n", + " return self.results" + ] + }, + { + "cell_type": "markdown", + "id": "5ef72d8c", + "metadata": {}, + "source": [ + "### Demo cell (Lemniscate of Gerono)\n", + " \n", + "The notebook includes a worked 2-variable example:\n", + " \n", + "- Curve: `f = x**4 - x**2 + y**2` (Lemniscate of Gerono).\n", + "- Two lines: `l = 2x - 7y - 1` (start slice) and `m = x - 3y` (end slice).\n", + "- Start witness points are computed independently and directly, by\n", + " solving `f ∩ l` with `nag_algorithm.ZeroDimSolver(..., mptype='adaptive')`\n", + " and taking `.finite_solutions()` — this does **not** use `SliceMover`\n", + " itself, it's just how the initial points for the demo are obtained.\n", + "- `SliceMover(sys, start_slice, end_slice, start_points_results)` is built\n", + " and `.move_slice()` is called, tracking each point on `f ∩ l` to its\n", + " counterpart on `f ∩ m`." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "fb660fdb", + "metadata": {}, + "outputs": [], + "source": [ + "# driver\n", + "\n", + "# uncomment below for demo\n", + "\"\"\"\n", + "import bertini as b2\n", + "import numpy as np\n", + "from bertini import nag_algorithm\n", + "from bertini import multiprec\n", + "\n", + "x,y = b2.Variable('x'), b2.Variable('y')\n", + "vars = b2.VariableGroup([x, y])\n", + "#Lemniscate of Gerono system\n", + "f = x**4 - x**2 + y**2\n", + "\n", + "#intersection lines from l to m\n", + "l = 2 * x - 7 * y - 1\n", + "m = x - 3 * y\n", + "\n", + "#system setup\n", + "sys = b2.System()\n", + "sys.add_variable_group(vars)\n", + "\n", + "start_slice = b2.system.clone(sys)\n", + "end_slice = b2.system.clone(sys)\n", + "\n", + "sys.add(f)\n", + "start_slice.add(l)\n", + "end_slice.add(m)\n", + "\n", + "start_points_sys = b2.system.clone(sys)\n", + "start_points_sys.add(l)\n", + "start_points = nag_algorithm.ZeroDimSolver(start_points_sys, mptype='adaptive')\n", + "start_points.solve()\n", + "start_points_results = start_points.finite_solutions()\n", + "\n", + "sm = SliceMover(sys, start_slice, end_slice, start_points_results)\n", + "results = sm.move_slice()\n", + "\n", + "\n", + "print(results)\n", + "print(results.solutions)\n", + "s = results.solutions[0]\n", + "print(s)\n", + "\"\"\"" + ] + }, + { + "cell_type": "markdown", + "id": "c0b4f14c", + "metadata": {}, + "source": [ + "### Plotting utilities\n", + " \n", + "- **`_eval_slice_y(slice_func, x_range)`** — parses the printed string form\n", + " of a linear slice equation with `sympy`, extracts the `x`/`y`\n", + " coefficients and constant term, and returns `y = -(a·x + c)/b` over\n", + " `x_range` (returns `NaN` if the line is ~vertical, i.e. `b ≈ 0`).\n", + "- **`plot_slice_move(sm, results, curve_plotter, slice1_expr=None, slice2_expr=None)`**\n", + " — draws the background curve (via `curve_plotter`), both slice lines\n", + " (dashed blue/red), the start points (blue dots) and result points (red\n", + " stars), with arrows from each start point to its tracked endpoint.\n", + "- **`lemniscate_plotter(ax)`**, **`circle_plotter(ax)`** — parametric plots\n", + " of the Lemniscate of Gerono and unit circle, used as backgrounds.\n", + "- The actual call, `plot_slice_move(sm, results, lemniscate_plotter)`, is\n", + " left commented out (\"uncomment to see graph\")." + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "a5993c6e", + "metadata": {}, + "outputs": [], + "source": [ + "#graph driver\n", + "\n", + "#uncomment to see graph (requires driver cell to run first)\n", + "\"\"\"\n", + "def _eval_slice_y(slice_func, x_range):\n", + " s = str(slice_func).split('=')[-1].strip().replace(' ', '')\n", + " x, y = sympy.symbols('x y')\n", + " expr = sympy.sympify(s.replace('^', '**'))\n", + " a = float(expr.coeff(x))\n", + " b = float(expr.coeff(y))\n", + " c = float(expr.subs([(x, 0), (y, 0)]))\n", + " if abs(b) < 1e-12:\n", + " return np.full_like(x_range, np.nan)\n", + " return -(a * x_range + c) / b\n", + "\n", + "\n", + "def plot_slice_move(sm, results, curve_plotter, slice1_expr=None, slice2_expr=None):\n", + " fig, ax = plt.subplots(figsize=(8, 8))\n", + " curve_plotter(ax)\n", + "\n", + " def real_xy(sol):\n", + " return float(complex(sol[0]).real), float(complex(sol[1]).real)\n", + "\n", + " start_points = sm.start_points\n", + " start_xys = [real_xy(s) for s in start_points]\n", + " result_xys = [real_xy(r) for r in results]\n", + "\n", + " all_xs = [x for x, y in start_xys + result_xys]\n", + " margin = 0.3\n", + " x_range = np.linspace(min(all_xs) - margin, max(all_xs) + margin, 200)\n", + "\n", + " label1 = slice1_expr or 'Start slice'\n", + " y_start = _eval_slice_y(sm.start_sys, x_range)\n", + " ax.plot(x_range, y_start, 'b--', linewidth=1.5, label=label1)\n", + "\n", + " label2 = slice2_expr or 'End slice'\n", + " y_end = _eval_slice_y(sm.end_sys, x_range)\n", + " ax.plot(x_range, y_end, 'r--', linewidth=1.5, label=label2)\n", + "\n", + " for sx, sy in start_xys:\n", + " ax.plot(sx, sy, 'bo', markersize=8)\n", + "\n", + " for (sx, sy), (rx, ry) in zip(start_xys, result_xys):\n", + " ax.plot(rx, ry, 'r*', markersize=10)\n", + " ax.annotate('', xy=(rx, ry), xytext=(sx, sy),\n", + " arrowprops=dict(arrowstyle='->', color='gray',\n", + " lw=1.2, connectionstyle='arc3,rad=0.3'))\n", + "\n", + " curve_line = next((l for l in ax.get_lines() if l.get_label() != '_nolegend_'), None)\n", + " legend_handles = []\n", + " if curve_line:\n", + " legend_handles.append(plt.Line2D([0], [0], color=curve_line.get_color(),\n", + " lw=2, label=curve_line.get_label()))\n", + " legend_handles += [\n", + " plt.Line2D([0], [0], color='b', lw=1.5, ls='--', label=label1),\n", + " plt.Line2D([0], [0], color='r', lw=1.5, ls='--', label=label2),\n", + " mpatches.Patch(color='blue', label='Start points'),\n", + " mpatches.Patch(color='red', label='End points'),\n", + " ]\n", + " ax.legend(handles=legend_handles)\n", + " ax.set_xlabel('x')\n", + " ax.set_ylabel('y')\n", + " ax.set_title('Slice Move')\n", + " ax.set_aspect('equal')\n", + " ax.axhline(0, color='k', linewidth=0.5)\n", + " ax.axvline(0, color='k', linewidth=0.5)\n", + " ax.grid(True, alpha=0.3)\n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "\n", + "def lemniscate_plotter(ax):\n", + " t = np.linspace(0, 2 * np.pi, 1000)\n", + " x = np.sin(t)\n", + " y = np.sin(t) * np.cos(t)\n", + " ax.plot(x, y, 'k-', linewidth=2, label='Lemniscate of Gerono')\n", + "\n", + "def circle_plotter(ax):\n", + " t = np.linspace(0, 2 * np.pi, 1000)\n", + " ax.plot(np.cos(t), np.sin(t), 'k-', linewidth=2, label='Unit circle')\n", + "\n", + "\n", + "# --- usage ---\n", + "plot_slice_move(sm, results, lemniscate_plotter)\n", + "\"\"\"" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}