diff --git a/Project.toml b/Project.toml index c987a09..211fdd8 100644 --- a/Project.toml +++ b/Project.toml @@ -5,9 +5,12 @@ version = "1.0.0-DEV" [deps] AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" GAP = "c863536a-3901-11e9-33e7-d5cd0df7b904" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] julia = "1.8" diff --git a/src/DecomposingPolynomialSystems.jl b/src/DecomposingPolynomialSystems.jl index d6d3ff1..bd3c8e7 100644 --- a/src/DecomposingPolynomialSystems.jl +++ b/src/DecomposingPolynomialSystems.jl @@ -8,12 +8,15 @@ export @var, Variable, Expression, System # Exports for testing # export degree +using SparseArrays: SparseVector, SparseMatrixCSC, spzeros, AbstractSparseVector, findnz +using Combinatorics: partitions, multiset_permutations, combinations using LinearAlgebra: nullspace +using UnPack: @unpack include("utils.jl") -include("monomials.jl") -include("expression_map.jl") include("sampled_system.jl") +include("monomials.jl") +# include("expression_map.jl") include("scalings.jl") # include("interpolation.jl") # include("deck_transformations.jl") diff --git a/src/deck_transformations.jl b/src/deck_transformations.jl index 9af60a0..3344a13 100644 --- a/src/deck_transformations.jl +++ b/src/deck_transformations.jl @@ -214,7 +214,6 @@ function symmetries_fixing_parameters_graded!( tols::Tolerances=Tolerances(), logging::Bool=false ) - n_unknowns, n_sols, _ = size(F.samples.solutions) # TODO: what if n_sols is huge? C = F.deck_permutations symmetries = _init_symmetries(length(C), unknowns(F)) @@ -236,9 +235,9 @@ function symmetries_fixing_parameters_graded!( g = gcd(vcat(num_mons, denom_mons)) if isone(g) && !only_param_dep(vcat(num_mons, denom_mons), Vector(1:n_unknowns)) if isnothing(eval_num_mons) - eval_num_mons = HC.evaluate(num_mons, F.samples) + eval_num_mons = evaluate(num_mons, F.samples) end - eval_denom_mons = HC.evaluate(denom_mons, F.samples) + eval_denom_mons = evaluate(denom_mons, F.samples) for (j, symmetry) in enumerate(symmetries) if ismissing(symmetry[i]) symmetry[i] = _interpolate_deck_function( @@ -300,7 +299,10 @@ function symmetries_fixing_parameters_dense!( sample_system!(F, n_instances) logging && println("Evaluating monomials...\n") - evaluated_mons = HC.evaluate(mons, F.samples) + # TODO: pick samples at which to evaluate monomials + # TODO: do it according to n_mons (influences n_samples) + # and n_param_dep_only_mons (influences n_instances) + evaluated_mons = evaluate(mons, F.samples) for (i, symmetry) in enumerate(symmetries) logging && printstyled("Interpolating the ", i, "-th symmetry map...\n"; color=:blue) diff --git a/src/interpolation.jl b/src/interpolation.jl index da7216d..44f74ee 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -4,12 +4,6 @@ export rational_function, using LinearAlgebra: norm, dot -function check_func_type(func_type::String) - if (func_type != "polynomial" && func_type != "rational") - error("func_type argument must be either \"polynomial\" or \"rational\"") - end -end - function rational_function( coeffs::AbstractVector{<:Number}, num_mons::MonomialVector, @@ -17,7 +11,7 @@ function rational_function( logging::Bool=false, tol::Float64=1e-5 ) - n_num_mons, n_denom_mons = length(num_mons.mds), length(denom_mons.mds) + n_num_mons, n_denom_mons = length(num_mons), length(denom_mons) @assert length(coeffs) == n_num_mons + n_denom_mons num, denom = coeffs[1:n_num_mons], coeffs[n_num_mons+1:end] if norm(num) < tol @@ -53,6 +47,12 @@ function polynomial_function( return p end +function check_func_type(func_type::String) + if (func_type != "polynomial" && func_type != "rational") + error("func_type argument must be either \"polynomial\" or \"rational\"") + end +end + function reconstruct_function( coeffs::AbstractVector{<:Number}, mons::MonomialVector; @@ -160,7 +160,7 @@ end # supposes each md in mds is a multidegree in both unknowns and parameters # TODO: The implementation below (with _) is more efficient (approx 2x), # TODO: since it exploits the sparsity of multidegrees. REMOVE THIS METHOD? -# function HomotopyContinuation.evaluate(mons::MonomialVector, samples::VarietySamples) +# function HomotopyContinuation.evaluate(mons::MonomialVector, samples::Samples) # solutions = samples.solutions # parameters = samples.parameters @@ -180,38 +180,6 @@ end # return evaluated_mons # end -# TODO: consider view for slices -function HC.evaluate( - mons::MonomialVector, - samples::VarietySamples; - sparse::Bool=false -) - solutions = samples.solutions - parameters = samples.parameters - - n_unknowns, n_sols, n_instances = size(solutions) - mds = mons.mds - n_mds = length(mds) - - nonzero_ids_unknowns = [findall(!iszero, md[1:n_unknowns]) for md in mds] - nonzero_ids_params = [findall(!iszero, md[n_unknowns+1:end]) for md in mds] - - evaluated_mons = zeros(CC, n_mds, n_sols, n_instances) - for i in 1:n_instances - params = view(parameters, :, i) - sols = view(solutions, :, :, i) - for (j, md) in enumerate(mds) - params_part = params[nonzero_ids_params[j]] - md_params_part = md[n_unknowns+1:end][nonzero_ids_params[j]] - params_eval = prod(params_part.^md_params_part) - sols_part = view(sols, nonzero_ids_unknowns[j], :) - md_sols_part = md[1:n_unknowns][nonzero_ids_unknowns[j]] - evaluated_mons[j, :, i] = vec(prod(sols_part.^md_sols_part, dims=1)).*params_eval - end - end - return evaluated_mons -end - function interpolate( A::AbstractMatrix{CC}, mons::MonomialVector; @@ -254,7 +222,7 @@ function interpolate( end function interpolate_vanishing_polynomials( - samples::VarietySamples, + samples::Samples, vars::Vector{Variable}, mons::MonomialVector; tol::Float64=1e-5 diff --git a/src/monomials.jl b/src/monomials.jl index 40cd19f..fb511f0 100644 --- a/src/monomials.jl +++ b/src/monomials.jl @@ -1,131 +1,353 @@ export Monomial, - MonomialVector, + DenseMonomialVector, + multiexponents, + iterate, + extend!, + evaluate, to_expressions, to_classes, only_param_dep, n_only_param_dep # TODO: remove? -# TODO: extend Number or nothing at all? -struct Monomial{T<:Integer} #<: Number - md::Vector{T} - vars::Vector{Variable} +struct Monomial{T<:Integer,S<:Integer} + unkn_mexp::SparseVector{T,S} + param_mexp::SparseVector{T,S} + unknowns::Vector{Variable} + parameters::Vector{Variable} - function Monomial{T}(md, vars) where {T<:Integer} + function Monomial{T,S}( + unkn_mexp, + param_mexp, + unknowns, + parameters + ) where {T,S<:Integer} # TODO: check if lengths are equal, if mds contains numbers >= 0 - return new(md, vars) + return new(unkn_mexp, param_mexp, unknowns, parameters) end end -Monomial(md::Vector{T}, vars::Vector{Variable}) where {T<:Integer} = Monomial{T}(md, vars) -function Monomial{T}(var::Variable, vars::Vector{Variable}) where {T<:Integer} - md = zeros(T, length(vars)) +Monomial( + md::SparseVector{T,S}, + vars::Vector{Variable} +) where {T,S<:Integer} = Monomial{T,S}(md, vars) + +function Monomial{T,S}(var::Variable, vars::Vector{Variable}) where {T,S<:Integer} + md = spzeros(T, S, length(vars)) md[findfirst(x->x==var, vars)] = 1 - return Monomial{T}(md, vars) + return Monomial{T,S}(md, vars) end -function Monomial{T}( +function Monomial{T,S}( mon::Expression, vars::Vector{Variable} -) where {T<:Integer} +) where {T,S<:Integer} es, cs = exponents_coefficients(mon, vars) if length(cs) > 1 || !isone(cs[1]) throw(ArgumentError("Input expression is not a monomial")) else - return Monomial{T}(es[:,1], vars) + return Monomial{T,S}(sparse(es[:,1]), vars) end end +prodpow(v::AbstractVector, e::AbstractSparseVector) = prod(v[e.nzind].^e.nzval) Base.isone(mon::Monomial) = iszero(mon.md) -Base.convert(::Type{Expression}, mon::Monomial) = prod(mon.vars.^mon.md) +Base.convert( + ::Type{Expression}, + mon::Monomial +) = prodpow(mon.unknowns, mon.unkn_mexp)*prodpow(mon.parameters, mon.param_mexp) Base.:(==)(m₁::Monomial, m₂::Monomial) = Expression(m₁) == Expression(m₂) Base.show(io::IO, mon::Monomial) = show(io, Expression(mon)) -mutable struct MonomialVector{T<:Integer} # <: AbstractVector{Expression} - mds::Vector{Vector{T}} - vars::Vector{Variable} - function MonomialVector{T}(mds, vars) where {T<:Integer} - # TODO: check if lengths are equal, if mds contains numbers >= 0 - return new(mds, vars) +""" + multiexponents(; degree::T, nvars::S) where {T<:Integer,S<:Integer} + +Generates vector of multiexponents of total degree `degree` each represented by `SparseVector{T,S}` of length `nvars`. +""" +function multiexponents(; degree::T, nvars::S) where {T<:Integer,S<:Integer} + mexps = [spzeros(T, S, nvars) for _ in 1:num_mons(nvars, degree)] + k = 1 + for n in 1:nvars + for part::Vector{T} in partitions(degree, n) + for vals in multiset_permutations(part, n) + for inds in combinations(S.(1:nvars), n) + mexps[k][inds] = vals + k += 1 + end + end + end end + return mexps end -MonomialVector( - mds::Vector{Vector{T}}, - vars::Vector{Variable} -) where {T<:Integer} = MonomialVector{T}(mds, vars) +abstract type AbstractMonomialVector end -Base.length(mons::MonomialVector) = length(mons.mds) -Base.getindex(mons::MonomialVector, i::Integer) = Monomial(mons.mds[i], mons.vars) -Base.getindex( - mons::MonomialVector, - inds... -) = MonomialVector(getindex(mons.mds, inds...), mons.vars) -# Base.findfirst(f::Function, mons::MonomialVector) = findfirst(f, mons.mds) -Base.findfirst(m::Monomial, mons::MonomialVector) = findfirst(x->x==m.md, mons.mds) -function Base.push!(mons::MonomialVector{T}, md::Vector{T}) where {T<:Integer} - push!(mons.mds, md) +# Structure representing vector of monomials +# Advantage: monomial evalution can be made more efficient +mutable struct DenseMonomialVector{T<:Integer,S<:Integer} <: AbstractMonomialVector + degree::T + unknowns_mexps::Dict{T, Vector{SparseVector{T,S}}} # keys are degrees + parameters_mexps::Dict{T, Vector{SparseVector{T,S}}} # keys are degrees + unknowns::Vector{Variable} + parameters::Vector{Variable} end -function Base.vcat(monVs::MonomialVector...) - # TODO: check if mons.vars are all equal - MonomialVector(vcat([mons.mds for mons in monVs]...), monVs[1].vars) +function DenseMonomialVector{T,S}(; + unknowns::Vector{Variable}, + parameters::Vector{Variable} +) where {T,S<:Integer} + + degree = zero(T) + unknowns_mexps = Dict(0 => [spzeros(T, S, length(unknowns))]) + parameters_mexps = Dict(0 => [spzeros(T, S, length(parameters))]) + return DenseMonomialVector{T,S}( + degree, + unknowns_mexps, + parameters_mexps, + unknowns, + parameters + ) end -# TODO -function Base.show(io::IO, mons::MonomialVector) - # println(io, "$(length(mons.mds))-element $(typeof(mons))") +unknowns(mons::DenseMonomialVector) = mons.unknowns +HC.parameters(mons::DenseMonomialVector) = mons.parameters +variables(mons::DenseMonomialVector) = vcat(unknowns(mons), parameters(mons)) + +nunknowns(mons::DenseMonomialVector{T,S}) where {T<:Integer,S<:Integer} = S(length(mons.unknowns)) +HC.nparameters(mons::DenseMonomialVector{T,S}) where {T<:Integer,S<:Integer} = S(length(mons.parameters)) +nvariables(mons::DenseMonomialVector{T,S}) where {T<:Integer,S<:Integer} = nunknowns(mons)+nparameters(mons) + +Base.length(mons::DenseMonomialVector) = num_mons_upto(nvariables(mons), mons.degree) +nparam_only(mons::DenseMonomialVector) = num_mons_upto(nparameters(mons), mons.degree) +nunkn_dep(mons::DenseMonomialVector) = length(mons) - nparam_only(mons) + +to_expression( + unknowns::Vector{Variable}, + unkn_mexp::AbstractSparseVector, + parameters::Vector{Variable}, + param_mexp::AbstractSparseVector +) = prodpow(unknowns, unkn_mexp)*prodpow(parameters, param_mexp) + +function to_expressions(mons::DenseMonomialVector) + return [to_expression( + mons.unknowns, + unkn_mexp, + mons.parameters, + param_mexp + ) for (unkn_mexp, param_mexp) in mons] +end + +function Base.show(io::IO, mons::DenseMonomialVector) + println(io, "$(length(mons))-element $(typeof(mons))") print(io, "[", join(to_expressions(mons), ", "), "]") end -function multidegrees_homogeneous(md::Vector{T}, n::Integer, d::T) where {T<:Integer} - n == 1 && return [vcat(md, d)] - return vcat([multidegrees_homogeneous(vcat(md, convert(T, i)), n-1, d-i) for i::T in 0:d]...) +""" + extend!(mons::DenseMonomialVector{T,S}; degree::Integer) where {T<:Integer,S<:Integer} + +Extends monomial vector `mons` by new multiexponents of total degree upto `degree`. +""" +function extend!(mons::DenseMonomialVector{T,S}; degree::Integer) where {T<:Integer,S<:Integer} + degree = convert(T, degree) + if degree > mons.degree + for d::T in mons.degree+1:degree + mons.unknowns_mexps[d] = multiexponents(degree=d, nvars=nunknowns(mons)) + mons.parameters_mexps[d] = multiexponents(degree=d, nvars=nparameters(mons)) + end + mons.degree = degree + end + return nothing end -function multidegrees_homogeneous(n::Integer, d::T) where {T<:Integer} - return multidegrees_homogeneous(Vector{T}([]), n, d) +struct DenseMonVecState + udeg::Int + uid::Int + pdeg::Int + pid::Int end -function multidegrees_affine(n::Integer, d::T) where {T<:Integer} - return vcat([multidegrees_homogeneous(n, i) for i::T in 0:d]...) +function pick_at_state( + mons::DenseMonomialVector, + state::DenseMonVecState +) + @unpack udeg, uid, pdeg, pid = state + uexps, pexps = mons.unknowns_mexps, mons.parameters_mexps + return uexps[udeg][uid], pexps[pdeg][pid] end -function MonomialVector{T}( - vars::Vector{Variable}; - degree::Integer, - homogeneous::Bool=false -) where {T<:Integer} +function next_state( + mons::DenseMonomialVector, + state::DenseMonVecState +) + @unpack udeg, uid, pdeg, pid = state + uexps, pexps = mons.unknowns_mexps, mons.parameters_mexps + if uid < length(uexps[udeg]) + return DenseMonVecState(udeg, uid+1, pdeg, pid) + elseif pid < length(pexps[pdeg]) + return DenseMonVecState(udeg, 1, pdeg, pid+1) + elseif pdeg > 0 + return DenseMonVecState(udeg, 1, pdeg-1, 1) + elseif udeg > 0 + return DenseMonVecState(udeg-1, 1, mons.degree-udeg+1, 1) + else + return nothing + end +end - degree = convert(T, degree) - if homogeneous - return MonomialVector(multidegrees_homogeneous(length(vars), degree), vars) +""" + iterate(mons::DenseMonomialVector, state::DenseMonVecState) + +Returns the next item (multiexponent) together with the next state. +The state is defined to be a tuple of (udeg, uid, pdeg, pid), where +`udeg` and `pdeg` are the keys for `mons.unknowns_mexps` and `mons.parameters_mexps`, +respectively, and `uid` and `pid` are the indices in the corresponding values of +these dictionaries. +""" +function Base.iterate( + mons::DenseMonomialVector, + state::DenseMonVecState +) + new_state = next_state(mons, state) + isnothing(new_state) && return nothing + return pick_at_state(mons, new_state), new_state +end + +""" + iterate(mons::DenseMonomialVector) + +Returns the first item (multiexponent) together with the state. +""" +function Base.iterate(mons::DenseMonomialVector) + state = DenseMonVecState(mons.degree, 1, 0, 1) + return pick_at_state(mons, state), state +end + +# Structure representing sparse monomial vector +# Advantages: +# 1) monomial evalution can be made more efficient +# 2) we know how many param_only monomials are there (needed for picking samples for evaluation) +struct SparseMonomialVector{T<:Integer,S<:Integer} <: AbstractMonomialVector + unkn_dep_mexps::NTuple{2, Vector{SparseVector{T,S}}} + param_only_mexps::Vector{SparseVector{T,S}} + unknowns::Vector{Variable} + parameters::Vector{Variable} +end + +SparseMonomialVector{Tv,Ti}( + unknowns::Vector{Variable}, + parameters::Vector{Variable} +) where {Tv<:Integer,Ti<:Integer} = SparseMonomialVector{Tv,Ti}(([], []), [], unknowns, parameters) + +Base.length(mons::SparseMonomialVector) = length(mons.unkn_dep_mexps[1]) + length(mons.param_only_mexps) +is_param_only(mons::SparseMonomialVector) = length(mons.unkn_dep_mexps[1]) == 0 +nunkn_dep(mons::SparseMonomialVector) = length(mons.unkn_dep_mexps[1]) +nparam_only(mons::SparseMonomialVector) = length(mons.param_only_mexps) + +function Base.push!( + mons::SparseMonomialVector{T,S}, + (uexp, pexp)::NTuple{2, SparseVector{T,S}} +) where {T<:Integer,S<:Integer} + + if iszero(uexp) + push!(mons.param_only_mexps, pexp) + else + push!(mons.unkn_dep_mexps[1], uexp) + push!(mons.unkn_dep_mexps[2], pexp) end - return MonomialVector(multidegrees_affine(length(vars), degree), vars) end -function to_expressions(mons::MonomialVector) - nonzero_ids = [findall(!iszero, md) for md in mons.mds] - return [prod(mons.vars[ids].^md[ids]) for (md, ids) in zip(mons.mds, nonzero_ids)] +function to_expressions(mons::SparseMonomialVector) + return vcat( + [prodpow(mons.unknowns, unkn_md)*prodpow(mons.parameters, param_md) for (unkn_md, param_md) in mons.mixed_mds], + [prodpow(mons.parameters, md) for md in mons.param_only_mexps] + ) +end + +# TODO +function Base.show(io::IO, mons::SparseMonomialVector) + # println(io, "$(length(mons.mds))-element $(typeof(mons))") + print(io, "[", join(to_expressions(mons), ", "), "]") +end + +function Base.gcd(mons::SparseMonomialVector) + return Monomial( + min.(mons.unkn_dep_mexps[1]...), + min.(mons.unkn_dep_mexps[2]..., mons.param_only_mexps...), + mons.unknowns, + mons.parameters + ) +end + +# Structure for evaluated monomials +# Advatage: requires less memory (no need to duplicate values for param_only monomials) +struct EvaluatedMonomials{T <: AbstractMonomialVector} + unkn_dep::Array{ComplexF64, 3} + param_only::Array{ComplexF64, 2} + mons::T end -function Base.gcd(mons::MonomialVector) - return Monomial(vec(minimum(hcat(mons.mds...); dims=2)), mons.vars) +# TODO: test timings w/ and w/o selectdim +prodpow( + v::AbstractArray, + e::AbstractSparseVector; +) = dropdims(prod(selectdim(v,1,e.nzind).^e.nzval; dims=1); dims=1) + +function HC.evaluate( + mexps_dict::Dict{T, Vector{SparseVector{T,S}}}, + samples::Array{P,N} +) where {T, S, P, N} + + evals = Dict{T, Array{P,N}}() + for (d, mexps) in mexps_dict + eval = zeros(P, length(mexps), size(samples)[2:end]...) + for (i, mexp) in enumerate(mexps) + eval[i, repeat([:], N-1)...] = prodpow(samples, mexp) + end + evals[d] = eval + end + return evals end -only_param_dep(md::Vector{<:Integer}, unknowns_ids::Vector{Int}) = iszero(md[unknowns_ids]) -only_param_dep(mon::Monomial, unknowns_ids::Vector{Int}) = only_param_dep(mon.md, unknowns_ids) -only_param_dep( - mds::Vector{Vector{T}}, - unknowns_ids::Vector{Int} -) where {T<:Integer} = all([only_param_dep(md, unknowns_ids) for md in mds]) -only_param_dep( - mons::MonomialVector, - unknowns_ids::Vector{Int} -) = only_param_dep(mons.mds, unknowns_ids) +function HC.evaluate( + mons::DenseMonomialVector, + samples::Samples +) + unkn_evals = evaluate(mons.unknowns_mexps, samples.solutions) + param_evals = evaluate(mons.parameters_mexps, samples.parameters) + unkn_dep = zeros(ComplexF64, nunkn_dep(mons), nsolutions(samples), ninstances(samples)) + param_only = zeros(ComplexF64, nparam_only(mons), ninstances(samples)) + k = 0 + for udeg in mons.degree:-1:1 + unkn_eval = unkn_evals[udeg] + for pdeg in (mons.degree-udeg):-1:0 + param_eval = param_evals[pdeg] + for pid in axes(param_eval, 1) + unkn_dep[(k+1):(k+size(unkn_eval, 1)), :, :] = unkn_eval.*reshape(param_eval[pid, :], 1, 1, :) + k += size(unkn_eval, 1) + end + end + end + param_only = vcat([param_evals[pdeg] for pdeg in mons.degree:-1:0]...) + return EvaluatedMonomials(unkn_dep, param_only, mons) +end + +function HC.evaluate( + mons::SparseMonomialVector, + samples::Samples +) + param_only = zeros(ComplexF64, nparam_only(mons), ninstances(samples)) + for (i, mexp) in enumerate(mons.param_only_mexps) + param_only[i, :] = prodpow(samples.parameters, mexp) + end + unkn_dep = zeros(ComplexF64, nmixed(mons), nsolutions(samples), ninstances(samples)) + for (i, (uexp, pexp)) in enumerate(zip(mons.unkn_dep_mexps...)) + unkn_dep[i, :, :] = prodpow(samples.solutions, uexp).*prodpow(samples.parameters, pexp)' + end + return EvaluatedMonomials(unkn_dep, param_only, mons) +end -n_only_param_dep( - mons::MonomialVector, - unknown_ids::Vector{Int} -) = sum([only_param_dep(md, unknown_ids) for md in mons.mds]) \ No newline at end of file +HC.evaluate( + mons::AbstractMonomialVector, + samples::Dict{Vector{Int}, Samples} +) = Dict(zip(keys(samples), [evaluate(mons, s) for (_, s) in samples])) diff --git a/src/sampled_system.jl b/src/sampled_system.jl index dcdeb3c..c6b29eb 100644 --- a/src/sampled_system.jl +++ b/src/sampled_system.jl @@ -1,5 +1,6 @@ export SampledSystem, - VarietySamples, + MonodromyInfo, + Samples, run_monodromy, sample_system!, unknowns, @@ -11,10 +12,7 @@ export SampledSystem, nsolutions, nsamples, ninstances, - monodromy_solutions, - monodromy_parameters, - monodromy_samples, - tracked_samples, + samples, monodromy_permutations, block_partitions, deck_permutations @@ -41,8 +39,9 @@ Samples( params::Vector{ComplexF64} ) = Samples(reshape(sols, size(sols)..., 1), reshape(params, :, 1)) -nsamples(samples::Samples) = size(samples.solutions, 2)*size(samples.solutions, 3) +HC.nsolutions(samples::Samples) = size(samples.solutions, 2) ninstances(samples::Samples) = size(samples.parameters, 2) +nsamples(samples::Samples) = nsolutions(samples)*ninstances(samples) mutable struct SampledSystem system::System @@ -129,7 +128,7 @@ function Base.show(io::IO, F::SampledSystem) print(io, " deck permutations: $(length(deck_permutations(F)))") end -function random_instance(samples::Samples) +function random_samples(samples::Samples) instance_id = rand(1:ninstances(samples)) return M2VV(samples.solutions[:, :, instance_id]), samples.parameters[:, instance_id] end @@ -230,7 +229,7 @@ function sample_system!( sols, params = extract_samples(res, F; resample=true) F.samples[path_ids] = Samples(sols, params) else - sols₀, p₀ = random_instance(samples) + sols₀, p₀ = random_samples(samples) res = HC.solve( F.system, sols₀, diff --git a/src/scalings.jl b/src/scalings.jl index 487e78c..180d7c2 100644 --- a/src/scalings.jl +++ b/src/scalings.jl @@ -5,25 +5,33 @@ export Grading, using AbstractAlgebra: ZZ, matrix, GF, lift, hnf, snf_with_transform using LinearAlgebra: diag -mutable struct Grading - free_part::Union{Nothing, Matrix{Int}} - mod_part::Vector{Tuple{Int, Matrix{Int}}} +mutable struct Grading{Tv<:Integer,Ti<:Integer} + nscalings::Int + free_part::Union{Nothing, SparseMatrixCSC{Tv,Ti}} + mod_part::Vector{Tuple{Tv, SparseMatrixCSC{Tv,Ti}}} end -Grading() = Grading(nothing, []) +Grading{Tv,Ti}() where {Tv<:Integer,Ti<:Integer} = Grading{Tv,Ti}(0, nothing, []) -function Grading(s::Vector{Int}, U::Vector{Matrix{Int}}) - grading = Grading() +function Grading{Tv,Ti}( + s::Vector{Int}, + U::Vector{Matrix{Int}} +) where {Tv<:Integer,Ti<:Integer} + + grading = Grading{Tv,Ti}() for (sᵢ, Uᵢ) in zip(s, U) if sᵢ == 0 grading.free_part = Uᵢ else push!(grading.mod_part, (sᵢ, Uᵢ)) end + grading.nscalings += size(Uᵢ, 1) end return grading end +nfree(grading::Grading) = isnothing(grading.free_part) ? 0 : size(grading.free_part, 1) +nscalings(grading::Grading) = grading.nscalings Base.isempty(grading::Grading) = isnothing(grading.free_part) && isempty(grading.mod_part) Base.copy(grading::Grading) = Grading(grading.free_part, grading.mod_part) @@ -50,16 +58,21 @@ end SparseAction = Vector{Tuple{Variable, Expression}} -struct ScalingGroup - grading::Grading +struct ScalingGroup{Tv<:Integer,Ti<:Integer} + grading::Grading{Tv,Ti} structure::String vars::Vector{Variable} - action::Tuple{Vector{SparseAction}, Vector{Tuple{Int, Vector{SparseAction}}}} + action::Tuple{Vector{SparseAction}, Vector{Tuple{Tv, Vector{SparseAction}}}} end -ScalingGroup(vars::Vector{Variable}) = ScalingGroup(Grading(), "N/A", vars, ([], [])) +ScalingGroup{Tv,Ti}( + vars::Vector{Variable} +) where {Tv<:Integer,Ti<:Integer} = ScalingGroup(Grading{Tv,Ti}(), "N/A", vars, ([], [])) -function ScalingGroup(grading::Grading, vars::Vector{Variable}) +function ScalingGroup( + grading::Grading{Tv,Ti}, + vars::Vector{Variable} +) where {Tv<:Integer,Ti<:Integer} free_action = Vector{SparseAction}([]) U₀ = grading.free_part if !isnothing(U₀) @@ -70,25 +83,25 @@ function ScalingGroup(grading::Grading, vars::Vector{Variable}) @var λ[1:size(U₀, 1)] end for j in axes(U₀, 1) - nonzero_ids = findall(!iszero, U₀[j, :]) - exprs = (λ[j].^U₀[j, :][nonzero_ids]).*vars[nonzero_ids] - push!(free_action, collect(zip(vars[nonzero_ids], exprs))) + nzind, nzval = findnz(U₀[j, :]) + exprs = (λ[j].^nzval).*vars[nzind] + push!(free_action, collect(zip(vars[nzind], exprs))) end end - mod_action = Vector{Tuple{Int, Vector{SparseAction}}}([]) + mod_action = Vector{Tuple{Tv, Vector{SparseAction}}}([]) for (sᵢ, Uᵢ) in grading.mod_part sᵢ_action = Vector{SparseAction}([]) if sᵢ == 2 for j in axes(Uᵢ, 1) - nonzero_ids = findall(!iszero, Uᵢ[j, :]) - push!(sᵢ_action, collect(zip(vars[nonzero_ids], -vars[nonzero_ids]))) + nzind, _ = findnz(Uᵢ[j, :]) + push!(sᵢ_action, collect(zip(vars[nzind], -vars[nzind]))) end else - @var ω[sᵢ] + @var ω[Int(sᵢ)] for j in axes(Uᵢ, 1) - nonzero_ids = findall(!iszero, Uᵢ[j, :]) - exprs = (ω[1].^Uᵢ[j, :][nonzero_ids]).*vars[nonzero_ids] - push!(sᵢ_action, collect(zip(vars[nonzero_ids], exprs))) + nzind, nzval = findnz(Uᵢ[j, :]) + exprs = (ω[1].^nzval).*vars[nzind] + push!(sᵢ_action, collect(zip(vars[nzind], exprs))) end end push!(mod_action, (sᵢ, sᵢ_action)) @@ -155,19 +168,41 @@ function _snf_scaling_symmetries(F::System)::Tuple{Vector{Int}, Vector{Matrix{In return s, Us end -function _hnf_reduce(grading::Grading) +function _hnf_reduce!( + s::Vector{Int}, + U::Vector{Matrix{Int}} +) + for (i, (sᵢ, Uᵢ)) in enumerate(zip(s, U)) + if sᵢ == 0 + U[i] = Matrix(hnf(matrix(ZZ, Uᵢ))) + else + try + U[i] = Int.(lift.(Matrix(hnf(matrix(GF(sᵢ), Uᵢ))))) + catch + U[i] = mod.(Uᵢ, sᵢ) + end + end + end + return s, U +end + +function _hnf_reduce(grading::Grading{Tv,Ti}) where {Tv<:Integer,Ti<:Integer} U₀ = grading.free_part - hnf_U₀ = isnothing(U₀) ? nothing : Matrix(hnf(matrix(ZZ, U₀))) - hnf_Uᵢs = Vector{Tuple{Int, Matrix{Int}}}([]) + red_grading = Grading{Tv,Ti}() + if !isnothing(U₀) + red_grading.free_part = Matrix(hnf(matrix(ZZ, U₀))) + red_grading.nscalings = size(U₀, 1) + end for (sᵢ, Uᵢ) in grading.mod_part try - Uᵢ = Int.(lift.(Matrix(hnf(matrix(GF(sᵢ), Uᵢ))))) + Uᵢ = lift.(Matrix(hnf(matrix(GF(sᵢ), Uᵢ)))) catch Uᵢ = mod.(Uᵢ, sᵢ) end - push!(hnf_Uᵢs, (sᵢ, Uᵢ)) + push!(red_grading.mod_part, (sᵢ, Uᵢ)) + red_grading.nscalings += size(Uᵢ, 1) end - return Grading(hnf_U₀, hnf_Uᵢs) + return red_grading end """ @@ -194,29 +229,31 @@ ScalingGroup isomorphic to ℤ² × ℤ₄ × ℤ₂ a ↦ -a, x ↦ -x, y ↦ -y, z ↦ -z ``` """ -function scaling_symmetries(F::System; in_hnf::Bool=true) +function scaling_symmetries(F::System) vars = vcat(F.variables, F.parameters) s, U = _snf_scaling_symmetries(F) - length(s) == 0 && return ScalingGroup(vars) - grading = Grading(s, U) - in_hnf && return ScalingGroup(_hnf_reduce(grading), vars) - return ScalingGroup(grading, vars) + length(s) == 0 && return ScalingGroup{Int8, Int16}(vars) + _hnf_reduce!(s, U) + return ScalingGroup(Grading{Int8, Int16}(s, U), vars) end -scaling_symmetries( - F::SampledSystem; - in_hnf::Bool=true -) = scaling_symmetries(F.system; in_hnf=in_hnf) +scaling_symmetries(F::SampledSystem) = scaling_symmetries(F.system) # TODO: extend to remove rows dependent on other blocks -function reduce(grading::Grading) +function reduce(grading::Grading{Tv,Ti}) where {Tv<:Integer,Ti<:Integer} hnf_grading = _hnf_reduce(grading) - red_grading = Grading() - U₀ = filter_rows(!iszero, hnf_grading.free_part) - red_grading.free_part = size(U₀, 1) == 0 ? nothing : U₀ + red_grading = Grading{Tv,Ti}() + U₀ = take_rows(!iszero, hnf_grading.free_part) + if size(U₀, 1) != 0 + red_grading.free_part = U₀ + red_grading.nscalings = size(U₀, 1) + end for (sᵢ, Uᵢ) in hnf_grading.mod_part - Uᵢ = filter_rows(!iszero, Uᵢ) - size(Uᵢ, 1) != 0 && push!(red_grading.mod_part, (sᵢ, Uᵢ)) + Uᵢ = take_rows(!iszero, Uᵢ) + if size(Uᵢ, 1) != 0 + push!(red_grading.mod_part, (sᵢ, Uᵢ)) + red_grading.nscalings += size(Uᵢ, 1) + end end return red_grading end @@ -231,7 +268,6 @@ function restrict_scalings(scalings::ScalingGroup, var_ids::Vector{Int}) return ScalingGroup(reduce(restr_grading), scalings.vars[var_ids]) end -# TODO: what if vars is not a subset of scalings.vars? function restrict_scalings(scalings::ScalingGroup, vars::Vector{Variable}) var_ids = [findfirst(v->v==var, scalings.vars) for var in vars] if nothing in var_ids @@ -249,38 +285,35 @@ scaling_symmetries( vars::Vector{Variable} ) = scaling_symmetries(F.system, vars) -# TODO: extend to scalings::ScalingGroup? -function HC.degree(md::Vector{<:Integer}, grading::Grading) +function HC.degree( + mexp::SparseVector{Tv,Ti}, + grading::Grading{Tv,Ti} +) where {Tv<:Integer,Ti<:Integer} + deg = spzeros(Tv, Ti, nscalings(grading)) U₀ = grading.free_part - deg_free = isnothing(U₀) ? Vector{Int}([]) : U₀*md - return vcat(deg_free, [mod.(Uᵢ*md, sᵢ) for (sᵢ, Uᵢ) in grading.mod_part]...) + if !isnothing(U₀) + deg[1:size(U₀,1)] = U₀*mexp + end + k = nfree(grading) + for (sᵢ, Uᵢ) in grading.mod_part + deg[(k+1):(k+size(Uᵢ,1))] = mod.(Uᵢ*mexp, sᵢ) + k += size(Uᵢ, 1) + end + return deg end -HC.degree(mon::Monomial, grading::Grading) = degree(mon.md, grading) +function to_classes( + mons::DenseMonomialVector{Tv,Ti}, + grading::Grading{Tv,Ti} +) where {Tv<:Integer,Ti<:Integer} -function to_classes(mons::MonomialVector{T}, grading::Grading) where {T<:Integer} - classes = Dict{Vector{Int}, MonomialVector{T}}() - for md in mons.mds - deg = HC.degree(md, grading) + classes = Dict{SparseVector{Tv,Ti}, SparseMonomialVector{Tv,Ti}}() + for (uexp, pexp) in mons + deg = HC.degree(vcat(uexp, pexp), grading) if isnothing(get(classes, deg, nothing)) # the key doesn't exist - classes[deg] = MonomialVector{T}([md], mons.vars) - else - push!(classes[deg], md) + classes[deg] = SparseMonomialVector{Tv,Ti}(mons.unknowns, mons.parameters) end + push!(classes[deg], (uexp, pexp)) end return classes end - -# TODO: Can we save multidegrees immediately? -# function to_classes(mons::MonomialVector, grading::Grading) -# classes = Dict{Vector{Int}, Vector{Int}}() -# for (i, md) in enumerate(mons.mds) -# deg = HC.degree(md, grading) -# if isnothing(get(classes, deg, nothing)) # the key doesn't exist -# classes[deg] = [i] -# else -# push!(classes[deg], i) -# end -# end -# return classes -# end diff --git a/src/utils/basic.jl b/src/utils/basic.jl index d0fccb9..9496233 100644 --- a/src/utils/basic.jl +++ b/src/utils/basic.jl @@ -1,6 +1,7 @@ export sparsify!, simplify_numbers, - eye, a2p, p2a + eye, a2p, p2a, + num_mons, num_mons_upto a2p(M::AbstractMatrix{<:Number}) = [M; ones(eltype(M), 1, size(M, 2))] p2a(M::AbstractMatrix{<:Number}) = (M./M[end:end,:])[1:end-1,:] @@ -18,8 +19,8 @@ xx(v) = [0 -v[3] v[2]; v[3] 0 -v[1]; -v[2] v[1] 0] xx2v(xx) = [-xx[2,3], xx[1,3], -xx[1,2]] eye(T, n::Integer) = Matrix{T}(I(n)) -num_mons(n::Integer, d::Integer) = n > 0 ? binomial(n - 1 + d, d) : 0 -num_mons_upto(n::Integer, d::Integer) = n > 0 ? binomial(n + d, d) : 0 +num_mons(n::Integer, d::Integer) = n > 0 ? binomial(Int(n - 1 + d), Int(d)) : 0 +num_mons_upto(n::Integer, d::Integer) = n > 0 ? binomial(Int(n + d), Int(d)) : 0 # TODO: test this function sparsify!(v::AbstractVector{<:Number}, tol::Real; digits::Integer=0) @@ -38,8 +39,8 @@ function sparsify!(v::AbstractVector{<:Number}, tol::Real; digits::Integer=0) end function sparsify!(M::AbstractMatrix{<:Number}, tol::Real; digits::Integer=0) - for i in axes(M, 1) - sparsify!(view(M, i, :), tol; digits=digits) + for r in eachrow(M) + sparsify!(r, tol; digits=digits) end end @@ -98,28 +99,17 @@ function superscript(n::Integer)::String return join(c) end +# TODO: eachcol? Base.findfirst( v::AbstractVector{<:Number}, M::AbstractMatrix{<:Number}; tol::Real=1e-5 ) = findfirst(i->norm(M[:,i]-v)