Converging in complex equations with absolute truth

[gm89uk]

Hi All,
I am working on a problem of measuring the distance between the apices of two triangles that share the same base, which can be a useful prognostic sign in ocular imaging.

I have derived the formula geometrically:

$\frac{\sqrt{\left(\sqrt{4a^2c^2 - (a^2 - b^2 + c^2)^2} - \sqrt{4c^2x^2 - (c^2 + x^2 - y^2)^2}\right)^2 + (a^2 - b^2 - x^2 + y^2)^2}}{2c}$
(complexity 74)

I can actually break it up into much simpler equations:

    xP = (a^2 - b^2 + c^2)/(2c)
    yP = sqrt(a^2 - xP^2)
    xQ = (x^2 - y^2 + c^2)/(2c)
    yQ = sqrt(x^2 - xQ^2)
    sqrt((xP - xQ)^2 + (yP - yQ)^2)

I was very interested to see how close SR would cope with this.
I managed to get an MSE of about 2e+03. I was wondering if anyone has solved more complicated expressions like this with SR and GP?

This code generates 50 datapoints for two triangles that share a base

using SymbolicRegression, IterTools, LoopVectorization, Bumper, MLJ, DataFrames, Random

is_valid_triangle(u,v,w) = (u+v>w) && (v+w>u) && (w+u>v)
function random_triangle(low,high)
    while true
        a,b,c = rand(low:high), rand(low:high), rand(low:high)
        is_valid_triangle(a,b,c) && return (a,b,c)
    end
end

#distance between PQ on shared base c using circle formula.
function dist_PQ(a,b,c,x,y)
    xP = (a^2 - b^2 + c^2)/(2c)
    yP = sqrt(a^2 - xP^2)
    xQ = (x^2 - y^2 + c^2)/(2c)
    yQ = sqrt(x^2 - xQ^2)
    sqrt((xP - xQ)^2 + (yP - yQ)^2) #all-in-one function: sqrt((sqrt(4*a^2*c^2 - (a^2 - b^2 + c^2)^2) - sqrt(4*c^2*x^2 - (c^2 + x^2 - y^2)^2))^2 + (a^2 - b^2 - x^2 + y^2)^2)/(2*c) complexity: 74
end
function generate_data(N=50; low=200, high=5000, seed=1) #larger triangles have higher MSE. 
    Random.seed!(seed)
    df = DataFrame(c=Float64[], a=Float64[], b=Float64[],
                   x=Float64[], y=Float64[], Distance=Float64[])
    for _ in 1:N
        a1,b1,c1 = random_triangle(low,high)
        while true
            x1,y1 = rand(low:high), rand(low:high)
            is_valid_triangle(x1,y1,c1) && begin
                push!(df, (c1,a1,b1,x1,y1, dist_PQ(a1,b1,c1,x1,y1)))
                break
            end
        end
    end
    df
end

df_triang = generate_data(50)

X = select(df_triang, :c, :a, :b, :x, :y)
y = df_triang[!, :Distance]

model = SRRegressor(
    niterations=1000000,
    binary_operators=[+, -, *, /],
    unary_operators=([square, sqrt]),
    nested_constraints=[
    sqrt => [sqrt => 5, square => 20],
    square => [sqrt => 5, square => 20]
    ],
    maxsize=100,
    bumper=true,
    turbo=true,
    populations=18,
    population_size=1000,
)
mach = machine(model, X, y)
fit!(mach)

Perhaps breaking the problem down into smaller functions with TemplateExpressions is the way to go, but it's difficult to do that if there was no priori, or much more compute might be necessary.

[gm89uk]

FYI @MilesCranmer

This relates to the discussion here #798
I am getting no convergence with default adaptive_parsimony_scaling and it only converges when this is set to 20 (v1.5.2 SR.jl).

No convergence:

using SymbolicRegression, IterTools, LoopVectorization, Bumper, MLJ, DataFrames, Random

is_valid_triangle(u,v,w) = (u+v>w) && (v+w>u) && (w+u>v)
function random_triangle(low,high)
    while true
        a,b,c = rand(low:high), rand(low:high), rand(low:high)
        is_valid_triangle(a,b,c) && return (a,b,c)
    end
end
function dist_PQ(a,b,c,x,y)
    xP = (a^2 - b^2 + c^2)/(2c)
    yP = sqrt(a^2 - xP^2)
    xQ = (x^2 - y^2 + c^2)/(2c)
    yQ = sqrt(x^2 - xQ^2)
    sqrt((xP - xQ)^2 + (yP - yQ)^2) #all-in-one function: sqrt((sqrt(4*a^2*c^2 - (a^2 - b^2 + c^2)^2) - sqrt(4*c^2*x^2 - (c^2 + x^2 - y^2)^2))^2 + (a^2 - b^2 - x^2 + y^2)^2)/(2*c) complexity: 74
end
function generate_data(N=20; low=200, high=5000, seed=1)
    Random.seed!(seed)
    df = DataFrame(c=Float64[], a=Float64[], b=Float64[],
                   x=Float64[], y=Float64[], Distance=Float64[])
    for _ in 1:N
        a1,b1,c1 = random_triangle(low,high)
        while true
            x1,y1 = rand(low:high), rand(low:high)
            is_valid_triangle(x1,y1,c1) && begin
                push!(df, (c1,a1,b1,x1,y1, dist_PQ(a1,b1,c1,x1,y1)))
                break
            end
        end
    end
    df
end

df_triang = generate_data(200)
X = select(df_triang, :c, :a, :b, :x, :y)
y = df_triang[!, :Distance]
double(x)=2*x
model = SRRegressor(
    niterations=1000000,
    binary_operators=[+, -, *, /],
    unary_operators=([square, sqrt, double]),
    nested_constraints=[
    sqrt => [sqrt => 5, square => 20],
    square => [sqrt => 5, square => 20],
    double => [double => 2],
    ],
    maxsize=120,
    bumper=true,
    turbo=true,
    populations=18,
    #adaptive_parsimony_scaling=20,
    population_size=100,
    complexity_of_constants = 1000,
)
mach = machine(model, X, y)
fit!(mach)

Changing adaptive_parsimony_scaling=20 works.