Skip to content

Commit

Permalink
Add functionality to automatically deal with rings with wrong order o…
Browse files Browse the repository at this point in the history 
…f variables
  • Loading branch information
@linus-md
linus-md committed Sep 2, 2024
1 parent 8beff34 commit e83cace
Show file tree
Hide file tree
Showing 11 changed files with 379 additions and 94 deletions.
260 changes: 210 additions & 50 deletions src/algorithms/dev.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 75,
"execution_count": 146,
"metadata": {},
"outputs": [
{
Expand All @@ -23,7 +23,7 @@
" R, variables = AlgebraicSolving.polynomial_ring(\n",
" AlgebraicSolving.GF(101),[\"x\",\"y\",\"u\",\"v\",\"l\",\"dl\"], \n",
" internal_ordering=:degrevlex)\n",
" (x,y,u,v,l,dl) = variables\n",
" (x,y,u,v,l,dl) = R_vars\n",
"\n",
" derivatives = Dict(\n",
" x => u,\n",
Expand All @@ -33,117 +33,277 @@
" l => dl)\n",
"\n",
" ideal = AlgebraicSolving.Ideal([x^2 + y^2 - 1])\n",
" return ideal, derivatives, R\n",
" return ideal, derivatives, R, R_vars\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 80,
"execution_count": 154,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"_manage_rings (generic function with 1 method)"
"_swap_vars"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"function _manage_rings(ideal, derivatives, R)\n",
"using AbstractAlgebra\n",
"using AlgebraicSolving\n",
"\n",
"\"\"\"\n",
" diff_op(q, derivatives)\n",
"\n",
" This function evaluates the linear differential operator from \n",
" Definition 7 and Notation 3 in the thesis for a polynomial and derivatives.\n",
" The inputs should not explicitly depend on time.\n",
"\n",
" # Arguments\n",
" - `q`: a polynomial\n",
" - `derivatives`: a dictionary of derivatives\n",
"\n",
" # Returns\n",
" - the linear differential operator applied to the polynomial\n",
"\"\"\"\n",
"function diff_op(q, derivatives) \n",
" n = length(q.parent.data.S)\n",
" @assert length(derivatives) <= n \"There can't be more derivatives than variables.\"\n",
" \n",
" result = 0\n",
" for (var, value) in derivatives\n",
" result += value * AbstractAlgebra.derivative(q, var)\n",
" end\n",
" return result\n",
"end\n",
"\n",
"\"\"\"\n",
" intersect(G, S_vars)\n",
"\n",
" This function computes the intersection of a Groebner basis \n",
" with a set of variables. It is important that a proper elimination\n",
" order is used to ensure that the variables that are to be eliminated\n",
" by intersection can actually be eliminated.\n",
"\n",
" An example would be lexico-graphical ordering where the variables that are\n",
" to be eliminated are the first ones. Or alternatively and superior in\n",
" practice is to use block ordering where the variables that are to be\n",
" eliminated are in the larger block w.r.t. to the ordering.\n",
"\n",
" # Arguments\n",
" - `G`: a Groebner basis w.r.t a monomial odering that eliminates all \n",
" variables except the ones in `S_vars`\n",
" - `S_vars`: the list of variables that are not to be eliminated\n",
"\n",
" # Returns\n",
" - the intersection of the Groebner basis with the subring S\n",
"\"\"\"\n",
"function intersect(G, S_vars)\n",
" sub_ideal = []\n",
" for generator in G\n",
" symbols = [Symbol(var) for var in AbstractAlgebra.vars(generator)]\n",
" if issubset(symbols, S_vars)\n",
" push!(sub_ideal, generator)\n",
" end\n",
" end\n",
" return sub_ideal\n",
"end\n",
"\n",
"\"\"\"\n",
" _manage_rings(ideal, derivatives, R)\n",
"\n",
" This function creates a new ring that permits the elimination of variables.\n",
"\n",
" # Arguments\n",
" - `ideal`: an ideal\n",
" - `derivatives`: a dictionary of derivatives\n",
" - `R`: a polynomial ring with two blocks of variables, where the first \n",
" block corresponds to the variables in `R` that are not in `derivatives` \n",
" and the second block corresponds to the variables in `derivatives`\n",
"\n",
" # Returns\n",
" - the new ring\n",
" - the new variables\n",
" - a map from the old variables to the new variables\n",
"\"\"\"\n",
"function _manage_rings(derivatives, R)\n",
" S_proper = []\n",
" R_elim = []\n",
" s = length(R.data.S) - length(derivatives)\n",
" map = Dict()\n",
" \n",
" index_n_el = s + 1\n",
" index_el = 1\n",
"\n",
" # Move the variables to the correct blocks and create a map\n",
" for var in R.data.S\n",
" if Symbol(var) in [Symbol(key) for key in keys(derivatives)]\n",
" push!(S_proper, Symbol(var))\n",
" map[Symbol(var)] = index_n_el\n",
" index_n_el += 1\n",
" else \n",
" push!(R_elim, Symbol(var))\n",
" map[Symbol(var)] = index_el\n",
" index_el += 1\n",
" end\n",
" end\n",
"\n",
" # Create new ring, where the first block is the eliminated variables\n",
" R_vars_proper = append!(R_elim, S_proper)\n",
" R_names = [string(var) for var in R_vars_proper]\n",
" R_new, R_names = AlgebraicSolving.polynomial_ring(\n",
" base_ring(R), R_names, internal_ordering=:degrevlex)\n",
" R_new_vars = [string(var) for var in R_vars_proper]\n",
" R_new, R_new_vars = AlgebraicSolving.polynomial_ring(\n",
" base_ring(R), R_new_vars, internal_ordering=:degrevlex)\n",
" \n",
" return R_new\n",
" return R_new, R_new_vars, map\n",
"end\n",
"\n",
"\n",
"\"\"\"\n",
" _swap_vars(poly, R_new, R_new_vars, map)\n",
"\n",
" This function swaps the variables in a polynomial to a new ring.\n",
"\n",
" # Arguments\n",
" - `poly`: a polynomial\n",
" - `R_new`: a new ring\n",
" - `R_new_vars`: the new variables\n",
" - `map`: a map from the old variables to the new variables\n",
"\n",
" # Returns\n",
" - the polynomial in the new ring\n",
"\"\"\"\n",
"function _swap_vars(poly, R_vars, R_new_vars, map)\n",
" vars_subst = [R_new_vars[map[Symbol(var)]] for var in R_vars]\n",
" poly_new = poly(vars_subst...)\n",
" return poly_new\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Multivariate polynomial ring in 6 variables over GF(101)[:x, :y, :u, :v, :l, :dl]\n",
"Multivariate polynomial ring in 6 variables over GF(101)[:dl, :x, :y, :u, :v, :l]\n"
]
}
],
"source": [
"ideal, derivatives, R = simple_pendulum()\n",
"println(R, R.data.S)\n",
"R_new = _manage_rings(ideal, derivatives, R)\n",
"println(R_new, R_new.data.S)"
]
},
{
"cell_type": "code",
"execution_count": 96,
"execution_count": 184,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Int64}:\n",
" 2\n",
" 0\n",
" 0\n",
" 0\n",
" 0\n",
" 0"
"differential_basis"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"s = ideal[1]\n",
"exponen = AbstractAlgebra.exponent_vector(s, 1)"
"\n",
"\n",
"\"\"\"\n",
" differential_basis(ideal, derivatives, R, nf=false, info_level=0)\n",
"\n",
" This function computes the differential Gröbner basis of an ideal or \n",
" an approximation of a full differential ideal if not all derivatives\n",
" are known. See Algorithm 4 in the thesis for reference.\n",
"\n",
" # Arguments\n",
" - `ideal`: an ideal\n",
" - `derivatives`: a dictionary of derivatives\n",
" - `R`: a polynomial ring with two blocks of variables, where the first \n",
" block corresponds to the variables in `R` that are not in `derivatives` \n",
" and the second block corresponds to the variables in `derivatives`\n",
" - `R_vars`: the variables in the ring, only for elimination, else use `[]``\n",
" - `nf`: a boolean indicating whether to compute the normal form\n",
" - `info_level`: an integer indicating the level of information to print\n",
"\n",
" # Returns\n",
" - the differential Gröbner basis\n",
"\"\"\"\n",
"function differential_basis(ideal, derivatives, R, R_vars, nf=false, info_level=0)\n",
" n = R.data.nvars\n",
" s = length(derivatives) \n",
" eliminate = n - s\n",
"\n",
" # If elimination is necessary reorganize the ring and substitute the variables\n",
" if eliminate > 0\n",
" R_new, R_new_vars, map = _manage_rings(derivatives, R)\n",
" \n",
" ideal_new_gens = [_swap_vars(elem, R_vars, R_new_vars, map) for elem in ideal.gens]\n",
" ideal = AlgebraicSolving.Ideal(ideal_new_gens)\n",
"\n",
" derivatives_new = Dict()\n",
" for (var, expr) in derivatives\n",
" derivatives_new[_swap_vars(var, R_vars, R_new_vars, map)] = _swap_vars(expr, R_vars, R_new_vars, map)\n",
" end\n",
" derivatives = derivatives_new\n",
" end \n",
"\n",
" # Infer and create the subring\n",
" S_vars = [Symbol(var) for var in R_new_vars[eliminate+1:end]]\n",
" \n",
" # Start computing the differential basis\n",
" G1 = groebner_basis(ideal)\n",
" pG1 = [diff_op(g, derivatives) for g in intersect(G1, S_vars)]\n",
" if nf\n",
" pG1 = [AlgebraicSolving.normal_form(pg, AlgebraicSolving.Ideal(G1)) for pg in pG1]\n",
" end\n",
" append!(pG1, G1)\n",
" G2 = groebner_basis(AlgebraicSolving.Ideal(pG1), eliminate=eliminate,\n",
" intersect=false, info_level=info_level)\n",
" if info_level > 0\n",
" i = 1\n",
" println(\"iteration \", i)\n",
" println(\"#G = \", length(G1))\n",
" end\n",
"\n",
" # Repeat until closed under diff_op\n",
" while G1 != G2\n",
" if info_level > 0\n",
" i += 1\n",
" println(\"iteration \", i)\n",
" println(\"#G = \", length(G2))\n",
" end\n",
" G1 = G2\n",
" pG1 = [diff_op(g, derivatives) for g in intersect(G1, S_vars)]\n",
" if nf == true\n",
" pG1 = [AlgebraicSolving.normal_form(pg, AlgebraicSolving.Ideal(G1)) for pg in pG1]\n",
" end\n",
" append!(pG1, G1)\n",
" G2 = groebner_basis(AlgebraicSolving.Ideal(pG1), eliminate=eliminate,\n",
" intersect=false, info_level=info_level)\n",
" end\n",
" return G1\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 97,
"execution_count": 186,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Symbol}:\n",
" :dl\n",
" :x\n",
" :y\n",
" :u\n",
" :v\n",
" :l"
"8-element Vector{FqMPolyRingElem}:\n",
" u^2 + v^2 + 100*y + l\n",
" x*u + y*v\n",
" x^2 + y^2 + 100\n",
" y*u*v + 100*x*v^2 + x*y + 100*x*l\n",
" y^2*u + 100*x*y*v + 100*u\n",
" y^3 + 100*y^2*l + v^2 + 100*y + l\n",
" x*y^2 + 100*x*y*l + u*v\n",
" dl + 98*v"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"R_new.data.S"
"ideal, derivatives, R, R_vars = simple_pendulum()\n",
"res = differential_basis(ideal, derivatives, R, R_vars)"
]
},
{
Expand Down
Loading

0 comments on commit e83cace

Please sign in to comment.