JuliaHomotopyContinuation / HomotopyContinuation.jl

A Julia package for solving systems of polynomials via homotopy continuation.
https://www.JuliaHomotopyContinuation.org
MIT License
175 stars 30 forks source link

`OP_POW` not defined #536

Closed oameye closed 1 year ago

oameye commented 1 year ago

Running some code through HarmonicBalance.jl I bump into the error UndefVarError: OP_POW not defined with stacktrace:

  [1] split_into_num_denom!(...)
    @ HomotopyContinuation.ModelKit HomotopyContinuation\src\model_kit\intermediate_representation.jl:391
  [2] process_mul!(...)
    @ HomotopyContinuation.ModelKit HomotopyContinuation\src\model_kit\intermediate_representation.jl:436
  [3] expr_to_ir_statements!(...)
    @ HomotopyContinuation.ModelKit HomotopyContinuation\src\model_kit\intermediate_representation.jl:479
  [4] expr_to_ir_statements!(...)
    @ HomotopyContinuation.ModelKit HomotopyContinuation\src\model_kit\intermediate_representation.jl:465

Indeed, investigating the source code in the model_kit module the definition of OP_POW is commented out. Is there a reason for that? Should OP_POW_INT be used instead at model_kit\intermediate_representation.jl:391?

saschatimme commented 1 year ago

The code path is only hit if the exponent is not an integer, so OP_POW_INT is most likely wrong. What type of exponents do you have?

When implementing the procedure, I wanted to deal with general exponents initially but I abandoned this approach since it would only be allowed in very specific circumstances. However, I think just supporting rational numbers or floating point values should be easily feasible.

oameye commented 1 year ago

Many thanks for responding so quickly. I have managed to make a minimum working example:

using HomotopyContinuation

@var x
@var a b

vars = HomotopyContinuation.Variable[x]
para = HomotopyContinuation.Variable[a, b]

pol = x*sqrt((a^2 - b^2))
eq = Expression[pol]

sys = System(eq, variables = vars, parameters = para)

solve(sys, start_system=:total_degree, target_parameters=Number[0.99, 0.98])

To problem seems to be the sqrt((a^2 - b^2)) in the expression, i.e., we have an exponent of 1/2.

saschatimme commented 1 year ago

Ah yes the square root is the problem, but this is actually exactly the case that I wanted to support in the first place. I will try to take a look soon :)

oameye commented 1 year ago

I have tested out the new PR and it solves my original problem I experienced in HarmonicBalance.jl. Thank you!

saschatimme commented 1 year ago

Happy to hear that :) Just tagged a new release (2.8.3) that contains the fix