GitHub

HypersurfaceRegions.jl

We present a Julia package for computing the regions (i.e., connected components) in the complement of an arrangement of real hypersurfaces.

Our input consists of $k$ polynomials in $n$ variables. The output is a list of all regions of the complement of the zero set of these polynomials.

For each region $C$ we report the Euler characteristic and well-chosen sample points. The list is grouped according to sign vectors $\sigma \in \{-1,+1 \}^k$, where $\sigma_i$ is the sign of the $i$th polynomial on $C$. For hyperplane arrangements, each region is uniquely characterized by its sign vector. However, in our situation, each sign vector $\sigma$ typically corresponds to multiple connected components.

Example: two concentric circles

Let us consider two concentric circles. For instance, we could take the two circles $f_1 = x^2 + y^2 - 1=0$ and $f_2=x^2 + y^2 - 4=0$ centered at the origin. To compute the regions of $\mathcal{U} = \{ u \in \mathbb{R}^2 \mid f_1(u) \cdot f_2(u) \not= 0 \}$ we can use the following code: 

julia> using HypersurfaceRegions
julia> @var x y;
julia> f_1 = x^2 + y^2 - 1;
julia> f_2 = x^2 + y^2 - 4;
julia> f = [f_1; f_2]
julia> C = regions(f)
egionsResult with 3 regions:
=============================
9 complex critical points
5 real critical points
╭──────────────┬───────────────────────╮
│ sign pattern │ regions               │
├──────────────┼───────────────────────┤
│ + +          │ number = 1            │
│              │  χ = 0, μ = [1, 1, 0] │
│ - -          │ number = 1            │
│              │  χ = 1, μ = [1, 0, 0] │
│ + -          │ number = 1            │
│              │  χ = 0, μ = [1, 1, 0] │
╰──────────────┴───────────────────────╯

The output shows that $\mathcal U$ has three regions. The first region has sign pattern $++$. This means that, on this region, both $f_1$ and $f_2$ are positive. On the second region, both $f_1$ and $f_2$ are negative, so it is the contractible region in the middle. The software correctly reports that this region has Euler characteristic 1. The other two regions each have one hole and thus have Euler characteristic 0.

Let us visualize the three regions. We use the Plots.jl package. For this, we define a function that spits out a color corresponding to a region given a point. 

function region_col(x, y)
    Ci = membership(C, [x; y]; warning = false)
    if !isnothing(Ci)
        return number(Ci)
    else
        return nothing
    end
end 

u = -4:0.05:4; v = -3:0.05:3;
X = repeat(reshape(u, 1, :), length(v), 1);
Y = repeat(v, 1, length(u));
Z = map(region_col, X, Y);

using Plots
contour(u, v, Z, fill = true, alpha = 0.25, 
                                    legend = false, 
                                    aspect_ratio = :equal)

This produces the following picture:

circles

Example: 3 random polynomials

We can set up a random example as follows.

using HypersurfaceRegions
@var v[1:3];
f = [rand_poly(Float64, v, d) for d in [2, 3, 3]];
C = regions(f)

Here, f consists of 3 random polynomials in v. The degrees of these polynomials are [2, 3, 3]. The coefficients are chosen from a Gaussian distribution.

At a certain number of rows, the output of C in terminal is cropped. To avoid this, we can use the following command.

show(C; crop = false)

Documentation: Output

Documentation: Main functions

HypersurfaceRegions.regions — Function
regions(C::RegionsResult)

Returns the vector of regions in C.

source
regions(f::Vector{Expression})
regions(f::System)

Input a list of hypersurfaces 'f = [f1,...fk]'. Outputs the regions in the complement of the hypersurface arrangement, their sign patterns, Euler characteristic and the indices of the critical points in each region.

Options:

  • bounded_check = false: if true, the algorithm also computes if regions are bounded or unbounded.
  • projective_fusion = false: if true, the algorithm computes which of the regions are fused at infinity.
  • target_parameters: Specify parameters of the System f (if its has any).
  • show_progress = true: if true, prints the progress of the computation to the terminal.
  • s: exponents of the Morse function f_1^(s_1) * ... * f_k^(s_k) * q^(s_k+1). Here, s is a list of integers [s_1, ..., s_k, s_{k+1}] such that s_1, ..., s_k>0, s_{k+1}<0 and 2 s_{k+1} > s_1 deg(f_1) + ... + s_k deg(f_k).
  • epsilon = 1e-6: how close from each critical point do we do the path tracking.
  • reltol = 1e-6, abstol = 1e-9: parameters for the accuracy of the ODE solver.
  • monodromy_options = MonodromyOptions(max_loops_no_progress = 25): pass options for monodromy.
  • start_pair_using_newton::Bool = false: if true, the algorithm tries to compute a start pair for monodromy by using Newton's methods. Can reduce the number of critical points, but is less stable.

Options for when bounded_check = true:

  • δ::Float64 = 1e-8: Parameter that defines the strip around infinity.

Example

using HypersurfaceRegions
@var x y
f = [x^2 + y^2 - 1; x^2 + y^2 - 4];
regions(f)

Example with information on boundedness

regions(f; bounded_check = true)
source
HypersurfaceRegions.projective_regions — Function
projective_regions(C::RegionsResult)

Returns a vector of vectors. The entries of the vectors are those regions in C, which are fused in projective space.

The following code will return a vector of vectors of type Region:

using HypersurfaceRegions
@var x y
f = [x^2 + y^2 - 1; x^2 + y^2 - 4];
C = regions(f)
p = projective_regions(C)
source

Documentation: Helper functions

Functions to call on RegionsResult.

HypersurfaceRegions.nbounded — Function
nbounded(C::RegionsResult)

Returns the number of regions in C that were labelled bounded (if that information was computed).

source
HypersurfaceRegions.nundecided — Function
nundecided(C::RegionsResult)

Returns the number of regions in C, where bounded or unbounded could not be decided (if that information was computed).

source

Functions to call on Region.