Optimizing Objective Function by Reducing Variable Complexity

[kmegahed]

Hello everyone,

I need the 4 trees in the following objective function to be functions of only variables 3, 4, and 8. However, my objective function requires 18 inputs. Searching between 18 variables requires more iterations. Is there an option for a custom objective function to be of the form my_custom_objective(tree, dataset::Dataset{T,L}, options, constant_dataset::Dataset{T,L}) where {T,L}, where dataset contains only variables 3, 4, and 8 that the trees use, and constant_dataset contains the remaining variables used for calculating the objective function? This will speed up the training.

objective = """
function my_custom_objective(tree, dataset::Dataset{T,L}, options) where {T,L}
    penalty_term = L(0)
    tot_features = 1:18
    concrete_contribution = L(100)
    web_reinf_contribution = L(100)
    tol = L(0.0000000000000001)

    r_EIs = dataset.X[1, :]
    r_EIc = dataset.X[2, :]
    eta = dataset.X[3, :]
    e_ratio = dataset.X[4, :]
    e_d = dataset.X[5, :]
    e1_d = dataset.X[6, :]
    lam = dataset.X[7, :]
    slender = dataset.X[8, :]
    ee = dataset.X[9, :]
    ee1 = dataset.X[10, :]
    zf_c = dataset.X[11, :]
    zf_s = dataset.X[12, :]
    EcIc = dataset.X[13, :]
    EsIs = dataset.X[14, :]
    pc1 = dataset.X[15, :]
    ps1 = dataset.X[16, :]
    Ll = dataset.X[17, :]
    f_f = dataset.X[18, :]

    TL1 = 0.94 .+ 0.003 .* f_f ./ lam
    fac = 0.1
    TL2 = 1.15
    TR1 = 1.1 .+ slender .* fac
    TR2 = 1.0
    prop = 0.0001 .* dataset.y

    if tree.degree != 2
        penalty_term += L(10000)
    else
        left = tree.l
        if left.degree != 2
            penalty_term += L(1000)
        else
            fn = left.l
            features = [3, 4, 8]
            remain_feat = setdiff(tot_features, features)
            should_nt = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in remain_feat
            ]
            should = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in features
            ]
            penalty_term +=
                L(100) * (sum(should_nt)) + L(100) * (length(features) - sum(should))
            TL1, flag = eval_tree_array(fn, dataset.X, options)
            if !flag
                return L(Inf)
            end

            fn = left.r
            features = [3, 4, 8]
            remain_feat = setdiff(tot_features, features)
            should_nt = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in remain_feat
            ]
            should = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in features
            ]
            penalty_term +=
                L(100) * (sum(should_nt)) + L(100) * (length(features) - sum(should))
            TL2, flag = eval_tree_array(fn, dataset.X, options)
            if !flag
                return L(Inf)
            end
        end
        right = tree.r
        if right.degree != 2
            penalty_term += L(1000)
        else
            fn = right.l
            features = [3, 4, 8]
            remain_feat = setdiff(tot_features, features)
            should_nt = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in remain_feat
            ]
            should = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in features
            ]
            penalty_term +=
                L(100) * (sum(should_nt)) + L(100) * (length(features) - sum(should))
            TR1, flag = eval_tree_array(fn, dataset.X, options)
            if !flag
                return L(Inf)
            end

            fn = right.r
            features = [3, 4, 8]
            remain_feat = setdiff(tot_features, features)
            should_nt = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in remain_feat
            ]
            should = [
                any(n -> n.degree == 0 && n.constant == false && n.feature == f, fn) for
                f in features
            ]
            penalty_term +=
                L(100) * (sum(should_nt)) + L(100) * (length(features) - sum(should))
            TR2, flag = eval_tree_array(fn, dataset.X, options)
            if !flag
                return L(Inf)
            end
        end
    end

    if penalty_term < tol
        f_fc_m = 1.0
        f_fy_m = 1.0
        MD = (zf_s .* f_fy_m .+ 0.5 .* zf_c .* f_fc_m)
        f_fc = min.(1.5, max.(TL1, 0.6))
        f_fy = min.(1.5, max.(TL2, 0.6))
        TR1 = min.(2.0, max.(TR1, 0.3))
        TR2 = min.(2.0, max.(TR2, 0.3))

        pa = ps1 .+ pc1
        pc = pc1
        #println("Size of the ps1: ", size(ps1))
        #println("Size of the f_fy: ", size(f_fy))
        #println("Size of the f_fc: ", size(f_fc))
        #println("Size of the pc1: ", size(pc1))
        #println("Size of the pa: ", size(pa))
        pa = max.(pa, pc1 .* f_fc .+ f_fy .* ps1)

        b1 = TR2 .* (0.6 .+ 0.4 .* ee1 ./ ee)
        EI2 = 0.9 .* (EsIs .+ EcIc .* 0.5)
        Ncr = (π^2 .* EI2) ./ (Ll .^ 2)

        ec = max.(Ll ./ 1000, ee .* b1)

        for j in 1:length(prop)
            ij = 0
            iter = 1000
            incr = 4
            while ij < iter
                ij += incr
                nn = ij / iter * pa[j]
                m = nn * max(ee[j], ec[j] / (1 - nn / Ncr[j] * TR1[j]))
                ratto_i = (nn + pa[j] - pc[j]) / (pa[j] - pc[j] / 2)

                th = ratto_i * π
                for k in 1:100
                    thy = sin(th) + ratto_i * π
                    if abs(thy - th) < 1e-3
                        break
                    end
                    th = thy
                end

                mm = MD[j] * sin(th / 2)^3
                net = m - mm
                if net > 0.0
                    ij -= incr
                    incr *= 0.5
                    if incr < 0.13
                        break
                    end
                end
                prop[j] = nn / 1000
            end
        end
    end

    r = dataset.y ./ prop .- 1.0
    MAPE = sum(abs.(r)) ./ length(r)
    return (penalty_term + MAPE)
end

my_custom_objective
"""
from pysr import PySRRegressor
import pandas as pd
import numpy as np

