Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fractional powers not supported by ModelKit #498

Open
jkosata opened this issue Dec 21, 2021 · 3 comments
Open

Fractional powers not supported by ModelKit #498

jkosata opened this issue Dec 21, 2021 · 3 comments

Comments

@jkosata
Copy link

jkosata commented Dec 21, 2021

julia> @var a
julia> a^(1/2)

ERROR: ^ not defined for Expression

Fractional powers are sometimes desirable for parameters in parameter homotopy.
One can fix ^ in src/model_kit/symengine.jl , but the fractional powers then still trigger AssertionError: class(n) == :Integer when solving.

@PBrdng
Copy link
Collaborator

PBrdng commented Dec 28, 2021

It's true: fractional powers are not supported by model_kit.

Can you send a short example?

Could it be an option to define a custom homotopy? See here:

abstract type AbstractHomotopy end
Base.size(H::AbstractHomotopy, i::Integer) = i == 1 ? first(size(H)) : last(size(H))
(H::AbstractHomotopy)(x, t, p = nothing) = evaluate(H, x, t, p)
function evaluate(H::AbstractHomotopy, x, t, p = nothing)
U = Vector{Any}(undef, first(size(H)))
to_smallest_eltype(evaluate!(U, H, x, t, p))
end
function jacobian(H::AbstractHomotopy, x, t, p = nothing)
n, m = size(H)
u = Vector{Any}(undef, size(H, 1))
U = Matrix{Any}(undef, size(H))
evaluate_and_jacobian!(u, U, H, x, t, p)
to_smallest_eltype(U)
end

@oameye
Copy link

oameye commented Oct 18, 2022

Could integration with Symbolics.jl, i.e., removal of the symengine dependency in favour of Symbolics.jl as suggested in #456, solve this issue?

using Symbolics
@variables a 
a^(1/2)

@saschatimme
Copy link
Member

No this is unfortunately not so straightforward. We internally have a fairly complex system for the evaluation of systems where the input are taylor polynomials. We need this for the predictor step to compute higher order derivatives of x(t). See Chapter 13 of Evaluating Derivatives from Griewank and Walther.

On the main branch we have already support for sqrt but not for general real powers. For this we need to first implement the taylor polynomial expansions of exp and log since you can rewrite x^a as exp(a * log(x)) and reuse the exp and log primitives. Here are the recursive relationships for exp and log (p. 306)

Screenshot 2022-10-19 at 21 30 11

If somebody is interested in giving this a shot I can also provide more guidance.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants