Utilising PySR on matrix input data

[daattavya98]

Hi! I have been trying to modify PySR for an application where my input features are pairs of matrices. That is, my input data is a list of pairs of matrices, say [(A1, B1), (A2, B2),....]. In general, the size of the matrices can vary. I am searching for symbolic expressions between A and B and want to utilize methods within Julia's LinearAlgebra package (there are predefined functions for determinant, trace etc.). As an example, my input data has the following form

x1 = Matrix{Float64}[]
x2 = Matrix{Float64}[]
for i in 1:10
    push!(x1, rand(1:10, 2, 2))
    push!(x2, rand(1:10, 2, 2))
end

X = DataFrame(x1=x1, x2=x2)
X = Tables.columntable(X)

and I would like to use the following operators for example

unary = [-, LinearAlgebra.det]
binary = [+, *, /]

and potential output symbolic expressions of the form: y = x1 + LinearAlgebra.det(x2)

I have been running into issues since the scitype of my input doesn't match the model specifications for SRRegressor. Is this possible to allow through a custom implementation? Or even passing the input data as a concatenation of vectors along with indices and then reconstructing my matrices while computing loss in a custom loss function?

My current error:

┌ Error: Problem fitting the machine machine(SRRegressor(defaults = nothing, …), …). 
└ @ MLJBase ~/.julia/packages/MLJBase/7nGJF/src/machines.jl:694
[ Info: Running type checks... 
┌ Warning: The number and/or types of data arguments do not match what the specified model
│ supports. Suppress this type check by specifying `scitype_check_level=0`.
│ 
│ Run `@doc SymbolicRegression.SRRegressor` to learn more about your model's requirements.
│ 
│ Commonly, but non exclusively, supervised models are constructed using the syntax
│ `machine(model, X, y)` or `machine(model, X, y, w)` while most other models are
│ constructed with `machine(model, X)`.  Here `X` are features, `y` a target, and `w`
│ sample or class weights.
│ 
│ In general, data in `machine(model, data...)` is expected to satisfy
│ 
│     scitype(data) <: MLJ.fit_data_scitype(model)
│ 
│ In the present case:
│ 
│ scitype(data) = Tuple{ScientificTypesBase.Table{AbstractVector{AbstractMatrix{ScientificTypesBase.Continuous}}}, AbstractVector{ScientificTypesBase.Continuous}}
│ 
│ fit_data_scitype(model) = Union{Tuple{Union{ScientificTypesBase.Table{<:Union{AbstractVector{<:ScientificTypesBase.Continuous}, AbstractVector{<:ScientificTypesBase.Count}}}, AbstractMatrix{<:ScientificTypesBase.Continuous}}, AbstractVector}, Tuple{Union{ScientificTypesBase.Table{<:Union{AbstractVector{<:ScientificTypesBase.Continuous}, AbstractVector{<:ScientificTypesBase.Count}}}, AbstractMatrix{<:ScientificTypesBase.Continuous}}, AbstractVector, AbstractVector{<:Union{ScientificTypesBase.Continuous, ScientificTypesBase.Count}}}}
└ @ MLJBase ~/.julia/packages/MLJBase/7nGJF/src/machines.jl:237
[ Info: It seems an upstream node in a learning network is providing data of incompatible scitype. See above. 
ERROR: MethodError: no method matching equation_search(::Matrix{…}, ::Vector{…}; niterations::Int64, weights::Nothing, variable_names::Vector{…}, display_variable_names::Vector{…}, options::Options{…}, parallelism::Symbol, numprocs::Nothing, procs::Nothing, addprocs_function::Nothing, heap_size_hint_in_bytes::Nothing, worker_imports::Nothing, runtests::Bool, saved_state::Nothing, return_state::Bool, run_id::Nothing, loss_type::DataType, X_units::Nothing, y_units::Nothing, verbosity::Int64, extra::@NamedTuple{}, logger::Nothing, v_dim_out::Val{…})
The function `equation_search` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  equation_search(::AbstractMatrix{T}, ::AbstractVector; kw...) where T<:Number
   @ SymbolicRegression ~/Desktop/research_projects_coding/matrix-op-dynamic-expressions/SymbolicRegression.jl/src/SymbolicRegression.jl:503
  equation_search(::Vector{D}; options, saved_state, runtime_options, runtime_options_kws...) where {T<:Number, L<:Real, D<:(Dataset{T, L, AX} where AX<:AbstractMatrix{T})}
   @ SymbolicRegression ~/Desktop/research_projects_coding/matrix-op-dynamic-expressions/SymbolicRegression.jl/src/SymbolicRegression.jl:513
  equation_search(::Dataset; kws...)
   @ SymbolicRegression ~/Desktop/research_projects_coding/matrix-op-dynamic-expressions/SymbolicRegression.jl/src/SymbolicRegression.jl:509
  ...

Reading the new version of DynamicExpressions.jl, it seemed to me that _eval_tree_array_generic should allow for input and output data of arbitrary type

Evaluate a generic binary tree (equation) over a given input data,
whatever that input data may be. The operators enum contains all
of the operators used. Unlike eval_tree_array with the normal
OperatorEnum, the array cX is sliced only along the first dimension.
i.e., if cX is a vector, then the output of a feature node
will be a scalar. If cX is a 3D tensor, then the output
of a feature node will be a 2D tensor.
Note also that tree.feature will index along the first axis of cX.

However, there is no requirement about input and output types in general.
You may set up your tree such that some operator nodes work on tensors, while
other operator nodes work on scalars. eval_tree_array will simply
return nothing if a given operator is not defined for the given input type."

Thank you again for the wonderful package!

[MilesCranmer]

So the first thing that needs to be changed is to use a GenericOperatorEnum instead of an OperatorEnum. You can see the examples at the bottom of DynamicExpressions.jl README which demonstrates some examples on arbitrary types: https://github.com/SymbolicML/DynamicExpressions.jl.

However, this still won't work out-of-the-box in SymbolicRegression.jl, because of the fact that https://github.com/MilesCranmer/SymbolicRegression.jl/blob/9fabc303dd33c624739e759d960374c7f85e56f7/src/ProgramConstants.jl#L5 is set:

const DATA_TYPE = Number

Meaning that all elements of X are constrained to be a Number. The reason this is set is because I simply haven't tried to get such a thing working yet, so I didn't want to have it working out-of-the-box for fear of people encountering some sort of bug. But in principle, it should be possible.

So the first to do is checkout a local version of SymbolicRegression.jl and make a change to:

const DATA_TYPE = Any

and then see what functions need to be updated to accommodate it. Hopefully not too many!

[daattavya98]

Thank you for the quick and helpful reply! I made the changes you suggested and am currently trying to debug further.

Here is my current code:

global x1 = Matrix{Float64}[]
global x2 = Matrix{Float64}[]
for i in 1:10
    push!(x1, rand(1:10, 2, 2))
    push!(x2, rand(1:10, 2, 2))
end

# global X = DataFrame(x1=x1, x2=x2)
# X = Tables.columntable(X)
global X = hcat(x1, x2)

global y = zeros(10)

@declare_expression_operator(LinearAlgebra.det, 1)
@declare_expression_operator(LinearAlgebra.rank, 1)

global operators = GenericOperatorEnum(; binary_operators=[+, *, /], unary_operators=[-, LinearAlgebra.det, LinearAlgebra.rank])

global model = SRRegressor(
    niterations=100,
    populations=30,
    binary_operators=operators.binops,
    unary_operators=operators.unaops,
    expression_type=Expression,
);

global mach = machine(model, X, y)
fit!(mach)
report(mach).equations[end]

The issue at the moment is constructing the Dataset (which is called on fit!(mach)).

ERROR: MethodError: no method matching one(::Type{Matrix{Float64}})
The function `one` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  one(::Type{Union{}}, Any...)
   @ Base number.jl:349
  one(::Type{Missing})
   @ Base missing.jl:104
  one(::Missing)
   @ Base missing.jl:101
  ...

Stacktrace:
  [1] Dataset(X::Matrix{…}, y::Vector{…}, loss_type::Type{…}; index::Int64, weights::Nothing, variable_names::Vector{…}, display_variable_names::Vector{…}, y_variable_name::String, extra::@NamedTuple{}, X_units::Nothing, y_units::Nothing, kws::@Kwargs{})

I am wondering if the right thing to do is to add functions for my argument type one by one as needed? This would be a bit concerning as it would require modifying core Julia files such as number.jl, missing.jl, bitarray,jl etc as that's where the function 'one' is defined.

[MilesCranmer]

This part will require some internal changes. I would just keep working through MethodErrors and changing the logic as required. I think this problem in paricular is about computing a baseline loss or something, so will need to edit Dataset(...) to change the logic somehow.

[daattavya98]

Thank you for your help with this! I managed to find an alternative solution to my original problem that doesn't require a full matrix implementation