PC  = pd.read_excel("da.xlsx",'M')[['D','t','L','Fy','fc','et','eb','P']]
D=PC['D'];t=PC['t'];L1=PC['L'];fy=PC['Fy'];fc=PC['fc'];et=PC['et'];eb=PC['eb'];P=PC['P']
ee = np.maximum(np.abs(et), np.abs(eb))
ee1=np.sign(et * eb) * np.minimum(np.abs(et), np.abs(eb))
Aco = (D - 2 * t) ** 2 / 4 * np.pi
Aso = D ** 2 / 4 * np.pi - Aco
zco = (D - 2 * t) ** 3 / 6; zf_c = zco * fc
zso = (D) ** 3 / 6 - zco;   zf_s = zso * fy
Ic = (D - 2 * t) ** 4 / 64 * np.pi
Is = (D) ** 4 / 64 * np.pi - Ic
Es = 210000
Ec = 22000 * ((fc + 8) / 10) ** 0.3  ;EsIs = Is * Es  ; EcIc = Ic * Ec;r_EIs=EsIs/(EsIs+EcIc);r_EIc=EcIc/(EsIs+EcIc)
pc1=Aco * fc;ps1=Aso * fy;eta=pc1/ps1;e_ratio=ee1/ee
lam = D / t / Es * fy

f_fc = 1.0
pc = pc1 * f_fc
pa = pc + ps1
EI1 = EsIs + EcIc * 0.6
Ncr = (np.pi ** 2 * EI1) / (L1 * L1)
slender = (pa / Ncr) ** 0.5

PC['r_EIs']=r_EIs;PC['r_EIc']=r_EIc;PC['eta']=eta;PC['e_ratio']=e_ratio;PC['ee']=ee;PC['ee1']=ee1;PC['e_d']=ee/D;PC['e1_d']=ee1/D
PC['Aco']=Aco;PC['Aso']=Aso;PC['EcIc']=EcIc;PC['EsIs']=EsIs;PC['lam']=lam;PC['slender']=slender
PC['pc1']=pc1;PC['ps1']=ps1;PC['Ll']=L1;PC['zf_s']=zf_s;PC['zf_c']=zf_c;PC['f_f']=fy/fc
features=['r_EIs','r_EIc','eta','e_ratio','e_d','e1_d','lam','slender','ee','ee1','zf_c','zf_s','EcIc','EsIs','pc1','ps1','Ll','f_f']

X=PC[features];y=PC['P']

modelior11=[]
for j in range(1):
    #population_size is the equation size
    model7 = PySRRegressor(niterations=200,populations=40,population_size=100,maxsize=80,
                      nested_constraints={"^":{"^":1}},parsimony= 0.02,#adaptive_parsimony_scaling=1000,
                      constraints={"^":(-1,10)},
                      binary_operators=["+", "*","^"],
                      denoise=True,loss_function=objective,model_selection ='accuracy',
                      )
    model7.fit(X,y)
    i=0
    modelior11.append(model7)
    while i <len(model7.equations_):
        print(model7.sympy(i))
        i=i+1

[MilesCranmer]

Initial tips:

I ran your loss function through @code_warntype in Julia to check whether there were any type instabilities, and it looks like there are a few:

julia> using SymbolicRegression

julia> X = randn(Float32, 5, 256); y = randn(Float32, 256);

julia> dataset = Dataset(X, y);

julia> tree = Node{Float32}(val=0.1);

julia> options = Options();

julia> @code_warntype my_custom_objective(tree, dataset, options)

This will highlight type instabilities in red which can hurt performance.

In particular it looks like it can't figure out the type of

  penalty_term::Any
  should::Vector
  should_nt::Vector

and there are also various types that are Union{Float32, Float64} which means Julia doesn't know if they are a 32-bit float or 64-bit float. To get rid of that you can use penalty_term += L(...) to convert to the right type (which in this case is Float32), because Dataset{T,L} will tell the compiler what L is.

It looks like this means the inferred return value of the custom objective is Any - which will slow things down a lot. Here is a page in the docs with more info: https://docs.julialang.org/en/v1/manual/performance-tips/#man-code-warntype.

[MilesCranmer]

Now, for passing a constant dataset, you could set a constant global variable in Julia, and then access those from your loss function? For example:

from pysr import jl
# ^jl is the Julia runtime

create_const = jl.seval("(const_name, ar) -> @eval const $(Symbol(const_name)) = convert(Array, $ar)")

This is a function in Julia, callable from Python, that creates constants with name const_name and value ar (converted to an array).

Make sure to convert the numpy array you pass this to the right type using e.g., .astype(np.float32). for example:

create_const("my_constant", np.random.randn(100).astype(np.float32))

jl.my_constant  # Now accessible within Julia

# Including within functions:
jl.seval("""
function my_fnc()
    return my_constant
end
""")
jl.my_fnc()

Therefore, you can use "my_constant" directly within the loss.

Just make sure your loss function creates a copy of the constant globals so it doesn't modify them. e.g., loc_r_EIs = copy(r_EIs).