diff --git a/Project.toml b/Project.toml index 57d393d..a24d2a0 100644 --- a/Project.toml +++ b/Project.toml @@ -3,32 +3,37 @@ uuid = "2f0a5bb0-5f4f-4f7f-b515-93a1e67611af" authors = ["Max Reissmann "] version = "0.3.4" - [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DynamicExpressions = "a40a106e-89c9-4ca8-8020-a735e8728b6b" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" - [compat] BenchmarkTools = "1" CSV = "0.10" diff --git a/README.md b/README.md index ad595ce..ab37d07 100644 --- a/README.md +++ b/README.md @@ -93,6 +93,66 @@ The repository contains an implementation of the Gene Expression Programming [1] - Remark: the tutorial folder contains notebook, that can be run with google-colab, while showing a step-by-step introduction +# How can I approximate functions involving vectors or matricies? +- To conduct a regression involving higher dimensional objects we swap the underlying evaluation from DynamicExpression.jl to Flux.jl +- Hint: By involving such objects, the performance deteriorates significantly + + ```julia +using GeneExpressionProgramming +using Random +using Tensors + +Random.seed!(1) + +#Define the iterations for the algorithm and the population size +epochs = 100 +population_size = 1000 + +#Number of features which needs to be inserted +number_features = 5 + +#define the +regressor = GepTensorRegressor(number_features, + gene_count=2, #2 works quite reliable + head_len=3) # 5 works quite reliable + +#create some testdata - testing simply on a few velocity vectors +size_test = 1000 +u1 = [randn(Tensor{1,3}) for _ in 1:size_test] +u2 = [randn(Tensor{1,3}) for _ in 1:size_test] +u3 = [randn(Tensor{1,3}) for _ in 1:size_test] + +x1 = [2.0 for _ in 1:size_test] + +x2 = [0.0 for _ in 1:size_test] + +a = 0.5 * u1 .+ x2 .* u2 + 2* u3 + +inputs = (x1,x2,u1,u2,u3) + + +@inline function loss_new(elem, validate::Bool) + if isnan(mean(elem.fitness)) || validate + model = elem.compiled_function + a_pred = model(inputs) + !isfinite(norm(a_pred)) && return (typemax(Float64),) + size(a_pred) != size(a) && return (typemax(Float64),) + size(a_pred[1]) != size(a[1]) && return (typemax(Float64),) + + loss = norm(a_pred .- a) + return (loss,) + else + return (elem.fitness,) + end +end +fit!(regressor, epochs, population_size, loss_new) +``` + +# Supported `Engines' for Symbolic Evaluation +- DynamicExpressions.jl +- Flux.jl --> in development + + # References - [1] Ferreira, C. (2001). Gene Expression Programming: a New Adaptive Algorithm for Solving Problems. Complex Systems, 13. - [2] Reissmann, M., Fang, Y., Ooi, A., & Sandberg, R. (2024). Constraining genetic symbolic regression via semantic backpropagation. arXiv. https://arxiv.org/abs/2409.07369 @@ -106,5 +166,5 @@ The repository contains an implementation of the Gene Expression Programming [1] - [x] Naming conventions! - [x] Improve usability for user interaction - [ ] Next operations: Tail flip, Connection symbol flip, wrapper class for easy usage, config class for predefinition, staggered exploration -- [ ] latest enhancements are provided in the branch 'feature/modular_error' -- [ ] Flexible underlying engine to evaluate the expressions -> Currently DynamicExpressions.jl, Flux in the future for GPU support +- [ ] nice print flux +- [ ] constant node needs to be fixed diff --git a/benchmark/Benchmark.md b/benchmark/Benchmark.md new file mode 100644 index 0000000..d7575ee --- /dev/null +++ b/benchmark/Benchmark.md @@ -0,0 +1,90 @@ +# Tensor Operation Benchmarks + +Comparing tensor computations using DynamicExpressions.jl vs custom Flux network. + +## Usage + +```bash +export JULIA_NUM_THREADS=1 +``` + +# Evaluating Performance of DynamicExpressions +- example from: https://github.com/SymbolicML/DynamicExpressions.jl (adapted for utilizing Tensors as Dtype) + +```julia +using DynamicExpressions +using DynamicExpressions: @declare_expression_operator +using BenchmarkTools +using LinearAlgebra + +# Operations +vec_add(x::Tensor, y::Tensor) = @fastmath x + y +vec_square(x::Tensor) = @fastmath dot(x,x) + +@declare_expression_operator(vec_add, 2) +@declare_expression_operator(vec_square, 1) + +# Build expression +operators = GenericOperatorEnum( + binary_operators=[vec_add], + unary_operators=[vec_square] +) +variable_names = ["x1"] +c1 = Expression(Node{T}(; val=ones(Tensor{2,3})); operators, variable_names); +expression = vec_add(vec_add(vec_square(c1), c1), c1); +X = ones(Tensor{2,3}); + +#Evalutate the expression: + +tests_n = 100000 +@show "Benchmark expression" +expression(X) # [[5.0 5.0 5.0], [5.0 5.0 5.0], [5.0 5.0 5.0]] +@btime for _ in 1:tests_n + expression(X) +end + +# 83.021 ms (1798979 allocations: 187.67 MiB) +``` + +# Evaluating performance when generated with Flux + +```julia +# create the inputs for Flux +c1_ = ones(Tensor{2,3}); +inputs = (c1_,); + +# create the arity map +arity_map = OrderedDict{Int8,Int}( + 1 => 2, # Addition + 2 => 2 # Multiplication +); + +#assign the callbacks +callbacks = Dict{Int8,Any}( + Int8(1) => AdditionNode, + Int8(2) => MultiplicationNode +); + +#define nodes +nodes = OrderedDict{Int8,Any}( + Int8(5) => InputSelector(1) +); + +rek_string = Int8[1, 1, 2, 5, 5, 5, 5]; +network = TensorRegUtils.compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0); +@show "Benchmark network" +result = network(inputs) # [[5.0 5.0 5.0], [5.0 5.0 5.0], [5.0 5.0 5.0]] +@btime for _ in 1:tests_n + result = network(inputs) +end + +#11.703 ms (998979 allocations: 59.49 MiB) + +``` + +## Conclusion +- Flux in handling higher dimension around 7.5 times faster +- Flux uses around 3.8 times less memory resources + + + diff --git a/benchmark/benchmark_djl_nn.jl b/benchmark/benchmark_djl_nn.jl new file mode 100644 index 0000000..daaca4c --- /dev/null +++ b/benchmark/benchmark_djl_nn.jl @@ -0,0 +1,81 @@ +using DynamicExpressions +using DynamicExpressions: @declare_expression_operator +using BenchmarkTools +using LinearAlgebra + +include("../src/TensorOps.jl") +using .TensorRegUtils +using Tensors +using OrderedCollections +using Flux + + +""" +Benchmark for comparing higher dim structures for tensor regression - run test with - export JULIA_NUM_THREADS=1 +example from: https://github.com/SymbolicML/DynamicExpressions.jl with changes according to utilize tensors +""" + + +T = Union{Float64,Vector{Float64},Tensor} +vec_add(x::Tensor, y::Tensor) = @fastmath x + y; +vec_square(x::Tensor) = @fastmath dot(x,x); + + +@declare_expression_operator(vec_add, 2); +@declare_expression_operator(vec_square, 1); + + +operators = GenericOperatorEnum(; binary_operators=[vec_add], unary_operators=[vec_square]); + +# Construct the expression: +variable_names = ["x1"] +c1 = Expression(Node{T}(; val=ones(Tensor{2,3})); operators, variable_names); +expression = vec_add(vec_add(vec_square(c1), c1), c1); + +X = ones(Tensor{2,3}); + + +# create the inputs for Flux +c1_ = ones(Tensor{2,3}); +inputs = (c1_,); + +# create the arity map +arity_map = OrderedDict{Int8,Int}( + 1 => 2, # Addition + 2 => 2 # Multiplication +); + +#assign the callbacks +callbacks = Dict{Int8,Any}( + Int8(1) => AdditionNode, + Int8(2) => MultiplicationNode +); + +#define nodes +nodes = OrderedDict{Int8,Any}( + Int8(5) => InputSelector(1) +); + +# Evaluate - expression +# Solution => [[5.0 5.0 5.0], [5.0 5.0 5.0], [5.0 5.0 5.0]] +tests_n = 100000 +@show "Benchmark expression" +expression(X) +@btime for _ in 1:tests_n + expression(X) +end + +#83.021 ms (1798979 allocations: 187.67 MiB) + +rek_string = Int8[1, 1, 2, 5, 5, 5, 5]; +network = TensorRegUtils.compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0); +@show "Benchmark network" +result = network(inputs) +@btime for _ in 1:tests_n + result = network(inputs) +end + +#11.703 ms (998979 allocations: 59.49 MiB) + + +#Conclusion ≈ 7 times faster than DynamicExpressions.jl for such structures \ No newline at end of file diff --git a/paper/ConstraintSymbolicRegression.md b/paper/ConstraintSymbolicRegression.md index a562296..a82a1eb 100644 --- a/paper/ConstraintSymbolicRegression.md +++ b/paper/ConstraintSymbolicRegression.md @@ -3,7 +3,8 @@ Evolutionary symbolic regression approaches are powerful tools that can approxim # Test Reproduction -The file ConstrainViaSBP.jl contains the test setup. Please follow the steps outlined there. +- The file ConstrainViaSBP.jl contains the test setup. Please follow the steps outlined there. +- Utilize release 0.4 # Source - [1] Reissmann, M., Fang, Y., Ooi, A., & Sandberg, R. (2024). Constraining genetic symbolic regression via semantic backpropagation. arXiv. https://arxiv.org/abs/2409.07369 \ No newline at end of file diff --git a/paper/ConstraintViaSBP.jl b/paper/ConstraintViaSBP.jl index 2160629..37af476 100644 --- a/paper/ConstraintViaSBP.jl +++ b/paper/ConstraintViaSBP.jl @@ -9,6 +9,21 @@ using Random using Logging using Dates using JSON +using Statistics + + +function break_condition(population, epoch) + return isclose(mean(population[1].fitness), 0.0) +end + +function loss_new(eqn::Node, operators::OperatorEnum, x_data::AbstractArray, y_data::AbstractArray) + try + y_pred = eqn(x_data, operators) + return get_loss_function("r2_score")(y_data, y_pred) + catch e + return zero(Float64) + end +end function setup_logger(log_file_path::String) @@ -88,7 +103,7 @@ function main() println(feature_names) println(case_name) phy_dims = get_feature_dims_json(case_data, feature_names, case_name) - phy_dims = Dict{Symbol, Vector{Float16}}( Symbol(x_n) => dim_n for (x_n, dim_n) in phy_dims) + phy_dims = Dict{Symbol,Vector{Float16}}(Symbol(x_n) => dim_n for (x_n, dim_n) in phy_dims) target_dim = get_target_dim_json(case_data, case_name) print(phy_dims) @@ -104,23 +119,42 @@ function main() start_time = time_ns() - regressor = GepRegressor(num_cols-1; + regressor = GepRegressor(num_cols - 1; considered_dimensions=phy_dims, entered_non_terminals=[:+, :-, :*, :/, :sqrt, :sin, :cos, :exp, :log], - max_permutations_lib=10000, rounds=7) + max_permutations_lib=10000, rounds=7, number_of_objectives=1) + + + @inline function loss_new_(elem, validate::Bool) + try + if isnan(mean(elem.fitness)) || validate + y_pred = elem.compiled_function(x_train', regressor.operators_) + return (get_loss_function("mse")(y_train, y_pred),) + else + return (elem.fitness, length(elem.expression_raw) * elem.fitness) + #return (elem.fitness,) + end + catch e + return (typemax(Float64),typemax(Float64)) + end + end + #perform the regression by entering epochs, population_size, the feature cols, the target col and the loss function fit!(regressor, epochs, population_size, x_train', y_train; x_test=x_test', y_test=y_test', - loss_fun="mse", target_dimension=target_dim) + loss_fun="mse", break_condition=break_condition) end_time = (time_ns() - start_time) / 1e9 elem = regressor.best_models_[1] + fitness_r2_train = loss_new(elem.compiled_function, regressor.operators_, x_train', y_train) + fitness_r2_test = loss_new(elem.compiled_function, regressor.operators_, x_test', y_test) + #log_results - push!(results, (seed, case_name, noise_level, elem.fitness, string(elem.compiled_function), - elem.fitness_r2_train, elem.fitness_r2_test, end_time, elem.dimension_homogene, target_dim)) + push!(results, (seed, case_name, noise_level, mean(elem.fitness), string(elem.compiled_function), + fitness_r2_train, fitness_r2_test, end_time, elem.dimension_homogene, target_dim)) - @show elem.fitness_r2_test + @show fitness_r2_test save_results_to_csv(file_name_save, results) end end diff --git a/src/Entities.jl b/src/Entities.jl index a621cc1..adb5faf 100644 --- a/src/Entities.jl +++ b/src/Entities.jl @@ -4,12 +4,6 @@ A module implementing core data structures and genetic operations for Gene Expression Programming (GEP). -# Core Types -## Symbol Types - depreacated -- `AbstractSymbol`: Base type for GEP symbols -- `BasicSymbol`: Terminal symbols (constants and variables) -- `FunctionalSymbol`: Function symbols with arithmetic operations -- `SymbolConfig`: Configuration container for all symbols ## Evolution Types - `Chromosome`: Individual solution representation @@ -52,13 +46,12 @@ Programming (GEP). # Exports ## Types -- `Chromosome`, `Toolbox`, `AbstractGepToolbox` -- `AbstractSymbol`, `FunctionalSymbol`, `BasicSymbol`, `SymbolConfig` +- `Chromosome`, `Toolbox`, `EvaluationStrategy``, `StandardRegressionStrategy`, `GenericRegressionStrategy` ## Functions ### Core Operations - `fitness`, `set_fitness!` -- `generate_gene`, `generate_preamle!`, `compile_expression!` +- `generate_gene`, `compile_expression!` - `generate_chromosome`, `generate_population` ### Genetic Operations @@ -79,18 +72,116 @@ Programming (GEP). module GepEntities -export Chromosome, Toolbox, AbstractGepToolbox +export Chromosome, Toolbox, EvaluationStrategy, StandardRegressionStrategy, GenericRegressionStrategy export fitness, set_fitness! -export generate_gene, generate_preamle!, compile_expression!, generate_chromosome, generate_population +export generate_gene, compile_expression!, generate_chromosome, generate_population export genetic_operations!, replicate, gene_inversion!, gene_mutation!, gene_one_point_cross_over!, gene_two_point_cross_over!, gene_fussion! - include("Util.jl") +include("TensorOps.jl") + using .GepUtils +using .TensorRegUtils using OrderedCollections +using DynamicExpressions + + +""" + Memory Buffer for reducing allocation during runtime! +""" +struct GeneticBuffers + alpha_operator::Vector{Int8} + beta_operator::Vector{Int8} + child_1_genes::Vector{Int8} + child_2_genes::Vector{Int8} + rolled_indices::Vector{Any} + arity_gene::Vector{Int8} +end + + +const THREAD_BUFFERS = let + default_size = 1000 + [GeneticBuffers( + zeros(Int8, default_size), + zeros(Int8, default_size), + Vector{Int8}(undef, default_size), + Vector{Int8}(undef, default_size), + Vector{Any}(undef, default_size), + Vector{Int8}(undef, default_size) + ) for _ in 1:Threads.nthreads()] +end + +function ensure_buffer_size!(head_len::Int, gene_count::Int) + gene_len = head_len * 2 + 1 + total_gene_length = gene_count - 1 + gene_count * gene_len + + for buffer in THREAD_BUFFERS + if length(buffer.alpha_operator) < total_gene_length + resize!(buffer.alpha_operator, total_gene_length) + resize!(buffer.beta_operator, total_gene_length) + resize!(buffer.child_1_genes, total_gene_length) + resize!(buffer.child_2_genes, total_gene_length) + resize!(buffer.rolled_indices, gene_count + 1) + resize!(buffer.arity_gene, gene_count * gene_len) + end + end +end + +abstract type EvaluationStrategy end + +struct StandardRegressionStrategy{T<:AbstractFloat} <: EvaluationStrategy + operators::OperatorEnum + number_of_objectives::Int + x_data::AbstractArray{T} + y_data::AbstractArray{T} + x_data_test::AbstractArray{T} + y_data_test::AbstractArray{T} + loss_function::Function + secOptimizer::Union{Function,Nothing} + break_condition::Union{Function,Nothing} + penalty::T + crash_value::T + + function StandardRegressionStrategy{T}(operators::OperatorEnum, + x_data::AbstractArray, + y_data::AbstractArray, + x_data_test::AbstractArray, + y_data_test::AbstractArray, + loss_function::Function; + secOptimizer::Union{Function,Nothing}=nothing, + break_condition::Union{Function,Nothing}=nothing, + penalty::T=zero(T), + crash_value::T=typemax(T)) where {T<:AbstractFloat} + new(operators, + 1, + x_data, + y_data, + x_data_test, + y_data_test, + loss_function, + secOptimizer, + break_condition, + penalty, + crash_value + ) + end +end + +struct GenericRegressionStrategy <: EvaluationStrategy + operators::Union{OperatorEnum,Nothing} + number_of_objectives::Int + loss_function::Function + secOptimizer::Union{Function,Nothing} + break_condition::Union{Function,Nothing} + + function GenericRegressionStrategy(operators::Union{OperatorEnum,Nothing}, number_of_objectives::Int, loss_function::Function; + secOptimizer::Union{Function,Nothing}, break_condition::Union{Function,Nothing}) + new(operators, number_of_objectives, loss_function, secOptimizer, break_condition) + end +end """ Toolbox @@ -138,19 +229,28 @@ struct Toolbox fitness_reset::Tuple preamble_syms::Vector{Int8} len_preamble::Int8 + operators_::Union{OperatorEnum,Nothing} + compile_function_::Function function Toolbox(gene_count::Int, head_len::Int, symbols::OrderedDict{Int8,Int8}, gene_connections::Vector{Int8}, callbacks::Dict, nodes::OrderedDict, gep_probs::Dict{String,AbstractFloat}; - unary_prob::Real=0.1, fitness_reset::Tuple=(Inf, NaN), preamble_syms=Int8[]) + unary_prob::Real=0.1, preamble_syms=Int8[], + number_of_objectives::Int=1, operators_::Union{OperatorEnum,Nothing}=nothing, function_complile::Function=compile_djl_datatype) + + fitness_reset = ( + ntuple(_ -> Inf, number_of_objectives), + ntuple(_ -> NaN, number_of_objectives) + ) gene_len = head_len * 2 + 1 headsyms = [key for (key, arity) in symbols if arity == 2] unary_syms = [key for (key, arity) in symbols if arity == 1] tailsyms = [key for (key, arity) in symbols if arity < 1 && !(key in preamble_syms)] - len_preamble = length(preamble_syms) == 0 ? 0 : gene_count - gen_start_indices = [gene_count + len_preamble + (gene_len * (i - 1)) for i in 1:gene_count] #depending on the usage should shift everthing + len_preamble = length(preamble_syms) + gen_start_indices = [gene_count + (gene_len * (i - 1)) for i in 1:gene_count] + ensure_buffer_size!(head_len, gene_count) new(gene_count, head_len, symbols, gene_connections, headsyms, unary_syms, tailsyms, symbols, - callbacks, nodes, gen_start_indices, gep_probs, unary_prob, fitness_reset, preamble_syms, len_preamble) + callbacks, nodes, gen_start_indices, gep_probs, unary_prob, fitness_reset, preamble_syms, len_preamble, operators_, function_complile) end end @@ -169,7 +269,6 @@ Represents an individual solution in GEP. - `fitness_r2_test::AbstractFloat`: R² score on testing - `expression_raw::Vector{Int8}`: Raw expression - `dimension_homogene::Bool`: Dimensional homogeneity -- `penalty::AbstractFloat`: Penalty value - `chromo_id::Int`: Chromosome identifier # Constructor @@ -177,29 +276,23 @@ Represents an individual solution in GEP. """ mutable struct Chromosome genes::Vector{Int8} - fitness::Union{AbstractFloat,Tuple} + fitness::Tuple toolbox::Toolbox compiled_function::Any compiled::Bool - fitness_r2_train::AbstractFloat - fitness_r2_test::AbstractFloat expression_raw::Vector{Int8} dimension_homogene::Bool - penalty::AbstractFloat chromo_id::Int function Chromosome(genes::Vector{Int8}, toolbox::Toolbox, compile::Bool=false) obj = new() obj.genes = genes - obj.fitness = toolbox.fitness_reset[2] obj.toolbox = toolbox - obj.fitness_r2_train = 0.0 - obj.fitness_r2_test = 0.0 obj.compiled = false obj.dimension_homogene = false - obj.penalty = 0.0 obj.chromo_id = -1 + obj.expression_raw = Int8[] if compile compile_expression!(obj) end @@ -229,17 +322,18 @@ Get chromosome's fitness value. # Returns Fitness value or tuple """ -function compile_expression!(chromosome::Chromosome; force_compile::Bool=false) +@inline function compile_expression!(chromosome::Chromosome; force_compile::Bool=false) if !chromosome.compiled || force_compile try expression = _karva_raw(chromosome) - expression_tree = compile_djl_datatype(expression, chromosome.toolbox.symbols, chromosome.toolbox.callbacks, - chromosome.toolbox.nodes) + expression_tree = chromosome.toolbox.compile_function_(expression, chromosome.toolbox.symbols, chromosome.toolbox.callbacks, + chromosome.toolbox.nodes, max(chromosome.toolbox.len_preamble, 1)) chromosome.compiled_function = expression_tree chromosome.expression_raw = expression chromosome.fitness = chromosome.toolbox.fitness_reset[2] chromosome.compiled = true - catch + catch e + @error "something went wrong" exception = (e, catch_backtrace()) chromosome.fitness = chromosome.toolbox.fitness_reset[1] end end @@ -267,7 +361,7 @@ Set chromosome's fitness value. - `chromosome::Chromosome`: Target chromosome - `value::AbstractFloat`: New fitness value """ -function set_fitness!(chromosome::Chromosome, value::AbstractFloat) +function set_fitness!(chromosome::Chromosome, value::Tuple) chromosome.fitness = value end @@ -316,13 +410,13 @@ Vector{Int8} representing the K-expression of the chromosome ``` """ -function _karva_raw(chromosome::Chromosome) + +@inline function _karva_raw(chromosome::Chromosome) gene_len = chromosome.toolbox.head_len * 2 + 1 gene_count = chromosome.toolbox.gene_count - len_preamble = chromosome.toolbox.len_preamble - connectionsym = @view chromosome.genes[1+len_preamble:gene_count+len_preamble-1] - genes = chromosome.genes[gene_count:end] + connectionsym = @view chromosome.genes[1:gene_count-1] + genes = chromosome.genes[gene_count:end] arity_gene_ = map(x -> chromosome.toolbox.arrity_by_id[x], genes) rolled_indices = Vector{Any}(undef, div(length(arity_gene_), gene_len) + 1) @@ -355,9 +449,10 @@ Generate a single gene for GEP. # Returns Vector{Int8} representing gene """ -function generate_gene(headsyms::Vector{Int8}, tailsyms::Vector{Int8}, headlen::Int; unarys::Vector{Int8}=[], unary_prob::Real=0.2) +@inline function generate_gene(headsyms::Vector{Int8}, tailsyms::Vector{Int8}, headlen::Int; + unarys::Vector{Int8}=[], unary_prob::Real=0.2, tensor_prob::Real=0.2) if !isempty(unarys) && rand() < unary_prob - heads = vcat(headsyms,tailsyms) + heads = vcat(headsyms, tailsyms) push!(heads, rand(unarys)) else heads = headsyms @@ -369,6 +464,7 @@ function generate_gene(headsyms::Vector{Int8}, tailsyms::Vector{Int8}, headlen:: end + """ generate_chromosome(toolbox::Toolbox) @@ -377,7 +473,7 @@ Generate a new chromosome using toolbox configuration. # Returns New Chromosome instance """ -function generate_chromosome(toolbox::Toolbox) +@inline function generate_chromosome(toolbox::Toolbox) connectors = rand(toolbox.gene_connections, toolbox.gene_count - 1) genes = vcat([generate_gene(toolbox.headsyms, toolbox.tailsyms, toolbox.head_len; unarys=toolbox.unary_syms, unary_prob=toolbox.unary_prob) for _ in 1:toolbox.gene_count]...) @@ -398,27 +494,39 @@ Generate initial population of chromosomes. # Returns Vector of Chromosomes """ -function generate_population(number::Int, toolbox::Toolbox) +@inline function generate_population(number::Int, toolbox::Toolbox) population = Vector{Chromosome}(undef, number) - for i in 1:number + + Threads.@threads for i in 1:number @inbounds population[i] = generate_chromosome(toolbox) end + return population end + @inline function create_operator_masks(gene_seq_alpha::Vector{Int8}, gene_seq_beta::Vector{Int8}, pb::Real=0.2) - alpha_operator = zeros(Int8, length(gene_seq_alpha)) - beta_operator = zeros(Int8, length(gene_seq_beta)) - indices_alpha = rand(1:length(gene_seq_alpha), min(round(Int, (pb * length(gene_seq_alpha))), length(gene_seq_alpha))) - indices_beta = rand(1:length(gene_seq_beta), min(round(Int, (pb * length(gene_seq_beta))), length(gene_seq_beta))) - alpha_operator[indices_alpha] .= Int8(1) - beta_operator[indices_beta] .= Int8(1) - return alpha_operator, beta_operator + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(gene_seq_alpha) + + buffer.alpha_operator[1:len_a] .= zeros(Int8) + buffer.beta_operator[1:len_a] .= zeros(Int8) + + indices_alpha = rand(1:len_a, min(round(Int, (pb * len_a)), len_a)) + indices_beta = rand(1:len_a, min(round(Int, (pb * len_a)), len_a)) + + buffer.alpha_operator[indices_alpha] .= Int8(1) + buffer.beta_operator[indices_beta] .= Int8(1) end + @inline function create_operator_point_one_masks(gene_seq_alpha::Vector{Int8}, gene_seq_beta::Vector{Int8}, toolbox::Toolbox) - alpha_operator = zeros(Int8, length(gene_seq_alpha)) - beta_operator = zeros(Int8, length(gene_seq_beta)) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(gene_seq_alpha) + + buffer.alpha_operator[1:len_a] .= zeros(Int8) + buffer.beta_operator[1:len_a] .= zeros(Int8) + head_len = toolbox.head_len gene_len = head_len * 2 + 1 @@ -428,20 +536,21 @@ end point1 = rand(ref:mid) point2 = rand((mid+1):(ref+gene_len-1)) - alpha_operator[point1:point2] .= Int8(1) + buffer.alpha_operator[point1:point2] .= Int8(1) point1 = rand(ref:mid) point2 = rand((mid+1):(ref+gene_len-1)) - beta_operator[point1:point2] .= Int8(1) + buffer.beta_operator[point1:point2] .= Int8(1) end - - return alpha_operator, beta_operator end @inline function create_operator_point_two_masks(gene_seq_alpha::Vector{Int8}, gene_seq_beta::Vector{Int8}, toolbox::Toolbox) - alpha_operator = zeros(Int8, length(gene_seq_alpha)) - beta_operator = zeros(Int8, length(gene_seq_beta)) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(gene_seq_alpha) + + buffer.alpha_operator[1:len_a] .= zeros(Int8) + buffer.beta_operator[1:len_a] .= zeros(Int8) head_len = toolbox.head_len gene_len = head_len * 2 + 1 @@ -455,17 +564,17 @@ end point1 = rand(start:quarter) point2 = rand(quarter+1:half) point3 = rand(half+1:end_gene) - alpha_operator[point1:point2] .= Int8(1) - alpha_operator[point3:end_gene] .= Int8(1) + buffer.alpha_operator[point1:point2] .= Int8(1) + buffer.alpha_operator[point3:end_gene] .= Int8(1) point1 = rand(start:end_gene) point2 = rand(point1:end_gene) - beta_operator[point1:point2] .= Int8(1) - beta_operator[point2+1:end_gene] .= Int8(1) + buffer.beta_operator[point1:point2] .= Int8(1) + buffer.beta_operator[point2+1:end_gene] .= Int8(1) end - return alpha_operator, beta_operator + end @inline function replicate(chromosome1::Chromosome, chromosome2::Chromosome, toolbox) @@ -474,98 +583,83 @@ end @inline function gene_dominant_fusion!(chromosome1::Chromosome, chromosome2::Chromosome, pb::Real=0.2) - gene_seq_alpha = chromosome1.genes - gene_seq_beta = chromosome2.genes - alpha_operator, beta_operator = create_operator_masks(gene_seq_alpha, gene_seq_beta, pb) - - child_1_genes = similar(gene_seq_alpha) - child_2_genes = similar(gene_seq_beta) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(chromosome1.genes) + create_operator_masks(chromosome1.genes, chromosome2.genes, pb) - @inbounds @simd for i in eachindex(gene_seq_alpha) - child_1_genes[i] = alpha_operator[i] == 1 ? max(gene_seq_alpha[i], gene_seq_beta[i]) : gene_seq_alpha[i] - child_2_genes[i] = beta_operator[i] == 1 ? max(gene_seq_alpha[i], gene_seq_beta[i]) : gene_seq_beta[i] + @inbounds @simd for i in eachindex(chromosome1.genes) + buffer.child_1_genes[i] = buffer.alpha_operator[i] == 1 ? max(chromosome1.genes[i], chromosome2.genes[i]) : chromosome1.genes[i] + buffer.child_2_genes[i] = buffer.beta_operator[i] == 1 ? max(chromosome1.genes[i], chromosome2.genes[i]) : chromosome2.genes[i] end - chromosome1.genes = child_1_genes - chromosome2.genes = child_2_genes + chromosome1.genes .= @view buffer.child_1_genes[1:len_a] + chromosome2.genes .= @view buffer.child_2_genes[1:len_a] end @inline function gen_rezessiv!(chromosome1::Chromosome, chromosome2::Chromosome, pb::Real=0.2) - gene_seq_alpha = chromosome1.genes - gene_seq_beta = chromosome2.genes - alpha_operator, beta_operator = create_operator_masks(gene_seq_alpha, gene_seq_beta, pb) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(chromosome1.genes) + create_operator_masks(chromosome1.genes, chromosome2.genes, pb) - child_1_genes = similar(gene_seq_alpha) - child_2_genes = similar(gene_seq_beta) - - @inbounds @simd for i in eachindex(gene_seq_alpha) - child_1_genes[i] = alpha_operator[i] == 1 ? min(gene_seq_alpha[i], gene_seq_beta[i]) : gene_seq_alpha[i] - child_2_genes[i] = beta_operator[i] == 1 ? min(gene_seq_alpha[i], gene_seq_beta[i]) : gene_seq_beta[i] + @inbounds @simd for i in eachindex(chromosome1.genes) + buffer.child_1_genes[i] = buffer.alpha_operator[i] == 1 ? min(chromosome1.genes[i], chromosome2.genes[i]) : chromosome1.genes[i] + buffer.child_2_genes[i] = buffer.beta_operator[i] == 1 ? min(chromosome1.genes[i], chromosome2.genes[i]) : chromosome2.genes[i] end - chromosome1.genes = child_1_genes - chromosome2.genes = child_2_genes + chromosome1.genes .= @view buffer.child_1_genes[1:len_a] + chromosome2.genes .= @view buffer.child_2_genes[1:len_a] end @inline function gene_fussion!(chromosome1::Chromosome, chromosome2::Chromosome, pb::Real=0.2) - gene_seq_alpha = chromosome1.genes - gene_seq_beta = chromosome2.genes - alpha_operator, beta_operator = create_operator_masks(gene_seq_alpha, gene_seq_beta, pb) - - child_1_genes = similar(gene_seq_alpha) - child_2_genes = similar(gene_seq_beta) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(chromosome1.genes) + create_operator_masks(chromosome1.genes, chromosome2.genes, pb) - @inbounds @simd for i in eachindex(gene_seq_alpha) - child_1_genes[i] = alpha_operator[i] == 1 ? Int8((gene_seq_alpha[i] + gene_seq_beta[i]) ÷ 2) : gene_seq_alpha[i] - child_2_genes[i] = beta_operator[i] == 1 ? Int8((gene_seq_alpha[i] + gene_seq_beta[i]) ÷ 2) : gene_seq_beta[i] + @inbounds @simd for i in eachindex(chromosome1.genes) + buffer.child_1_genes[i] = buffer.alpha_operator[i] == 1 ? Int8((chromosome1.genes[i] + chromosome2.genes[i]) ÷ 2) : chromosome1.genes[i] + buffer.child_2_genes[i] = buffer.beta_operator[i] == 1 ? Int8((chromosome1.genes[i] + chromosome2.genes[i]) ÷ 2) : chromosome2.genes[i] end - chromosome1.genes = child_1_genes - chromosome2.genes = child_2_genes + chromosome1.genes .= @view buffer.child_1_genes[1:len_a] + chromosome2.genes .= @view buffer.child_2_genes[1:len_a] end @inline function gene_one_point_cross_over!(chromosome1::Chromosome, chromosome2::Chromosome) - gene_seq_alpha = chromosome1.genes - gene_seq_beta = chromosome2.genes - alpha_operator, beta_operator = create_operator_point_one_masks(gene_seq_alpha, gene_seq_beta, chromosome1.toolbox) - - child_1_genes = similar(gene_seq_alpha) - child_2_genes = similar(gene_seq_beta) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(chromosome1.genes) + create_operator_point_one_masks(chromosome1.genes, chromosome2.genes, chromosome1.toolbox) - @inbounds @simd for i in eachindex(gene_seq_alpha) - child_1_genes[i] = alpha_operator[i] == 1 ? gene_seq_alpha[i] : gene_seq_beta[i] - child_2_genes[i] = beta_operator[i] == 1 ? gene_seq_beta[i] : gene_seq_alpha[i] + @inbounds @simd for i in eachindex(chromosome1.genes) + buffer.child_1_genes[i] = buffer.alpha_operator[i] == 1 ? chromosome1.genes[i] : chromosome2.genes[i] + buffer.child_2_genes[i] = buffer.beta_operator[i] == 1 ? chromosome2.genes[i] : chromosome1.genes[i] end - chromosome1.genes = child_1_genes - chromosome2.genes = child_2_genes + chromosome1.genes .= @view buffer.child_1_genes[1:len_a] + chromosome2.genes .= @view buffer.child_2_genes[1:len_a] end @inline function gene_two_point_cross_over!(chromosome1::Chromosome, chromosome2::Chromosome) - gene_seq_alpha = chromosome1.genes - gene_seq_beta = chromosome2.genes - alpha_operator, beta_operator = create_operator_point_two_masks(gene_seq_alpha, gene_seq_beta, chromosome1.toolbox) + buffer = THREAD_BUFFERS[Threads.threadid()] + len_a = length(chromosome1.genes) + create_operator_point_two_masks(chromosome1.genes, chromosome2.genes, chromosome1.toolbox) - child_1_genes = similar(gene_seq_alpha) - child_2_genes = similar(gene_seq_beta) - - @inbounds @simd for i in eachindex(gene_seq_alpha) - child_1_genes[i] = alpha_operator[i] == 1 ? gene_seq_alpha[i] : gene_seq_beta[i] - child_2_genes[i] = beta_operator[i] == 1 ? gene_seq_beta[i] : gene_seq_alpha[i] + @inbounds @simd for i in eachindex(chromosome1.genes) + buffer.child_1_genes[i] = buffer.alpha_operator[i] == 1 ? chromosome1.genes[i] : chromosome2.genes[i] + buffer.child_2_genes[i] = buffer.beta_operator[i] == 1 ? chromosome2.genes[i] : chromosome1.genes[i] end - chromosome1.genes = child_1_genes - chromosome2.genes = child_2_genes + chromosome1.genes .= @view buffer.child_1_genes[1:len_a] + chromosome2.genes .= @view buffer.child_2_genes[1:len_a] end @inline function gene_mutation!(chromosome1::Chromosome, pb::Real=0.25) - gene_seq_alpha = chromosome1.genes - alpha_operator, _ = create_operator_masks(gene_seq_alpha, gene_seq_alpha, pb) - mutation_seq_1 = generate_chromosome(chromosome1.toolbox) + buffer = THREAD_BUFFERS[Threads.threadid()] + create_operator_masks(chromosome1.genes, chromosome1.genes, pb) + buffer.child_1_genes[1:length(chromosome1.genes)] .= generate_chromosome(chromosome1.toolbox).genes[1:length(chromosome1.genes)] - @inbounds @simd for i in eachindex(gene_seq_alpha) - gene_seq_alpha[i] = alpha_operator[i] == 1 ? mutation_seq_1.genes[i] : gene_seq_alpha[i] - end + @inbounds @simd for i in eachindex(chromosome1.genes) + chromosome1.genes[i] = buffer.alpha_operator[i] == 1 ? buffer.child_1_genes[i] : chromosome1.genes[i] + end end @inline function gene_inversion!(chromosome1::Chromosome) @@ -587,7 +681,7 @@ end end @inline function reverse_insertion_tail!(chromosome::Chromosome) - start_1 = rand(chromosome.toolbox.gen_start_indices)+chromosome.toolbox.head_len+1 + start_1 = rand(chromosome.toolbox.gen_start_indices) + chromosome.toolbox.head_len + 1 rolled_array = circshift(chromosome.genes[start_1:start_1+chromosome.toolbox.head_len-1], rand(1:chromosome.toolbox.head_len-1)) chromosome.genes[start_1:start_1+chromosome.toolbox.head_len-1] = rolled_array end diff --git a/src/GeneExpressionProgramming.jl b/src/GeneExpressionProgramming.jl index 87a6418..8b01a4c 100644 --- a/src/GeneExpressionProgramming.jl +++ b/src/GeneExpressionProgramming.jl @@ -38,7 +38,6 @@ providing tools for symbolic regression and evolutionary computation. ## Utility Functions - `find_indices_with_sum`, `compile_djl_datatype`: Data processing utilities - `optimize_constants!`: Constants optimization -- `minmax_scale`, `float16_scale`: Scaling functions - `isclose`: Numerical comparison - `save_state`, `load_state`: State persistence - `train_test_split`: Data splitting utility @@ -53,7 +52,7 @@ providing tools for symbolic regression and evolutionary computation. ## History Recording - `HistoryRecorder`, `OptimizationHistory`: History tracking structures -- `create_history_recorder`, `record_history!`, `record!`: Recording functions +- `record_history!`, `record!`: Recording functions - `close_recorder!`, `get_history_arrays`: History management # Example Usage @@ -88,51 +87,137 @@ include("Selection.jl") include("Util.jl") include("RegressionWrapper.jl") -export GepEntities, LossFunction, PhysicalConstants, GepUtils, GepRegression, RegressionWrapper - -using .GepRegression -export runGep - -using .LossFunction +# First export the submodules themselves +export GepEntities, LossFunction, PhysicalConstants, + GepUtils, GepRegression, RegressionWrapper, GPSurrogate + +# Import GEP core functionality +import .GepRegression: + runGep + +# Import loss functions +import .LossFunction: + get_loss_function + +# Import utilities +import .GepUtils: + find_indices_with_sum, + compile_djl_datatype, + optimize_constants!, + minmax_scale, + isclose, + save_state, + load_state, + record_history!, + record!, + close_recorder!, + HistoryRecorder, + OptimizationHistory, + get_history_arrays, + train_test_split, + ARITY_LIB_COMMON, + FUNCTION_LIB_COMMON + +# Import selection mechanisms +import .EvoSelection: + tournament_selection, + nsga_selection, + dominates_, + fast_non_dominated_sort, + calculate_fronts, + determine_ranks, + assign_crowding_distance + +# Import physical constants functionality +import .PhysicalConstants: + physical_constants, + physical_constants_all, + get_constant, + get_constant_value, + get_constant_dims + +# Import symbolic computation utilities +import .SBPUtils: + TokenLib, + TokenDto, + LibEntry, + TempComputeTree, + create_lib, + create_compute_tree, + propagate_necessary_changes!, + calculate_vector_dimension!, + flush!, + flatten_dependents, + correct_genes!, + equal_unit_forward, + mul_unit_forward, + div_unit_forward, + zero_unit_backward, + zero_unit_forward, + sqr_unit_backward, + sqr_unit_forward, + mul_unit_backward, + div_unit_backward, + equal_unit_backward, + get_feature_dims_json, + get_target_dim_json, + retrieve_coeffs_based_on_similarity + +# Import core GEP entities +import .GepEntities: + Chromosome, + Toolbox, + EvaluationStrategy, + StandardRegressionStrategy, + GenericRegressionStrategy, + fitness, + set_fitness!, + generate_gene, + compile_expression!, + generate_chromosome, + generate_population, + genetic_operations! + +# Import regression wrapper functionality +import .RegressionWrapper: + GepRegressor, + GepTensorRegressor, + fit!, + list_all_functions, + list_all_arity, + list_all_forward_handlers, + list_all_backward_handlers, + list_all_genetic_params, + set_function!, + set_arity!, + set_forward_handler!, + set_backward_handler!, + update_function! + + + +export runGep, EvaluationStrategy, StandardRegressionStrategy, GenericRegressionStrategy export get_loss_function - - -using .GepUtils -export find_indices_with_sum, compile_djl_datatype, optimize_constants!, minmax_scale, float16_scale, isclose -export save_state, load_state, create_history_recorder, record_history!, record!, close_recorder! +export find_indices_with_sum, compile_djl_datatype, optimize_constants!, minmax_scale, isclose +export save_state, load_state, record_history!, record!, close_recorder! export HistoryRecorder, OptimizationHistory, get_history_arrays export train_test_split - -using .EvoSelection -export selection_NSGA, basic_tournament_selection, dominates_, fast_non_dominated_sort, calculate_fronts, determine_ranks, assign_crowding_distance - - -using .PhysicalConstants -export physical_constants, physical_constants_all -export get_constant, get_constant_value, get_constant_dims - - -using .SBPUtils +export tournament_selection, nsga_selection, dominates_, fast_non_dominated_sort, calculate_fronts, determine_ranks, assign_crowding_distance +export physical_constants, physical_constants_all, get_constant, get_constant_value, get_constant_dims export TokenLib, TokenDto, LibEntry, TempComputeTree -export create_lib, create_compute_tree, propagate_necessary_changes!, calculate_vector_dimension!, flush!, calculate_vector_dimension!, flatten_dependents -export propagate_necessary_changes!, correct_genes! -export equal_unit_forward, mul_unit_forward, div_unit_forward, zero_unit_backward, zero_unit_forward, sqr_unit_backward, sqr_unit_forward, mul_unit_backward, div_unit_backward, equal_unit_backward +export create_lib, create_compute_tree, propagate_necessary_changes!, calculate_vector_dimension!, flush!, flatten_dependents +export correct_genes!, equal_unit_forward, mul_unit_forward, div_unit_forward +export zero_unit_backward, zero_unit_forward, sqr_unit_backward, sqr_unit_forward, mul_unit_backward, div_unit_backward, equal_unit_backward export get_feature_dims_json, get_target_dim_json, retrieve_coeffs_based_on_similarity - - -using .GepEntities -export Chromosome, Toolbox -export AbstractSymbol, FunctionalSymbol, BasicSymbol, SymbolConfig -export fitness, set_fitness! -export generate_gene, generate_preamle!, compile_expression!, generate_chromosome, generate_population +export Chromosome, Toolbox, fitness, set_fitness! +export generate_gene, compile_expression!, generate_chromosome, generate_population export genetic_operations! +export GepRegressor,GepTensorRegressor, fit! +export list_all_functions, list_all_arity, list_all_forward_handlers +export list_all_backward_handlers, list_all_genetic_params +export set_function!, set_arity!, set_forward_handler!, set_backward_handler! +export update_function! +export ARITY_LIB_COMMON, FUNCTION_LIB_COMMON -using .RegressionWrapper -export GepRegressor, fit! -export list_all_functions, list_all_arity, list_all_forward_handlers, - list_all_backward_handlers, list_all_genetic_params, - set_function!, set_arity!, set_forward_handler!, set_backward_handler!, - update_function! - end diff --git a/src/Gep.jl b/src/Gep.jl index 26fd865..198fb14 100644 --- a/src/Gep.jl +++ b/src/Gep.jl @@ -57,7 +57,7 @@ See also: - [`GepEntities.fitness`](@ref): Fitness value access - [`GepEntities.set_fitness!`](@ref): Fitness value modification - +#TODO need to create different strategies for wrapping costum functions """ module GepRegression @@ -66,14 +66,13 @@ module GepRegression include("Losses.jl") include("Util.jl") include("Selection.jl") -include("Entities.jl") -using .LossFunction -export runGep +using .LossFunction using .GepUtils using .EvoSelection -using .GepEntities +using ..GepEntities + using Random using Statistics @@ -83,17 +82,21 @@ using OrderedCollections using DynamicExpressions using Logging using Printf +using Base.Threads: SpinLock +using .Threads +export runGep -const Chromosome = GepEntities.Chromosome -const Toolbox = GepEntities.Toolbox - +#const Chromosome = GepEntities.Chromosome +#const Toolbox = GepEntities.Toolbox +#const EvaluationStrategy = GepEntities.EvaluationStrategy +#redesign -> compute fitness should return fitness and crash, we just need to insert the chromosome """ compute_fitness(elem::Chromosome, operators::OperatorEnum, x_data::AbstractArray{T}, y_data::AbstractArray{T}, loss_function::Function, crash_value::T; - validate::Bool=false, penalty_consideration::Real=0.0) where {T<:AbstractFloat} + validate::Bool=false) where {T<:AbstractFloat} Computes the fitness score for a chromosome using the specified loss function. @@ -105,32 +108,33 @@ Computes the fitness score for a chromosome using the specified loss function. - `loss_function::Function`: The loss function used to compute fitness - `crash_value::T`: Default value returned if computation fails - `validate::Bool=false`: If true, forces recomputation of fitness even if already calculated -- `penalty_consideration::Real=0.0`: Additional penalty term added to the fitness score # Returns -Returns the computed fitness value (loss + penalty) or crash_value if computation fails +Returns the computed fitness value (loss) or crash_value if computation fails # Details - Checks if fitness needs to be computed (if NaN or validate=true) - Evaluates the chromosome's compiled function on input data -- Applies loss function and adds any penalty consideration - Returns crash_value if any errors occur during computation """ -@inline function compute_fitness(elem::Chromosome, operators::OperatorEnum, x_data::AbstractArray{T}, - y_data::AbstractArray{T}, loss_function::Function, - crash_value::T; validate::Bool=false, penalty_consideration::Real=0.0) where {T<:AbstractFloat} +@inline function compute_fitness(elem::Chromosome, evalArgs::StandardRegressionStrategy; validate::Bool=false) try - if isnan(elem.fitness) || validate - y_pred = elem.compiled_function(x_data, operators) - return loss_function(y_data, y_pred) + penalty_consideration + if isnan(mean(elem.fitness)) || validate + y_pred = elem.compiled_function(evalArgs.x_data, evalArgs.operators) + return (evalArgs.loss_function(evalArgs.y_data, y_pred),) else - return elem.fitness + return (elem.fitness,) end catch e - return crash_value + return (evalArgs.crash_value,) end end +@inline function compute_fitness(elem::Chromosome, evalArgs::GenericRegressionStrategy; validate::Bool=false) + return evalArgs.loss_function(elem, validate) +end + + """ perform_step!(population::Vector{Chromosome}, parents::Vector{Chromosome}, next_gen::Vector{Chromosome}, toolbox::Toolbox, mating_size::Int) @@ -176,6 +180,10 @@ Performs one evolutionary step in the GEP algorithm, creating and evaluating new end +@inline function update_surrogate!(::EvaluationStrategy) + nothing +end + """ perform_correction_callback!(population::Vector{Chromosome}, epoch::Int, correction_epochs::Int, correction_amount::Real, @@ -201,172 +209,155 @@ Applies correction operations to ensure dimensional homogeneity in chromosomes. if !isnothing(correction_callback) && epoch % correction_epochs == 0 pop_amount = Int(ceil(length(population) * correction_amount)) - Threads.@threads for i in 1:pop_amount + Threads.@threads for i in 1:pop_amount if !(population[i].dimension_homogene) && population[i].compiled distance, correction = correction_callback(population[i].genes, population[i].toolbox.gen_start_indices, population[i].expression_raw) if correction compile_expression!(population[i]; force_compile=true) - population[i].dimension_homogene = true + population[i].dimension_homogene = true else - population[i].fitness += distance + #population[i].fitness += distance end end end end end + """ - runGep(epochs::Int, population_size::Int, operators::OperatorEnum, - x_data::AbstractArray{T}, y_data::AbstractArray{T}, toolbox::Toolbox; - hof::Int=3, x_data_test::Union{AbstractArray{T},Nothing}=nothing, - y_data_test::Union{AbstractArray{T},Nothing}=nothing, - loss_fun_::Union{String,Function}="mae", - correction_callback::Union{Function,Nothing}=nothing, + runGep(epochs::Int, population_size::Int, toolbox::Toolbox, evalStrategy::EvaluationStrategy; + hof::Int=3, correction_callback::Union{Function,Nothing}=nothing, correction_epochs::Int=1, correction_amount::Real=0.6, - tourni_size::Int=3, opt_method_const::Symbol=:cg, - optimisation_epochs::Int=500) where {T<:AbstractFloat} + tourni_size::Int=3) -Main function that executes the GEP algorithm for regression problems. +Main function that executes the GEP algorithm using a specified evaluation strategy. # Arguments - `epochs::Int`: Number of evolutionary epochs to run - `population_size::Int`: Size of the chromosome population -- `operators::OperatorEnum`: Mathematical operators for the DynamicExpressions -- `x_data::AbstractArray{T}`: Training input features -- `y_data::AbstractArray{T}`: Training target values - `toolbox::Toolbox`: Contains genetic operators and algorithm parameters +- `evalStrategy::EvaluationStrategy`: Strategy for evaluating chromosomes, handling fitness computation, and optimization + +# Optional Arguments - `hof::Int=3`: Number of best solutions to return (Hall of Fame size) -- `x_data_test::Union{AbstractArray{T},Nothing}`: Optional test input features -- `y_data_test::Union{AbstractArray{T},Nothing}`: Optional test target values -- `loss_fun_::Union{String,Function}="mae"`: Loss function for fitness computation -- `correction_callback::Union{Function,Nothing}`: Function for dimensional homogeneity correction +- `correction_callback::Union{Function,Nothing}=nothing`: Function for dimensional homogeneity correction - `correction_epochs::Int=1`: Frequency of correction operations - `correction_amount::Real=0.6`: Proportion of population for correction - `tourni_size::Int=3`: Tournament selection size -- `opt_method_const::Symbol=:cg`: Optimization method for constant optimization -- `optimisation_epochs::Int=500`: Frequency of constant optimization # Returns -Tuple{Vector{Chromosome}, Any}: Returns best solutions and training history -Best solutions include both training and test R² scores when test data is provided +`Tuple{Vector{Chromosome}, Any}`: Returns best solutions and training history # Details 1. Initializes population and evolution parameters 2. For each epoch: - - Applies dimensional homogeneity corrections - - Computes fitness for all chromosomes - - Optimizes constants for best solution periodically + - Applies dimensional homogeneity corrections if provided + - Computes fitness for all chromosomes using evaluation strategy + - Sorts population based on fitness + - Applies secondary optimization if specified in strategy - Records training progress - - Performs tournament selection - - Creates new generation through genetic operations -3. Computes final R² scores for best solutions -4. Returns best solutions and training history + - Checks break condition from evaluation strategy + - Performs selection and creates new generation +3. Returns best solutions and training history Progress is monitored through a progress bar showing: - Current epoch - Training loss - Validation loss -Early stopping occurs if perfect R² score is achieved + +The evolution process stops when either: +- Maximum epochs is reached +- Break condition specified in evaluation strategy is met => needs to be informed as break_condition(population, epoch) """ -function runGep(epochs::Int, +@inline function runGep(epochs::Int, population_size::Int, - operators::OperatorEnum, - x_data::AbstractArray{T}, - y_data::AbstractArray{T}, - toolbox::Toolbox; + toolbox::Toolbox, + evalStrategy::EvaluationStrategy; hof::Int=3, - x_data_test::Union{AbstractArray{T},Nothing}=nothing, - y_data_test::Union{AbstractArray{T},Nothing}=nothing, - loss_fun_::Union{String,Function}="mae", correction_callback::Union{Function,Nothing}=nothing, correction_epochs::Int=1, correction_amount::Real=0.6, tourni_size::Int=3, - opt_method_const::Symbol=:cg, - optimisation_epochs::Int=500) where {T<:AbstractFloat} - - loss_fun = typeof(loss_fun_) == String ? get_loss_function(loss_fun_) : loss_fun_ - - recorder = HistoryRecorder(epochs, Float64) - - function optimizer_function(sub_tree::Node) - y_pred, flag = eval_tree_array(sub_tree, x_data, operators) - return get_loss_function("mse")(y_pred, y_data) - end + optimization_epochs::Int=500) - penalty_consideration = convert(Real,toolbox.gep_probs["penalty_consideration"]) + recorder = HistoryRecorder(epochs, Tuple) mating_ = toolbox.gep_probs["mating_size"] mating_size = Int(ceil(population_size * mating_)) mating_size = mating_size % 2 == 0 ? mating_size : mating_size - 1 - fits_representation = Vector{T}(undef, population_size) + fits_representation = Vector{Tuple}(undef, population_size) + fit_cache = Dict{Vector{Int8},Tuple}() + cache_lock = SpinLock() population = generate_population(population_size, toolbox) next_gen = Vector{eltype(population)}(undef, mating_size) progBar = Progress(epochs; showspeed=true, desc="Training: ") - - prev_best = -1 - + prev_best = (typemax(Float64),) for epoch in 1:epochs + same = Atomic{Int}(0) perform_correction_callback!(population, epoch, correction_epochs, correction_amount, correction_callback) + Threads.@threads for i in eachindex(population) - if isnan(population[i].fitness) - population[i].fitness = compute_fitness(population[i], operators, x_data, y_data, loss_fun, typemax(T); - penalty_consideration=population[i].penalty * penalty_consideration) + if isnan(mean(population[i].fitness)) + cache_value = get(fit_cache, population[i].expression_raw, nothing) + if isnothing(cache_value) + + population[i].fitness = compute_fitness(population[i], evalStrategy) + lock(cache_lock) + fit_cache[population[i].expression_raw] = population[i].fitness + unlock(cache_lock) + else + atomic_add!(same, 1) + population[i].fitness = cache_value + end end end - - sort!(population, by=x -> x.fitness) - - try - if (prev_best == -1 || prev_best > population[1].fitness) && epoch % optimisation_epochs == 0 - eqn, result = optimize_constants!(population[1].compiled_function, optimizer_function; - opt_method=opt_method_const, max_iterations=150, n_restarts=5) - population[1].fitness = result - population[1].compiled_function = eqn - prev_best = result - end - catch - @show "Ignored constant opt." + sort!(population, by=x -> mean(x.fitness)) + Threads.@threads for index in eachindex(population) + fits_representation[index] = population[index].fitness end - Threads.@threads for index in eachindex(population) - fits_representation[index] = population[index].fitness + if !isnothing(evalStrategy.secOptimizer) && epoch % optimization_epochs == 0 && population[1].fitness < prev_best + evalStrategy.secOptimizer(population) + fits_representation[1] = population[1].fitness + prev_best = fits_representation[1] end - val_loss = compute_fitness(population[1], operators, x_data_test, y_data_test, loss_fun, typemax(T); validate=true) + val_loss = compute_fitness(population[1], evalStrategy; validate=true) record!(recorder, epoch, fits_representation[1], val_loss, fits_representation) ProgressMeter.update!(progBar, epoch, showvalues=[ (:epoch_, @sprintf("%.0f", epoch)), - (:train_loss, @sprintf("%.6e", fits_representation[1])), - (:validation_loss, @sprintf("%.6e", val_loss)) + (:duplicates_per_epoch, @sprintf("%.0f", same[])), + (:train_loss, @sprintf("%.6e", mean(fits_representation[1]))), + (:validation_loss, @sprintf("%.6e", mean(val_loss))) ]) + update_surrogate!(evalStrategy) + - if isclose(fits_representation[1], zero(T)) + if !isnothing(evalStrategy.break_condition) && evalStrategy.break_condition(population, epoch) break end if epoch < epochs - indices = basic_tournament_selection(fits_representation[1:mating_size], tourni_size, mating_size) - parents = population[indices] + if length(fits_representation[1]) == 1 + selectedMembers = tournament_selection(fits_representation, mating_size, tourni_size) + else + selectedMembers = nsga_selection(fits_representation) + end + + parents = population[selectedMembers.indices] perform_step!(population, parents, next_gen, toolbox, mating_size) end end best = population[1:hof] - for elem in best - elem.fitness_r2_train = compute_fitness(elem, operators, x_data, y_data, get_loss_function("r2_score"), zero(T); validate=true) - if !isnothing(x_data_test) - elem.fitness_r2_test = compute_fitness(elem, operators, x_data_test, y_data_test, get_loss_function("r2_score"), zero(T); validate=true) - end - end - close_recorder!(recorder) return best, recorder.history end + end diff --git a/src/Losses.jl b/src/Losses.jl index 6e9f20e..3a8da0d 100644 --- a/src/Losses.jl +++ b/src/Losses.jl @@ -72,6 +72,7 @@ module LossFunction export get_loss_function using Statistics +using LoopVectorization function floor_to_n10p(x::T) where T<:AbstractFloat abs_x = abs(x) @@ -159,15 +160,48 @@ function r2_score_floor(y_true::AbstractArray{T}, y_pred::AbstractArray{T}) wher return r2_score(y_true_scaled, y_pred_scaled) end -function mean_squared_error(y_true::AbstractArray{T}, y_pred::AbstractArray{T}) where T<:AbstractFloat - d::T = zero(T) - @assert length(y_true) == length(y_pred) - @fastmath @inbounds @simd for i in eachindex(y_true, y_pred) - temp = (y_true[i]-y_pred[i]) - d += temp*temp + +@inline function mean_squared_error_(y_true::AbstractArray{T}, + y_pred::AbstractArray{T}) where {T<:AbstractFloat} + sum = zero(T) + len = length(y_true) + + @fastmath @turbo for i in eachindex(y_true, y_pred) + diff = y_true[i] - y_pred[i] + sum += diff * diff + end + + return sum / len +end + + +@inline function mean_squared_error(y_true::AbstractArray{T}, + y_pred::AbstractArray{T}) where {T<:AbstractFloat} + len = length(y_true) + if len < 100_000 + return mean_squared_error_(y_true, y_pred) + end + + n_chunks = Threads.nthreads() + chunk_size = div(len, n_chunks) + + + partial_sums = zeros(T, n_chunks) + + Threads.@threads for chunk in 1:n_chunks + start_idx = (chunk - 1) * chunk_size + 1 + end_idx = chunk == n_chunks ? len : chunk * chunk_size + + sum = zero(T) + @fastmath @turbo for i in start_idx:end_idx + diff = y_true[i] - y_pred[i] + sum += diff * diff end - return d/length(y_true) -end + partial_sums[chunk] = sum + end + + return sum(partial_sums) / len +end function root_mean_squared_error(y_true::AbstractArray{T}, y_pred::AbstractArray{T}) where T<:AbstractFloat d::T = zero(T) diff --git a/src/RegressionWrapper.jl b/src/RegressionWrapper.jl index c31567b..219b0be 100644 --- a/src/RegressionWrapper.jl +++ b/src/RegressionWrapper.jl @@ -27,17 +27,6 @@ simplified interfaces and utilities for symbolic regression tasks with physical - `GENE_COMMON_PROBS`: Default genetic operation probabilities - `FUNCTION_LIB_BACKWARD_COMMON`: Backward dimension propagation functions - `FUNCTION_LIB_FORWARD_COMMON`: Forward dimension propagation functions -- `FUNCTION_LIB_COMMON`: Available mathematical functions - -# Function Library -Includes extensive mathematical operations: -- Basic arithmetic: +, -, *, /, ^ -- Comparisons: min, max -- Rounding: floor, ceil, round -- Exponential: exp, log, log10, log2 -- Trigonometric: sin, cos, tan, asin, acos, atan -- Hyperbolic: sinh, cosh, tanh, asinh, acosh, atanh -- Special: sqr, sqrt, sign, abs # Usage Example ```julia @@ -132,15 +121,15 @@ See also: module RegressionWrapper -export GepRegressor +export GepRegressor, GepTensorRegressor export create_function_entries, create_feature_entries, create_constants_entries, create_physical_operations -export GENE_COMMON_PROBS, FUNCTION_LIB_BACKWARD_COMMON, FUNCTION_LIB_FORWARD_COMMON, FUNCTION_LIB_COMMON +export GENE_COMMON_PROBS, FUNCTION_LIB_BACKWARD_COMMON, FUNCTION_LIB_FORWARD_COMMON export fit! -export list_all_functions, list_all_arity, list_all_forward_handlers, - list_all_backward_handlers, list_all_genetic_params, - set_function!, set_arity!, set_forward_handler!, set_backward_handler!, - update_function! +export list_all_functions, list_all_arity, list_all_forward_handlers, + list_all_backward_handlers, list_all_genetic_params, + set_function!, set_arity!, set_forward_handler!, set_backward_handler!, + update_function!, vec_add, vec_mul include("Entities.jl") include("Gep.jl") @@ -149,115 +138,22 @@ include("PhyConstants.jl") include("Sbp.jl") include("Selection.jl") include("Util.jl") +include("TensorOps.jl") using .GepEntities using .LossFunction -using .GepEntities using .EvoSelection using .GepRegression using .SBPUtils using .GepUtils +using .TensorRegUtils using DynamicExpressions using OrderedCollections +using LinearAlgebra -const Toolbox = GepRegression.GepEntities.Toolbox -const TokenDto = SBPUtils.TokenDto - -function sqr(x::Vector{T}) where {T<:AbstractFloat} - return x .* x -end - -function sqr(x::T) where {T<:Union{AbstractFloat,Node{<:AbstractFloat}}} - return x * x -end - - -""" - FUNCTION_LIB_COMMON::Dict{Symbol,Function} - -Dictionary mapping function symbols to their corresponding functions. -Contains basic mathematical operations, trigonometric, and other common functions. - -# Available Functions -- Basic arithmetic: `+`, `-`, `*`, `/`, `^` -- Comparison: `min`, `max` -- Rounding: `floor`, `ceil`, `round` -- Exponential & Logarithmic: `exp`, `log`, `log10`, `log2` -- Trigonometric: `sin`, `cos`, `tan`, `asin`, `acos`, `atan` -- Hyperbolic: `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh` -- Other: `abs`, `sqr`, `sqrt`, `sign` - -To add a new function, ensure you also add corresponding entries in `ARITY_LIB_COMMON`, -`FUNCTION_LIB_FORWARD_COMMON`, and `FUNCTION_LIB_BACKWARD_COMMON`. -""" -const FUNCTION_LIB_COMMON = Dict{Symbol,Function}( - :+ => +, - :- => -, - :* => *, - :/ => /, - :^ => ^, - :min => min, - :max => max, :abs => abs, - :floor => floor, - :ceil => ceil, - :round => round, :exp => exp, - :log => log, - :log10 => log10, - :log2 => log2, :sin => sin, - :cos => cos, - :tan => tan, - :asin => asin, - :acos => acos, - :atan => atan, :sinh => sinh, - :cosh => cosh, - :tanh => tanh, - :asinh => asinh, - :acosh => acosh, - :atanh => atanh, :sqr => sqr, - :sqrt => sqrt, :sign => sign -) - -""" - ARITY_LIB_COMMON::Dict{Symbol,Int8} - -Dictionary specifying the number of arguments (arity) for each function in the library. -- Value of 1 indicates unary functions (e.g., `sin`, `cos`, `abs`) -- Value of 2 indicates binary functions (e.g., `+`, `-`, `*`, `/`) - -When adding new functions to `FUNCTION_LIB_COMMON`, ensure to specify their arity here. -""" -const ARITY_LIB_COMMON = Dict{Symbol,Int8}( - :+ => 2, - :- => 2, - :* => 2, - :/ => 2, - :^ => 2, - :min => 2, - :max => 2, :abs => 1, - :floor => 1, - :ceil => 1, - :round => 1, - :exp => 1, - :log => 1, - :log10 => 1, - :log2 => 1, - :sin => 1, - :cos => 1, - :tan => 1, - :asin => 1, - :acos => 1, - :atan => 1, - :sinh => 1, - :cosh => 1, - :tanh => 1, - :asinh => 1, - :acosh => 1, - :atanh => 1, - :sqrt => 1, - :sqr => 1 -) +const InputSelector = TensorRegUtils.InputSelector """ FUNCTION_LIB_FORWARD_COMMON::Dict{Symbol,Function} @@ -363,7 +259,6 @@ Dictionary containing default probabilities and parameters for genetic algorithm - `fusion_rate`: Rate of general fusion (0.0) - `inversion_prob`: Probability of inversion (0.1) - `mating_size`: Relative size of mating pool (0.5) -- `penalty_consideration`: Weight of penalty in fitness evaluation (0.2) These values can be adjusted to fine-tune the genetic algorithm's behavior. """ @@ -381,8 +276,7 @@ const GENE_COMMON_PROBS = Dict{String,AbstractFloat}( "inversion_prob" => 0.1, "reverse_insertion" => 0.05, "reverse_insertion_tail" => 0.05, - "mating_size" => 0.5, - "penalty_consideration" => 0.2) + "mating_size" => 0.5) const SymbolDict = OrderedDict{Int8,Int8} const CallbackDict = Dict{Int8,Function} @@ -487,7 +381,7 @@ function create_constants_entries( for elem in entered_terminal_nums utilized_symbols[cur_idx] = 0 - nodes[cur_idx] = parse(node_type, string(elem)) + nodes[cur_idx] = Node{node_type}(; val=parse(node_type, string(elem))) dimension_information[cur_idx] = get(dimensions_to_consider, elem, ZERO_DIM) cur_idx += 1 end @@ -495,7 +389,7 @@ function create_constants_entries( for _ in 1:rnd_count utilized_symbols[cur_idx] = 0 - nodes[cur_idx] = rand() + nodes[cur_idx] = Node{node_type}(; val=rand()) dimension_information[cur_idx] = ZERO_DIM cur_idx += 1 end @@ -504,6 +398,7 @@ function create_constants_entries( end +#preamble syms are just dummies function create_preamble_entries( preamble_syms_raw::Vector{Symbol}, dimensions_to_consider::Dict{Symbol,Vector{Float16}}, @@ -519,7 +414,7 @@ function create_preamble_entries( for elem in preamble_syms_raw utilized_symbols[cur_idx] = 0 - nodes[cur_idx] = Node{AbstractArray}(feature=cur_idx) + nodes[cur_idx] = Node{node_type}(feature=cur_idx) dimension_information[cur_idx] = get(dimensions_to_consider, elem, ZERO_DIM) push!(preamble_syms, cur_idx) cur_idx += 1 @@ -563,12 +458,13 @@ Create a Gene Expression Programming regressor for symbolic regression. - `preamble_syms::Vector{Symbol}=Symbol[]`: Preamble symbols - `max_permutations_lib::Int=10000`: Maximum permutations for dimension library - `rounds::Int=4`: Rounds for dimension library creation +- `number_of_objectives::Int=1`: Defines the number of objectives considered by the search """ mutable struct GepRegressor toolbox_::Toolbox operators_::OperatorEnum dimension_information_::OrderedDict{Int8,Vector{Float16}} - best_models_::Union{Nothing,Vector{GepRegression.GepEntities.Chromosome}} + best_models_::Union{Nothing,Vector{Chromosome}} fitness_history_::Any token_dto_::Union{TokenDto,Nothing} @@ -581,10 +477,11 @@ mutable struct GepRegressor considered_dimensions::Dict{Symbol,Vector{Float16}}=Dict{Symbol,Vector{Float16}}(), rnd_count::Int=1, node_type::Type=Float64, - gene_count::Int=2, - head_len::Int=10, + gene_count::Int=3, + head_len::Int=6, preamble_syms::Vector{Symbol}=Symbol[], - max_permutations_lib::Int=10000, rounds::Int=4 + max_permutations_lib::Int=10000, rounds::Int=4, + number_of_objectives::Int=1 ) entered_features_ = isempty(entered_features) ? @@ -637,7 +534,8 @@ mutable struct GepRegressor end toolbox = GepRegression.GepEntities.Toolbox(gene_count, head_len, utilized_symbols, gene_connections_, - callbacks, nodes, GENE_COMMON_PROBS; preamble_syms=preamble_syms_) + callbacks, nodes, GENE_COMMON_PROBS; preamble_syms=preamble_syms_, number_of_objectives=number_of_objectives, + operators_=operators) obj = new() obj.toolbox_ = toolbox @@ -648,6 +546,115 @@ mutable struct GepRegressor end end +""" + GepTensorRegressor + +A Gene Expression Programming (GEP) regressor that evolves higher order mathematical expressions (e.g tensor-based). + +# Fields +- `toolbox_::Toolbox`: Contains configuration and operators for GEP evolution +- `best_models_::Union{Nothing,Vector{Chromosome}}`: Best models found during evolution +- `fitness_history_::Any`: History of fitness values during training + +# Constructor + GepTensorRegressor(feature_amount::Int; kwargs...) + +Create a new GEP tensor regressor with specified number of input features. + +# Arguments +- `scalar_feature_amount::Int`: Number of input features representing scalar quantities +- `higher_dim_feature_amount::Int`: Number of input features representing hihger quantities + +# Keyword Arguments +- `entered_non_terminals::Vector{Symbol}=[:+, :-, :*, :/]`: Available mathematical operators +- `entered_terminal_nums::Vector{<:AbstractFloat}=[0.0, 0.5]`: Constants available as terminals +- `gene_connections::Vector{Symbol}=[:+, :-, :*, :/]`: Operators for connecting genes +- `rnd_count::Int=1`: Number of random constant terminals to generate +- `gene_count::Int=3`: Number of genes in each chromosome +- `head_len::Int=6`: Length of the head section in each gene +- `number_of_objectives::Int=1`: Number of optimization objectives + +The regressor uses GEP to evolve tensor-based mathematical expressions that map input features +to target values. It supports multiple genes connected by operators and can optimize for +multiple objectives. + +# Implementation Details +- Uses InputSelector nodes for features +- Combines fixed and random constant terminals +- Maps operators to TENSOR_NODES callbacks +- Uses TENSOR_NODES_ARITY for operator arity +- Compiles expressions to Flux networks via compile_to_flux_network +""" +#TODO => adapt probs for occurance of ! +mutable struct GepTensorRegressor + toolbox_::Toolbox + best_models_::Union{Nothing,Vector{Chromosome}} + fitness_history_::Any + + + function GepTensorRegressor(scalar_feature_amount::Int; + higher_dim_feature_amount::Int=0, + entered_non_terminals::Vector{Symbol}=[:+, :-, :*], + entered_terminal_nums::Vector{<:AbstractFloat}=Float64[], + gene_connections::Vector{Symbol}=[:+, :*], + rnd_count::Int=0, + gene_count::Int=2, + head_len::Int=3, + number_of_objectives::Int=1 + ) + #Creating the feature Nodes -> asuming a data dict pointing to + cur_idx = Int8(1) + nodes = OrderedDict{Int8,Any}() + utilized_symbols = SymbolDict() + callbacks = Dict{Int8,Any}() + gene_connections_ = Int8[] + tensor_syms_idx = Int8[] + tensor_function_idx = Int8[] + for _ in 1:scalar_feature_amount + nodes[cur_idx] = InputSelector(cur_idx) + utilized_symbols[cur_idx] = Int8(0) + cur_idx += 1 + end + + for _ in 1:higher_dim_feature_amount + nodes[cur_idx] = InputSelector(cur_idx) + utilized_symbols[cur_idx] = Int8(0) + cur_idx += 1 + end + + #Creating the const_nodes + for elem in entered_terminal_nums + nodes[cur_idx] = elem + utilized_symbols[cur_idx] = Int8(0) + cur_idx += 1 + end + + for _ in 1:rnd_count + nodes[cur_idx] = rand() + utilized_symbols[cur_idx] = Int8(0) + cur_idx += 1 + end + + #callback - index => function + for elem in entered_non_terminals + callbacks[cur_idx] = TENSOR_NODES[elem] + utilized_symbols[cur_idx] = TENSOR_NODES_ARITY[elem] + if elem in gene_connections + push!(gene_connections_,cur_idx) + end + cur_idx += 1 + end + + toolbox = GepRegression.GepEntities.Toolbox(gene_count, head_len, utilized_symbols, gene_connections_, + callbacks, nodes, GENE_COMMON_PROBS; number_of_objectives=number_of_objectives, + operators_=nothing, function_complile=compile_to_flux_network) + + obj = new() + obj.toolbox_ = toolbox + return obj + end +end + """ fit!(regressor::GepRegressor, epochs::Int, population_size::Int, x_train::AbstractArray, @@ -674,14 +681,15 @@ Train the GEP regressor model. - `opt_method_const::Symbol=:cg`: Optimization method for constants - `target_dimension::Union{Vector{Float16},Nothing}=nothing`: Target physical dimension """ -function fit!(regressor::GepRegressor, epochs::Int, population_size, x_train::AbstractArray, +function fit!(regressor::GepRegressor, epochs::Int, population_size::Int, x_train::AbstractArray, y_train::AbstractArray; x_test::Union{AbstractArray,Nothing}=nothing, y_test::Union{AbstractArray,Nothing}=nothing, - optimization_epochs::Int=500, + optimization_epochs::Int=100, hof::Int=3, loss_fun::Union{String,Function}="mse", correction_epochs::Int=1, correction_amount::Real=0.05, - tourni_size::Int=3, opt_method_const::Symbol=:cg, + opt_method_const::Symbol=:cg, target_dimension::Union{Vector{Float16},Nothing}=nothing, - cycles::Int=10 + cycles::Int=10, max_iterations::Int=150, n_starts::Int=5, + break_condition::Union{Function,Nothing}=nothing ) correction_callback = if !isnothing(target_dimension) @@ -697,28 +705,134 @@ function fit!(regressor::GepRegressor, epochs::Int, population_size, x_train::Ab nothing end - - best, history = runGep(epochs, - population_size, + @inline function optimizer_function_(sub_tree::Node) + y_pred, flag = eval_tree_array(sub_tree, x_train, regressor.operators_) + return get_loss_function("mse")(y_pred, y_train) + end + + function optimizer_wrapper(population::Vector{Chromosome}) + try + eqn, result = optimize_constants!(population[1].compiled_function, optimizer_function_; + opt_method=opt_method_const, max_iterations=max_iterations, n_restarts=n_starts) + population[1].fitness = (result,) + population[1].compiled_function = eqn + catch + @show "Ignored constant opt." + end + end + + evalStrat = StandardRegressionStrategy{typeof(first(x_train))}( regressor.operators_, x_train, y_train, - regressor.toolbox_; + !isnothing(x_test) ? x_test : x_train, + !isnothing(y_test) ? y_test : y_train, + get_loss_function(loss_fun); + secOptimizer=optimizer_wrapper, + break_condition=break_condition + ) + + best, history = runGep(epochs, + population_size, + regressor.toolbox_, + evalStrat; + hof=hof, + correction_callback=correction_callback, + correction_epochs=correction_epochs, + correction_amount=correction_amount, + tourni_size=max(Int(ceil(population_size * 0.03)), 3), + optimization_epochs=optimization_epochs + ) + + regressor.best_models_ = best + regressor.fitness_history_ = history +end + +function fit!(regressor::GepRegressor, epochs::Int, population_size::Int, loss_function::Function; + optimizer_function_::Union{Function,Nothing}=nothing, + optimization_epochs::Int=100, + hof::Int=3, + correction_epochs::Int=1, + correction_amount::Real=0.05, + opt_method_const::Symbol=:cg, + target_dimension::Union{Vector{Float16},Nothing}=nothing, + cycles::Int=10, max_iterations::Int=150, n_starts::Int=5, + break_condition::Union{Function,Nothing}=nothing +) + + correction_callback = if !isnothing(target_dimension) + (genes, start_indices, expression) -> correct_genes!( + genes, + start_indices, + expression, + target_dimension, + regressor.token_dto_; + cycles=cycles + ) + else + nothing + end + + + function optimizer_wrapper(population::Vector{Chromosome}) + try + eqn, result = optimize_constants!(population[1].compiled_function, optimizer_function_; + opt_method=opt_method_const, max_iterations=max_iterations, n_restarts=n_starts) + population[1].fitness = result + population[1].compiled_function = eqn + catch e + @show "Ignored constant opt." + end + end + + evalStrat = GenericRegressionStrategy( + regressor.operators_, + length(regressor.toolbox_.fitness_reset[1]), + loss_function; + secOptimizer=nothing, + break_condition=break_condition + ) + + best, history = runGep(epochs, + population_size, + regressor.toolbox_, + evalStrat; hof=hof, - x_data_test=!isnothing(x_test) ? x_test : x_train, - y_data_test=!isnothing(y_test) ? y_test : y_train, - loss_fun_=loss_fun, correction_callback=correction_callback, correction_epochs=correction_epochs, correction_amount=correction_amount, - tourni_size=tourni_size, - opt_method_const=opt_method_const, - optimisation_epochs=optimization_epochs) + tourni_size=max(Int(ceil(population_size * 0.03)), 3), + optimization_epochs=optimization_epochs + ) regressor.best_models_ = best regressor.fitness_history_ = history end +function fit!(regressor::GepTensorRegressor, epochs::Int, population_size::Int, loss_function::Function; + hof::Int=3, + break_condition::Union{Function,Nothing}=nothing +) + + evalStrat = GenericRegressionStrategy( + nothing, + length(regressor.toolbox_.fitness_reset[1]), + loss_function; + secOptimizer=nothing, + break_condition=break_condition + ) + + best, history = runGep(epochs, + population_size, + regressor.toolbox_, + evalStrat; + hof=hof, + tourni_size=max(Int(ceil(population_size * 0.03)), 3) + ) + + regressor.best_models_ = best + regressor.fitness_history_ = history +end """ (regressor::GepRegressor)(x_data::AbstractArray; ensemble::Bool=false) @@ -732,7 +846,7 @@ Make predictions using the trained regressor. - `ensemble::Bool=false`: Whether to use ensemble predictions # Returns -- Predicted values for the input features +- Predicted values for the input features -> only employable when using compile_djl_datatype """ function (regressor::GepRegressor)(x_data::AbstractArray; ensemble::Bool=false) return regressor.best_models_[1].compiled_function(x_data, regressor.operators_) @@ -761,11 +875,11 @@ println(sin_info.arity) # 1 """ function list_all_functions() return Dict(sym => ( - function_ = _FUNCTION_LIB_COMMON[sym], - arity = _ARITY_LIB_COMMON[sym], - forward_handler = _FUNCTION_LIB_FORWARD_COMMON[sym], - backward_handler = _FUNCTION_LIB_BACKWARD_COMMON[sym] - ) for sym in keys(_FUNCTION_LIB_COMMON)) + function_=FUNCTION_LIB_COMMON[sym], + arity=ARITY_LIB_COMMON[sym], + forward_handler=FUNCTION_LIB_FORWARD_COMMON[sym], + backward_handler=FUNCTION_LIB_BACKWARD_COMMON[sym] + ) for sym in keys(FUNCTION_LIB_COMMON)) end """ @@ -860,8 +974,8 @@ set_function!(:sin, new_sin_implementation) - ArgumentError if the function symbol doesn't exist in the library """ function set_function!(sym::Symbol, func::Function) - haskey(_FUNCTION_LIB_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) - _FUNCTION_LIB_COMMON[sym] = func + haskey(_UNCTION_LIB_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) + FUNCTION_LIB_COMMON[sym] = func return nothing end @@ -883,9 +997,9 @@ set_arity!(:custom_func, 2) - ArgumentError if the function symbol doesn't exist or arity is invalid """ function set_arity!(sym::Symbol, arity::Int8) - haskey(_ARITY_LIB_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) + haskey(ARITY_LIB_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) arity in (1, 2) || throw(ArgumentError("Arity must be 1 or 2")) - _ARITY_LIB_COMMON[sym] = arity + ARITY_LIB_COMMON[sym] = arity return nothing end @@ -907,8 +1021,8 @@ set_forward_handler!(:custom_func, zero_unit_forward) - ArgumentError if the function symbol doesn't exist """ function set_forward_handler!(sym::Symbol, handler::Function) - haskey(_FUNCTION_LIB_FORWARD_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) - _FUNCTION_LIB_FORWARD_COMMON[sym] = handler + haskey(FUNCTION_LIB_FORWARD_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) + FUNCTION_LIB_FORWARD_COMMON[sym] = handler return nothing end @@ -930,8 +1044,8 @@ set_backward_handler!(:custom_func, zero_unit_backward) - ArgumentError if the function symbol doesn't exist """ function set_backward_handler!(sym::Symbol, handler::Function) - haskey(_FUNCTION_LIB_BACKWARD_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) - _FUNCTION_LIB_BACKWARD_COMMON[sym] = handler + haskey(FUNCTION_LIB_BACKWARD_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) + FUNCTION_LIB_BACKWARD_COMMON[sym] = handler return nothing end @@ -962,13 +1076,13 @@ update_function!(:custom_func, # Throws - ArgumentError if the function symbol doesn't exist or parameters are invalid """ -function update_function!(sym::Symbol; - func::Union{Function,Nothing}=nothing, - arity::Union{Int8,Nothing}=nothing, - forward_handler::Union{Function,Nothing}=nothing, - backward_handler::Union{Function,Nothing}=nothing) +function update_function!(sym::Symbol; + func::Union{Function,Nothing}=nothing, + arity::Union{Int8,Nothing}=nothing, + forward_handler::Union{Function,Nothing}=nothing, + backward_handler::Union{Function,Nothing}=nothing) haskey(_FUNCTION_LIB_COMMON, sym) || throw(ArgumentError("Function $sym not found in library")) - + if !isnothing(func) set_function!(sym, func) end @@ -981,7 +1095,7 @@ function update_function!(sym::Symbol; if !isnothing(backward_handler) set_backward_handler!(sym, backward_handler) end - + return nothing end diff --git a/src/Sbp.jl b/src/Sbp.jl index 2e8b490..e68b852 100644 --- a/src/Sbp.jl +++ b/src/Sbp.jl @@ -1014,7 +1014,7 @@ function propagate_necessary_changes!( return true end - if check_crit_up!(tree.depend_on_total_number+1, expected_dim, tree) && distance_to_change <= 0 && rand() > 0.25 + if check_crit_up!(tree.depend_on_total_number+1, expected_dim, tree) && distance_to_change <= 0 && rand() > 0.05 return enforce_changes!(tree, expected_dim) end diff --git a/src/Selection.jl b/src/Selection.jl index 82895f7..2d36850 100644 --- a/src/Selection.jl +++ b/src/Selection.jl @@ -1,15 +1,18 @@ module EvoSelection using LinearAlgebra -export selection_NSGA, basic_tournament_selection, dominates_, fast_non_dominated_sort, calculate_fronts, determine_ranks, assign_crowding_distance +export tournament_selection, nsga_selection, dominates_, fast_non_dominated_sort, calculate_fronts, determine_ranks, assign_crowding_distance +struct SelectedMembers + indices::Vector{Int} + fronts::Dict{Int,Vector{Int}} +end #Note: selection is constructed to allways return a list of indices => {just care about the data not the objects} - -@inline function basic_tournament_selection(population::AbstractArray{T}, tournament_size::Int, number_of_winners::Int) where {T<:Number} +@inline function tournament_selection(population::AbstractArray{Tuple}, number_of_winners::Int, tournament_size::Int) selected_indices = Vector{Int}(undef, number_of_winners) - valid_indices_ = findall(isfinite, population) + valid_indices_ = findall(x -> isfinite(x[1]), population) valid_indices = [] doubles = Set() for elem in valid_indices_ @@ -25,7 +28,7 @@ export selection_NSGA, basic_tournament_selection, dominates_, fast_non_dominate selected_indices[index] = winner end selected_indices[end] = 1 - return selected_indices + return SelectedMembers(selected_indices,Dict{Int,Vector{Int}}()) end @@ -161,38 +164,35 @@ end return distances end -@inline function selection_NSGA(population::Vector{T}, num_to_select::Int) where {T<:Tuple} + +@inline function nsga_selection(population::Vector{T}) where {T<:Tuple} fronts = calculate_fronts(population) n_fronts = length(fronts) selected_indices = Int[] pareto_fronts = Dict{Int,Vector{Int}}() - for front_idx in 1:n_fronts - front = fronts[front_idx] - pareto_fronts[front_idx] = front - - if length(selected_indices) + length(front) <= num_to_select - append!(selected_indices, front) - else - remaining_slots = num_to_select - length(selected_indices) - crowding_distances = assign_crowding_distance(front, population) - - sorted_front = sort(front, by=i -> crowding_distances[i], rev=true) + estimated_size = length(population) + selected_indices = Vector{Int}(undef, estimated_size) + current_idx = 1 - append!(selected_indices, sorted_front[1:remaining_slots]) - break - end - if length(selected_indices) >= num_to_select - break - end + @inbounds for front_idx in 1:n_fronts + front = fronts[front_idx] + crowding_distances = assign_crowding_distance(front, population) + + sorted_front = sort(front, by=i -> crowding_distances[i], rev=true) + front_size = length(sorted_front) + copyto!(selected_indices, current_idx, sorted_front, 1, front_size) + + pareto_fronts[front_idx] = sorted_front + current_idx +=front_size end + resize!(selected_indices, current_idx - 1) - if length(selected_indices) < num_to_select - @warn "Not enough individuals in the population to select $(num_to_select). Returning $(length(selected_indices)) individuals." - end + return SelectedMembers(selected_indices, pareto_fronts) - return selected_indices, pareto_fronts end + + end diff --git a/src/TensorOps.jl b/src/TensorOps.jl new file mode 100644 index 0000000..ddfd735 --- /dev/null +++ b/src/TensorOps.jl @@ -0,0 +1,409 @@ +module TensorRegUtils + +using Flux, LinearAlgebra, OrderedCollections, ChainRulesCore, Tensors, PrecompileTools + + +struct ThreadBuffer + vector::Vector{Union{Number,Tensor,SymmetricTensor}} +end + +const THREAD_BUFFERS = Vector{ThreadBuffer}(undef, Threads.nthreads()) + +function __init__() + for i in 1:Threads.nthreads() + THREAD_BUFFERS[i] = ThreadBuffer(Vector{Union{Number,Tensor,SymmetricTensor}}(undef, 0)) + end +end + +@inline function get_thread_buffer(n::Integer) + buffer = THREAD_BUFFERS[Threads.threadid()].vector + if length(buffer) < n + resize!(buffer, n) + end + buffer +end + + +# Abstract base type with parametric types for improved type stability +abstract type AbstractOperationNode{T} end + +# Input selector with strict typing +struct InputSelector{T<:Integer} + idx::T +end + +@inline function (l::InputSelector{T})(x::Tuple) where {T} + @inbounds x[l.idx] +end + +@inline function (l::InputSelector{T})(x::Any) where {T} + @inbounds x +end + +# Macro to generate specialized operation nodes with functors +macro generate_operation_node(name) + return quote + struct $(esc(name)){T<:Union{Nothing,Chain}} <: AbstractOperationNode{T} + chain::T + $(esc(name))(chain=nothing) = new{typeof(chain)}(chain) + end + Flux.Functors.@functor $(esc(name)) + end +end + +# Generate concrete operation nodes +@generate_operation_node AdditionNode +@generate_operation_node SubtractionNode +@generate_operation_node MultiplicationNode +@generate_operation_node DivisionNode +@generate_operation_node PowerNode +@generate_operation_node MinNode +@generate_operation_node MaxNode +@generate_operation_node InversionNode +@generate_operation_node TraceNode +@generate_operation_node DeterminantNode +@generate_operation_node SymmetricNode +@generate_operation_node SkewNode +@generate_operation_node VolumetricNode +@generate_operation_node DeviatricNode +@generate_operation_node TdotNode +@generate_operation_node DottNode +@generate_operation_node DoubleContractionNode +@generate_operation_node DeviatoricNode + +# Specialized nodes with their functors +struct ConstantNode{T<:Number,C<:Union{Nothing,Chain}} <: AbstractOperationNode{C} + value::T + chain::C + ConstantNode(value::T, chain=nothing) where {T<:Number} = new{T,typeof(chain)}(value, chain) +end +Flux.Functors.@functor ConstantNode + +struct UnaryNode{F<:Function,C<:Union{Nothing,Chain}} <: AbstractOperationNode{C} + operation::F + chain::C + UnaryNode(operation::F, chain=nothing) where {F<:Function} = new{F,typeof(chain)}(operation, chain) +end +Flux.Functors.@functor UnaryNode + +# Operation implementations +@inline function (l::AdditionNode{T})(x::Union{Tensor,SymmetricTensor}, + y::Union{Tensor,SymmetricTensor}) where {T} + @fastmath (x + y)::Union{Tensor,SymmetricTensor} +end + +@inline function (l::AdditionNode{T})(x::Number, y::Number) where {T} + @fastmath (x + y)::Number +end + + +@inline function (l::AdditionNode{T})(x::Any, y::Any) where {T} + Inf::Number +end + +@inline function (l::MultiplicationNode{T})(x::Any, y::Any) where {T} + Inf::Number +end + +@inline function (l::SubtractionNode{T})(x::Union{Tensor,SymmetricTensor}, y::Union{Tensor,SymmetricTensor}) where {T} + @fastmath (x - y)::Union{Tensor,SymmetricTensor} +end + + +@inline function (l::SubtractionNode{T})(x::Number, y::Number) where {T} + @fastmath (x - y)::Number +end + +@inline function (l::SubtractionNode{T})(x::Any, y::Any) where {T} + Inf::Number +end + +@inline function (l::MultiplicationNode{T})(x::Number, y::Number) where {T} + @fastmath (x * y)::Number +end + +@inline function (l::MultiplicationNode{T})(x::Union{Tensor,SymmetricTensor}, y::Number) where {T} + @fastmath (x * y)::Union{Tensor,SymmetricTensor} +end + +@inline function (l::MultiplicationNode{T})(x::Number, y::Union{Tensor,SymmetricTensor}) where {T} + @fastmath (x * y)::Union{Tensor,SymmetricTensor} +end + + + +@inline function (l::MultiplicationNode{T})(x::Union{Tensor,SymmetricTensor}, y::Union{Tensor,SymmetricTensor}) where {T} + @fastmath dot(x, y)::Union{Tensor,SymmetricTensor,Number} +end + +@inline function (l::DivisionNode{T})(x::Union{Tensor,SymmetricTensor,Number}, y::Number) where {T} + @fastmath (x / y)::Union{Tensor,SymmetricTensor,Number} +end + +@inline function (l::DivisionNode{T})(x::Any, y::Any) where {T} + Inf::Number +end + +@inline function (l::PowerNode{T})(x::Union{Tensor,SymmetricTensor,Number}, y::Number) where {T} + @fastmath (x^y)::Union{Tensor,SymmetricTensor} +end + +@inline function (l::DoubleContractionNode{T})(x, y) where {T} + @fastmath dcontract(x, y) +end + +@inline function (l::DeviatoricNode{T})(x) where {T} + @fastmath dev(x) +end + +@inline function (l::MinNode{T})(x, y) where {T} + @fastmath min(x, y) +end + +@inline function (l::MaxNode{T})(x, y) where {T} + @fastmath max(x, y) +end + +@inline function (l::InversionNode{T})(x) where {T} + @fastmath inv(x) +end + +@inline function (l::TraceNode{T})(x) where {T} + @fastmath tr(x) +end + +@inline function (l::DeterminantNode{T})(x) where {T} + @fastmath det(x)::Number +end + +@inline function (l::SymmetricNode{T})(x::Union{Tensor,SymmetricTensor}) where {T} + @fastmath symmetric(x)::Union{Tensor,SymmetricTensor} +end + +@inline function (l::SkewNode{T})(x::Union{Tensor,SymmetricTensor}) where {T} + @fastmath skew(x)::Union{Tensor,SymmetricTensor} +end + +@inline function (l::VolumetricNode{T})(x) where {T} + @fastmath vol(x) +end + +@inline function (l::DeviatricNode{T})(x) where {T} + @fastmath dev(x) +end + +@inline function (l::TdotNode{T})(x) where {T} + @fastmath tdot(x) +end + +@inline function (l::DottNode{T})(x) where {T} + @fastmath dott(x) +end + +@inline function (l::ConstantNode{V,T})() where {V,T} + l.value +end + +@inline function (l::UnaryNode{F,T})(x) where {F,T} + @fastmath l.operation.(x) +end + + +@inline function (l::AdditionNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::AdditionNode{T})(x::AbstractVector{Number}, y::AbstractVector{Number}) where {T} + (x .+ y)::AbstractVector{Number} +end + +@inline function (l::SubtractionNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::MultiplicationNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::DivisionNode{T})(x::AbstractVector, y::Number) where {T} + map(a -> l(a, y), x)::AbstractVector +end + +@inline function (l::DivisionNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::PowerNode{T})(x::AbstractVector, y::Number) where {T} + map(a -> l(a, y), x)::AbstractVector +end + +@inline function (l::PowerNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::MinNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::MaxNode{T})(x::AbstractVector, y::AbstractVector) where {T} + map((a, b) -> l(a, b), x, y)::AbstractVector +end + +@inline function (l::TraceNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::DeterminantNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::SymmetricNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::SkewNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::VolumetricNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::DeviatricNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::InversionNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::TdotNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::DottNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + +@inline function (l::DeviatoricNode{T})(x::AbstractVector) where {T} + map(l, x)::AbstractVector +end + + + +function compile_to_flux_network(rek_string::Vector, arity_map::OrderedDict, callbacks::Dict, nodes::OrderedDict, pre_len::Int) + stack = [] + for elem in reverse(rek_string) + if get(arity_map, elem, 0) == 2 + op1 = pop!(stack) + op2 = pop!(stack) + node = callbacks[elem]() + push!(stack, (inputs) -> node(op1(inputs), op2(inputs))) + elseif get(arity_map, elem, 0) == 1 + op = pop!(stack) + node = callbacks[elem]() + push!(stack, (inputs) -> node(op(inputs))) + else + if nodes[elem] isa Number + num = nodes[elem] + push!(stack, (inputs) -> ConstantNode(num)()) + else + idx = nodes[elem].idx + push!(stack, (inputs) -> InputSelector(idx)(inputs)) + end + end + + end + return Chain(pop!(stack)) +end + +# Constant mappings +const TENSOR_NODES = Dict{Symbol,Type}( + :+ => AdditionNode, + :- => SubtractionNode, + :* => MultiplicationNode, + :/ => DivisionNode, + :^ => PowerNode, + :min => MinNode, + :max => MaxNode, + :inv => InversionNode, + :tr => TraceNode, + :det => DeterminantNode, + :symmetric => SymmetricNode, + :skew => SkewNode, + :vol => VolumetricNode, + :dev => DeviatricNode, + :tdot => TdotNode, + :dott => DottNode, + :dcontract => DoubleContractionNode, + :deviator => DeviatoricNode +) + +const TENSOR_NODES_ARITY = Dict{Symbol,Int8}( + :+ => 2, :- => 2, :* => 2, :/ => 2, :^ => 2, + :min => 2, :max => 2, + :inv => 1, :tr => 1, :det => 1, + :symmetric => 1, :skew => 1, :vol => 1, :dev => 1, + :tdot => 1, :dott => 1, :dcontract => 2, :deviator => 1 +) + +# Exports +export OperationNode, InputSelector +export AdditionNode, SubtractionNode, MultiplicationNode, DivisionNode, PowerNode +export MinNode, MaxNode, InversionNode +export TraceNode, DeterminantNode, SymmetricNode, SkewNode +export VolumetricNode, DeviatricNode, TdotNode, DottNode +export DoubleContractionNode, DeviatoricNode +export ConstantNode, UnaryNode +export compile_to_flux_network +export TENSOR_NODES, TENSOR_NODES_ARITY + + +@setup_workload begin + dim = 3 + t2 = rand(Tensor{2,dim}) + vec3 = rand(Tensor{1,dim}) + const_val = 2.0 + nodes = OrderedDict{Int8,Any}( + 5 => InputSelector(1), + 6 => InputSelector(2), + 7 => const_val + ) + arity_map = OrderedDict{Int8,Int}( + 1 => 2, # + (AdditionNode) + 2 => 2, # * (MultiplicationNode) + 3 => 2, # dcontract + 4 => 1 # tr + ) + callbacks = Dict{Int8,Any}( + 1 => AdditionNode, + 2 => MultiplicationNode, + 3 => DoubleContractionNode, + 4 => TraceNode + ) + + expressions = [ + Int8[2, 5, 5], + Int8[2, 5, 5], + Int8[1, 5, 5], + Int8[2, 6, 7], + Int8[1, 5, 6], + Int8[3, 5, 5], + Int8[4, 5] + ] + inputs = (t2, vec3, const_val) + inputs2 = ([t2 for _ in 1:10],[vec3 for _ in 1:10],[const_val for _ in 1:10]) + @compile_workload begin + for expr in expressions + net = compile_to_flux_network(expr, arity_map, callbacks, nodes, 0) + try + net(inputs) + net(inputs2) + catch + # Ignore runtime dimension mismatch for precompile + end + end + end +end + +end \ No newline at end of file diff --git a/src/Util.jl b/src/Util.jl index d5eeaeb..59a2dc9 100644 --- a/src/Util.jl +++ b/src/Util.jl @@ -48,6 +48,20 @@ operations, including optimization, history tracking, data manipulation, and sta - `train_test_split`: Split dataset for training/testing - `save_state`, `load_state`: State persistence +## Constants +- `FUNCTION_LIB_COMMON`: Available mathematical functions + +# Function Library +Includes extensive mathematical operations: +- Basic arithmetic: +, -, *, /, ^ +- Comparisons: min, max +- Rounding: floor, ceil, round +- Exponential: exp, log, log10, log2 +- Trigonometric: sin, cos, tan, asin, acos, atan +- Hyperbolic: sinh, cosh, tanh, asinh, acosh, atanh +- Special: sqr, sqrt, sign, abs +- Tensorfunctions: + # Usage Example ```julia # History recording @@ -86,7 +100,10 @@ optimized_node, final_loss = optimize_constants!( - `Zygote`: Automatic differentiation - `Serialization`: State persistence - `Statistics`: Statistical computations +- `Flux`: ML-Package +- `Tensors`: mathematical objects of higher order - `Random`: Random number generation +- `CUDA`: Extension to run on CUDA-cores """ module GepUtils @@ -95,6 +112,8 @@ export save_state, load_state export create_history_recorder, record_history!, record!, close_recorder! export HistoryRecorder, OptimizationHistory, get_history_arrays export train_test_split +export FUNCTION_LIB_COMMON, ARITY_LIB_COMMON +export TensorNode, compile_network using OrderedCollections using DynamicExpressions @@ -105,16 +124,114 @@ using Zygote using Serialization using Statistics using Random +using Tensors +using Flux using Base.Threads: @spawn -struct OptimizationHistory{T<:AbstractFloat} + +function sqr(x::Vector{T}) where {T<:AbstractFloat} + return x .* x +end + +function sqr(x::T) where {T<:Union{AbstractFloat,Node{<:AbstractFloat}}} + return x * x +end + +""" + FUNCTION_LIB_COMMON::Dict{Symbol,Function} + +Dictionary mapping function symbols to their corresponding functions. +Contains basic mathematical operations, trigonometric, and other common functions. + +# Available Functions +- Basic arithmetic: `+`, `-`, `*`, `/`, `^` +- Comparison: `min`, `max` +- Rounding: `floor`, `ceil`, `round` +- Exponential & Logarithmic: `exp`, `log`, `log10`, `log2` +- Trigonometric: `sin`, `cos`, `tan`, `asin`, `acos`, `atan` +- Hyperbolic: `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh` +- Other: `abs`, `sqr`, `sqrt`, `sign` +- Tensor Functions: + +To add a new function, ensure you also add corresponding entries in `ARITY_LIB_COMMON`, +`FUNCTION_LIB_FORWARD_COMMON`, and `FUNCTION_LIB_BACKWARD_COMMON`. +""" +const FUNCTION_LIB_COMMON = Dict{Symbol,Function}( + :+ => +, + :- => -, + :* => *, + :/ => /, + :^ => ^, + :min => min, + :max => max, :abs => abs, + :floor => floor, + :ceil => ceil, + :round => round, :exp => exp, + :log => log, + :log10 => log10, + :log2 => log2, :sin => sin, + :cos => cos, + :tan => tan, + :asin => asin, + :acos => acos, + :atan => atan, :sinh => sinh, + :cosh => cosh, + :tanh => tanh, + :asinh => asinh, + :acosh => acosh, + :atanh => atanh, :sqr => sqr, + :sqrt => sqrt, :sign => sign +) + +""" + ARITY_LIB_COMMON::Dict{Symbol,Int8} + +Dictionary specifying the number of arguments (arity) for each function in the library. +- Value of 1 indicates unary functions (e.g., `sin`, `cos`, `abs`) +- Value of 2 indicates binary functions (e.g., `+`, `-`, `*`, `/`) + +When adding new functions to `FUNCTION_LIB_COMMON`, ensure to specify their arity here. +""" +const ARITY_LIB_COMMON = Dict{Symbol,Int8}( + :+ => 2, + :- => 2, + :* => 2, + :/ => 2, + :^ => 2, + :min => 2, + :max => 2, + :abs => 1, + :floor => 1, + :ceil => 1, + :round => 1, + :exp => 1, + :log => 1, + :log10 => 1, + :log2 => 1, + :sin => 1, + :cos => 1, + :tan => 1, + :asin => 1, + :acos => 1, + :atan => 1, + :sinh => 1, + :cosh => 1, + :tanh => 1, + :asinh => 1, + :acosh => 1, + :atanh => 1, + :sqrt => 1, + :sqr => 1 +) + +struct OptimizationHistory{T<:Union{AbstractFloat,Tuple}} train_loss::Vector{T} val_loss::Vector{T} train_mean::Vector{T} train_std::Vector{T} - function OptimizationHistory(epochs::Int, ::Type{T}) where {T<:AbstractFloat} + function OptimizationHistory(epochs::Int, ::Type{T}) where {T<:Union{AbstractFloat,Tuple}} return new{T}( Vector{T}(undef, epochs), Vector{T}(undef, epochs), @@ -242,12 +359,12 @@ final_history = recorder.history - Supports any AbstractFloat type - Channel depth can be adjusted for different recording patterns """ -struct HistoryRecorder{T<:AbstractFloat} +struct HistoryRecorder{T<:Union{AbstractFloat,Tuple}} channel::Channel{Tuple{Int,T,T,Vector{T}}} task::Task history::OptimizationHistory{T} - function HistoryRecorder(epochs::Int, ::Type{T}; buffer_size::Int=32) where {T<:AbstractFloat} + function HistoryRecorder(epochs::Int, ::Type{T}; buffer_size::Int=32) where {T<:Union{AbstractFloat,Tuple}} history = OptimizationHistory(epochs, T) channel = Channel{Tuple{Int,T,T,Vector{T}}}(buffer_size) task = @spawn record_history!(channel, history) @@ -255,17 +372,32 @@ struct HistoryRecorder{T<:AbstractFloat} end end -# Usage in record_history! + +@inline function tuple_agg(entries::Vector{T}, fun::Function) where {T<:Tuple} + isempty(entries) && return entries[1] + N = length(first(entries)) + L = length(entries) + + vectors = ntuple(i -> Vector{Float64}(undef, L), N) + + for (j, entry) in enumerate(entries) + for i in 1:length(entry) + vectors[i][j] = entry[i] + end + end + return tuple(i -> fun(vectors[i]), N) +end + @inline function record_history!( channel::Channel{Tuple{Int,T,T,Vector{T}}}, history::OptimizationHistory{T} -) where {T<:AbstractFloat} +) where {T<:Union{AbstractFloat,Tuple}} for (epoch, train_loss, val_loss, fit_vector) in channel @inbounds begin history.train_loss[epoch] = train_loss history.val_loss[epoch] = val_loss - history.train_mean[epoch] = mean(fit_vector) - history.train_std[epoch] = std(fit_vector) + #history.train_mean[epoch] = tuple_agg(fit_vector,mean) + #history.train_std[epoch] = tuple_agg(fit_vector, std) end end end @@ -276,11 +408,10 @@ end train_loss::T, val_loss::T, fit_vector::Vector{T} -) where {T<:AbstractFloat} +) where {T<:Union{AbstractFloat,Tuple}} put!(recorder.channel, (epoch, train_loss, val_loss, fit_vector)) end - @inline function close_recorder!(recorder::HistoryRecorder) close(recorder.channel) wait(recorder.task) @@ -403,27 +534,26 @@ result = compile_djl_datatype(rek_string, arity_map, callbacks, nodes) See also: [`DynamicExpressions.Node`](@ref) """ -function compile_djl_datatype(rek_string::Vector, arity_map::OrderedDict, callbacks::Dict, nodes::OrderedDict) - #just let it fail when it becomes invalid, because then the equation is not that use ful +function compile_djl_datatype(rek_string::Vector, arity_map::OrderedDict, callbacks::Dict, nodes::OrderedDict, pre_len::Int) stack = [] - for elem in reverse(rek_string) + for elem in reverse(rek_string[pre_len:end]) if get(arity_map, elem, 0) == 2 - op1 = (temp = pop!(stack); temp isa Int8 ? nodes[temp] : temp) - op2 = (temp = pop!(stack); temp isa Int8 ? nodes[temp] : temp) + op1 = pop!(stack) + op2 = pop!(stack) ops = callbacks[elem] push!(stack, ops(op1, op2)) elseif get(arity_map, elem, 0) == 1 - op1 = (temp = pop!(stack); temp isa Int8 ? nodes[temp] : temp) + op1 = pop!(stack) ops = callbacks[elem] push!(stack, ops(op1)) else - push!(stack, elem) + push!(stack, elem isa Int8 ? nodes[elem] : elem) end end - return last(stack) + return pre_len == 1 ? last(stack) : stack end -function retrieve_constants_from_node(node::Node) +@inline function retrieve_constants_from_node(node::Node) constants = AbstractFloat[] for op in node if op isa AbstractNode && op.degree == 0 && op.constant @@ -433,6 +563,7 @@ function retrieve_constants_from_node(node::Node) constants end + """ optimize_constants!( node::Node, @@ -528,7 +659,7 @@ See also: [`DynamicExpressions.Node`](@ref), [`Optim.optimize`](@ref), [`LineSea n_restarts::Int=3 ) - nconst = count_constants(node) + nconst = count_constant_nodes(node) if nconst == 0 return node, 0.0 @@ -609,32 +740,32 @@ end function train_test_split( - X::AbstractMatrix{T}, - y::AbstractVector{T}; - train_ratio::T=0.9, + X::AbstractMatrix{T}, + y::AbstractVector{T}; + train_ratio::T=0.9, consider::Int=1 ) where {T<:AbstractFloat} data = hcat(X, y) - + data = data[shuffle(1:size(data, 1)), :] - + split_point = floor(Int, size(data, 1) * train_ratio) - + data_train = data[1:split_point, :] - data_test = data[(split_point + 1):end, :] - + data_test = data[(split_point+1):end, :] + x_train = T.(data_train[1:consider:end, 1:(end-1)]) y_train = T.(data_train[1:consider:end, end]) - + x_test = T.(data_test[1:consider:end, 1:(end-1)]) y_test = T.(data_test[1:consider:end, end]) - - return x_train, y_train, x_test, y_test + + return x_train, y_train, x_test, y_test end diff --git a/test/Entity_test.jl b/test/Entity_test.jl index 7cafb33..5ebc581 100644 --- a/test/Entity_test.jl +++ b/test/Entity_test.jl @@ -67,7 +67,7 @@ end toolbox.gene_count * (2 * toolbox.head_len + 1)) @test chromosome.compiled == true @test chromosome.dimension_homogene == false - @test isnan(chromosome.fitness) + @test isnan(chromosome.fitness[1]) end @testset "Function Compilation" begin @@ -85,53 +85,6 @@ end @test typeof(result[1]) <: Real end end - - @testset "Genetic Operators" begin - toolbox = create_test_toolbox() - Random.seed!(42) - chromosome1 = generate_chromosome(toolbox) - chromosome2 = generate_chromosome(toolbox) - - genes1_original = copy(chromosome1.genes) - genes2_original = copy(chromosome2.genes) - - @testset "One Point Crossover" begin - c1, c2 = replicate(chromosome1, chromosome2, toolbox) - gene_one_point_cross_over!(c1, c2) - @test c1.genes != genes1_original || c2.genes != genes2_original - end - - @testset "Mutation" begin - c1 = Chromosome(copy(genes1_original), toolbox) - gene_mutation!(c1, 1.0) # Force mutation - @test c1.genes != genes1_original - end - - @testset "Gene Fusion" begin - c1, c2 = replicate(chromosome1, chromosome2, toolbox) - gene_fussion!(c1, c2, 1.0) # Force fusion - @test c1.genes != genes1_original || c2.genes != genes2_original - end - - @testset "Inversion" begin - pos_11 = toolbox.gene_count - pos_12 = toolbox.gene_count+toolbox.head_len-1 - pos_21 = pos_12 + toolbox.head_len+2 - pos_22 = pos_21 + toolbox.head_len-1 - - c1 = Chromosome(copy(genes1_original), toolbox) - - - initial_head = copy(c1.genes[pos_11:pos_12]) - initial_head2 = copy(c1.genes[pos_21:pos_22]) - - @show initial_head - @show initial_head2 - - gene_inversion!(c1) - @test c1.genes[toolbox.gene_count+1:toolbox.head_len+toolbox.gene_count] != initial_head || initial_head2 != c1.genes[9:11] - end - end @testset "Population Generation" begin toolbox = create_test_toolbox() @@ -142,18 +95,4 @@ end @test all(x -> x isa Chromosome, population) @test length(unique([p.genes for p in population])) == population_size end - - @testset "History Recording" begin - toolbox = create_test_toolbox() - population_size = 10 - population = generate_population(population_size, toolbox) - - for (i, chromo) in enumerate(population) - set_fitness!(chromo, Float64(i)) - @test fitness(chromo) == Float64(i) - end - - sorted_pop = sort(population, by=fitness) - @test issorted([fitness(c) for c in sorted_pop]) - end end \ No newline at end of file diff --git a/test/Main_min_bench.jl b/test/Main_min_bench.jl new file mode 100644 index 0000000..478c2e1 --- /dev/null +++ b/test/Main_min_bench.jl @@ -0,0 +1,22 @@ +include("../src/GeneExpressionProgramming.jl") + +using .GeneExpressionProgramming +using Random +using Plots +using BenchmarkTools + +Random.seed!(1) + +#Define the iterations for the algorithm and the population size +epochs = 10 +population_size = 20 + +#Number of features which needs to be inserted +number_features = 2 + +x_data = randn(Float64, 10000, number_features) +y_data = @. x_data[:,1] * x_data[:,1] + x_data[:,1] * x_data[:,2] - 2 * x_data[:,2] * x_data[:,2] + +#define the +regressor = GepRegressor(number_features) +@btime fit!(regressor, epochs, population_size, x_data', y_data; loss_fun="mse") \ No newline at end of file diff --git a/test/Main_min_example.jl b/test/Main_min_example.jl index 34017a3..a1d89f2 100644 --- a/test/Main_min_example.jl +++ b/test/Main_min_example.jl @@ -7,19 +7,20 @@ using Plots Random.seed!(1) #Define the iterations for the algorithm and the population size -epochs = 1000 -population_size = 1000 +epochs = 100 +population_size = 100 #Number of features which needs to be inserted number_features = 2 x_data = randn(Float64, 100, number_features) -y_data = @. x_data[:,1] * x_data[:,1] + x_data[:,1] * x_data[:,2] - 2 * x_data[:,1] * x_data[:,2] +y_data = @. x_data[:,1] * x_data[:,1] + x_data[:,1] * x_data[:,2] - 2 * x_data[:,2] * x_data[:,2] #define the regressor = GepRegressor(number_features) fit!(regressor, epochs, population_size, x_data', y_data; loss_fun="mse") + pred = regressor(x_data') @show regressor.best_models_[1].compiled_function @@ -40,7 +41,7 @@ color=:red) #train loss vs validation loss train_validation = plot( - regressor.fitness_history_.train_loss, + [x[1] for x in regressor.fitness_history_.train_loss], label="Training Loss", ylabel="Loss", xlabel="Epoch", @@ -49,7 +50,7 @@ train_validation = plot( plot!( train_validation, - regressor.fitness_history_.val_loss, + [x[1] for x in regressor.fitness_history_.val_loss], label="Validation Loss", linewidth=2 ) diff --git a/test/Main_min_example_mo.jl b/test/Main_min_example_mo.jl new file mode 100644 index 0000000..960e211 --- /dev/null +++ b/test/Main_min_example_mo.jl @@ -0,0 +1,46 @@ +include("../src/GeneExpressionProgramming.jl") + +using .GeneExpressionProgramming +using Random +using Plots +using DynamicExpressions +using Statistics + +const Chromosome = GepEntities.Chromosome + + +Random.seed!(1) + +#Define the iterations for the algorithm and the population size +epochs = 1000 +population_size = 1000 + +#Number of features which needs to be inserted +number_features = 2 + +x_data = randn(Float64, 100, number_features) +y_data = @. x_data[:,1] * x_data[:,1] + x_data[:,1] * x_data[:,2] - 2 * x_data[:,2] * x_data[:,2] + +#define the +regressor = GepRegressor(number_features; number_of_objectives=2) + +@inline function loss_new(elem,validate::Bool) + try + if isnan(mean(elem.fitness)) || validate + y_pred = elem.compiled_function(x_data', regressor.operators_) + return (get_loss_function("mse")(y_data, y_pred), length(elem.expression_raw)*0.01) + else + return (elem.fitness,length(elem.expression_raw)*elem.fitness) + end + catch e + return (typemax(Float64),typemax(Float64)) + end +end + +fit!(regressor, epochs, population_size, loss_new) + + +pred = regressor(x_data') + +@show regressor.best_models_[1].compiled_function +@show regressor.best_models_[1].fitness \ No newline at end of file diff --git a/test/Main_min_with_csv.jl b/test/Main_min_with_csv.jl index b0f2487..f77b55f 100644 --- a/test/Main_min_with_csv.jl +++ b/test/Main_min_with_csv.jl @@ -53,7 +53,7 @@ plot!(pred_vs_actual, vec(y_test), vec(y_test), #train loss vs validation loss train_validation = plot( - regressor.fitness_history_.train_loss, + [x[1] for x in regressor.fitness_history_.train_loss], label="Training Loss", ylabel="Loss", xlabel="Epoch", @@ -62,7 +62,7 @@ train_validation = plot( plot!( train_validation, - regressor.fitness_history_.val_loss, + [x[1] for x in regressor.fitness_history_.val_loss], label="Validation Loss", linewidth=2 ) diff --git a/test/Spb_core_test.jl b/test/Spb_core_test.jl index 8538a12..f044909 100644 --- a/test/Spb_core_test.jl +++ b/test/Spb_core_test.jl @@ -6,7 +6,7 @@ include("../src/Util.jl") using .SBPUtils using Random using OrderedCollections - +Random.seed!(1) function create_token_lib_test() physical_dimension_dict = OrderedDict{Int8, Vector{Float16}}( @@ -52,7 +52,7 @@ end features = Int8[7, 8] # x1, x2 functions = Int8[1, 2, 3, 4, 5, 6] # mul, div, add, sub, sqr, sin constants = Int8[] - lib = create_lib(token_lib, features, functions, constants; rounds=2, max_permutations=10000) + lib = create_lib(token_lib, features, functions, constants; rounds=8, max_permutations=10000) total_len_lib = sum(length(entry) for entry in values(lib)) @show ("Lib Entries:" , total_len_lib) @test !isempty(lib) @@ -106,7 +106,7 @@ end features = Int8[7, 8] # x1, x2 functions = Int8[1, 2, 3, 4, 5, 6] # mul, div, add, sub, sqr, sin constants = Int8[] - lib = create_lib(token_lib, features, functions, constants; rounds=6, max_permutations=10000) + lib = create_lib(token_lib, features, functions, constants; rounds=8, max_permutations=100000) point_operations = Int8[1,2] inverse_operation = Dict{Int8, Function}( 1 => (mul_unit_backward), diff --git a/test/non_dom_sort_test.jl b/test/non_dom_sort_test.jl index 09a425b..9bfbba8 100644 --- a/test/non_dom_sort_test.jl +++ b/test/non_dom_sort_test.jl @@ -55,7 +55,6 @@ end ]) fronts3 = calculate_fronts(pop3) - @show fronts3 @test length(fronts3) == 3 @test Set(fronts3[1]) == Set([1]) @test Set(fronts3[2]) == Set([2, 4, 5, 6, 7, 8, 9]) @@ -90,10 +89,10 @@ end (1, 2), (2, 1), (1, 1) ]) - selected2, fronts2 = selection_NSGA(pop2, 10) - @test length(selected2) == 10 - @test 15 in selected2 # Preserve the best!! alllllways - @test length(fronts2) == 4 + selected2 = nsga_selection(pop2) + @test length(selected2.indices) == 15 + @test 15 in selected2.indices # Preserve the best!! alllllways + @test length(selected2.fronts) == 5 # 3 objectives pop3 = create_population([ @@ -101,8 +100,8 @@ end (1, 2, 3), (2, 3, 1), (3, 1, 2), (1, 3, 2), (2, 1, 3), (3, 2, 1) ]) - selected3, fronts3 = selection_NSGA(pop3, 5) - @test length(selected3) == 5 - @test 1 in selected3 # The best individual should always be selected - @test length(fronts3) == 2 -end \ No newline at end of file + selected3 = nsga_selection(pop3) + @test length(selected3.indices) == 9 + @test 1 in selected3.indices # The best individual should always be selected + @test length(selected3.fronts) == 3 +end diff --git a/test/regression_wrapper_test.jl b/test/regression_wrapper_test.jl index 508d150..a3f4ebf 100644 --- a/test/regression_wrapper_test.jl +++ b/test/regression_wrapper_test.jl @@ -51,8 +51,8 @@ Random.seed!(1) @test length(syms) == 4 # 2 constants + 2 random @test all(v -> v == 0, values(syms)) @test length(nodes) == 4 - @test nodes[1] == 1.0 - @test nodes[2] == 2.5 + @test nodes[1] == Node{Float64}(;val=1.0) + @test nodes[2] == Node{Float64}(;val=2.5) end @testset "Physical Operations" begin diff --git a/test/runtests.jl b/test/runtests.jl index 3c1208b..ab370c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ include("Entity_test.jl") include("Spb_core_test.jl") include("non_dom_sort_test.jl") -include("regression_wrapper_test.jl") \ No newline at end of file +include("regression_wrapper_test.jl") +include("tensor_ops_test.jl") \ No newline at end of file diff --git a/test/tensor_ops_test.jl b/test/tensor_ops_test.jl new file mode 100644 index 0000000..5c48ab0 --- /dev/null +++ b/test/tensor_ops_test.jl @@ -0,0 +1,211 @@ +include("../src/TensorOps.jl") +using BenchmarkTools +using Test +using .TensorRegUtils +using Tensors +using OrderedCollections +using Flux + +@testset "TensorRegUtils" begin + # Test setup + dim = 3 + t2 = rand(Tensor{2,dim}) + vec3 = rand(Tensor{1,dim}) + const_val = 2.0 + + data_ = Dict( + 5 => t2, + 6 => vec3, + 7 => const_val + ) + inputs = (t2, vec3, const_val) + + + arity_map = OrderedDict{Int8,Int}( + 1 => 2, # Addition + 2 => 2, # Multiplication + 3 => 2, # Double Contraction + 4 => 1 # Trace + ) + + callbacks = Dict{Int8,Any}( + Int8(1) => AdditionNode, + Int8(2) => MultiplicationNode, + Int8(3) => DoubleContractionNode, + Int8(4) => TraceNode + ) + + @testset "Basic Operations" begin + nodes = OrderedDict{Int8,Any}( + Int8(5) => InputSelector(1), + Int8(6) => InputSelector(2), + Int8(7) => const_val + ) + + # Scalar multiplication + rek_string = Int8[2, 6, 7] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + result = network(inputs) + @test result ≈ vec3 * const_val + + # Addition with dimension mismatch + rek_string = Int8[1, 5, 6] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test_throws DimensionMismatch network(inputs) + end + + @testset "Tensor Operations" begin + nodes = OrderedDict{Int8,Any}( + Int8(5) => InputSelector(1), + Int8(6) => InputSelector(2), + Int8(7) => const_val + ) + + # Double contraction + rek_string = Int8[3, 5, 5] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + result = network(inputs) + @test result ≈ dcontract(t2, t2) + + # Trace + rek_string = Int8[4, 5] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + result = network(inputs) + @test result ≈ tr(t2) + end + + @testset "Complex Expressions" begin + nodes = OrderedDict{Int8,Any}( + Int8(5) => InputSelector(1), + Int8(6) => InputSelector(2), + Int8(7) => const_val + ) + + # (vec3 * const_val) + vec3 + rek_string = Int8[1, 2, 6, 7, 6] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + + result = network(inputs) + @test result ≈ (vec3 * const_val) + vec3 + + # (t2 * vec3) + const_val - should fail + rek_string = Int8[1, 2, 5, 6, 7] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test !isfinite(network(inputs)) + + # tr(t2) * const_val + tr(t2) + rek_string = Int8[1, 2, 4, 5, 7, 4, 5] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + result = network(inputs) + @test result ≈ tr(t2) * const_val + tr(t2) + end + +end + +@testset "All Node Operations" begin + dim = 3 + t2 = rand(Tensor{2,dim}) + t2_2 = rand(Tensor{2,dim}) + vec3 = rand(Vec{dim}) + const_val = 2.0 + + data_ = Dict(19 => t2, 20 => t2_2, 21 => vec3, 22 => const_val) + + inputs = (t2, t2_2, vec3, const_val) + + arity_map = OrderedDict{Int8,Int}( + 1 => 2, 2 => 2, 3 => 2, 4 => 2, 5 => 2, 6 => 2, 7 => 2, # Binary ops + 8 => 1, 9 => 1, 10 => 1, 11 => 1, 12 => 1, 13 => 1, # Unary ops part 1 + 14 => 1, 15 => 1, 16 => 1, 17 => 2, 18 => 1 # Unary ops part 2 + DC + ) + + callbacks = Dict{Int8,Any}( + Int8(1) => AdditionNode, + Int8(2) => SubtractionNode, + Int8(3) => MultiplicationNode, + Int8(4) => DivisionNode, + Int8(5) => PowerNode, + Int8(6) => MinNode, + Int8(7) => MaxNode, + Int8(8) => InversionNode, + Int8(9) => TraceNode, + Int8(10) => DeterminantNode, + Int8(11) => SymmetricNode, + Int8(12) => SkewNode, + Int8(13) => VolumetricNode, + Int8(14) => DeviatricNode, + Int8(15) => TdotNode, + Int8(16) => DottNode, + Int8(17) => DoubleContractionNode, + Int8(18) => DeviatoricNode + ) + + @testset "Basic Node Operations" begin + nodes = OrderedDict{Int8,Any}( + Int8(19) => InputSelector(1), + Int8(20) => InputSelector(2), + Int8(21) => InputSelector(3), + Int8(22) => const_val + ) + + # Test 1: Addition + rek_string = Int8[1, 19, 20] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ t2 + t2_2 + + # Test 2: Subtraction + rek_string = Int8[2, 19, 20] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ t2 - t2_2 + + # Test 3: Multiplication with constant + rek_string = Int8[3, 19, 22] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ t2 * const_val + end + + @testset "Tensor Operations" begin + nodes = OrderedDict{Int8,Any}( + Int8(19) => InputSelector(1), + Int8(20) => InputSelector(2) + ) + + # Test 1: Double Contraction + rek_string = Int8[17, 19, 20] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ dcontract(t2, t2_2) + + # Test 2: Deviatoric + rek_string = Int8[14, 19] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ dev(t2) + + # Test 3: Trace + Symmetric + rek_string = Int8[9, 11, 19] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ tr(symmetric(t2)) + end + + @testset "Complex Compositions" begin + nodes = OrderedDict{Int8,Any}( + Int8(19) => InputSelector(1), + Int8(20) => InputSelector(2), + Int8(22) => const_val + ) + + # Test 1: (t2 ⊡ t2_2) * const_val + rek_string = Int8[3, 17, 19, 20, 22] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ dcontract(t2, t2_2) * const_val + + # Test 2: dev(symmetric(t2)) + t2_2 + rek_string = Int8[1, 14, 11, 19, 20] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ dev(symmetric(t2)) + t2_2 + + # Test 3: tr(t2) * tr(t2_2) + rek_string = Int8[3, 9, 19, 9, 20] + network = compile_to_flux_network(rek_string, arity_map, callbacks, nodes, 0) + @test network(inputs) ≈ tr(t2) * tr(t2_2) + end +end \ No newline at end of file