From 0ac81c6284feddadf3dbc8e26843d77aa5f26375 Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Mon, 13 Jul 2026 21:51:51 +0100 Subject: [PATCH 1/7] Add directed hypergraph prototype examples --- .../01_directed_hypergraph_represtation.jl | 220 ++++ .../02_incidence_matrix_variations.jl | 198 ++++ .../03_masking_operations.jl | 333 ++++++ .../04_gather_operations.jl | 405 +++++++ .../05_basic_transformations.jl | 410 +++++++ .../06_preprocessing_pipeline.jl | 280 +++++ .../07_toy_ml_feature_model.jl | 219 ++++ .../08_crn_parser_to_dihypergraph.jl | 399 +++++++ .../09_package_directed_hypergraph.jl | 320 +++++ .../10_lux_directed_message_passing_layer.jl | 452 +++++++ .../shonali_prototypes/11bz_crn_case_study.jl | 513 ++++++++ .../12_formose_crn_case_study.jl | 514 ++++++++ .../13_species_feature_engineering.jl | 417 +++++++ .../14_bidirectional message passing.jl | 446 +++++++ .../15_reaction level regression.jl | 1043 +++++++++++++++++ .../16_architecture_comparison.jl | 919 +++++++++++++++ .../17 lux directed hypergraph layer.jl | 259 ++++ 17 files changed, 7347 insertions(+) create mode 100644 examples/shonali_prototypes/01_directed_hypergraph_represtation.jl create mode 100644 examples/shonali_prototypes/02_incidence_matrix_variations.jl create mode 100644 examples/shonali_prototypes/03_masking_operations.jl create mode 100644 examples/shonali_prototypes/04_gather_operations.jl create mode 100644 examples/shonali_prototypes/05_basic_transformations.jl create mode 100644 examples/shonali_prototypes/06_preprocessing_pipeline.jl create mode 100644 examples/shonali_prototypes/07_toy_ml_feature_model.jl create mode 100644 examples/shonali_prototypes/08_crn_parser_to_dihypergraph.jl create mode 100644 examples/shonali_prototypes/09_package_directed_hypergraph.jl create mode 100644 examples/shonali_prototypes/10_lux_directed_message_passing_layer.jl create mode 100644 examples/shonali_prototypes/11bz_crn_case_study.jl create mode 100644 examples/shonali_prototypes/12_formose_crn_case_study.jl create mode 100644 examples/shonali_prototypes/13_species_feature_engineering.jl create mode 100644 examples/shonali_prototypes/14_bidirectional message passing.jl create mode 100644 examples/shonali_prototypes/15_reaction level regression.jl create mode 100644 examples/shonali_prototypes/16_architecture_comparison.jl create mode 100644 examples/shonali_prototypes/17 lux directed hypergraph layer.jl diff --git a/examples/shonali_prototypes/01_directed_hypergraph_represtation.jl b/examples/shonali_prototypes/01_directed_hypergraph_represtation.jl new file mode 100644 index 0000000..601e01d --- /dev/null +++ b/examples/shonali_prototypes/01_directed_hypergraph_represtation.jl @@ -0,0 +1,220 @@ +#directed hypergraph representation +#starting with a small toy reaction network + +#packages I am keeping in mind for the project: +#SimpleHypergraphs.jl +#SimpleDirectedHypergraphs.jl +#HyperGraphNeuralNetworks.jl + + +#toy species/nodes + +nodes = ["A", "B", "C", "D", "E", "F"] + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nNodes:") +println(nodes) + +println("\nNode to ID mapping:") +println(node_to_id) + + +#toy reactions as directed hyperedges + +hyperedges = [ + ( + id = 1, + reactants = ["A", "B"], + products = ["C", "D"] + ), + ( + id = 2, + reactants = ["C"], + products = ["E"] + ), + ( + id = 3, + reactants = ["D", "E"], + products = ["F"] + ) +] + +println("\nDirected hyperedges / toy reactions:") +for edge in hyperedges + println("e", edge.id, ": ", + join(edge.reactants, " + "), + " -> ", + join(edge.products, " + ")) +end + + +#just pulling out source and target sides + +function source_set(edge) + return edge.reactants +end + +function target_set(edge) + return edge.products +end + +println("\nSource and target sets:") +for edge in hyperedges + println("Hyperedge e", edge.id) + println(" Source/reactants: ", source_set(edge)) + println(" Target/products: ", target_set(edge)) +end + + +#signed incidence style matrix +#-1 for reactants, +1 for products, 0 otherwise + +function build_directed_incidence(nodes, hyperedges, node_to_id) + incidence = zeros(Int, length(nodes), length(hyperedges)) + + for (j, edge) in enumerate(hyperedges) + for r in edge.reactants + incidence[node_to_id[r], j] = -1 + end + + for p in edge.products + incidence[node_to_id[p], j] = 1 + end + end + + return incidence +end + +incidence = build_directed_incidence(nodes, hyperedges, node_to_id) + +println("\nDirected incidence-style matrix:") +println("Rows = nodes: ", nodes) +println("Columns = hyperedges: e1, e2, e3") +println(incidence) + + +#separate reactant/product membership matrices + +function build_source_target_incidence(nodes, hyperedges, node_to_id) + source_incidence = zeros(Int, length(nodes), length(hyperedges)) + target_incidence = zeros(Int, length(nodes), length(hyperedges)) + + for (j, edge) in enumerate(hyperedges) + for r in edge.reactants + source_incidence[node_to_id[r], j] = 1 + end + + for p in edge.products + target_incidence[node_to_id[p], j] = 1 + end + end + + return source_incidence, target_incidence +end + +source_incidence, target_incidence = + build_source_target_incidence(nodes, hyperedges, node_to_id) + +println("\nSource/reactant incidence matrix:") +println(source_incidence) + +println("\nTarget/product incidence matrix:") +println(target_incidence) + + +#basic size checks for each reaction + +println("\nHyperedge statistics:") + +for edge in hyperedges + source_size = length(edge.reactants) + target_size = length(edge.products) + total_size = source_size + target_size + + println("Hyperedge e", edge.id) + println(" Source size: ", source_size) + println(" Target size: ", target_size) + println(" Total size: ", total_size) +end + + +#how often each species appears + +node_participation = Dict(node => 0 for node in nodes) + +for edge in hyperedges + involved_nodes = vcat(edge.reactants, edge.products) + + for node in involved_nodes + node_participation[node] += 1 + end +end + +println("\nNode participation counts:") +for node in nodes + println(node, " participates in ", node_participation[node], " hyperedge(s)") +end + + +#source vs target counts for each node + +out_degree = Dict(node => 0 for node in nodes) +in_degree = Dict(node => 0 for node in nodes) + +for edge in hyperedges + for r in edge.reactants + out_degree[r] += 1 + end + + for p in edge.products + in_degree[p] += 1 + end +end + +println("\nDirected node degree-style counts:") +for node in nodes + println( + node, + " | out/source count = ", out_degree[node], + " | in/target count = ", in_degree[node] + ) +end + + +#same reactions but using ids instead of names + +function convert_edges_to_ids(hyperedges, node_to_id) + id_edges = [] + + for edge in hyperedges + push!( + id_edges, + ( + id = edge.id, + reactant_ids = [node_to_id[r] for r in edge.reactants], + product_ids = [node_to_id[p] for p in edge.products] + ) + ) + end + + return id_edges +end + +id_hyperedges = convert_edges_to_ids(hyperedges, node_to_id) + +println("\nID-based directed hyperedge representation:") +for edge in id_hyperedges + println("e", edge.id) + println(" Reactant IDs: ", edge.reactant_ids) + println(" Product IDs: ", edge.product_ids) +end + + +#quick summary + +println("\nSummary:") +println("Number of nodes: ", length(nodes)) +println("Number of directed hyperedges: ", length(hyperedges)) +println("Incidence matrix size: ", size(incidence)) \ No newline at end of file diff --git a/examples/shonali_prototypes/02_incidence_matrix_variations.jl b/examples/shonali_prototypes/02_incidence_matrix_variations.jl new file mode 100644 index 0000000..7495628 --- /dev/null +++ b/examples/shonali_prototypes/02_incidence_matrix_variations.jl @@ -0,0 +1,198 @@ +#incidence matrix variations +#playing around with a few ways to store the same reaction network + +#toy reaction network + +nodes = ["A", "B", "C", "D", "E", "F", "G"] + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +reactions = [ + ( + id = 1, + reactants = ["A", "B"], + products = ["C"] + ), + ( + id = 2, + reactants = ["C"], + products = ["D", "E"] + ), + ( + id = 3, + reactants = ["E", "F"], + products = ["G"] + ) +] + +println("\nToy reactions:") +for r in reactions + println("e", r.id, ": ", join(r.reactants, " + "), " -> ", join(r.products, " + ")) +end + + +#signed incidence matrix +#-1 for reactants, +1 for products + +function signed_incidence_matrix(nodes, reactions, node_to_id) + H = zeros(Int, length(nodes), length(reactions)) + + for (j, r) in enumerate(reactions) + for reactant in r.reactants + H[node_to_id[reactant], j] = -1 + end + + for product in r.products + H[node_to_id[product], j] = 1 + end + end + + return H +end + +H_signed = signed_incidence_matrix(nodes, reactions, node_to_id) + +println("\nSigned incidence matrix:") +println("Rows = nodes: ", nodes) +println("Columns = reactions: e1, e2, e3") +println(H_signed) + + +#separate source and target matrices + +function source_incidence_matrix(nodes, reactions, node_to_id) + H_source = zeros(Int, length(nodes), length(reactions)) + + for (j, r) in enumerate(reactions) + for reactant in r.reactants + H_source[node_to_id[reactant], j] = 1 + end + end + + return H_source +end + +function target_incidence_matrix(nodes, reactions, node_to_id) + H_target = zeros(Int, length(nodes), length(reactions)) + + for (j, r) in enumerate(reactions) + for product in r.products + H_target[node_to_id[product], j] = 1 + end + end + + return H_target +end + +H_source = source_incidence_matrix(nodes, reactions, node_to_id) +H_target = target_incidence_matrix(nodes, reactions, node_to_id) + +println("\nSource/reactant incidence matrix:") +println(H_source) + +println("\nTarget/product incidence matrix:") +println(H_target) + + +#unsigned version, ignoring direction for now + +function unsigned_membership_matrix(H_signed) + return abs.(H_signed) +end + +H_unsigned = unsigned_membership_matrix(H_signed) + +println("\nUnsigned hypergraph membership matrix:") +println(H_unsigned) + + +#checking if I can recover reactants/products from one matrix column + +function extract_reactants_from_column(H_signed, edge_index, id_to_node) + reactants = String[] + + for i in 1:size(H_signed, 1) + if H_signed[i, edge_index] == -1 + push!(reactants, id_to_node[i]) + end + end + + return reactants +end + +function extract_products_from_column(H_signed, edge_index, id_to_node) + products = String[] + + for i in 1:size(H_signed, 1) + if H_signed[i, edge_index] == 1 + push!(products, id_to_node[i]) + end + end + + return products +end + +println("\nReconstruct reactions from signed incidence matrix:") +for j in 1:size(H_signed, 2) + reactants = extract_reactants_from_column(H_signed, j, id_to_node) + products = extract_products_from_column(H_signed, j, id_to_node) + + println("e", j, ": ", join(reactants, " + "), " -> ", join(products, " + ")) +end + + +#edge list version of the same incidence information + +function incidence_to_edge_list(H_signed) + edge_list = [] + + for j in 1:size(H_signed, 2) + for i in 1:size(H_signed, 1) + value = H_signed[i, j] + + if value != 0 + push!( + edge_list, + ( + node_id = i, + hyperedge_id = j, + role = value + ) + ) + end + end + end + + return edge_list +end + +edge_list = incidence_to_edge_list(H_signed) + +println("\nEdge-list style representation:") +println("(role = -1 means source/reactant, role = +1 means target/product)") +for item in edge_list + println(item) +end + + +#quick stats from the matrices + +source_sizes = vec(sum(H_source, dims = 1)) +target_sizes = vec(sum(H_target, dims = 1)) +total_sizes = vec(sum(H_unsigned, dims = 1)) + +println("\nHyperedge sizes computed from matrices:") +for j in 1:length(reactions) + println("e", j, + " | source size = ", source_sizes[j], + " | target size = ", target_sizes[j], + " | total size = ", total_sizes[j]) +end + +node_membership_counts = vec(sum(H_unsigned, dims = 2)) + +println("\nNode membership counts computed from unsigned matrix:") +for i in 1:length(nodes) + println(nodes[i], " appears in ", node_membership_counts[i], " hyperedge(s)") +end \ No newline at end of file diff --git a/examples/shonali_prototypes/03_masking_operations.jl b/examples/shonali_prototypes/03_masking_operations.jl new file mode 100644 index 0000000..242de0d --- /dev/null +++ b/examples/shonali_prototypes/03_masking_operations.jl @@ -0,0 +1,333 @@ +#masking operations +#trying out masks on a small reaction network + +using Statistics +using SimpleHypergraphs +using SimpleDirectedHypergraphs + + +#toy reactions + +nodes = ["A", "B", "C", "D", "E", "F", "G", "H"] + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +reactions = [ + ( + id = 1, + reactants = ["A", "B"], + products = ["C", "D"] + ), + ( + id = 2, + reactants = ["C"], + products = ["E"] + ), + ( + id = 3, + reactants = ["D", "E"], + products = ["F"] + ), + ( + id = 4, + reactants = ["F"], + products = ["G", "H"] + ) +] + +println("\nToy directed hypergraph reactions:") +for r in reactions + println("e", r.id, ": ", + join(r.reactants, " + "), + " -> ", + join(r.products, " + ")) +end + + +#small helper functions + +function empty_bool_mask(n) + return falses(n) +end + +function mask_from_nodes(selected_nodes, node_to_id, total_nodes) + mask = falses(total_nodes) + + for node in selected_nodes + mask[node_to_id[node]] = true + end + + return mask +end + +function print_mask_with_nodes(mask, nodes, label) + println("\n", label) + println("Mask: ", mask) + println("Selected nodes: ", nodes[mask]) +end + + +#reactant/source masks + +println("\nReactant/source masks per hyperedge:") + +reactant_masks = Dict{Int, Vector{Bool}}() + +for r in reactions + mask = mask_from_nodes(r.reactants, node_to_id, length(nodes)) + reactant_masks[r.id] = mask + + print_mask_with_nodes(mask, nodes, "e$(r.id) reactant/source mask") +end + + +#product/target masks + +println("\nProduct/target masks per hyperedge:") + +product_masks = Dict{Int, Vector{Bool}}() + +for r in reactions + mask = mask_from_nodes(r.products, node_to_id, length(nodes)) + product_masks[r.id] = mask + + print_mask_with_nodes(mask, nodes, "e$(r.id) product/target mask") +end + + +#active nodes means reactants and products together + +println("\nActive node masks per hyperedge:") + +active_masks = Dict{Int, Vector{Bool}}() + +for r in reactions + involved_nodes = vcat(r.reactants, r.products) + mask = mask_from_nodes(involved_nodes, node_to_id, length(nodes)) + active_masks[r.id] = mask + + print_mask_with_nodes(mask, nodes, "e$(r.id) active node mask") +end + + +#opposite of active nodes + +println("\nInactive node masks per hyperedge:") + +inactive_masks = Dict{Int, Vector{Bool}}() + +for r in reactions + inactive_mask = .!active_masks[r.id] + inactive_masks[r.id] = inactive_mask + + print_mask_with_nodes(inactive_mask, nodes, "e$(r.id) inactive node mask") +end + + +#same masks but as integers + +println("\nInteger masks for reactants/products:") + +for r in reactions + reactant_int_mask = Int.(reactant_masks[r.id]) + product_int_mask = Int.(product_masks[r.id]) + active_int_mask = Int.(active_masks[r.id]) + + println("\ne", r.id) + println("Reactant integer mask: ", reactant_int_mask) + println("Product integer mask: ", product_int_mask) + println("Active integer mask: ", active_int_mask) +end + + +#global masks across the whole reaction network + +global_reactant_mask = falses(length(nodes)) +global_product_mask = falses(length(nodes)) + +for r in reactions + global_reactant_mask .|= reactant_masks[r.id] + global_product_mask .|= product_masks[r.id] +end + +print_mask_with_nodes(global_reactant_mask, nodes, "Global reactant/source mask") +print_mask_with_nodes(global_product_mask, nodes, "Global product/target mask") + +global_active_mask = global_reactant_mask .| global_product_mask +print_mask_with_nodes(global_active_mask, nodes, "Global active node mask") + + +#building masks from an incidence matrix too + +function build_signed_incidence(nodes, reactions, node_to_id) + H = zeros(Int, length(nodes), length(reactions)) + + for (j, r) in enumerate(reactions) + for reactant in r.reactants + H[node_to_id[reactant], j] = -1 + end + + for product in r.products + H[node_to_id[product], j] = 1 + end + end + + return H +end + +H = build_signed_incidence(nodes, reactions, node_to_id) + +println("\nSigned incidence matrix:") +println(H) + +println("\nMasks derived from incidence matrix:") + +for j in 1:size(H, 2) + source_mask = H[:, j] .== -1 + target_mask = H[:, j] .== 1 + active_mask = H[:, j] .!= 0 + + println("\nHyperedge e", j) + println("Source mask from incidence: ", source_mask) + println("Target mask from incidence: ", target_mask) + println("Active mask from incidence: ", active_mask) + println("Source nodes: ", nodes[source_mask]) + println("Target nodes: ", nodes[target_mask]) + println("Active nodes: ", nodes[active_mask]) +end + + +#trying masks on toy node features + +node_features = [ + 1.0 0.2 0.5; + 0.8 0.1 0.4; + 0.3 0.9 0.7; + 0.4 0.6 0.8; + 0.9 0.2 0.1; + 0.5 0.5 0.5; + 0.7 0.3 0.6; + 0.2 0.8 0.9 +] + +println("\nNode feature matrix:") +println(node_features) + +chosen_edge = reactions[1] + +chosen_reactant_mask = reactant_masks[chosen_edge.id] +chosen_product_mask = product_masks[chosen_edge.id] +chosen_active_mask = active_masks[chosen_edge.id] + +println("\nFeature masking for hyperedge e", chosen_edge.id) + +reactant_features = node_features[chosen_reactant_mask, :] +product_features = node_features[chosen_product_mask, :] +active_features = node_features[chosen_active_mask, :] + +println("Reactant features:") +println(reactant_features) + +println("Product features:") +println(product_features) + +println("Active node features:") +println(active_features) + + +#filtering reactions with a simple hyperedge mask + +hyperedge_mask = falses(length(reactions)) + +for (i, r) in enumerate(reactions) + if length(r.reactants) >= 2 + hyperedge_mask[i] = true + end +end + +println("\nHyperedge mask where source size >= 2:") +println(hyperedge_mask) + +selected_reactions = reactions[hyperedge_mask] + +println("Selected reactions:") +for r in selected_reactions + println("e", r.id, ": ", + join(r.reactants, " + "), + " -> ", + join(r.products, " + ")) +end + + +#padding ids and keeping masks for the real entries + +max_source_size = maximum(length(r.reactants) for r in reactions) +max_target_size = maximum(length(r.products) for r in reactions) + +source_id_batch = zeros(Int, length(reactions), max_source_size) +source_padding_mask = falses(length(reactions), max_source_size) + +target_id_batch = zeros(Int, length(reactions), max_target_size) +target_padding_mask = falses(length(reactions), max_target_size) + +for (i, r) in enumerate(reactions) + reactant_ids = [node_to_id[node] for node in r.reactants] + product_ids = [node_to_id[node] for node in r.products] + + for (j, id) in enumerate(reactant_ids) + source_id_batch[i, j] = id + source_padding_mask[i, j] = true + end + + for (j, id) in enumerate(product_ids) + target_id_batch[i, j] = id + target_padding_mask[i, j] = true + end +end + +println("\nPadded source/reactant ID batch:") +println(source_id_batch) + +println("\nSource padding mask:") +println(source_padding_mask) + +println("\nPadded target/product ID batch:") +println(target_id_batch) + +println("\nTarget padding mask:") +println(target_padding_mask) + + +#masked aggregation example + +println("\nMasked feature aggregation per hyperedge:") + +for r in reactions + src_mask = reactant_masks[r.id] + tgt_mask = product_masks[r.id] + + src_features = node_features[src_mask, :] + tgt_features = node_features[tgt_mask, :] + + src_mean = vec(mean(src_features, dims = 1)) + tgt_mean = vec(mean(tgt_features, dims = 1)) + + println("\ne", r.id, ": ", + join(r.reactants, " + "), + " -> ", + join(r.products, " + ")) + + println("Source mean feature: ", src_mean) + println("Target mean feature: ", tgt_mean) +end + + +#quick summary + +println("\nSummary:") +println("Number of nodes: ", length(nodes)) +println("Number of reactions/hyperedges: ", length(reactions)) +println("Signed incidence size: ", size(H)) +println("Maximum source size: ", max_source_size) +println("Maximum target size: ", max_target_size) \ No newline at end of file diff --git a/examples/shonali_prototypes/04_gather_operations.jl b/examples/shonali_prototypes/04_gather_operations.jl new file mode 100644 index 0000000..5bf0f80 --- /dev/null +++ b/examples/shonali_prototypes/04_gather_operations.jl @@ -0,0 +1,405 @@ +#gather operations +#trying out different ways of collecting node information from reactions + +using Statistics +using SimpleHypergraphs +using SimpleDirectedHypergraphs + + +#toy reaction network + +nodes = ["A", "B", "C", "D", "E", "F", "G", "H"] + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +reactions = [ + ( + id = 1, + reactants = ["A", "B"], + products = ["C", "D"] + ), + ( + id = 2, + reactants = ["C"], + products = ["E"] + ), + ( + id = 3, + reactants = ["D", "E"], + products = ["F"] + ), + ( + id = 4, + reactants = ["F"], + products = ["G", "H"] + ) +] + +println("\nToy reactions:") +for r in reactions + println("e", r.id, ": ", join(r.reactants, " + "), " -> ", join(r.products, " + ")) +end + + +#node features for testing + +#rows = nodes +#columns = feature values + +node_features = [ + 1.0 0.2 0.5; + 0.8 0.1 0.4; + 0.3 0.9 0.7; + 0.4 0.6 0.8; + 0.9 0.2 0.1; + 0.5 0.5 0.5; + 0.7 0.3 0.6; + 0.2 0.8 0.9 +] + +println("\nNode feature matrix:") +println(node_features) + +println("\nNode feature lookup:") +for i in 1:length(nodes) + println(nodes[i], " => ", node_features[i, :]) +end + + +#small helper functions + +function node_ids(node_names, node_to_id) + return [node_to_id[node] for node in node_names] +end + +function gather_features(feature_matrix, ids) + return feature_matrix[ids, :] +end + +function gather_reactant_features(reaction, node_to_id, feature_matrix) + ids = node_ids(reaction.reactants, node_to_id) + return gather_features(feature_matrix, ids) +end + +function gather_product_features(reaction, node_to_id, feature_matrix) + ids = node_ids(reaction.products, node_to_id) + return gather_features(feature_matrix, ids) +end + +function gather_active_features(reaction, node_to_id, feature_matrix) + active_nodes = vcat(reaction.reactants, reaction.products) + ids = node_ids(active_nodes, node_to_id) + return gather_features(feature_matrix, ids) +end + + +#gather features for each reaction + +println("\nGathered features per reaction:") + +for r in reactions + reactant_features = gather_reactant_features(r, node_to_id, node_features) + product_features = gather_product_features(r, node_to_id, node_features) + active_features = gather_active_features(r, node_to_id, node_features) + + println("\ne", r.id, ": ", join(r.reactants, " + "), " -> ", join(r.products, " + ")) + + println("Reactant/source features:") + println(reactant_features) + + println("Product/target features:") + println(product_features) + + println("All active node features:") + println(active_features) +end + + +#some simple aggregation ideas + +#mean, sum, max etc. + +function mean_feature(features) + return vec(mean(features, dims = 1)) +end + +function sum_feature(features) + return vec(sum(features, dims = 1)) +end + +function max_feature(features) + return vec(maximum(features, dims = 1)) +end + +println("\nAggregated hyperedge features:") + +for r in reactions + src = gather_reactant_features(r, node_to_id, node_features) + tgt = gather_product_features(r, node_to_id, node_features) + + src_mean = mean_feature(src) + tgt_mean = mean_feature(tgt) + + src_sum = sum_feature(src) + tgt_sum = sum_feature(tgt) + + direction_difference = tgt_mean .- src_mean + concatenated = vcat(src_mean, tgt_mean) + + println("\ne", r.id) + println("Source mean: ", src_mean) + println("Target mean: ", tgt_mean) + println("Source sum: ", src_sum) + println("Target sum: ", tgt_sum) + println("Target - source mean: ", direction_difference) + println("Concatenated source/target mean: ", concatenated) +end + + +#trying the same thing from the incidence matrix + +function build_signed_incidence(nodes, reactions, node_to_id) + H = zeros(Int, length(nodes), length(reactions)) + + for (j, r) in enumerate(reactions) + for reactant in r.reactants + H[node_to_id[reactant], j] = -1 + end + + for product in r.products + H[node_to_id[product], j] = 1 + end + end + + return H +end + +H = build_signed_incidence(nodes, reactions, node_to_id) + +println("\nSigned incidence matrix:") +println(H) + +println("\nGathering features using incidence columns:") + +for j in 1:size(H, 2) + source_ids = findall(H[:, j] .== -1) + target_ids = findall(H[:, j] .== 1) + active_ids = findall(H[:, j] .!= 0) + + source_features = gather_features(node_features, source_ids) + target_features = gather_features(node_features, target_ids) + active_features = gather_features(node_features, active_ids) + + println("\nHyperedge e", j) + println("Source IDs: ", source_ids) + println("Target IDs: ", target_ids) + println("Active IDs: ", active_ids) + + println("Source features:") + println(source_features) + + println("Target features:") + println(target_features) + + println("Active features:") + println(active_features) +end + + +#padding source/target ids so reactions can be processed together + +max_source_size = maximum(length(r.reactants) for r in reactions) +max_target_size = maximum(length(r.products) for r in reactions) + +source_id_batch = zeros(Int, length(reactions), max_source_size) +target_id_batch = zeros(Int, length(reactions), max_target_size) + +source_mask = falses(length(reactions), max_source_size) +target_mask = falses(length(reactions), max_target_size) + +for (i, r) in enumerate(reactions) + src_ids = node_ids(r.reactants, node_to_id) + tgt_ids = node_ids(r.products, node_to_id) + + for (j, id) in enumerate(src_ids) + source_id_batch[i, j] = id + source_mask[i, j] = true + end + + for (j, id) in enumerate(tgt_ids) + target_id_batch[i, j] = id + target_mask[i, j] = true + end +end + +println("\nPadded source ID batch:") +println(source_id_batch) + +println("\nSource mask:") +println(source_mask) + +println("\nPadded target ID batch:") +println(target_id_batch) + +println("\nTarget mask:") +println(target_mask) + + +#turn padded ids into padded feature batches + +num_reactions = length(reactions) +num_features = size(node_features, 2) + +source_feature_batch = zeros(Float64, num_reactions, max_source_size, num_features) +target_feature_batch = zeros(Float64, num_reactions, max_target_size, num_features) + +for i in 1:num_reactions + for j in 1:max_source_size + if source_mask[i, j] + node_id = source_id_batch[i, j] + source_feature_batch[i, j, :] = node_features[node_id, :] + end + end + + for j in 1:max_target_size + if target_mask[i, j] + node_id = target_id_batch[i, j] + target_feature_batch[i, j, :] = node_features[node_id, :] + end + end +end + +println("\nSource feature batch:") +println(source_feature_batch) + +println("\nTarget feature batch:") +println(target_feature_batch) + + +#masked aggregation on padded batches + +function masked_mean_feature_batch(feature_batch, mask) + batch_size = size(feature_batch, 1) + num_features = size(feature_batch, 3) + + output = zeros(Float64, batch_size, num_features) + + for i in 1:batch_size + valid_count = sum(mask[i, :]) + + if valid_count > 0 + for f in 1:num_features + total = 0.0 + + for j in 1:size(feature_batch, 2) + if mask[i, j] + total += feature_batch[i, j, f] + end + end + + output[i, f] = total / valid_count + end + end + end + + return output +end + +source_mean_batch = masked_mean_feature_batch(source_feature_batch, source_mask) +target_mean_batch = masked_mean_feature_batch(target_feature_batch, target_mask) + +println("\nSource mean feature batch:") +println(source_mean_batch) + +println("\nTarget mean feature batch:") +println(target_mean_batch) + + +#building a simple reaction/hyperedge representation + +hyperedge_embeddings = zeros(Float64, num_reactions, num_features * 3) + +for i in 1:num_reactions + src_mean = source_mean_batch[i, :] + tgt_mean = target_mean_batch[i, :] + diff = tgt_mean .- src_mean + + hyperedge_embeddings[i, :] = vcat(src_mean, tgt_mean, diff) +end + +println("\nToy hyperedge embeddings:") +println(hyperedge_embeddings) + + +#checking which reactions each node belongs to + +node_to_hyperedges = Dict(node => Int[] for node in nodes) + +for r in reactions + involved_nodes = vcat(r.reactants, r.products) + + for node in involved_nodes + push!(node_to_hyperedges[node], r.id) + end +end + +println("\nNode-to-hyperedge gather mapping:") + +for node in nodes + println(node, " participates in hyperedges ", node_to_hyperedges[node]) +end + + +#gather reaction embeddings for each node + +println("\nGather hyperedge embeddings per node:") + +for node in nodes + connected_edges = node_to_hyperedges[node] + + if isempty(connected_edges) + println(node, " has no connected hyperedges") + else + gathered = hyperedge_embeddings[connected_edges, :] + println("\nNode ", node) + println("Connected hyperedges: ", connected_edges) + println("Gathered hyperedge embeddings:") + println(gathered) + end +end + + +#very simple node update experiment + +#take the average of connected reaction embeddings + +node_updated_features = Dict{String, Vector{Float64}}() + +for node in nodes + connected_edges = node_to_hyperedges[node] + + if isempty(connected_edges) + node_updated_features[node] = zeros(Float64, size(hyperedge_embeddings, 2)) + else + gathered = hyperedge_embeddings[connected_edges, :] + node_updated_features[node] = vec(mean(gathered, dims = 1)) + end +end + +println("\nToy node updates from gathered hyperedge embeddings:") + +for node in nodes + println(node, " updated representation: ", node_updated_features[node]) +end + + +#quick summary + +println("\nSummary:") +println("Number of nodes: ", length(nodes)) +println("Number of hyperedges/reactions: ", length(reactions)) +println("Node feature size: ", size(node_features)) +println("Source feature batch size: ", size(source_feature_batch)) +println("Target feature batch size: ", size(target_feature_batch)) +println("Hyperedge embedding size: ", size(hyperedge_embeddings)) \ No newline at end of file diff --git a/examples/shonali_prototypes/05_basic_transformations.jl b/examples/shonali_prototypes/05_basic_transformations.jl new file mode 100644 index 0000000..bdad3c8 --- /dev/null +++ b/examples/shonali_prototypes/05_basic_transformations.jl @@ -0,0 +1,410 @@ +#basic transformations +#just trying to turn toy reaction strings into the pieces I might need later + +using Statistics +using SimpleHypergraphs +using SimpleDirectedHypergraphs + + +#toy reactions to play around with + +raw_reactions = [ + "A + B -> C + D", + "C -> E", + "D + E -> F", + "F -> G + H", + "H + A -> I" +] + +println("\nRaw reaction strings:") +for r in raw_reactions + println(r) +end + + +#parsing the reaction strings into reactants and products + +function clean_species_name(x) + return strip(x) +end + +function parse_side(side_string) + species = split(side_string, "+") + return [clean_species_name(s) for s in species] +end + +function parse_reaction_string(reaction_string, reaction_id) + sides = split(reaction_string, "->") + + reactants = parse_side(sides[1]) + products = parse_side(sides[2]) + + return ( + id = reaction_id, + raw = reaction_string, + reactants = reactants, + products = products + ) +end + +reactions = [ + parse_reaction_string(raw_reactions[i], i) + for i in 1:length(raw_reactions) +] + +println("\nParsed reactions:") +for r in reactions + println("e", r.id, ": ", join(r.reactants, " + "), " -> ", join(r.products, " + ")) +end + + +#get all unique species/nodes + +function extract_nodes(reactions) + node_set = Set{String}() + + for r in reactions + for node in vcat(r.reactants, r.products) + push!(node_set, node) + end + end + + return sort(collect(node_set)) +end + +nodes = extract_nodes(reactions) + +println("\nExtracted node/species set:") +println(nodes) + + +#assigning integer ids to species + +function build_node_mappings(nodes) + node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) + id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + return node_to_id, id_to_node +end + +node_to_id, id_to_node = build_node_mappings(nodes) + +println("\nNode to ID mapping:") +println(node_to_id) + +println("\nID to node mapping:") +println(id_to_node) + + +#converting each reaction to an id-based hyperedge + +function reaction_to_id_hyperedge(reaction, node_to_id) + reactant_ids = [node_to_id[x] for x in reaction.reactants] + product_ids = [node_to_id[x] for x in reaction.products] + + return ( + id = reaction.id, + reactants = reaction.reactants, + products = reaction.products, + reactant_ids = reactant_ids, + product_ids = product_ids + ) +end + +id_hyperedges = [ + reaction_to_id_hyperedge(r, node_to_id) + for r in reactions +] + +println("\nID-based hyperedge representation:") +for e in id_hyperedges + println("\ne", e.id) + println("Reactants: ", e.reactants, " => ", e.reactant_ids) + println("Products: ", e.products, " => ", e.product_ids) +end + + +#building a signed incidence matrix +#using -1 for reactants and +1 for products + +function build_signed_incidence(nodes, id_hyperedges) + H = zeros(Int, length(nodes), length(id_hyperedges)) + + for (j, e) in enumerate(id_hyperedges) + for id in e.reactant_ids + H[id, j] = -1 + end + + for id in e.product_ids + H[id, j] = 1 + end + end + + return H +end + +H_signed = build_signed_incidence(nodes, id_hyperedges) + +println("\nSigned incidence matrix:") +println("Rows = nodes/species") +println("Columns = reactions/hyperedges") +println(H_signed) + + +#splitting the signed matrix into source, target, and membership versions + +function source_matrix(H_signed) + return Int.(H_signed .== -1) +end + +function target_matrix(H_signed) + return Int.(H_signed .== 1) +end + +function unsigned_matrix(H_signed) + return abs.(H_signed) +end + +H_source = source_matrix(H_signed) +H_target = target_matrix(H_signed) +H_unsigned = unsigned_matrix(H_signed) + +println("\nSource/reactant matrix:") +println(H_source) + +println("\nTarget/product matrix:") +println(H_target) + +println("\nUnsigned membership matrix:") +println(H_unsigned) + + +#another representation: node-hyperedge-role triples + +function incidence_to_edge_list(H_signed) + edge_list = [] + + for j in 1:size(H_signed, 2) + for i in 1:size(H_signed, 1) + role = H_signed[i, j] + + if role != 0 + push!( + edge_list, + ( + node_id = i, + hyperedge_id = j, + role = role + ) + ) + end + end + end + + return edge_list +end + +edge_list = incidence_to_edge_list(H_signed) + +println("\nEdge-list representation:") +println("(role = -1 source/reactant, role = +1 target/product)") +for item in edge_list + println(item) +end + + +#quick reaction size checks + +function hyperedge_sizes(id_hyperedges) + source_sizes = [length(e.reactant_ids) for e in id_hyperedges] + target_sizes = [length(e.product_ids) for e in id_hyperedges] + total_sizes = source_sizes .+ target_sizes + + return source_sizes, target_sizes, total_sizes +end + +source_sizes, target_sizes, total_sizes = hyperedge_sizes(id_hyperedges) + +println("\nHyperedge size statistics:") +for i in 1:length(id_hyperedges) + println( + "e", i, + " | source size = ", source_sizes[i], + " | target size = ", target_sizes[i], + " | total size = ", total_sizes[i] + ) +end + +println("\nAverage source size: ", mean(source_sizes)) +println("Average target size: ", mean(target_sizes)) +println("Average total hyperedge size: ", mean(total_sizes)) + + +#counting how often each node appears anywhere + +function node_participation_counts(H_unsigned, nodes) + counts = vec(sum(H_unsigned, dims = 2)) + + result = Dict{String, Int}() + + for i in 1:length(nodes) + result[nodes[i]] = counts[i] + end + + return result +end + +participation = node_participation_counts(H_unsigned, nodes) + +println("\nNode participation counts:") +for node in nodes + println(node, " appears in ", participation[node], " hyperedge(s)") +end + + +#source vs target participation counts + +function source_target_degrees(H_source, H_target, nodes) + source_counts = vec(sum(H_source, dims = 2)) + target_counts = vec(sum(H_target, dims = 2)) + + source_degree = Dict{String, Int}() + target_degree = Dict{String, Int}() + + for i in 1:length(nodes) + source_degree[nodes[i]] = source_counts[i] + target_degree[nodes[i]] = target_counts[i] + end + + return source_degree, target_degree +end + +source_degree, target_degree = + source_target_degrees(H_source, H_target, nodes) + +println("\nSource/target degree-style counts:") +for node in nodes + println( + node, + " | source count = ", source_degree[node], + " | target count = ", target_degree[node] + ) +end + + +#simple structural node features for now + +node_feature_names = [ + "source_count", + "target_count", + "participation_count" +] + +node_features = zeros(Float64, length(nodes), length(node_feature_names)) + +for i in 1:length(nodes) + node = nodes[i] + + node_features[i, 1] = source_degree[node] + node_features[i, 2] = target_degree[node] + node_features[i, 3] = participation[node] +end + +println("\nNode feature names:") +println(node_feature_names) + +println("\nNode feature matrix:") +println(node_features) + + +#same idea for reaction/hyperedge features + +hyperedge_feature_names = [ + "source_size", + "target_size", + "total_size" +] + +hyperedge_features = zeros(Float64, length(id_hyperedges), length(hyperedge_feature_names)) + +for i in 1:length(id_hyperedges) + hyperedge_features[i, 1] = source_sizes[i] + hyperedge_features[i, 2] = target_sizes[i] + hyperedge_features[i, 3] = total_sizes[i] +end + +println("\nHyperedge feature names:") +println(hyperedge_feature_names) + +println("\nHyperedge feature matrix:") +println(hyperedge_features) + + +#checking source, target, and active masks + +source_masks = Dict{Int, Vector{Bool}}() +target_masks = Dict{Int, Vector{Bool}}() +active_masks = Dict{Int, Vector{Bool}}() + +for e in id_hyperedges + source_mask = falses(length(nodes)) + target_mask = falses(length(nodes)) + + for id in e.reactant_ids + source_mask[id] = true + end + + for id in e.product_ids + target_mask[id] = true + end + + active_mask = source_mask .| target_mask + + source_masks[e.id] = source_mask + target_masks[e.id] = target_mask + active_masks[e.id] = active_mask +end + +println("\nMasks from transformed representation:") +for e in id_hyperedges + println("\ne", e.id) + println("Source mask: ", source_masks[e.id]) + println("Target mask: ", target_masks[e.id]) + println("Active mask: ", active_masks[e.id]) +end + + +#storing everything together so it is easier to inspect later + +transformed_data = Dict( + "nodes" => nodes, + "node_to_id" => node_to_id, + "id_to_node" => id_to_node, + "raw_reactions" => raw_reactions, + "parsed_reactions" => reactions, + "id_hyperedges" => id_hyperedges, + "signed_incidence" => H_signed, + "source_matrix" => H_source, + "target_matrix" => H_target, + "unsigned_matrix" => H_unsigned, + "edge_list" => edge_list, + "node_features" => node_features, + "hyperedge_features" => hyperedge_features, + "source_masks" => source_masks, + "target_masks" => target_masks, + "active_masks" => active_masks +) + +println("\nTransformed data keys:") +println(collect(keys(transformed_data))) + + +#quick summary + +println("\nSummary:") +println("Raw reactions: ", length(raw_reactions)) +println("Nodes/species: ", length(nodes)) +println("Directed hyperedges/reactions: ", length(id_hyperedges)) +println("Signed incidence size: ", size(H_signed)) +println("Node feature matrix size: ", size(node_features)) +println("Hyperedge feature matrix size: ", size(hyperedge_features)) +println("Edge-list length: ", length(edge_list)) \ No newline at end of file diff --git a/examples/shonali_prototypes/06_preprocessing_pipeline.jl b/examples/shonali_prototypes/06_preprocessing_pipeline.jl new file mode 100644 index 0000000..18b9491 --- /dev/null +++ b/examples/shonali_prototypes/06_preprocessing_pipeline.jl @@ -0,0 +1,280 @@ +#preprocessing pipeline +#trying to connect the small pieces from the earlier files: +#parsing reactions, ids, incidence matrix, masks, and feature prep + +using Statistics +using SimpleHypergraphs +using SimpleDirectedHypergraphs + + +#toy style reactions +raw_reactions = [ + "A + B -> C", + "C + D -> E", + "E -> F + G", + "G -> H", + "H + A -> I" +] + +println("\nraw reactions") +for r in raw_reactions + println(r) +end + +#small helper to split one side of a reaction +function parse_side(x) + return [strip(s) for s in split(x, "+")] +end + +function parse_reaction(r, id) + sides = split(r, "->") + + reactants = parse_side(sides[1]) + products = parse_side(sides[2]) + + return ( + id = id, + raw = r, + reactants = reactants, + products = products + ) +end + +reactions = [parse_reaction(raw_reactions[i], i) for i in 1:length(raw_reactions)] + +println("\nparsed reactions") +for r in reactions + println("r", r.id, ": ", join(r.reactants, " + "), " -> ", join(r.products, " + ")) +end + +#collectig all species/nodes +species_set = Set{String}() + +for r in reactions + for s in vcat(r.reactants, r.products) + push!(species_set, s) + end +end + +nodes = sort(collect(species_set)) + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nnodes/species") +println(nodes) + +println("\nnode ids") +println(node_to_id) + +#converting reaction names into ids +id_reactions = [] + +for r in reactions + reactant_ids = [node_to_id[x] for x in r.reactants] + product_ids = [node_to_id[x] for x in r.products] + + push!( + id_reactions, + ( + id = r.id, + reactants = r.reactants, + products = r.products, + reactant_ids = reactant_ids, + product_ids = product_ids + ) + ) +end + +println("\nid-based reactions") +for r in id_reactions + println("r", r.id) + println(" reactants: ", r.reactants, " -> ", r.reactant_ids) + println(" products: ", r.products, " -> ", r.product_ids) +end + +#build signed incidence matrix +#-1 = reactant/source side +#+1 = product/target side +#0 = not involved + +H = zeros(Int, length(nodes), length(id_reactions)) + +for (j, r) in enumerate(id_reactions) + for id in r.reactant_ids + H[id, j] = -1 + end + + for id in r.product_ids + H[id, j] = 1 + end +end + +println("\nsigned incidence matrix") +println(H) + +H_source = Int.(H .== -1) +H_target = Int.(H .== 1) +H_membership = abs.(H) + +println("\nsource matrix") +println(H_source) + +println("\ntarget matrix") +println(H_target) + +println("\nunsigned membership matrix") +println(H_membership) + +#masks from the incidence matrix +source_masks = Dict{Int, Vector{Bool}}() +target_masks = Dict{Int, Vector{Bool}}() +active_masks = Dict{Int, Vector{Bool}}() + +for j in 1:size(H, 2) + source_masks[j] = H[:, j] .== -1 + target_masks[j] = H[:, j] .== 1 + active_masks[j] = H[:, j] .!= 0 +end + +println("\nchecking masks") +for j in 1:length(id_reactions) + println("\nr", j) + println("source nodes: ", nodes[source_masks[j]]) + println("target nodes: ", nodes[target_masks[j]]) + println("active nodes: ", nodes[active_masks[j]]) +end + +#simple structural features for nodes +#this is not chemical descriptor data yet, just useful toy features + +source_count = vec(sum(H_source, dims = 2)) +target_count = vec(sum(H_target, dims = 2)) +participation_count = vec(sum(H_membership, dims = 2)) + +node_features = hcat(source_count, target_count, participation_count) + +println("\nnode features") +println("columns = source_count, target_count, participation_count") +println(node_features) + +#hyperedge/reaction features +source_size = vec(sum(H_source, dims = 1)) +target_size = vec(sum(H_target, dims = 1)) +total_size = vec(sum(H_membership, dims = 1)) + +hyperedge_features = hcat(source_size, target_size, total_size) + +println("\nhyperedge features") +println("columns = source_size, target_size, total_size") +println(hyperedge_features) + +#edge list version +#useful because some ML/message passing code works with index lists + +edge_list = [] + +for j in 1:size(H, 2) + for i in 1:size(H, 1) + if H[i, j] != 0 + push!( + edge_list, + ( + node_id = i, + hyperedge_id = j, + role = H[i, j] + ) + ) + end + end +end + +println("\nedge list") +for e in edge_list + println(e) +end + +#gather prep: for each reaction, storing source and target ids separately +source_id_lists = Dict{Int, Vector{Int}}() +target_id_lists = Dict{Int, Vector{Int}}() + +for r in id_reactions + source_id_lists[r.id] = r.reactant_ids + target_id_lists[r.id] = r.product_ids +end + +println("\nsource/target id lists") +for r in id_reactions + println("r", r.id, " source ids = ", source_id_lists[r.id], + " target ids = ", target_id_lists[r.id]) +end + +#small batching prep +#since reactions have variable numbers of reactants/products, +#pad the ids and keep masks so padded zeros are ignored + +max_source_len = maximum(length(r.reactant_ids) for r in id_reactions) +max_target_len = maximum(length(r.product_ids) for r in id_reactions) + +source_id_batch = zeros(Int, length(id_reactions), max_source_len) +target_id_batch = zeros(Int, length(id_reactions), max_target_len) + +source_batch_mask = falses(length(id_reactions), max_source_len) +target_batch_mask = falses(length(id_reactions), max_target_len) + +for (i, r) in enumerate(id_reactions) + for (j, id) in enumerate(r.reactant_ids) + source_id_batch[i, j] = id + source_batch_mask[i, j] = true + end + + for (j, id) in enumerate(r.product_ids) + target_id_batch[i, j] = id + target_batch_mask[i, j] = true + end +end + +println("\npadded source id batch") +println(source_id_batch) + +println("\nsource batch mask") +println(source_batch_mask) + +println("\npadded target id batch") +println(target_id_batch) + +println("\ntarget batch mask") +println(target_batch_mask) + +#packaging the output like a preprocessing function might return +processed = Dict( + "nodes" => nodes, + "node_to_id" => node_to_id, + "id_to_node" => id_to_node, + "reactions" => id_reactions, + "signed_incidence" => H, + "source_matrix" => H_source, + "target_matrix" => H_target, + "membership_matrix" => H_membership, + "node_features" => node_features, + "hyperedge_features" => hyperedge_features, + "edge_list" => edge_list, + "source_masks" => source_masks, + "target_masks" => target_masks, + "active_masks" => active_masks, + "source_id_batch" => source_id_batch, + "target_id_batch" => target_id_batch, + "source_batch_mask" => source_batch_mask, + "target_batch_mask" => target_batch_mask +) + +println("\nprocessed object keys") +println(collect(keys(processed))) + +println("\nsummary") +println("number of nodes: ", length(nodes)) +println("number of reactions/hyperedges: ", length(id_reactions)) +println("incidence matrix size: ", size(H)) +println("node feature matrix size: ", size(node_features)) +println("hyperedge feature matrix size: ", size(hyperedge_features)) + diff --git a/examples/shonali_prototypes/07_toy_ml_feature_model.jl b/examples/shonali_prototypes/07_toy_ml_feature_model.jl new file mode 100644 index 0000000..97e9677 --- /dev/null +++ b/examples/shonali_prototypes/07_toy_ml_feature_model.jl @@ -0,0 +1,219 @@ +#toy ml feature model +#small experiment to see how reaction/hyperedge features could become model input + +using Statistics +using LinearAlgebra +using SimpleHypergraphs +using SimpleDirectedHypergraphs + + +#toy reaction network + +nodes = ["A", "B", "C", "D", "E", "F", "G", "H", "I"] + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) + +reactions = [ + (id = 1, reactants = ["A", "B"], products = ["C"]), + (id = 2, reactants = ["C", "D"], products = ["E"]), + (id = 3, reactants = ["E"], products = ["F", "G"]), + (id = 4, reactants = ["G"], products = ["H"]), + (id = 5, reactants = ["H", "A"], products = ["I"]), + (id = 6, reactants = ["B"], products = ["D"]), + (id = 7, reactants = ["D", "F"], products = ["G"]), + (id = 8, reactants = ["I"], products = ["A", "E"]) +] + +println("\nToy reactions:") +for r in reactions + println("e", r.id, ": ", join(r.reactants, " + "), " -> ", join(r.products, " + ")) +end + + +#building incidence matrix first + +function build_signed_incidence(nodes, reactions, node_to_id) + H = zeros(Int, length(nodes), length(reactions)) + + for (j, r) in enumerate(reactions) + for reactant in r.reactants + H[node_to_id[reactant], j] = -1 + end + + for product in r.products + H[node_to_id[product], j] = 1 + end + end + + return H +end + +H = build_signed_incidence(nodes, reactions, node_to_id) + +println("\nSigned incidence matrix:") +println(H) + +H_source = Int.(H .== -1) +H_target = Int.(H .== 1) +H_membership = abs.(H) + + +#node features from structure + +source_count = vec(sum(H_source, dims = 2)) +target_count = vec(sum(H_target, dims = 2)) +participation_count = vec(sum(H_membership, dims = 2)) + +node_features = hcat(source_count, target_count, participation_count) + +println("\nNode features:") +println("columns = source_count, target_count, participation_count") +println(node_features) + + +#reaction/hyperedge features + +source_size = vec(sum(H_source, dims = 1)) +target_size = vec(sum(H_target, dims = 1)) +total_size = vec(sum(H_membership, dims = 1)) + +hyperedge_features = hcat(source_size, target_size, total_size) + +println("\nHyperedge features:") +println("columns = source_size, target_size, total_size") +println(hyperedge_features) + + +#toy label +#for now, label = 1 if the reaction has more than one reactant + +labels = Int.(source_size .> 1) + +println("\nToy labels:") +println(labels) + + +#normalising features a little bit + +function normalize_columns(X) + X_norm = zeros(Float64, size(X)) + + for j in 1:size(X, 2) + col = X[:, j] + μ = mean(col) + σ = std(col) + + if σ == 0 + X_norm[:, j] .= 0.0 + else + X_norm[:, j] = (col .- μ) ./ σ + end + end + + return X_norm +end + +X = normalize_columns(Float64.(hyperedge_features)) +y = labels + +println("\nNormalised hyperedge features:") +println(X) + + +#simple score-based classifier +#not a full ML model yet, just a first toy experiment + +weights = [1.0, -0.5, 0.25] +bias = 0.0 + +function sigmoid(z) + return 1.0 / (1.0 + exp(-z)) +end + +function predict_scores(X, weights, bias) + scores = zeros(Float64, size(X, 1)) + + for i in 1:size(X, 1) + scores[i] = dot(X[i, :], weights) + bias + end + + return scores +end + +scores = predict_scores(X, weights, bias) +probabilities = [sigmoid(s) for s in scores] +predictions = Int.(probabilities .>= 0.5) + +println("\nScores:") +println(scores) + +println("\nProbabilities:") +println(probabilities) + +println("\nPredictions:") +println(predictions) + +println("\nActual labels:") +println(y) + + +#small accuracy check + +function accuracy(y_true, y_pred) + correct = 0 + + for i in 1:length(y_true) + if y_true[i] == y_pred[i] + correct += 1 + end + end + + return correct / length(y_true) +end + +acc = accuracy(y, predictions) + +println("\nToy accuracy:") +println(acc) + + +#trying a few different weight settings manually + +weight_trials = [ + [1.0, -0.5, 0.25], + [0.8, -0.3, 0.1], + [1.5, -1.0, 0.2], + [0.5, 0.0, 0.5] +] + +println("\nTrying a few simple weight settings:") + +for (trial_id, w) in enumerate(weight_trials) + trial_scores = predict_scores(X, w, bias) + trial_probs = [sigmoid(s) for s in trial_scores] + trial_preds = Int.(trial_probs .>= 0.5) + trial_acc = accuracy(y, trial_preds) + + println("\ntrial ", trial_id) + println("weights: ", w) + println("predictions: ", trial_preds) + println("accuracy: ", trial_acc) +end + + +#small interpretation + +println("\nInterpretation:") +println("This is only a tiny ML-style experiment.") +println("The main point is to test how hyperedge features can be created from a directed reaction network.") +println("A later version could replace this with Lux.jl layers or a proper HGNN model.") + + +#quick summary + +println("\nSummary:") +println("Number of nodes: ", length(nodes)) +println("Number of reactions/hyperedges: ", length(reactions)) +println("Node feature matrix size: ", size(node_features)) +println("Hyperedge feature matrix size: ", size(hyperedge_features)) +println("Number of labels: ", length(labels)) \ No newline at end of file diff --git a/examples/shonali_prototypes/08_crn_parser_to_dihypergraph.jl b/examples/shonali_prototypes/08_crn_parser_to_dihypergraph.jl new file mode 100644 index 0000000..8216f2c --- /dev/null +++ b/examples/shonali_prototypes/08_crn_parser_to_dihypergraph.jl @@ -0,0 +1,399 @@ +#08_crn_parser_to_dihypergraph.jl + +#this file is for converting CRN reaction strings into a directed hypergraph style data structure. + +#the main idea is: +#species = nodes +#reactions = directed hyperedges +#reactants = source/tail side +#products = target/head side + +#i am starting with a small toy CRN first, because it makes the parsing logic easier to check before applying the same workflow to CRNs from literature. + + +#1.small CRN example + +raw_reactions = [ + "A + B -> C", + "C -> D + E", + "E + F -> G", + "G -> H", + "H + A -> I" +] + +println("\nRaw CRN reactions:") + +for reaction in raw_reactions + println(reaction) +end + +#2.cleaning helper + +#this just removes extra spaces from species names +#for example, " A " becomes "A" + +function clean_species_name(name) + return strip(name) +end + + +#3.parsing one reaction + + +#this function takes a reaction string like: +#A + B -> C + D +#and converts it into: +#reactants = ["A", "B"] +#products = ["C", "D"] + +function parse_reaction(reaction_string) + if !occursin("->", reaction_string) + error("Reaction is missing -> : $reaction_string") + end + + left_side, right_side = split(reaction_string, "->") + + reactants = clean_species_name.(split(strip(left_side), "+")) + products = clean_species_name.(split(strip(right_side), "+")) + + reactants = filter(x -> x != "", reactants) + products = filter(x -> x != "", products) + + return ( + raw = reaction_string, + reactants = reactants, + products = products + ) +end + + +#4.parsing the full CRN + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nParsed reactions:") + +for (i, reaction) in enumerate(parsed_reactions) + println( + "r", i, ": ", + join(reaction.reactants, " + "), + " -> ", + join(reaction.products, " + ") + ) +end + + +#5.extracting species / nodes + +#here i collect every species that appears either as a reactant or product +#these become the nodes of the directed hypergraph + +species_set = Set{String}() + +for reaction in parsed_reactions + for species in reaction.reactants + push!(species_set, species) + end + + for species in reaction.products + push!(species_set, species) + end +end + +nodes = sort(collect(species_set)) + +println("\nExtracted species / nodes:") +println(nodes) + + +#6.creating node ID mappings + + +#it is easier to use integer IDs for matrix construction and ML code +#so each species gets an integer ID + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nNode to ID mapping:") +println(node_to_id) + +println("\nID to node mapping:") +println(id_to_node) + + +#7.converting reactions into directed hyperedges + +#each reaction becomes one directed hyperedge +#the source side contains reactant IDs +#the target side contains product IDs + +directed_hyperedges = [] + +for (i, reaction) in enumerate(parsed_reactions) + source_ids = [node_to_id[x] for x in reaction.reactants] + target_ids = [node_to_id[x] for x in reaction.products] + + push!( + directed_hyperedges, + ( + id = i, + raw = reaction.raw, + reactants = reaction.reactants, + products = reaction.products, + source_ids = source_ids, + target_ids = target_ids + ) + ) +end + +println("\nDirected hyperedge representation:") + +for edge in directed_hyperedges + println("\nr", edge.id) + println(" raw reaction: ", edge.raw) + println(" source/reactants: ", edge.reactants, " -> ", edge.source_ids) + println(" target/products: ", edge.products, " -> ", edge.target_ids) +end + + +#8.signed incidence matrix + + +#this matrix stores direction information + +#-1 means the species is on the source/reactant side +#+1 means the species is on the target/product side +#0 means the species is not involved in that reaction + +function build_signed_incidence(nodes, directed_hyperedges) + H = zeros(Int, length(nodes), length(directed_hyperedges)) + + for edge in directed_hyperedges + for node_id in edge.source_ids + H[node_id, edge.id] = -1 + end + + for node_id in edge.target_ids + H[node_id, edge.id] = 1 + end + end + + return H +end + +signed_incidence = build_signed_incidence(nodes, directed_hyperedges) + +println("\nSigned incidence matrix:") +println("Rows = species / nodes") +println("Columns = reactions / directed hyperedges") +println(signed_incidence) + + +#9.source and target incidence matrices + + +#these are split versions of the signed incidence matrix +#they make it easier to work separately with reactants and products + +function build_source_matrix(nodes, directed_hyperedges) + S = zeros(Int, length(nodes), length(directed_hyperedges)) + + for edge in directed_hyperedges + for node_id in edge.source_ids + S[node_id, edge.id] = 1 + end + end + + return S +end + +function build_target_matrix(nodes, directed_hyperedges) + T = zeros(Int, length(nodes), length(directed_hyperedges)) + + for edge in directed_hyperedges + for node_id in edge.target_ids + T[node_id, edge.id] = 1 + end + end + + return T +end + +source_matrix = build_source_matrix(nodes, directed_hyperedges) +target_matrix = build_target_matrix(nodes, directed_hyperedges) + +println("\nSource / reactant matrix:") +println(source_matrix) + +println("\nTarget / product matrix:") +println(target_matrix) + + + +#10.unsigned membership matrix + + +#this ignores direction and only checks whether a species participates in a reaction. + +membership_matrix = abs.(signed_incidence) + +println("\nUnsigned membership matrix:") +println(membership_matrix) + + +#11.reconstructing reactions from the matrix + + +#this is a sanity check +#if the incidence matrix is correct, the original reactions should be recoverable from it + +function reconstruct_reaction(H, reaction_index, id_to_node) + reactants = String[] + products = String[] + + for i in 1:size(H, 1) + if H[i, reaction_index] == -1 + push!(reactants, id_to_node[i]) + elseif H[i, reaction_index] == 1 + push!(products, id_to_node[i]) + end + end + + return reactants, products +end + +println("\nReconstructed reactions from signed incidence matrix:") + +for j in 1:size(signed_incidence, 2) + reactants, products = reconstruct_reaction(signed_incidence, j, id_to_node) + + println( + "r", j, ": ", + join(reactants, " + "), + " -> ", + join(products, " + ") + ) +end + + +#12.edge-list representation + + +#this is another way of storing the same information +#it may be useful later for batching and indexing + +#role = -1 means source/reactant +#role = +1 means target/product + +edge_list = [] + +for edge in directed_hyperedges + for node_id in edge.source_ids + push!( + edge_list, + ( + node_id = node_id, + hyperedge_id = edge.id, + role = -1 + ) + ) + end + + for node_id in edge.target_ids + push!( + edge_list, + ( + node_id = node_id, + hyperedge_id = edge.id, + role = 1 + ) + ) + end +end + +println("\nEdge-list representation:") +println("(role = -1 source/reactant, role = +1 target/product)") + +for item in edge_list + println(item) +end + + +#13.simple node features + +#these are basic structural features +#they are not chemistry-specific yet, but they are useful as a first ML input + +source_count = vec(sum(source_matrix, dims = 2)) +target_count = vec(sum(target_matrix, dims = 2)) +participation_count = vec(sum(membership_matrix, dims = 2)) + +node_features = hcat(source_count, target_count, participation_count) + +println("\nNode feature names:") +println(["source_count", "target_count", "participation_count"]) + +println("\nNode feature matrix:") +println(node_features) + + +#14.simple hyperedge features + + +#these describe the size of each reaction. +#again, this is a simple starting point before using richer information. + +source_size = vec(sum(source_matrix, dims = 1)) +target_size = vec(sum(target_matrix, dims = 1)) +total_size = vec(sum(membership_matrix, dims = 1)) + +hyperedge_features = hcat(source_size, target_size, total_size) + +println("\nHyperedge feature names:") +println(["source_size", "target_size", "total_size"]) + +println("\nHyperedge feature matrix:") +println(hyperedge_features) + + +#15.storing the processed CRN data together + + +#this keeps all the processed pieces in one object +#later this can be adapted to use package-specific DirectedHypergraph or HGNNHypergraph objects + +crn_dihypergraph_data = ( + raw_reactions = raw_reactions, + parsed_reactions = parsed_reactions, + nodes = nodes, + node_to_id = node_to_id, + id_to_node = id_to_node, + directed_hyperedges = directed_hyperedges, + signed_incidence = signed_incidence, + source_matrix = source_matrix, + target_matrix = target_matrix, + membership_matrix = membership_matrix, + edge_list = edge_list, + node_features = node_features, + hyperedge_features = hyperedge_features +) + +println("\nProcessed CRN directed hypergraph data keys:") +println(keys(crn_dihypergraph_data)) + + +#16.printing the summary + + +println("\nSummary:") +println("Number of raw reactions: ", length(raw_reactions)) +println("Number of species / nodes: ", length(nodes)) +println("Number of directed hyperedges / reactions: ", length(directed_hyperedges)) +println("Signed incidence matrix size: ", size(signed_incidence)) +println("Source matrix size: ", size(source_matrix)) +println("Target matrix size: ", size(target_matrix)) +println("Membership matrix size: ", size(membership_matrix)) +println("Node feature matrix size: ", size(node_features)) +println("Hyperedge feature matrix size: ", size(hyperedge_features)) +println("Edge-list length: ", length(edge_list)) \ No newline at end of file diff --git a/examples/shonali_prototypes/09_package_directed_hypergraph.jl b/examples/shonali_prototypes/09_package_directed_hypergraph.jl new file mode 100644 index 0000000..436599f --- /dev/null +++ b/examples/shonali_prototypes/09_package_directed_hypergraph.jl @@ -0,0 +1,320 @@ +#09_package_directed_hypergraph.jl + +#this file connects the crn parsing work to the actual package that evan mentioned: SimpleDirectedHypergraphs.jl + +#in file 8, i built my own directed hypergraph-style structures using tuples and matrices +#here i am taking the same idea and converting the source/target matrices into a DirectedHypergraph object from the package. + +using SimpleDirectedHypergraphs +using Statistics + +raw_reactions = [ + "A + B -> C", + "C -> D + E", + "E + F -> G", + "G -> H", + "H + A -> I" +] + +println("\nraw crn reactions:") + +for reaction in raw_reactions + println(reaction) +end + +#small helper for removing extra spaces + +function clean_species_name(name) + return strip(name) +end + +#parsing one reaction string into reactants and products + +function parse_reaction(reaction_string) + if !occursin("->", reaction_string) + error("reaction is missing -> : $reaction_string") + end + + left_side, right_side = split(reaction_string, "->") + + reactants = clean_species_name.(split(strip(left_side), "+")) + products = clean_species_name.(split(strip(right_side), "+")) + + reactants = filter(x -> x != "", reactants) + products = filter(x -> x != "", products) + + return ( + raw = reaction_string, + reactants = reactants, + products = products + ) +end + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nparsed reactions:") + +for (i, reaction) in enumerate(parsed_reactions) + println( + "r", i, ": ", + join(reaction.reactants, " + "), + " -> ", + join(reaction.products, " + ") + ) +end + +#collecting all species in the crn +#these will become the nodes of the directed hypergraph + +species_set = Set{String}() + +for reaction in parsed_reactions + for species in reaction.reactants + push!(species_set, species) + end + + for species in reaction.products + push!(species_set, species) + end +end + +nodes = sort(collect(species_set)) + +println("\nextracted species / nodes:") +println(nodes) + +#mapping species names to integer ids +#the package constructor works with matrices, so integer indexing is useful + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nnode to id mapping:") +println(node_to_id) + +println("\nid to node mapping:") +println(id_to_node) + +#converting parsed reactions into source and target id lists + +directed_hyperedges = [] + +for (i, reaction) in enumerate(parsed_reactions) + source_ids = [node_to_id[x] for x in reaction.reactants] + target_ids = [node_to_id[x] for x in reaction.products] + + push!( + directed_hyperedges, + ( + id = i, + raw = reaction.raw, + reactants = reaction.reactants, + products = reaction.products, + source_ids = source_ids, + target_ids = target_ids + ) + ) +end + +println("\nid-based directed hyperedges:") + +for edge in directed_hyperedges + println("\nr", edge.id) + println(" source/reactants: ", edge.reactants, " -> ", edge.source_ids) + println(" target/products: ", edge.products, " -> ", edge.target_ids) +end + +#building normal integer source and target matrices first +#rows = species +#columns = reactions + +source_matrix = zeros(Int, length(nodes), length(directed_hyperedges)) +target_matrix = zeros(Int, length(nodes), length(directed_hyperedges)) + +for edge in directed_hyperedges + for node_id in edge.source_ids + source_matrix[node_id, edge.id] = 1 + end + + for node_id in edge.target_ids + target_matrix[node_id, edge.id] = 1 + end +end + +println("\nsource / reactant matrix:") +println(source_matrix) + +println("\ntarget / product matrix:") +println(target_matrix) + +#the DirectedHypergraph constructor accepts matrices with either numbers or nothing +#i am converting 1 values to 1.0 and 0 values to nothing +#this keeps only actual membership entries in the package object + +function matrix_to_package_format(M) + converted = Matrix{Union{Nothing, Float64}}(nothing, size(M, 1), size(M, 2)) + + for i in 1:size(M, 1) + for j in 1:size(M, 2) + if M[i, j] != 0 + converted[i, j] = Float64(M[i, j]) + end + end + end + + return converted +end + +tail_matrix = matrix_to_package_format(source_matrix) +head_matrix = matrix_to_package_format(target_matrix) + +println("\ntail matrix for DirectedHypergraph:") +println(tail_matrix) + +println("\nhead matrix for DirectedHypergraph:") +println(head_matrix) + +#creating the package-level directed hypergraph +#this is the main package connection in this file + +dh = DirectedHypergraph(tail_matrix, head_matrix) + +println("\nDirectedHypergraph object:") +println(dh) + +println("\nobject type:") +println(typeof(dh)) + +#inspecting the object a little bit +#this is useful while learning the package structure + +println("\nfields stored in the DirectedHypergraph object:") +println(fieldnames(typeof(dh))) + +#signed incidence is still useful for checking and for later ml code + +signed_incidence = target_matrix .- source_matrix + +println("\nsigned incidence matrix:") +println(signed_incidence) + +membership_matrix = abs.(signed_incidence) + +println("\nunsigned membership matrix:") +println(membership_matrix) + +#reconstructing reactions from the signed incidence matrix as a check + +function reconstruct_reaction(H, reaction_index, id_to_node) + reactants = String[] + products = String[] + + for i in 1:size(H, 1) + if H[i, reaction_index] == -1 + push!(reactants, id_to_node[i]) + elseif H[i, reaction_index] == 1 + push!(products, id_to_node[i]) + end + end + + return reactants, products +end + +println("\nreconstructed reactions from signed incidence matrix:") + +for j in 1:size(signed_incidence, 2) + reactants, products = reconstruct_reaction(signed_incidence, j, id_to_node) + + println( + "r", j, ": ", + join(reactants, " + "), + " -> ", + join(products, " + ") + ) +end + +#trying package helper functions +#i am wrapping these in try/catch because i am still exploring the package api + +println("\ntrying package helper: to_undirected") + +try + undirected_version = to_undirected(dh) + println(undirected_version) +catch err + println("to_undirected did not run here:") + println(err) +end + +println("\ntrying package helper: get_weakly_connected_components") + +try + components = get_weakly_connected_components(dh) + println(components) +catch err + println("get_weakly_connected_components did not run here:") + println(err) +end + +#simple structural features +#these are the same type of features used in earlier files + +source_count = vec(sum(source_matrix, dims = 2)) +target_count = vec(sum(target_matrix, dims = 2)) +participation_count = vec(sum(membership_matrix, dims = 2)) + +node_features = hcat(source_count, target_count, participation_count) + +println("\nnode feature names:") +println(["source_count", "target_count", "participation_count"]) + +println("\nnode feature matrix:") +println(node_features) + +source_size = vec(sum(source_matrix, dims = 1)) +target_size = vec(sum(target_matrix, dims = 1)) +total_size = vec(sum(membership_matrix, dims = 1)) + +hyperedge_features = hcat(source_size, target_size, total_size) + +println("\nhyperedge feature names:") +println(["source_size", "target_size", "total_size"]) + +println("\nhyperedge feature matrix:") +println(hyperedge_features) + +#keeping everything together +#this object contains both my preprocessing output and the package object + +package_crn_data = ( + raw_reactions = raw_reactions, + parsed_reactions = parsed_reactions, + nodes = nodes, + node_to_id = node_to_id, + id_to_node = id_to_node, + directed_hyperedges = directed_hyperedges, + source_matrix = source_matrix, + target_matrix = target_matrix, + tail_matrix = tail_matrix, + head_matrix = head_matrix, + directed_hypergraph = dh, + signed_incidence = signed_incidence, + membership_matrix = membership_matrix, + node_features = node_features, + hyperedge_features = hyperedge_features +) + +println("\nprocessed package crn data keys:") +println(keys(package_crn_data)) + +println("\nsummary:") +println("number of reactions: ", length(raw_reactions)) +println("number of species / nodes: ", length(nodes)) +println("number of directed hyperedges: ", length(directed_hyperedges)) +println("source matrix size: ", size(source_matrix)) +println("target matrix size: ", size(target_matrix)) +println("tail matrix size: ", size(tail_matrix)) +println("head matrix size: ", size(head_matrix)) +println("node feature matrix size: ", size(node_features)) +println("hyperedge feature matrix size: ", size(hyperedge_features)) +println("package object type: ", typeof(dh)) \ No newline at end of file diff --git a/examples/shonali_prototypes/10_lux_directed_message_passing_layer.jl b/examples/shonali_prototypes/10_lux_directed_message_passing_layer.jl new file mode 100644 index 0000000..16c7281 --- /dev/null +++ b/examples/shonali_prototypes/10_lux_directed_message_passing_layer.jl @@ -0,0 +1,452 @@ +#10_lux_directed_message_passing_layer.jl +#this file is a first draft of a lux-based directed message passing layer +#the main flow is source nodes -> hyperedge embeddings -> target node updates +#i am keeping this version simple so the logic is easy to check before making it more complex + +using Lux +using Random +using Statistics +using SimpleDirectedHypergraphs + +raw_reactions = [ + "A + B -> C", + "C -> D + E", + "E + F -> G", + "G -> H", + "H + A -> I" +] + +println("\nraw crn reactions:") +for reaction in raw_reactions + println(reaction) +end + +#cleaning species names + +function clean_species_name(name) + return String(strip(name)) +end + +#parsing one reaction string into reactants and products + +function parse_reaction(reaction_string) + if !occursin("->", reaction_string) + error("reaction is missing -> : $reaction_string") + end + + left_side, right_side = split(reaction_string, "->") + + reactants = clean_species_name.(split(strip(left_side), "+")) + products = clean_species_name.(split(strip(right_side), "+")) + + reactants = filter(x -> x != "", reactants) + products = filter(x -> x != "", products) + + return ( + raw = reaction_string, + reactants = reactants, + products = products + ) +end + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nparsed reactions:") +for (i, reaction) in enumerate(parsed_reactions) + println( + "r", i, ": ", + join(reaction.reactants, " + "), + " -> ", + join(reaction.products, " + ") + ) +end + +#collecting species as nodes + +species_set = Set{String}() + +for reaction in parsed_reactions + for species in reaction.reactants + push!(species_set, species) + end + + for species in reaction.products + push!(species_set, species) + end +end + +nodes = sort(collect(species_set)) + +println("\nspecies / nodes:") +println(nodes) + +#creating node ids + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nnode to id mapping:") +println(node_to_id) + +#creating directed hyperedges from parsed reactions + +directed_hyperedges = [] + +for (i, reaction) in enumerate(parsed_reactions) + source_ids = [node_to_id[x] for x in reaction.reactants] + target_ids = [node_to_id[x] for x in reaction.products] + + push!( + directed_hyperedges, + ( + id = i, + raw = reaction.raw, + reactants = reaction.reactants, + products = reaction.products, + source_ids = source_ids, + target_ids = target_ids + ) + ) +end + +println("\nid-based directed hyperedges:") +for edge in directed_hyperedges + println("\nr", edge.id) + println("source ids: ", edge.source_ids) + println("target ids: ", edge.target_ids) +end + +#building source and target matrices + +source_matrix = zeros(Float32, length(nodes), length(directed_hyperedges)) +target_matrix = zeros(Float32, length(nodes), length(directed_hyperedges)) + +for edge in directed_hyperedges + for node_id in edge.source_ids + source_matrix[node_id, edge.id] = 1.0f0 + end + + for node_id in edge.target_ids + target_matrix[node_id, edge.id] = 1.0f0 + end +end + +println("\nsource matrix:") +println(source_matrix) + +println("\ntarget matrix:") +println(target_matrix) + +#creating a DirectedHypergraph package object too +#this keeps this file connected to SimpleDirectedHypergraphs.jl + +function matrix_to_package_format(M) + converted = Matrix{Union{Nothing, Float64}}(nothing, size(M, 1), size(M, 2)) + + for i in 1:size(M, 1) + for j in 1:size(M, 2) + if M[i, j] != 0 + converted[i, j] = Float64(M[i, j]) + end + end + end + + return converted +end + +tail_matrix = matrix_to_package_format(source_matrix) +head_matrix = matrix_to_package_format(target_matrix) + +dh = DirectedHypergraph(tail_matrix, head_matrix) + +println("\nDirectedHypergraph object type:") +println(typeof(dh)) + +#building simple node features +#these are just structural features for the first message passing test + +signed_incidence = target_matrix .- source_matrix +membership_matrix = abs.(signed_incidence) + +source_count = vec(sum(source_matrix, dims = 2)) +target_count = vec(sum(target_matrix, dims = 2)) +participation_count = vec(sum(membership_matrix, dims = 2)) + +node_features = Float32.(hcat(source_count, target_count, participation_count)) + +println("\nnode feature names:") +println(["source_count", "target_count", "participation_count"]) + +println("\nnode features:") +println(node_features) + +#normalising features before sending them through dense layers + +function normalise_columns(X) + X_norm = copy(X) + + for j in 1:size(X, 2) + col_mean = mean(X[:, j]) + col_std = std(X[:, j]) + + if col_std == 0 + X_norm[:, j] .= 0 + else + X_norm[:, j] .= (X[:, j] .- col_mean) ./ col_std + end + end + + return Float32.(X_norm) +end + +node_features_norm = normalise_columns(node_features) + +println("\nnormalised node features:") +println(node_features_norm) + +#taking mean features over selected source nodes + +function masked_mean_node_features(node_features, mask_vector) + selected = findall(mask_vector .> 0) + + if length(selected) == 0 + return zeros(Float32, size(node_features, 2)) + end + + return vec(mean(node_features[selected, :], dims = 1)) +end + +#building source aggregated features for every hyperedge + +function source_to_hyperedge_features(node_features, source_matrix) + n_edges = size(source_matrix, 2) + feature_dim = size(node_features, 2) + + edge_features = zeros(Float32, n_edges, feature_dim) + + for e in 1:n_edges + edge_features[e, :] .= masked_mean_node_features(node_features, source_matrix[:, e]) + end + + return edge_features +end + +initial_hyperedge_features = source_to_hyperedge_features(node_features_norm, source_matrix) + +println("\ninitial hyperedge features from source aggregation:") +println(initial_hyperedge_features) + +#setting up lux layers +#source_transform turns source aggregated features into hyperedge embeddings +#target_transform turns hyperedge embeddings into messages for target nodes + +in_dim = size(node_features_norm, 2) +hidden_dim = 4 + +source_transform = Dense(in_dim => hidden_dim, tanh) +target_transform = Dense(hidden_dim => hidden_dim, tanh) + +rng = Random.default_rng() + +ps_source, st_source = Lux.setup(rng, source_transform) +ps_target, st_target = Lux.setup(rng, target_transform) + +println("\nlux source transform:") +println(source_transform) + +println("\nlux target transform:") +println(target_transform) + +println("\nsource transform parameters:") +println(ps_source) + +println("\ntarget transform parameters:") +println(ps_target) + +#running directed message passing +#this is the first draft of the layer behaviour + +function directed_message_passing( + node_features, + source_matrix, + target_matrix, + source_transform, + target_transform, + ps_source, + st_source, + ps_target, + st_target +) + n_nodes = size(node_features, 1) + n_edges = size(source_matrix, 2) + + edge_input = source_to_hyperedge_features(node_features, source_matrix) + + edge_hidden, new_st_source = + source_transform(edge_input', ps_source, st_source) + + edge_hidden = edge_hidden' + + edge_message, new_st_target = + target_transform(edge_hidden', ps_target, st_target) + + edge_message = edge_message' + + node_updates = zeros(Float32, n_nodes, size(edge_message, 2)) + node_update_counts = zeros(Float32, n_nodes) + + for e in 1:n_edges + target_nodes = findall(target_matrix[:, e] .> 0) + + for node_id in target_nodes + node_updates[node_id, :] .+= edge_message[e, :] + node_update_counts[node_id] += 1.0f0 + end + end + + for node_id in 1:n_nodes + if node_update_counts[node_id] > 0 + node_updates[node_id, :] ./= node_update_counts[node_id] + end + end + + return ( + hyperedge_embeddings = edge_hidden, + node_updates = node_updates, + st_source = new_st_source, + st_target = new_st_target + ) +end + +message_passing_output = directed_message_passing( + node_features_norm, + source_matrix, + target_matrix, + source_transform, + target_transform, + ps_source, + st_source, + ps_target, + st_target +) + +hyperedge_embeddings = message_passing_output.hyperedge_embeddings +node_updates = message_passing_output.node_updates + +println("\nhyperedge embeddings:") +println(hyperedge_embeddings) + +println("\nnode updates from directed message passing:") +println(node_updates) + +println("\nhyperedge embedding size:") +println(size(hyperedge_embeddings)) + +println("\nnode update size:") +println(size(node_updates)) + +#showing message passing by reaction + +println("\nmessage passing by reaction:") + +for edge in directed_hyperedges + println("\nr", edge.id, ": ", edge.raw) + println("source nodes: ", edge.reactants) + println("target nodes: ", edge.products) + println("hyperedge embedding: ", hyperedge_embeddings[edge.id, :]) +end + +#showing node updates by species + +println("\nnode updates by species:") + +for i in 1:length(nodes) + println(nodes[i], " update: ", node_updates[i, :]) +end + +#adding a small prediction head +#this is only to check that hyperedge embeddings can feed into another lux layer + +prediction_head = Dense(hidden_dim => 1) + +ps_head, st_head = Lux.setup(rng, prediction_head) + +scores, st_head_new = prediction_head(hyperedge_embeddings', ps_head, st_head) +scores = vec(scores) + +println("\ntoy reaction scores from prediction head:") +println(scores) + +#creating toy labels for checking the forward pass +#label is 1 if the reaction has more than one reactant + +toy_labels = Float32[ + length(edge.source_ids) > 1 ? 1.0f0 : 0.0f0 + for edge in directed_hyperedges +] + +println("\ntoy labels:") +println(toy_labels) + +#turning scores into probabilities + +function sigmoid(x) + return 1.0f0 / (1.0f0 + exp(-x)) +end + +probabilities = sigmoid.(scores) + +println("\ntoy probabilities:") +println(probabilities) + +predictions = Int.(probabilities .>= 0.5f0) + +println("\ntoy predictions:") +println(predictions) + +println("\nactual toy labels:") +println(Int.(toy_labels)) + +accuracy = mean(predictions .== Int.(toy_labels)) + +println("\ntoy accuracy:") +println(accuracy) + +#storing outputs from this draft layer + +lux_message_passing_data = ( + raw_reactions = raw_reactions, + parsed_reactions = parsed_reactions, + nodes = nodes, + node_to_id = node_to_id, + id_to_node = id_to_node, + directed_hyperedges = directed_hyperedges, + directed_hypergraph = dh, + source_matrix = source_matrix, + target_matrix = target_matrix, + signed_incidence = signed_incidence, + membership_matrix = membership_matrix, + node_features = node_features, + node_features_norm = node_features_norm, + initial_hyperedge_features = initial_hyperedge_features, + hyperedge_embeddings = hyperedge_embeddings, + node_updates = node_updates, + toy_scores = scores, + toy_probabilities = probabilities, + toy_predictions = predictions, + toy_labels = toy_labels +) + +println("\nprocessed lux message passing data keys:") +println(keys(lux_message_passing_data)) + +println("\nsummary:") +println("number of nodes: ", length(nodes)) +println("number of reactions / hyperedges: ", length(directed_hyperedges)) +println("node feature matrix size: ", size(node_features)) +println("normalised node feature matrix size: ", size(node_features_norm)) +println("source matrix size: ", size(source_matrix)) +println("target matrix size: ", size(target_matrix)) +println("hyperedge embedding size: ", size(hyperedge_embeddings)) +println("node update size: ", size(node_updates)) +println("package directed hypergraph type: ", typeof(dh)) \ No newline at end of file diff --git a/examples/shonali_prototypes/11bz_crn_case_study.jl b/examples/shonali_prototypes/11bz_crn_case_study.jl new file mode 100644 index 0000000..cf52410 --- /dev/null +++ b/examples/shonali_prototypes/11bz_crn_case_study.jl @@ -0,0 +1,513 @@ +#11_bz_crn_case_study.jl +#this file applies the crn parser and lux message passing pipeline to a small bz-inspired crn case study +#the goal is to move from toy reactions to a more chemistry-style reaction network +#this is still a simplified case study, but it checks whether the full workflow works on a bz-style crn + +using Lux +using Random +using Statistics +using SimpleDirectedHypergraphs + +raw_reactions = [ + "BrO3 + Br -> HBrO2", + "HBrO2 + Br -> HOBr", + "BrO3 + HBrO2 -> BrO2 + HOBr", + "BrO2 + Ce3 -> Ce4 + HBrO2", + "Ce4 + MalonicAcid -> Ce3 + Br", + "HBrO2 + HBrO2 -> BrO3 + HOBr", + "HOBr + Br -> Br2", + "Br2 + MalonicAcid -> BrMalonicAcid + Br" +] + +println("\nbz-inspired crn reactions:") +for reaction in raw_reactions + println(reaction) +end + +#cleaning species names + +function clean_species_name(name) + return String(strip(name)) +end + +#parsing one reaction string into reactants and products + +function parse_reaction(reaction_string) + if !occursin("->", reaction_string) + error("reaction is missing -> : $reaction_string") + end + + left_side, right_side = split(reaction_string, "->") + + reactants = clean_species_name.(split(strip(left_side), "+")) + products = clean_species_name.(split(strip(right_side), "+")) + + reactants = filter(x -> x != "", reactants) + products = filter(x -> x != "", products) + + return ( + raw = reaction_string, + reactants = reactants, + products = products + ) +end + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nparsed reactions:") +for (i, reaction) in enumerate(parsed_reactions) + println( + "r", i, ": ", + join(reaction.reactants, " + "), + " -> ", + join(reaction.products, " + ") + ) +end + +#collecting all species as nodes + +species_set = Set{String}() + +for reaction in parsed_reactions + for species in reaction.reactants + push!(species_set, species) + end + + for species in reaction.products + push!(species_set, species) + end +end + +nodes = sort(collect(species_set)) + +println("\nspecies / nodes:") +println(nodes) + +#creating node ids + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nnode to id mapping:") +println(node_to_id) + +println("\nid to node mapping:") +println(id_to_node) + +#creating directed hyperedges from parsed reactions + +directed_hyperedges = [] + +for (i, reaction) in enumerate(parsed_reactions) + source_ids = [node_to_id[x] for x in reaction.reactants] + target_ids = [node_to_id[x] for x in reaction.products] + + push!( + directed_hyperedges, + ( + id = i, + raw = reaction.raw, + reactants = reaction.reactants, + products = reaction.products, + source_ids = source_ids, + target_ids = target_ids + ) + ) +end + +println("\nid-based directed hyperedges:") +for edge in directed_hyperedges + println("\nr", edge.id) + println("source/reactants: ", edge.reactants, " -> ", edge.source_ids) + println("target/products: ", edge.products, " -> ", edge.target_ids) +end + +#building source and target matrices +#rows are species and columns are reactions + +source_matrix = zeros(Float32, length(nodes), length(directed_hyperedges)) +target_matrix = zeros(Float32, length(nodes), length(directed_hyperedges)) + +for edge in directed_hyperedges + for node_id in edge.source_ids + source_matrix[node_id, edge.id] = 1.0f0 + end + + for node_id in edge.target_ids + target_matrix[node_id, edge.id] = 1.0f0 + end +end + +println("\nsource / reactant matrix:") +println(source_matrix) + +println("\ntarget / product matrix:") +println(target_matrix) + +#converting matrices into the format used by SimpleDirectedHypergraphs.jl + +function matrix_to_package_format(M) + converted = Matrix{Union{Nothing, Float64}}(nothing, size(M, 1), size(M, 2)) + + for i in 1:size(M, 1) + for j in 1:size(M, 2) + if M[i, j] != 0 + converted[i, j] = Float64(M[i, j]) + end + end + end + + return converted +end + +tail_matrix = matrix_to_package_format(source_matrix) +head_matrix = matrix_to_package_format(target_matrix) + +bz_dh = DirectedHypergraph(tail_matrix, head_matrix) + +println("\nDirectedHypergraph object type:") +println(typeof(bz_dh)) + +#creating incidence and membership matrices + +signed_incidence = target_matrix .- source_matrix +membership_matrix = abs.(signed_incidence) + +println("\nsigned incidence matrix:") +println(signed_incidence) + +println("\nunsigned membership matrix:") +println(membership_matrix) + +#reconstructing reactions from the signed incidence matrix +#this is a check that the matrix representation still matches the reaction list + +function reconstruct_reaction(H, reaction_index, id_to_node) + reactants = String[] + products = String[] + + for i in 1:size(H, 1) + if H[i, reaction_index] == -1 + push!(reactants, id_to_node[i]) + elseif H[i, reaction_index] == 1 + push!(products, id_to_node[i]) + end + end + + return reactants, products +end + +println("\nreconstructed reactions from signed incidence matrix:") +for j in 1:size(signed_incidence, 2) + reactants, products = reconstruct_reaction(signed_incidence, j, id_to_node) + + println( + "r", j, ": ", + join(reactants, " + "), + " -> ", + join(products, " + ") + ) +end + +#trying a package helper function +#this checks that the package object can be used for graph-style analysis + +println("\nweakly connected components:") + +try + components = get_weakly_connected_components(bz_dh) + println(components) +catch err + println("could not compute weakly connected components:") + println(err) +end + +#building structural node features +#these are simple features for the first message passing test + +source_count = vec(sum(source_matrix, dims = 2)) +target_count = vec(sum(target_matrix, dims = 2)) +participation_count = vec(sum(membership_matrix, dims = 2)) + +node_features = Float32.(hcat(source_count, target_count, participation_count)) + +println("\nnode feature names:") +println(["source_count", "target_count", "participation_count"]) + +println("\nnode features:") +println(node_features) + +#building hyperedge features + +source_size = vec(sum(source_matrix, dims = 1)) +target_size = vec(sum(target_matrix, dims = 1)) +total_size = vec(sum(membership_matrix, dims = 1)) + +hyperedge_features = Float32.(hcat(source_size, target_size, total_size)) + +println("\nhyperedge feature names:") +println(["source_size", "target_size", "total_size"]) + +println("\nhyperedge features:") +println(hyperedge_features) + +#normalising node features before using them in lux layers + +function normalise_columns(X) + X_norm = copy(X) + + for j in 1:size(X, 2) + col_mean = mean(X[:, j]) + col_std = std(X[:, j]) + + if col_std == 0 + X_norm[:, j] .= 0 + else + X_norm[:, j] .= (X[:, j] .- col_mean) ./ col_std + end + end + + return Float32.(X_norm) +end + +node_features_norm = normalise_columns(node_features) + +println("\nnormalised node features:") +println(node_features_norm) + +#taking the mean of selected source node features + +function masked_mean_node_features(node_features, mask_vector) + selected = findall(mask_vector .> 0) + + if length(selected) == 0 + return zeros(Float32, size(node_features, 2)) + end + + return vec(mean(node_features[selected, :], dims = 1)) +end + +#creating source aggregated features for every reaction + +function source_to_hyperedge_features(node_features, source_matrix) + n_edges = size(source_matrix, 2) + feature_dim = size(node_features, 2) + + edge_features = zeros(Float32, n_edges, feature_dim) + + for e in 1:n_edges + edge_features[e, :] .= masked_mean_node_features(node_features, source_matrix[:, e]) + end + + return edge_features +end + +initial_hyperedge_features = source_to_hyperedge_features(node_features_norm, source_matrix) + +println("\ninitial hyperedge features from source aggregation:") +println(initial_hyperedge_features) + +#setting up a simple lux message passing draft +#source_transform creates reaction embeddings from source node information +#target_transform creates messages from reaction embeddings + +in_dim = size(node_features_norm, 2) +hidden_dim = 4 + +source_transform = Dense(in_dim => hidden_dim, tanh) +target_transform = Dense(hidden_dim => hidden_dim, tanh) + +rng = Random.default_rng() + +ps_source, st_source = Lux.setup(rng, source_transform) +ps_target, st_target = Lux.setup(rng, target_transform) + +println("\nlux source transform:") +println(source_transform) + +println("\nlux target transform:") +println(target_transform) + +#running directed message passing + +function directed_message_passing( + node_features, + source_matrix, + target_matrix, + source_transform, + target_transform, + ps_source, + st_source, + ps_target, + st_target +) + n_nodes = size(node_features, 1) + n_edges = size(source_matrix, 2) + + edge_input = source_to_hyperedge_features(node_features, source_matrix) + + edge_hidden, new_st_source = + source_transform(edge_input', ps_source, st_source) + + edge_hidden = edge_hidden' + + edge_message, new_st_target = + target_transform(edge_hidden', ps_target, st_target) + + edge_message = edge_message' + + node_updates = zeros(Float32, n_nodes, size(edge_message, 2)) + node_update_counts = zeros(Float32, n_nodes) + + for e in 1:n_edges + target_nodes = findall(target_matrix[:, e] .> 0) + + for node_id in target_nodes + node_updates[node_id, :] .+= edge_message[e, :] + node_update_counts[node_id] += 1.0f0 + end + end + + for node_id in 1:n_nodes + if node_update_counts[node_id] > 0 + node_updates[node_id, :] ./= node_update_counts[node_id] + end + end + + return ( + hyperedge_embeddings = edge_hidden, + node_updates = node_updates, + st_source = new_st_source, + st_target = new_st_target + ) +end + +message_passing_output = directed_message_passing( + node_features_norm, + source_matrix, + target_matrix, + source_transform, + target_transform, + ps_source, + st_source, + ps_target, + st_target +) + +hyperedge_embeddings = message_passing_output.hyperedge_embeddings +node_updates = message_passing_output.node_updates + +println("\nhyperedge embeddings:") +println(hyperedge_embeddings) + +println("\nnode updates from directed message passing:") +println(node_updates) + +#showing the message passing result reaction by reaction + +println("\nmessage passing by reaction:") + +for edge in directed_hyperedges + println("\nr", edge.id, ": ", edge.raw) + println("source nodes: ", edge.reactants) + println("target nodes: ", edge.products) + println("hyperedge embedding: ", hyperedge_embeddings[edge.id, :]) +end + +#showing the node updates species by species + +println("\nnode updates by species:") + +for i in 1:length(nodes) + println(nodes[i], " update: ", node_updates[i, :]) +end + +#adding a tiny prediction head +#this is not the final learning task +#it only checks that reaction embeddings can feed into another lux layer + +prediction_head = Dense(hidden_dim => 1) + +ps_head, st_head = Lux.setup(rng, prediction_head) + +scores, st_head_new = prediction_head(hyperedge_embeddings', ps_head, st_head) +scores = vec(scores) + +println("\ntoy reaction scores:") +println(scores) + +#creating simple toy labels for testing +#label is 1 if a reaction has more than one reactant + +toy_labels = Float32[ + length(edge.source_ids) > 1 ? 1.0f0 : 0.0f0 + for edge in directed_hyperedges +] + +println("\ntoy labels:") +println(toy_labels) + +function sigmoid(x) + return 1.0f0 / (1.0f0 + exp(-x)) +end + +probabilities = sigmoid.(scores) +predictions = Int.(probabilities .>= 0.5f0) + +println("\ntoy probabilities:") +println(probabilities) + +println("\ntoy predictions:") +println(predictions) + +println("\nactual toy labels:") +println(Int.(toy_labels)) + +accuracy = mean(predictions .== Int.(toy_labels)) + +println("\ntoy accuracy:") +println(accuracy) + +#storing the bz case study data + +bz_case_study_data = ( + case_study_name = "bz-inspired crn", + raw_reactions = raw_reactions, + parsed_reactions = parsed_reactions, + nodes = nodes, + node_to_id = node_to_id, + id_to_node = id_to_node, + directed_hyperedges = directed_hyperedges, + directed_hypergraph = bz_dh, + source_matrix = source_matrix, + target_matrix = target_matrix, + signed_incidence = signed_incidence, + membership_matrix = membership_matrix, + node_features = node_features, + hyperedge_features = hyperedge_features, + node_features_norm = node_features_norm, + initial_hyperedge_features = initial_hyperedge_features, + hyperedge_embeddings = hyperedge_embeddings, + node_updates = node_updates, + toy_scores = scores, + toy_probabilities = probabilities, + toy_predictions = predictions, + toy_labels = toy_labels +) + +println("\nprocessed bz case study data keys:") +println(keys(bz_case_study_data)) + +println("\nsummary:") +println("case study: ", bz_case_study_data.case_study_name) +println("number of reactions: ", length(raw_reactions)) +println("number of species / nodes: ", length(nodes)) +println("number of directed hyperedges: ", length(directed_hyperedges)) +println("source matrix size: ", size(source_matrix)) +println("target matrix size: ", size(target_matrix)) +println("signed incidence matrix size: ", size(signed_incidence)) +println("node feature matrix size: ", size(node_features)) +println("hyperedge feature matrix size: ", size(hyperedge_features)) +println("hyperedge embedding size: ", size(hyperedge_embeddings)) +println("node update size: ", size(node_updates)) +println("package directed hypergraph type: ", typeof(bz_dh)) \ No newline at end of file diff --git a/examples/shonali_prototypes/12_formose_crn_case_study.jl b/examples/shonali_prototypes/12_formose_crn_case_study.jl new file mode 100644 index 0000000..7532d01 --- /dev/null +++ b/examples/shonali_prototypes/12_formose_crn_case_study.jl @@ -0,0 +1,514 @@ +#12_formose_crn_case_study.jl +#this file applies the crn parser and lux message passing pipeline to a small formose-inspired crn case study +#the goal is to test the same workflow on another chemistry-style reaction network +#this is still a simplified formose-inspired example, not the final full literature network + +using Lux +using Random +using Statistics +using SimpleDirectedHypergraphs + +raw_reactions = [ + "CH2O + CH2O -> Glycolaldehyde", + "Glycolaldehyde + CH2O -> Glyceraldehyde", + "Glyceraldehyde -> Dihydroxyacetone", + "Dihydroxyacetone + CH2O -> Tetrose", + "Tetrose -> Glycolaldehyde + Glycolaldehyde", + "Glyceraldehyde + Glycolaldehyde -> Pentose", + "Pentose -> Dihydroxyacetone + Glycolaldehyde", + "Tetrose + CH2O -> Hexose", + "Hexose -> Glyceraldehyde + Glyceraldehyde" +] + +println("\nformose-inspired crn reactions:") +for reaction in raw_reactions + println(reaction) +end + +#cleaning species names + +function clean_species_name(name) + return String(strip(name)) +end + +#parsing one reaction string into reactants and products + +function parse_reaction(reaction_string) + if !occursin("->", reaction_string) + error("reaction is missing -> : $reaction_string") + end + + left_side, right_side = split(reaction_string, "->") + + reactants = clean_species_name.(split(strip(left_side), "+")) + products = clean_species_name.(split(strip(right_side), "+")) + + reactants = filter(x -> x != "", reactants) + products = filter(x -> x != "", products) + + return ( + raw = reaction_string, + reactants = reactants, + products = products + ) +end + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nparsed reactions:") +for (i, reaction) in enumerate(parsed_reactions) + println( + "r", i, ": ", + join(reaction.reactants, " + "), + " -> ", + join(reaction.products, " + ") + ) +end + +#collecting all species as nodes + +species_set = Set{String}() + +for reaction in parsed_reactions + for species in reaction.reactants + push!(species_set, species) + end + + for species in reaction.products + push!(species_set, species) + end +end + +nodes = sort(collect(species_set)) + +println("\nspecies / nodes:") +println(nodes) + +#creating node ids + +node_to_id = Dict(node => i for (i, node) in enumerate(nodes)) +id_to_node = Dict(i => node for (i, node) in enumerate(nodes)) + +println("\nnode to id mapping:") +println(node_to_id) + +println("\nid to node mapping:") +println(id_to_node) + +#creating directed hyperedges from parsed reactions + +directed_hyperedges = [] + +for (i, reaction) in enumerate(parsed_reactions) + source_ids = [node_to_id[x] for x in reaction.reactants] + target_ids = [node_to_id[x] for x in reaction.products] + + push!( + directed_hyperedges, + ( + id = i, + raw = reaction.raw, + reactants = reaction.reactants, + products = reaction.products, + source_ids = source_ids, + target_ids = target_ids + ) + ) +end + +println("\nid-based directed hyperedges:") +for edge in directed_hyperedges + println("\nr", edge.id) + println("source/reactants: ", edge.reactants, " -> ", edge.source_ids) + println("target/products: ", edge.products, " -> ", edge.target_ids) +end + +#building source and target matrices +#rows are species and columns are reactions + +source_matrix = zeros(Float32, length(nodes), length(directed_hyperedges)) +target_matrix = zeros(Float32, length(nodes), length(directed_hyperedges)) + +for edge in directed_hyperedges + for node_id in edge.source_ids + source_matrix[node_id, edge.id] = 1.0f0 + end + + for node_id in edge.target_ids + target_matrix[node_id, edge.id] = 1.0f0 + end +end + +println("\nsource / reactant matrix:") +println(source_matrix) + +println("\ntarget / product matrix:") +println(target_matrix) + +#converting matrices into the format used by SimpleDirectedHypergraphs.jl + +function matrix_to_package_format(M) + converted = Matrix{Union{Nothing, Float64}}(nothing, size(M, 1), size(M, 2)) + + for i in 1:size(M, 1) + for j in 1:size(M, 2) + if M[i, j] != 0 + converted[i, j] = Float64(M[i, j]) + end + end + end + + return converted +end + +tail_matrix = matrix_to_package_format(source_matrix) +head_matrix = matrix_to_package_format(target_matrix) + +formose_dh = DirectedHypergraph(tail_matrix, head_matrix) + +println("\nDirectedHypergraph object type:") +println(typeof(formose_dh)) + +#creating incidence and membership matrices + +signed_incidence = target_matrix .- source_matrix +membership_matrix = abs.(signed_incidence) + +println("\nsigned incidence matrix:") +println(signed_incidence) + +println("\nunsigned membership matrix:") +println(membership_matrix) + +#reconstructing reactions from the signed incidence matrix +#this is a check that the matrix representation still matches the reaction list + +function reconstruct_reaction(H, reaction_index, id_to_node) + reactants = String[] + products = String[] + + for i in 1:size(H, 1) + if H[i, reaction_index] == -1 + push!(reactants, id_to_node[i]) + elseif H[i, reaction_index] == 1 + push!(products, id_to_node[i]) + end + end + + return reactants, products +end + +println("\nreconstructed reactions from signed incidence matrix:") +for j in 1:size(signed_incidence, 2) + reactants, products = reconstruct_reaction(signed_incidence, j, id_to_node) + + println( + "r", j, ": ", + join(reactants, " + "), + " -> ", + join(products, " + ") + ) +end + +#trying a package helper function +#this checks that the package object can be used for graph-style analysis + +println("\nweakly connected components:") + +try + components = get_weakly_connected_components(formose_dh) + println(components) +catch err + println("could not compute weakly connected components:") + println(err) +end + +#building structural node features +#these are simple features for the first message passing test + +source_count = vec(sum(source_matrix, dims = 2)) +target_count = vec(sum(target_matrix, dims = 2)) +participation_count = vec(sum(membership_matrix, dims = 2)) + +node_features = Float32.(hcat(source_count, target_count, participation_count)) + +println("\nnode feature names:") +println(["source_count", "target_count", "participation_count"]) + +println("\nnode features:") +println(node_features) + +#building hyperedge features + +source_size = vec(sum(source_matrix, dims = 1)) +target_size = vec(sum(target_matrix, dims = 1)) +total_size = vec(sum(membership_matrix, dims = 1)) + +hyperedge_features = Float32.(hcat(source_size, target_size, total_size)) + +println("\nhyperedge feature names:") +println(["source_size", "target_size", "total_size"]) + +println("\nhyperedge features:") +println(hyperedge_features) + +#normalising node features before using them in lux layers + +function normalise_columns(X) + X_norm = copy(X) + + for j in 1:size(X, 2) + col_mean = mean(X[:, j]) + col_std = std(X[:, j]) + + if col_std == 0 + X_norm[:, j] .= 0 + else + X_norm[:, j] .= (X[:, j] .- col_mean) ./ col_std + end + end + + return Float32.(X_norm) +end + +node_features_norm = normalise_columns(node_features) + +println("\nnormalised node features:") +println(node_features_norm) + +#taking the mean of selected source node features + +function masked_mean_node_features(node_features, mask_vector) + selected = findall(mask_vector .> 0) + + if length(selected) == 0 + return zeros(Float32, size(node_features, 2)) + end + + return vec(mean(node_features[selected, :], dims = 1)) +end + +#creating source aggregated features for every reaction + +function source_to_hyperedge_features(node_features, source_matrix) + n_edges = size(source_matrix, 2) + feature_dim = size(node_features, 2) + + edge_features = zeros(Float32, n_edges, feature_dim) + + for e in 1:n_edges + edge_features[e, :] .= masked_mean_node_features(node_features, source_matrix[:, e]) + end + + return edge_features +end + +initial_hyperedge_features = source_to_hyperedge_features(node_features_norm, source_matrix) + +println("\ninitial hyperedge features from source aggregation:") +println(initial_hyperedge_features) + +#setting up a simple lux message passing draft +#source_transform creates reaction embeddings from source node information +#target_transform creates messages from reaction embeddings + +in_dim = size(node_features_norm, 2) +hidden_dim = 4 + +source_transform = Dense(in_dim => hidden_dim, tanh) +target_transform = Dense(hidden_dim => hidden_dim, tanh) + +rng = Random.default_rng() + +ps_source, st_source = Lux.setup(rng, source_transform) +ps_target, st_target = Lux.setup(rng, target_transform) + +println("\nlux source transform:") +println(source_transform) + +println("\nlux target transform:") +println(target_transform) + +#running directed message passing + +function directed_message_passing( + node_features, + source_matrix, + target_matrix, + source_transform, + target_transform, + ps_source, + st_source, + ps_target, + st_target +) + n_nodes = size(node_features, 1) + n_edges = size(source_matrix, 2) + + edge_input = source_to_hyperedge_features(node_features, source_matrix) + + edge_hidden, new_st_source = + source_transform(edge_input', ps_source, st_source) + + edge_hidden = edge_hidden' + + edge_message, new_st_target = + target_transform(edge_hidden', ps_target, st_target) + + edge_message = edge_message' + + node_updates = zeros(Float32, n_nodes, size(edge_message, 2)) + node_update_counts = zeros(Float32, n_nodes) + + for e in 1:n_edges + target_nodes = findall(target_matrix[:, e] .> 0) + + for node_id in target_nodes + node_updates[node_id, :] .+= edge_message[e, :] + node_update_counts[node_id] += 1.0f0 + end + end + + for node_id in 1:n_nodes + if node_update_counts[node_id] > 0 + node_updates[node_id, :] ./= node_update_counts[node_id] + end + end + + return ( + hyperedge_embeddings = edge_hidden, + node_updates = node_updates, + st_source = new_st_source, + st_target = new_st_target + ) +end + +message_passing_output = directed_message_passing( + node_features_norm, + source_matrix, + target_matrix, + source_transform, + target_transform, + ps_source, + st_source, + ps_target, + st_target +) + +hyperedge_embeddings = message_passing_output.hyperedge_embeddings +node_updates = message_passing_output.node_updates + +println("\nhyperedge embeddings:") +println(hyperedge_embeddings) + +println("\nnode updates from directed message passing:") +println(node_updates) + +#showing the message passing result reaction by reaction + +println("\nmessage passing by reaction:") + +for edge in directed_hyperedges + println("\nr", edge.id, ": ", edge.raw) + println("source nodes: ", edge.reactants) + println("target nodes: ", edge.products) + println("hyperedge embedding: ", hyperedge_embeddings[edge.id, :]) +end + +#showing the node updates species by species + +println("\nnode updates by species:") + +for i in 1:length(nodes) + println(nodes[i], " update: ", node_updates[i, :]) +end + +#adding a tiny prediction head +#this is not the final learning task +#it only checks that reaction embeddings can feed into another lux layer + +prediction_head = Dense(hidden_dim => 1) + +ps_head, st_head = Lux.setup(rng, prediction_head) + +scores, st_head_new = prediction_head(hyperedge_embeddings', ps_head, st_head) +scores = vec(scores) + +println("\ntoy reaction scores:") +println(scores) + +#creating simple toy labels for testing +#label is 1 if a reaction has more than one reactant + +toy_labels = Float32[ + length(edge.source_ids) > 1 ? 1.0f0 : 0.0f0 + for edge in directed_hyperedges +] + +println("\ntoy labels:") +println(toy_labels) + +function sigmoid(x) + return 1.0f0 / (1.0f0 + exp(-x)) +end + +probabilities = sigmoid.(scores) +predictions = Int.(probabilities .>= 0.5f0) + +println("\ntoy probabilities:") +println(probabilities) + +println("\ntoy predictions:") +println(predictions) + +println("\nactual toy labels:") +println(Int.(toy_labels)) + +accuracy = mean(predictions .== Int.(toy_labels)) + +println("\ntoy accuracy:") +println(accuracy) + +#storing the formose case study data + +formose_case_study_data = ( + case_study_name = "formose-inspired crn", + raw_reactions = raw_reactions, + parsed_reactions = parsed_reactions, + nodes = nodes, + node_to_id = node_to_id, + id_to_node = id_to_node, + directed_hyperedges = directed_hyperedges, + directed_hypergraph = formose_dh, + source_matrix = source_matrix, + target_matrix = target_matrix, + signed_incidence = signed_incidence, + membership_matrix = membership_matrix, + node_features = node_features, + hyperedge_features = hyperedge_features, + node_features_norm = node_features_norm, + initial_hyperedge_features = initial_hyperedge_features, + hyperedge_embeddings = hyperedge_embeddings, + node_updates = node_updates, + toy_scores = scores, + toy_probabilities = probabilities, + toy_predictions = predictions, + toy_labels = toy_labels +) + +println("\nprocessed formose case study data keys:") +println(keys(formose_case_study_data)) + +println("\nsummary:") +println("case study: ", formose_case_study_data.case_study_name) +println("number of reactions: ", length(raw_reactions)) +println("number of species / nodes: ", length(nodes)) +println("number of directed hyperedges: ", length(directed_hyperedges)) +println("source matrix size: ", size(source_matrix)) +println("target matrix size: ", size(target_matrix)) +println("signed incidence matrix size: ", size(signed_incidence)) +println("node feature matrix size: ", size(node_features)) +println("hyperedge feature matrix size: ", size(hyperedge_features)) +println("hyperedge embedding size: ", size(hyperedge_embeddings)) +println("node update size: ", size(node_updates)) +println("package directed hypergraph type: ", typeof(formose_dh)) \ No newline at end of file diff --git a/examples/shonali_prototypes/13_species_feature_engineering.jl b/examples/shonali_prototypes/13_species_feature_engineering.jl new file mode 100644 index 0000000..6bc630a --- /dev/null +++ b/examples/shonali_prototypes/13_species_feature_engineering.jl @@ -0,0 +1,417 @@ +#13 species feature engineering + +#purpose: +#build species-level feature engineering for chemical reaction networks represented as directed hypergraphs. + +#this file supports the thesis question: +#which structural and chemistry-informed features improve reaction-level regression performance? + +using Statistics +using LinearAlgebra + +println("Species Feature Engineering") + + +#1. defining a small Formose-inspired CRN + + +raw_reactions = [ + "CH2O + CH2O -> Glycolaldehyde", + "Glycolaldehyde + CH2O -> Glyceraldehyde", + "Glyceraldehyde -> Dihydroxyacetone", + "Dihydroxyacetone + CH2O -> Tetrose", + "Tetrose -> Glycolaldehyde + Glycolaldehyde", + "Glyceraldehyde + Glycolaldehyde -> Pentose", + "Pentose -> Dihydroxyacetone + Glycolaldehyde", + "Tetrose + CH2O -> Hexose", + "Hexose -> Glyceraldehyde + Glyceraldehyde" +] + +println("Raw reactions:") +for r in raw_reactions + println(r) +end + +#2.basic reaction parser + +#abstractString avoids errors with SubString values returned by split(). + +function parse_side(side::AbstractString) + species = split(strip(side), "+") + return strip.(species) +end + +function parse_reaction(reaction::AbstractString) + sides = split(reaction, "->") + + if length(sides) != 2 + error("Reaction must contain exactly one -> symbol: $reaction") + end + + reactants = parse_side(sides[1]) + products = parse_side(sides[2]) + + return reactants, products +end + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nParsed reactions:") +for (i, (reactants, products)) in enumerate(parsed_reactions) + println("r$i: ", join(reactants, " + "), " -> ", join(products, " + ")) +end + + +#3.species / node mappings + + +all_species = String[] + +for (reactants, products) in parsed_reactions + append!(all_species, String.(reactants)) + append!(all_species, String.(products)) +end + +species = sort(unique(all_species)) + +node_to_id = Dict(s => i for (i, s) in enumerate(species)) +id_to_node = Dict(i => s for (s, i) in node_to_id) + +num_species = length(species) +num_reactions = length(parsed_reactions) + +println("\nSpecies / nodes:") +println(species) + +println("\nNode to ID mapping:") +println(node_to_id) + + +#4.build source and target matrices + + +#source_matrix[i, j] = number of times species i appears as a reactant in reaction j +#target_matrix[i, j] = number of times species i appears as a product in reaction j + +#this preserves repeated species such as: +#CH2O + CH2O -> Glycolaldehyde + +#here source_matrix[CH2O, r1] = 2.0 + +source_matrix = zeros(Float32, num_species, num_reactions) +target_matrix = zeros(Float32, num_species, num_reactions) + +for (j, (reactants, products)) in enumerate(parsed_reactions) + for r in reactants + source_matrix[node_to_id[String(r)], j] += 1.0f0 + end + + for p in products + target_matrix[node_to_id[String(p)], j] += 1.0f0 + end +end + +println("\nSource / reactant matrix:") +println(source_matrix) + +println("\nTarget / product matrix:") +println(target_matrix) + + +#5.structural species features + + +#these features describe the role of each species in the reaction network +#they are simple baseline features that can later be compared against chemistry-informed descriptors + +source_count = vec(sum(source_matrix .> 0, dims = 2)) +target_count = vec(sum(target_matrix .> 0, dims = 2)) +participation_count = source_count .+ target_count + +#stoichiometry-aware counts +source_stoich_sum = vec(sum(source_matrix, dims = 2)) +target_stoich_sum = vec(sum(target_matrix, dims = 2)) +total_stoich_sum = source_stoich_sum .+ target_stoich_sum + +function neighbouring_species(species_id::Int, parsed_reactions, node_to_id) + neighbours = Set{Int}() + + for (reactants, products) in parsed_reactions + reaction_species = vcat(String.(reactants), String.(products)) + ids = [node_to_id[s] for s in reaction_species] + + if species_id in ids + for id in ids + if id != species_id + push!(neighbours, id) + end + end + end + end + + return length(neighbours) +end + +neighbour_count = Float32[ + neighbouring_species(i, parsed_reactions, node_to_id) + for i in 1:num_species +] + +reaction_count = participation_count + +structural_feature_names = [ + "source_count", + "target_count", + "participation_count", + "source_stoich_sum", + "target_stoich_sum", + "total_stoich_sum", + "neighbour_count", + "reaction_count" +] + +structural_features = hcat( + source_count, + target_count, + participation_count, + source_stoich_sum, + target_stoich_sum, + total_stoich_sum, + neighbour_count, + reaction_count +) + +structural_features = Float32.(structural_features) + +println("\nStructural feature names:") +println(structural_feature_names) + +println("\nStructural species features:") +println(structural_features) + + +#6.chemistry-informed species features + + +#in a complete version, these could come from molecular descriptors, molecular fingerprints, or chemistry ML tools + +#here they are manually added to demonstrate how chemical information can be attached to species nodes. + +chemical_feature_names = [ + "approx_molecular_weight", + "num_carbon", + "num_hydrogen", + "num_oxygen", + "num_atoms", + "formal_charge" +] + +chemical_features_dict = Dict( + "CH2O" => Float32[30.03, 1, 2, 1, 4, 0], + "Glycolaldehyde" => Float32[60.05, 2, 4, 2, 8, 0], + "Glyceraldehyde" => Float32[90.08, 3, 6, 3, 12, 0], + "Dihydroxyacetone" => Float32[90.08, 3, 6, 3, 12, 0], + "Tetrose" => Float32[120.10, 4, 8, 4, 16, 0], + "Pentose" => Float32[150.13, 5, 10, 5, 20, 0], + "Hexose" => Float32[180.16, 6, 12, 6, 24, 0] +) + +chemical_features = zeros(Float32, num_species, length(chemical_feature_names)) + +for (i, s) in enumerate(species) + if haskey(chemical_features_dict, s) + chemical_features[i, :] .= chemical_features_dict[s] + else + @warn "No chemistry-informed features found for species $s. Using zeros." + end +end + +println("\nChemistry-informed feature names:") +println(chemical_feature_names) + +println("\nChemistry-informed species features:") +println(chemical_features) + + +#7.combine feature sets + + +combined_feature_names = vcat(structural_feature_names, chemical_feature_names) +combined_features = hcat(structural_features, chemical_features) + +println("\nCombined feature names:") +println(combined_feature_names) + +println("\nCombined species features:") +println(combined_features) + + +#8.feature normalisation + + +#standardise each feature column: +#x_norm = (x - mean) / standard deviation + +#if a column has zero variance, use standard deviation = 1 to avoid division by zero. + +function standardise_features(X::AbstractMatrix) + μ = mean(X, dims = 1) + σ = std(X, dims = 1) + + σ_safe = similar(σ) + + for i in eachindex(σ) + σ_safe[i] = σ[i] == 0 ? 1 : σ[i] + end + + X_norm = (X .- μ) ./ σ_safe + + return Float32.(X_norm), Float32.(μ), Float32.(σ_safe) +end + +structural_features_norm, structural_mean, structural_std = + standardise_features(structural_features) + +chemical_features_norm, chemical_mean, chemical_std = + standardise_features(chemical_features) + +combined_features_norm, combined_mean, combined_std = + standardise_features(combined_features) + +println("\nNormalised structural features:") +println(structural_features_norm) + +println("\nNormalised chemistry-informed features:") +println(chemical_features_norm) + +println("\nNormalised combined features:") +println(combined_features_norm) + + +#9.store feature sets for later experiments + +#later regression experiments can compare: +#1. structural features only +#2. chemistry-informed features only +#3. combined features + +feature_sets = Dict( + :structural => structural_features_norm, + :chemical => chemical_features_norm, + :combined => combined_features_norm +) + +feature_names = Dict( + :structural => structural_feature_names, + :chemical => chemical_feature_names, + :combined => combined_feature_names +) + +println("\nAvailable feature sets:") +for key in sort(collect(keys(feature_sets))) + println(key, " => size ", size(feature_sets[key])) +end + + +#10.toy reaction level regression targets + +#these are placeholder values only. +#In the final project, they would be replaced by real: +#- energy barriers +#- reaction rates +#- reaction yields + +toy_energy_barriers = Float32[ + 12.5, + 18.2, + 7.4, + 21.0, + 15.6, + 25.3, + 10.2, + 30.1, + 13.8 +] + +toy_reaction_rates = Float32[ + 1.2, + 0.8, + 2.1, + 0.5, + 1.4, + 0.3, + 1.9, + 0.2, + 1.1 +] + +toy_reaction_yields = Float32[ + 65.0, + 58.0, + 72.0, + 49.0, + 61.0, + 38.0, + 70.0, + 31.0, + 67.0 +] + +reaction_targets = Dict( + :energy_barrier => toy_energy_barriers, + :reaction_rate => toy_reaction_rates, + :reaction_yield => toy_reaction_yields +) + +println("\nToy reaction-level regression targets:") +for key in sort(collect(keys(reaction_targets))) + println(key, " => ", reaction_targets[key]) +end + + +#11.package processed output + +#this makes it easier to reuse the data in later files + +processed_feature_data = ( + raw_reactions = raw_reactions, + parsed_reactions = parsed_reactions, + species = species, + node_to_id = node_to_id, + id_to_node = id_to_node, + source_matrix = source_matrix, + target_matrix = target_matrix, + structural_feature_names = structural_feature_names, + chemical_feature_names = chemical_feature_names, + combined_feature_names = combined_feature_names, + structural_features = structural_features, + chemical_features = chemical_features, + combined_features = combined_features, + structural_features_norm = structural_features_norm, + chemical_features_norm = chemical_features_norm, + combined_features_norm = combined_features_norm, + feature_sets = feature_sets, + feature_names = feature_names, + reaction_targets = reaction_targets +) + +println("\nProcessed feature data keys:") +println(keys(processed_feature_data)) + + +#12.final summary + + + +println("Summary") +println("Number of reactions: ", num_reactions) +println("Number of species: ", num_species) +println("Source matrix size: ", size(source_matrix)) +println("Target matrix size: ", size(target_matrix)) +println("Structural feature matrix size: ", size(structural_features)) +println("Chemical feature matrix size: ", size(chemical_features)) +println("Combined feature matrix size: ", size(combined_features)) +println("Normalised combined feature matrix size: ", size(combined_features_norm)) + +println("\nThis file prepares species-level structural and chemistry-informed") +println("features for later reaction-level regression experiments.") \ No newline at end of file diff --git a/examples/shonali_prototypes/14_bidirectional message passing.jl b/examples/shonali_prototypes/14_bidirectional message passing.jl new file mode 100644 index 0000000..3910d89 --- /dev/null +++ b/examples/shonali_prototypes/14_bidirectional message passing.jl @@ -0,0 +1,446 @@ +#bidirectional message passing + +#purpose: +#comparing simple directed hypergraph message-passing strategies for CRNs. + +#this file connects to the thesis question: +#"What message-passing strategy is most appropriate for directed hypergraphs representing chemical reaction networks?" + +#tt compares: +#1. Species -> Reaction message passing only +#2. Species -> Reaction -> Species -> Reaction bidirectional message passing + +using Statistics +using LinearAlgebra +using Random + +Random.seed!(42) + + +println("Bidirectional Message Passing") + + + +#1.defineing a small Formose-inspired CRN + + +raw_reactions = [ + "CH2O + CH2O -> Glycolaldehyde", + "Glycolaldehyde + CH2O -> Glyceraldehyde", + "Glyceraldehyde -> Dihydroxyacetone", + "Dihydroxyacetone + CH2O -> Tetrose", + "Tetrose -> Glycolaldehyde + Glycolaldehyde", + "Glyceraldehyde + Glycolaldehyde -> Pentose", + "Pentose -> Dihydroxyacetone + Glycolaldehyde", + "Tetrose + CH2O -> Hexose", + "Hexose -> Glyceraldehyde + Glyceraldehyde" +] + +println("Raw reactions:") +for r in raw_reactions + println(r) +end + +#2.parser + + +function parse_side(side::AbstractString) + species = split(strip(side), "+") + return strip.(species) +end + +function parse_reaction(reaction::AbstractString) + sides = split(reaction, "->") + + if length(sides) != 2 + error("Reaction must contain exactly one -> symbol: $reaction") + end + + reactants = parse_side(sides[1]) + products = parse_side(sides[2]) + + return reactants, products +end + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nParsed reactions:") +for (i, (reactants, products)) in enumerate(parsed_reactions) + println("r$i: ", join(reactants, " + "), " -> ", join(products, " + ")) +end + +#3.species mappings + + +all_species = String[] + +for (reactants, products) in parsed_reactions + append!(all_species, String.(reactants)) + append!(all_species, String.(products)) +end + +species = sort(unique(all_species)) + +node_to_id = Dict(s => i for (i, s) in enumerate(species)) +id_to_node = Dict(i => s for (s, i) in node_to_id) + +num_species = length(species) +num_reactions = length(parsed_reactions) + +println("\nSpecies:") +println(species) + +println("\nNode to ID mapping:") +println(node_to_id) + + +#4.source and target matrices + +#source_matrix[i, j] gives the stoichiometric count of species i as a reactant in reaction j. + +#target_matrix[i, j] gives the stoichiometric count of species i as a product in reaction j. + +source_matrix = zeros(Float32, num_species, num_reactions) +target_matrix = zeros(Float32, num_species, num_reactions) + +for (j, (reactants, products)) in enumerate(parsed_reactions) + for r in reactants + source_matrix[node_to_id[String(r)], j] += 1.0f0 + end + + for p in products + target_matrix[node_to_id[String(p)], j] += 1.0f0 + end +end + +membership_matrix = source_matrix .+ target_matrix + +println("\nSource matrix:") +println(source_matrix) + +println("\nTarget matrix:") +println(target_matrix) + +println("\nMembership matrix:") +println(membership_matrix) + + +#5.building species features + +#these are structural baseline features. +#they can later be replaced or extended using chemistry-informed features. + +source_count = vec(sum(source_matrix .> 0, dims = 2)) +target_count = vec(sum(target_matrix .> 0, dims = 2)) +participation_count = source_count .+ target_count + +source_stoich_sum = vec(sum(source_matrix, dims = 2)) +target_stoich_sum = vec(sum(target_matrix, dims = 2)) +total_stoich_sum = source_stoich_sum .+ target_stoich_sum + +species_features = hcat( + source_count, + target_count, + participation_count, + source_stoich_sum, + target_stoich_sum, + total_stoich_sum +) + +species_features = Float32.(species_features) + +species_feature_names = [ + "source_count", + "target_count", + "participation_count", + "source_stoich_sum", + "target_stoich_sum", + "total_stoich_sum" +] + +println("\nSpecies feature names:") +println(species_feature_names) + +println("\nSpecies features:") +println(species_features) + + +#6.normalise features + + +function standardise_features(X::AbstractMatrix) + μ = mean(X, dims = 1) + σ = std(X, dims = 1) + + σ_safe = similar(σ) + + for i in eachindex(σ) + σ_safe[i] = σ[i] == 0 ? 1 : σ[i] + end + + X_norm = (X .- μ) ./ σ_safe + + return Float32.(X_norm) +end + +species_features_norm = standardise_features(species_features) + +println("\nNormalised species features:") +println(species_features_norm) + + +#7.helper functions for aggregation + + +function safe_column_normalise(M::AbstractMatrix) + col_sums = sum(M, dims = 1) + + safe_sums = similar(col_sums) + + for i in eachindex(col_sums) + safe_sums[i] = col_sums[i] == 0 ? 1 : col_sums[i] + end + + return Float32.(M ./ safe_sums) +end + +function safe_row_normalise(M::AbstractMatrix) + row_sums = sum(M, dims = 2) + + safe_sums = similar(row_sums) + + for i in eachindex(row_sums) + safe_sums[i] = row_sums[i] == 0 ? 1 : row_sums[i] + end + + return Float32.(M ./ safe_sums) +end + + +#8.species -> reaction message passing + +#this aggregates species features into reaction embeddings. + +#reactants and products are treated separately, then concatenated. +#this allows the model representation to preserve directionality. + +function species_to_reaction( + X_species::AbstractMatrix, + source_matrix::AbstractMatrix, + target_matrix::AbstractMatrix +) + source_norm = safe_column_normalise(source_matrix) + target_norm = safe_column_normalise(target_matrix) + + source_reaction_features = transpose(source_norm) * X_species + target_reaction_features = transpose(target_norm) * X_species + + reaction_features = hcat(source_reaction_features, target_reaction_features) + + return Float32.(reaction_features) +end + +reaction_embeddings_oneway = + species_to_reaction(species_features_norm, source_matrix, target_matrix) + + +println("Architecture A: Species -> Reaction") + + +println("\nReaction embeddings from one-way message passing:") +println(reaction_embeddings_oneway) + +println("\nReaction embedding size:") +println(size(reaction_embeddings_oneway)) + + +#9.reaction -> species message passing + + +#this propagates reaction information back to species. + +#it allows species representations to be updated using information from the reactions they participate in. + +function reaction_to_species( + X_reaction::AbstractMatrix, + membership_matrix::AbstractMatrix +) + membership_norm = safe_row_normalise(membership_matrix) + + updated_species = membership_norm * X_reaction + + return Float32.(updated_species) +end + +updated_species_from_reactions = + reaction_to_species(reaction_embeddings_oneway, membership_matrix) + +println("\nUpdated species embeddings from Reaction -> Species:") +println(updated_species_from_reactions) + +println("\nUpdated species embedding size:") +println(size(updated_species_from_reactions)) + + +#10.bidirectional message passing + + +# Architecture B: +# +# Species features +# ↓ +# Species -> Reaction +# ↓ +# Reaction embeddings +# ↓ +# Reaction -> Species +# ↓ +# Updated species embeddings +# ↓ +# Species -> Reaction again +# +#this gives a second reaction representation after information has moved in both directions. + +reaction_embeddings_bidirectional = + species_to_reaction( + updated_species_from_reactions, + source_matrix, + target_matrix + ) + + +println("Architecture B: Species -> Reaction -> Species -> Reaction") + + +println("\nReaction embeddings from bidirectional message passing:") +println(reaction_embeddings_bidirectional) + +println("\nBidirectional reaction embedding size:") +println(size(reaction_embeddings_bidirectional)) + + +#11.toy regression target + + +#placeholder energy barriers. +#in the final project, these would come from a real chemical dataset. + +toy_energy_barriers = Float32[ + 12.5, + 18.2, + 7.4, + 21.0, + 15.6, + 25.3, + 10.2, + 30.1, + 13.8 +] + +println("\nToy energy barrier targets:") +println(toy_energy_barriers) + + +#12.simple ridge regression helper + + +#this is not the final neural network model. +#it is a simple baseline to check whether reaction embeddings can be used for reaction-level regression. + +# β = (X'X + λI)^(-1) X'y + +function add_intercept(X::AbstractMatrix) + ones_col = ones(Float32, size(X, 1), 1) + return hcat(ones_col, Float32.(X)) +end + +function fit_ridge_regression(X::AbstractMatrix, y::AbstractVector; λ = 1.0f-3) + X_aug = add_intercept(X) + + I_reg = Matrix{Float32}(I, size(X_aug, 2), size(X_aug, 2)) + I_reg[1, 1] = 0.0f0 # do not regularise intercept + + β = (transpose(X_aug) * X_aug + λ * I_reg) \ (transpose(X_aug) * y) + + return Float32.(β) +end + +function predict_ridge(X::AbstractMatrix, β::AbstractVector) + X_aug = add_intercept(X) + return Float32.(X_aug * β) +end + +function mae(y_true::AbstractVector, y_pred::AbstractVector) + return mean(abs.(y_true .- y_pred)) +end + +function rmse(y_true::AbstractVector, y_pred::AbstractVector) + return sqrt(mean((y_true .- y_pred).^2)) +end + + +#13.comparing one-way vs bidirectional embeddings + + +β_oneway = fit_ridge_regression(reaction_embeddings_oneway, toy_energy_barriers) +pred_oneway = predict_ridge(reaction_embeddings_oneway, β_oneway) + +β_bidirectional = fit_ridge_regression(reaction_embeddings_bidirectional, toy_energy_barriers) +pred_bidirectional = predict_ridge(reaction_embeddings_bidirectional, β_bidirectional) + + +println("Toy Regression Comparison") + + +println("\nOne-way predictions:") +println(pred_oneway) + +println("\nBidirectional predictions:") +println(pred_bidirectional) + +mae_oneway = mae(toy_energy_barriers, pred_oneway) +rmse_oneway = rmse(toy_energy_barriers, pred_oneway) + +mae_bidirectional = mae(toy_energy_barriers, pred_bidirectional) +rmse_bidirectional = rmse(toy_energy_barriers, pred_bidirectional) + +println("\nOne-way message passing MAE: ", mae_oneway) +println("One-way message passing RMSE: ", rmse_oneway) + +println("\nBidirectional message passing MAE: ", mae_bidirectional) +println("Bidirectional message passing RMSE: ", rmse_bidirectional) + + +#14.reaction-level output by reaction + + + +println("Reaction-level comparison") + + +for i in 1:num_reactions + reactants, products = parsed_reactions[i] + + println("\nr$i: ", join(reactants, " + "), " -> ", join(products, " + ")) + println("true energy barrier: ", toy_energy_barriers[i]) + println("one-way prediction: ", pred_oneway[i]) + println("bidirectional prediction: ", pred_bidirectional[i]) +end + + +#15.summary + + +println("\n====================================") +println("Summary") +println("====================================") + +println("Number of species: ", num_species) +println("Number of reactions: ", num_reactions) +println("Original species feature size: ", size(species_features_norm)) +println("One-way reaction embedding size: ", size(reaction_embeddings_oneway)) +println("Updated species embedding size: ", size(updated_species_from_reactions)) +println("Bidirectional reaction embedding size: ", size(reaction_embeddings_bidirectional)) + +println("\nThis file compares one-way and bidirectional message passing") +println("for reaction-level regression using toy energy barrier targets.") \ No newline at end of file diff --git a/examples/shonali_prototypes/15_reaction level regression.jl b/examples/shonali_prototypes/15_reaction level regression.jl new file mode 100644 index 0000000..d8a230a --- /dev/null +++ b/examples/shonali_prototypes/15_reaction level regression.jl @@ -0,0 +1,1043 @@ +#reaction level regression + +#purpose: +#this file evaluates whether reaction embeddings learned from the directed hypergraph neural network can predict continuous reaction +#properties such as energy barriers, reaction rates and reaction yields. + +#this file connects directly to the thesis question: +#"Which directed hypergraph neural network architecture is most suitable +#for reaction-level regression tasks involving chemical reaction networks?" + +#overall workflow: +# +#Chemical Reaction Network +# ↓ +#Species Features +# ↓ +#Species → Reaction Message Passing +# ↓ +#Reaction Embeddings +# ↓ +#Regression Model +# ↓ +#Predicted Reaction Property +# ↓ +#Model Evaluation + +using Statistics +using LinearAlgebra +using Random + +Random.seed!(42) + +println("Reaction Level Regression") + + +#1.defining a small Formose-inspired CRN + + +raw_reactions = [ + "CH2O + CH2O -> Glycolaldehyde", + "Glycolaldehyde + CH2O -> Glyceraldehyde", + "Glyceraldehyde -> Dihydroxyacetone", + "Dihydroxyacetone + CH2O -> Tetrose", + "Tetrose -> Glycolaldehyde + Glycolaldehyde", + "Glyceraldehyde + Glycolaldehyde -> Pentose", + "Pentose -> Dihydroxyacetone + Glycolaldehyde", + "Tetrose + CH2O -> Hexose", + "Hexose -> Glyceraldehyde + Glyceraldehyde" +] + +println("\nRaw reactions:") + +for reaction in raw_reactions + println(reaction) +end + + +#2.parser + + +function parse_side(side::AbstractString) + + species = split(strip(side), "+") + + return strip.(species) + +end + + +function parse_reaction(reaction::AbstractString) + + sides = split(reaction, "->") + + if length(sides) != 2 + error("Reaction must contain exactly one -> symbol.") + end + + reactants = parse_side(sides[1]) + products = parse_side(sides[2]) + + return reactants, products + +end + + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + + +println("\nParsed reactions:") + +for (i, (reactants, products)) in enumerate(parsed_reactions) + + println( + "r$i: ", + join(reactants, " + "), + " -> ", + join(products, " + ") + ) + +end + + +#3.species mappings + + +all_species = String[] + +for (reactants, products) in parsed_reactions + + append!(all_species, String.(reactants)) + append!(all_species, String.(products)) + +end + +species = sort(unique(all_species)) + +node_to_id = Dict(s => i for (i, s) in enumerate(species)) +id_to_node = Dict(i => s for (s, i) in node_to_id) + +num_species = length(species) +num_reactions = length(parsed_reactions) + +println("\nSpecies:") +println(species) + +println("\nNode to ID mapping:") +println(node_to_id) +#4.source and target matrices + +#source_matrix[i, j] gives the stoichiometric count of species i +#appearing as a reactant in reaction j. + +#target_matrix[i, j] gives the stoichiometric count of species i +#appearing as a product in reaction j. + +source_matrix = zeros(Float32, num_species, num_reactions) +target_matrix = zeros(Float32, num_species, num_reactions) + +for (j, (reactants, products)) in enumerate(parsed_reactions) + + for reactant in reactants + source_matrix[node_to_id[String(reactant)], j] += 1.0f0 + end + + for product in products + target_matrix[node_to_id[String(product)], j] += 1.0f0 + end + +end + +membership_matrix = source_matrix .+ target_matrix + +println("\nSource matrix:") +println(source_matrix) + +println("\nTarget matrix:") +println(target_matrix) + +println("\nMembership matrix:") +println(membership_matrix) + + +#5.building structural species features + + +#these are simple topology-based features. + +source_count = vec(sum(source_matrix .> 0, dims = 2)) +target_count = vec(sum(target_matrix .> 0, dims = 2)) +participation_count = source_count .+ target_count + +source_stoich_sum = vec(sum(source_matrix, dims = 2)) +target_stoich_sum = vec(sum(target_matrix, dims = 2)) +total_stoich_sum = source_stoich_sum .+ target_stoich_sum + +structural_features = hcat( + + source_count, + target_count, + participation_count, + source_stoich_sum, + target_stoich_sum, + total_stoich_sum + +) + +structural_features = Float32.(structural_features) + +println("\nStructural species features:") +println(structural_features) + + +#6.building chemistry-informed species features + + +#toy chemistry descriptors. + +#these can later be replaced by real molecular descriptors. + +carbon_atoms = Float32[ + 1, + 2, + 3, + 3, + 4, + 5, + 6 +] + +oxygen_atoms = Float32[ + 1, + 2, + 3, + 3, + 4, + 5, + 6 +] + +molecular_weight = Float32[ + 30.0, + 60.0, + 90.0, + 90.0, + 120.0, + 150.0, + 180.0 +] + +chemistry_features = hcat( + + carbon_atoms, + oxygen_atoms, + molecular_weight + +) + +println("\nChemistry-informed features:") +println(chemistry_features) + + +#7.combining structural and chemistry features + + +combined_features = hcat( + + structural_features, + chemistry_features + +) + +combined_features = Float32.(combined_features) + +println("\nCombined species features:") +println(combined_features) + + +#8.standardising feature matrices + + +function standardise_features(X::AbstractMatrix) + + μ = mean(X, dims = 1) + σ = std(X, dims = 1) + + σ_safe = similar(σ) + + for i in eachindex(σ) + + σ_safe[i] = σ[i] == 0 ? 1 : σ[i] + + end + + X_norm = (X .- μ) ./ σ_safe + + return Float32.(X_norm) + +end + +structural_features_norm = standardise_features(structural_features) + +chemistry_features_norm = standardise_features(chemistry_features) + +combined_features_norm = standardise_features(combined_features) + +println("\nNormalised structural features:") +println(structural_features_norm) + +println("\nNormalised chemistry features:") +println(chemistry_features_norm) + +println("\nNormalised combined features:") +println(combined_features_norm) +#9.helper functions for message passing + + +function safe_column_normalise(M::AbstractMatrix) + + col_sums = sum(M, dims = 1) + + safe_sums = similar(col_sums) + + for i in eachindex(col_sums) + + safe_sums[i] = col_sums[i] == 0 ? 1 : col_sums[i] + + end + + return Float32.(M ./ safe_sums) + +end + + +function safe_row_normalise(M::AbstractMatrix) + + row_sums = sum(M, dims = 2) + + safe_sums = similar(row_sums) + + for i in eachindex(row_sums) + + safe_sums[i] = row_sums[i] == 0 ? 1 : row_sums[i] + + end + + return Float32.(M ./ safe_sums) + +end + + +#10.species -> reaction message passing + + +#this aggregates species features into reaction embeddings. + +#reactant features and product features are aggregated +#separately before being concatenated. + +function species_to_reaction( + + X_species::AbstractMatrix, + source_matrix::AbstractMatrix, + target_matrix::AbstractMatrix + +) + + source_norm = safe_column_normalise(source_matrix) + target_norm = safe_column_normalise(target_matrix) + + reactant_embeddings = transpose(source_norm) * X_species + product_embeddings = transpose(target_norm) * X_species + + reaction_embeddings = hcat( + + reactant_embeddings, + product_embeddings + + ) + + return Float32.(reaction_embeddings) + +end + + +#11.select feature representation + + +#this allows easy comparison of different +#species feature representations. + +feature_set = "combined" + +species_features = combined_features_norm + +if feature_set == "structural" + + species_features = structural_features_norm + +elseif feature_set == "chemistry" + + species_features = chemistry_features_norm + +elseif feature_set == "combined" + + species_features = combined_features_norm + +else + + error("Unknown feature set.") + +end + + +println("\nSelected feature set:") +println(feature_set) + + +#12.construct reaction embeddings + + +reaction_embeddings = species_to_reaction( + + species_features, + source_matrix, + target_matrix + +) + +println("\nReaction embeddings:") +println(reaction_embeddings) + +println("\nReaction embedding size:") +println(size(reaction_embeddings)) + + +#13.select prediction task + + +#currently toy targets are used. + +#these can later be replaced by +#real reaction datasets. + +prediction_task = "Energy Barrier" + +reaction_targets = Float32[ + + 12.5, + 18.2, + 7.4, + 21.0, + 15.6, + 25.3, + 10.2, + 30.1, + 13.8 + +] + +println("\nPrediction task:") +println(prediction_task) + +println("\nReaction targets:") +println(reaction_targets) + + +#14.build regression dataset + + +#X contains the reaction embeddings. + +#rows represent reactions. + +#columns represent embedding dimensions. + +X = reaction_embeddings + +y = reaction_targets + +println("\nRegression dataset") + +println("Feature matrix size:") +println(size(X)) + +println("Target vector size:") +println(size(y)) + + +#15.shuffle dataset + + +indices = collect(1:num_reactions) + +Random.shuffle!(indices) + +X = X[indices, :] + +y = y[indices] + +println("\nShuffled reaction indices:") +println(indices) + + +#16.training and testing split + + +train_ratio = 0.80 + +num_train = floor(Int, train_ratio * num_reactions) + +train_indices = 1:num_train +test_indices = (num_train + 1):num_reactions + +X_train = X[train_indices, :] +X_test = X[test_indices, :] + +y_train = y[train_indices] +y_test = y[test_indices] + +println("\nTraining samples:") +println(length(y_train)) + +println("Testing samples:") +println(length(y_test)) + +#17.helper functions for ridge regression + + +#an intercept term is added to the regression model. + +function add_intercept(X::AbstractMatrix) + + ones_column = ones(Float32, size(X, 1), 1) + + return hcat(ones_column, Float32.(X)) + +end + + +#ridge regression solution: +# +#β = (X'X + λI)^(-1)X'y + +function fit_ridge_regression( + + X::AbstractMatrix, + y::AbstractVector; + + λ = 1.0f-3 + +) + + X_aug = add_intercept(X) + + I_reg = Matrix{Float32}(I, size(X_aug, 2), size(X_aug, 2)) + + I_reg[1,1] = 0.0f0 + + β = + + (transpose(X_aug) * X_aug + λ * I_reg) \ + + (transpose(X_aug) * y) + + return Float32.(β) + +end + + +function predict_ridge( + + X::AbstractMatrix, + β::AbstractVector + +) + + X_aug = add_intercept(X) + + return Float32.(X_aug * β) + +end + + +#18.regression evaluation metrics + + +function mse( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + return mean((y_true .- y_pred).^2) + +end + + +function rmse( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + return sqrt( + + mean( + + (y_true .- y_pred).^2 + + ) + + ) + +end + + +function mae( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + return mean( + + abs.(y_true .- y_pred) + + ) + +end + + +function r_squared( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + ss_res = + + sum( + + (y_true .- y_pred).^2 + + ) + + ss_tot = + + sum( + + (y_true .- mean(y_true)).^2 + + ) + + return 1.0 - ss_res / ss_tot + +end + + +#19.training the regression model + + +println("\nTraining ridge regression model...") + +β = fit_ridge_regression( + + X_train, + y_train + +) + +println("\nRegression coefficients:") + +println(β) + + +#20.generate predictions + + +train_predictions = + + predict_ridge( + + X_train, + β + + ) + +test_predictions = + + predict_ridge( + + X_test, + β + + ) + + +println("\nTraining predictions:") + +println(train_predictions) + +println("\nTesting predictions:") + +println(test_predictions) + + +#21.evaluate training performance + + +train_mse = + + mse( + + y_train, + train_predictions + + ) + +train_rmse = + + rmse( + + y_train, + train_predictions + + ) + +train_mae = + + mae( + + y_train, + train_predictions + + ) + +train_r2 = + + r_squared( + + y_train, + train_predictions + + ) + + +println("\nTraining Performance") + +println("MSE: ", train_mse) + +println("RMSE: ", train_rmse) + +println("MAE: ", train_mae) + +println("R²: ", train_r2) + + +#22.evaluate testing performance + + +test_mse = + + mse( + + y_test, + test_predictions + + ) + +test_rmse = + + rmse( + + y_test, + test_predictions + + ) + +test_mae = + + mae( + + y_test, + test_predictions + + ) + +test_r2 = + + r_squared( + + y_test, + test_predictions + + ) + + +println("\nTesting Performance") + +println("MSE: ", test_mse) + +println("RMSE: ", test_rmse) + +println("MAE: ", test_mae) + +println("R²: ", test_r2) + +#23.display reaction-level predictions + + +println("\n========================================") +println("Reaction-Level Predictions") +println("========================================") + +for i in 1:length(y_test) + + absolute_error = abs(y_test[i] - test_predictions[i]) + + println("\nReaction ", i) + + println("True ", prediction_task, ": ", y_test[i]) + + println("Predicted: ", test_predictions[i]) + + println("Absolute Error: ", absolute_error) + +end + + +#24.compare different feature sets + + +println("\n========================================") +println("Feature Set Comparison") +println("========================================") + + +feature_sets = Dict( + + "Structural" => structural_features_norm, + + "Chemistry" => chemistry_features_norm, + + "Combined" => combined_features_norm + +) + + +comparison_results = Dict{String, Dict{String, Float32}}() + + +for (feature_name, feature_matrix) in feature_sets + + println("\nRunning experiment using ", feature_name, " features...") + + + reaction_embeddings = species_to_reaction( + + feature_matrix, + + source_matrix, + + target_matrix + + ) + + + X = reaction_embeddings + + y = reaction_targets + + + X = X[indices, :] + y = y[indices] + + + X_train = X[train_indices, :] + X_test = X[test_indices, :] + + y_train = y[train_indices] + y_test = y[test_indices] + + + β = fit_ridge_regression( + + X_train, + + y_train + + ) + + + predictions = predict_ridge( + + X_test, + + β + + ) + + + current_mse = mse( + + y_test, + + predictions + + ) + + current_rmse = rmse( + + y_test, + + predictions + + ) + + current_mae = mae( + + y_test, + + predictions + + ) + + current_r2 = r_squared( + + y_test, + + predictions + + ) + + + comparison_results[feature_name] = Dict( + + "MSE" => Float32(current_mse), + + "RMSE" => Float32(current_rmse), + + "MAE" => Float32(current_mae), + + "R2" => Float32(current_r2) + + ) + + + println("MSE : ", current_mse) + println("RMSE: ", current_rmse) + println("MAE : ", current_mae) + println("R² : ", current_r2) + +end + + +#25.determine best feature representation + + +best_feature_set = "" + +best_rmse = Inf + +for (feature_name, metrics) in comparison_results + + global best_rmse + global best_feature_set + + if metrics["RMSE"] < best_rmse + + best_rmse = metrics["RMSE"] + best_feature_set = feature_name + + end + +end + + +println("\nBest feature representation:") + +println(best_feature_set) + + +#26.save experiment results + + +experiment_results = Dict( + + "prediction_task" => prediction_task, + + "selected_feature_set" => feature_set, + + "training_samples" => length(y_train), + + "testing_samples" => length(y_test), + + "mse" => test_mse, + + "rmse" => test_rmse, + + "mae" => test_mae, + + "r2" => test_r2, + + "predictions" => test_predictions, + + "ground_truth" => y_test, + + "best_feature_set" => best_feature_set, + + "comparison_results" => comparison_results + +) + + +println("\nExperiment results saved.") + + +#27.final summary + + +println("\n====================================") +println("Summary") +println("====================================") + +println("Prediction Task: ", prediction_task) + +println("Feature Set: ", feature_set) + +println("Training Samples: ", length(y_train)) + +println("Testing Samples: ", length(y_test)) + +println("MSE : ", test_mse) + +println("RMSE: ", test_rmse) + +println("MAE : ", test_mae) + +println("R² : ", test_r2) + +println("Best Feature Set: ", best_feature_set) + +println("\nThis file demonstrates reaction-level regression") + +println("using directed hypergraph reaction embeddings.") + +println("The learned reaction embeddings are evaluated") + +println("using a ridge regression baseline and compared") + +println("across structural, chemistry-informed and") + +println("combined feature representations.") + +println("\nThe outputs from this file will be used") + +println("in File 16 for architecture comparison.") + +println("\nReaction-level regression completed successfully.") \ No newline at end of file diff --git a/examples/shonali_prototypes/16_architecture_comparison.jl b/examples/shonali_prototypes/16_architecture_comparison.jl new file mode 100644 index 0000000..8afe1ff --- /dev/null +++ b/examples/shonali_prototypes/16_architecture_comparison.jl @@ -0,0 +1,919 @@ +#16_architecture_comparison + +#purpose: +#this file compares different directed hypergraph message-passing +#architectures for reaction-level regression. + +#this file directly answers the thesis question: +#"Which directed hypergraph neural network architecture is most suitable +#for reaction-level regression tasks involving chemical reaction networks?" + +#overall workflow: +# +#Chemical Reaction Network +# ↓ +#Species Features +# ↓ +#Different Message Passing Architectures +# ↓ +#Reaction Embeddings +# ↓ +#Regression Model +# ↓ +#Performance Evaluation +# ↓ +#Architecture Comparison + +using Statistics +using LinearAlgebra +using Random + +Random.seed!(42) + +println("Architecture Comparison") + + +#1.defining a small Formose-inspired CRN + + +raw_reactions = [ + + "CH2O + CH2O -> Glycolaldehyde", + "Glycolaldehyde + CH2O -> Glyceraldehyde", + "Glyceraldehyde -> Dihydroxyacetone", + "Dihydroxyacetone + CH2O -> Tetrose", + "Tetrose -> Glycolaldehyde + Glycolaldehyde", + "Glyceraldehyde + Glycolaldehyde -> Pentose", + "Pentose -> Dihydroxyacetone + Glycolaldehyde", + "Tetrose + CH2O -> Hexose", + "Hexose -> Glyceraldehyde + Glyceraldehyde" + +] + +println("\nRaw reactions:") + +for reaction in raw_reactions + + println(reaction) + +end + + +#2.parser + + +function parse_side(side::AbstractString) + + species = split(strip(side), "+") + + return strip.(species) + +end + + +function parse_reaction(reaction::AbstractString) + + sides = split(reaction, "->") + + if length(sides) != 2 + + error("Reaction must contain exactly one -> symbol.") + + end + + reactants = parse_side(sides[1]) + + products = parse_side(sides[2]) + + return reactants, products + +end + + +parsed_reactions = [parse_reaction(r) for r in raw_reactions] + +println("\nParsed reactions:") + +for (i,(reactants,products)) in enumerate(parsed_reactions) + + println( + + "r$i: ", + join(reactants," + "), + " -> ", + join(products," + ") + + ) + +end + + +#3.species mappings + + +all_species = String[] + +for (reactants,products) in parsed_reactions + + append!(all_species,String.(reactants)) + append!(all_species,String.(products)) + +end + + +species = sort(unique(all_species)) + +node_to_id = Dict(s => i for (i,s) in enumerate(species)) + +id_to_node = Dict(i => s for (s,i) in node_to_id) + +num_species = length(species) + +num_reactions = length(parsed_reactions) + + +println("\nSpecies:") + +println(species) + +println("\nNode to ID mapping:") + +println(node_to_id) + + +#4.source and target matrices + + +source_matrix = zeros(Float32,num_species,num_reactions) + +target_matrix = zeros(Float32,num_species,num_reactions) + + +for (j,(reactants,products)) in enumerate(parsed_reactions) + + for reactant in reactants + + source_matrix[node_to_id[String(reactant)],j] += 1.0f0 + + end + + for product in products + + target_matrix[node_to_id[String(product)],j] += 1.0f0 + + end + +end + + +membership_matrix = source_matrix .+ target_matrix + + +println("\nSource matrix:") + +println(source_matrix) + +println("\nTarget matrix:") + +println(target_matrix) + +println("\nMembership matrix:") + +println(membership_matrix) +#5.building structural species features + + +#these are simple topology-based features. + +source_count = vec(sum(source_matrix .> 0, dims = 2)) +target_count = vec(sum(target_matrix .> 0, dims = 2)) +participation_count = source_count .+ target_count + +source_stoich_sum = vec(sum(source_matrix, dims = 2)) +target_stoich_sum = vec(sum(target_matrix, dims = 2)) +total_stoich_sum = source_stoich_sum .+ target_stoich_sum + +structural_features = hcat( + + source_count, + target_count, + participation_count, + source_stoich_sum, + target_stoich_sum, + total_stoich_sum + +) + +structural_features = Float32.(structural_features) + +println("\nStructural species features:") +println(structural_features) + + +#6.building chemistry-informed species features + + +#toy chemistry descriptors. + +carbon_atoms = Float32[ + 1, + 2, + 3, + 3, + 4, + 5, + 6 +] + +oxygen_atoms = Float32[ + 1, + 2, + 3, + 3, + 4, + 5, + 6 +] + +molecular_weight = Float32[ + 30.0, + 60.0, + 90.0, + 90.0, + 120.0, + 150.0, + 180.0 +] + +chemistry_features = hcat( + + carbon_atoms, + oxygen_atoms, + molecular_weight + +) + +println("\nChemistry-informed features:") +println(chemistry_features) + + +#7.combined species features + + +combined_features = hcat( + + structural_features, + chemistry_features + +) + +combined_features = Float32.(combined_features) + +println("\nCombined species features:") +println(combined_features) + + +#8.standardise features + + +function standardise_features(X::AbstractMatrix) + + μ = mean(X, dims = 1) + σ = std(X, dims = 1) + + σ_safe = similar(σ) + + for i in eachindex(σ) + + σ_safe[i] = σ[i] == 0 ? 1 : σ[i] + + end + + X_norm = (X .- μ) ./ σ_safe + + return Float32.(X_norm) + +end + + +combined_features_norm = + standardise_features(combined_features) + +println("\nNormalised combined features:") +println(combined_features_norm) + + +#9.helper functions + + +function safe_column_normalise(M::AbstractMatrix) + + col_sums = sum(M, dims = 1) + + safe_sums = similar(col_sums) + + for i in eachindex(col_sums) + + safe_sums[i] = + col_sums[i] == 0 ? 1 : col_sums[i] + + end + + return Float32.(M ./ safe_sums) + +end + + +function safe_row_normalise(M::AbstractMatrix) + + row_sums = sum(M, dims = 2) + + safe_sums = similar(row_sums) + + for i in eachindex(row_sums) + + safe_sums[i] = + row_sums[i] == 0 ? 1 : row_sums[i] + + end + + return Float32.(M ./ safe_sums) + +end + + +#10.architecture A + + +#Species -> Reaction + +function species_to_reaction( + + X_species::AbstractMatrix, + source_matrix::AbstractMatrix, + target_matrix::AbstractMatrix + +) + + source_norm = + safe_column_normalise(source_matrix) + + target_norm = + safe_column_normalise(target_matrix) + + reactant_embeddings = + transpose(source_norm) * X_species + + product_embeddings = + transpose(target_norm) * X_species + + reaction_embeddings = + hcat( + reactant_embeddings, + product_embeddings + ) + + return Float32.(reaction_embeddings) + +end + + +#11.architecture B + + +#Species -> Reaction +# +#Reaction -> Species +# +#Species -> Reaction + +function reaction_to_species( + + X_reaction::AbstractMatrix, + membership_matrix::AbstractMatrix + +) + + membership_norm = + safe_row_normalise(membership_matrix) + + updated_species = + membership_norm * X_reaction + + return Float32.(updated_species) + +end + + +reaction_embeddings_A = + species_to_reaction( + + combined_features_norm, + + source_matrix, + + target_matrix + + ) + + +updated_species = + reaction_to_species( + + reaction_embeddings_A, + + membership_matrix + + ) + + +reaction_embeddings_B = + species_to_reaction( + + updated_species, + + source_matrix, + + target_matrix + + ) + + +println("\nArchitecture A embedding size:") +println(size(reaction_embeddings_A)) + +println("\nArchitecture B embedding size:") +println(size(reaction_embeddings_B)) +#12.toy regression target + + +#placeholder energy barrier values. +#these can later be replaced with +#real chemical reaction datasets. + +reaction_targets = Float32[ + + 12.5, + 18.2, + 7.4, + 21.0, + 15.6, + 25.3, + 10.2, + 30.1, + 13.8 + +] + +println("\nReaction targets:") +println(reaction_targets) + + +#13.build regression datasets + + +X_A = reaction_embeddings_A + +X_B = reaction_embeddings_B + +y = reaction_targets + +println("\nArchitecture A dataset size:") +println(size(X_A)) + +println("\nArchitecture B dataset size:") +println(size(X_B)) + + +#14.shuffle dataset + + +indices = collect(1:num_reactions) + +Random.shuffle!(indices) + +X_A = X_A[indices,:] +X_B = X_B[indices,:] + +y = y[indices] + +println("\nShuffled reaction indices:") +println(indices) + + +#15.train-test split + + +train_ratio = 0.80 + +num_train = floor(Int, train_ratio * num_reactions) + +train_indices = 1:num_train + +test_indices = (num_train + 1):num_reactions + + +X_A_train = X_A[train_indices,:] +X_A_test = X_A[test_indices,:] + +X_B_train = X_B[train_indices,:] +X_B_test = X_B[test_indices,:] + +y_train = y[train_indices] +y_test = y[test_indices] + + +println("\nTraining samples:") +println(length(y_train)) + +println("Testing samples:") +println(length(y_test)) + + +#16.ridge regression helper functions + + +function add_intercept(X::AbstractMatrix) + + ones_column = ones(Float32,size(X,1),1) + + return hcat(ones_column,Float32.(X)) + +end + + +function fit_ridge_regression( + + X::AbstractMatrix, + y::AbstractVector; + + λ = 1.0f-3 + +) + + X_aug = add_intercept(X) + + I_reg = Matrix{Float32}(I,size(X_aug,2),size(X_aug,2)) + + I_reg[1,1] = 0.0f0 + + β = + + (transpose(X_aug) * X_aug + λ * I_reg) \ + + (transpose(X_aug) * y) + + return Float32.(β) + +end + + +function predict_ridge( + + X::AbstractMatrix, + β::AbstractVector + +) + + X_aug = add_intercept(X) + + return Float32.(X_aug * β) + +end + + +#17.evaluation metrics + + +function mse( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + mean((y_true .- y_pred).^2) + +end + + +function rmse( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + sqrt(mean((y_true .- y_pred).^2)) + +end + + +function mae( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + mean(abs.(y_true .- y_pred)) + +end + + +function r_squared( + + y_true::AbstractVector, + y_pred::AbstractVector + +) + + ss_res = sum((y_true .- y_pred).^2) + + ss_tot = sum((y_true .- mean(y_true)).^2) + + return 1.0 - ss_res / ss_tot + +end + + +#18.train Architecture A + + +println("\nTraining Architecture A") + +beta_A = fit_ridge_regression( + + X_A_train, + y_train + +) + +predictions_A = predict_ridge( + + X_A_test, + beta_A + +) + + +#19.train Architecture B + + +println("\nTraining Architecture B") + +beta_B = fit_ridge_regression( + + X_B_train, + y_train + +) + +predictions_B = predict_ridge( + + X_B_test, + beta_B + +) + + +#20.evaluate Architecture A + + +mse_A = mse( + + y_test, + predictions_A + +) + +rmse_A = rmse( + + y_test, + predictions_A + +) + +mae_A = mae( + + y_test, + predictions_A + +) + +r2_A = r_squared( + + y_test, + predictions_A + +) + + +#21.evaluate Architecture B + + +mse_B = mse( + + y_test, + predictions_B + +) + +rmse_B = rmse( + + y_test, + predictions_B + +) + +mae_B = mae( + + y_test, + predictions_B + +) + +r2_B = r_squared( + + y_test, + predictions_B + +) + + +println("\nArchitecture A metrics") + +println("MSE : ",mse_A) +println("RMSE: ",rmse_A) +println("MAE : ",mae_A) +println("R² : ",r2_A) + + +println("\nArchitecture B metrics") + +println("MSE : ",mse_B) +println("RMSE: ",rmse_B) +println("MAE : ",mae_B) +println("R² : ",r2_B) +#22.architecture comparison table + + +println("\n========================================") +println("Architecture Comparison") +println("========================================") + +println(rpad("Architecture",45), + rpad("RMSE",12), + rpad("MAE",12), + "R²") + +println("-"^80) + +println(rpad("Species -> Reaction",45), + rpad(string(round(rmse_A,digits=4)),12), + rpad(string(round(mae_A,digits=4)),12), + round(r2_A,digits=4)) + +println(rpad("Species -> Reaction -> Species -> Reaction",45), + rpad(string(round(rmse_B,digits=4)),12), + rpad(string(round(mae_B,digits=4)),12), + round(r2_B,digits=4)) + + +#23.determine best architecture + + +best_architecture = "" + +lowest_rmse = Inf + +if rmse_A < lowest_rmse + + lowest_rmse = rmse_A + + best_architecture = "Species -> Reaction" + +end + +if rmse_B < lowest_rmse + + lowest_rmse = rmse_B + + best_architecture = + "Species -> Reaction -> Species -> Reaction" + +end + + +println("\nBest architecture:") + +println(best_architecture) + + +#24.save comparison results + + +architecture_results = Dict( + + "Architecture A" => Dict( + + "Name" => "Species -> Reaction", + + "MSE" => mse_A, + + "RMSE" => rmse_A, + + "MAE" => mae_A, + + "R2" => r2_A + + ), + + "Architecture B" => Dict( + + "Name" => "Species -> Reaction -> Species -> Reaction", + + "MSE" => mse_B, + + "RMSE" => rmse_B, + + "MAE" => mae_B, + + "R2" => r2_B + + ), + + "Best Architecture" => best_architecture + +) + +println("\nArchitecture comparison results saved.") + + +#25.final summary + + +println("\n====================================") +println("Summary") +println("====================================") + +println("Research Question:") + +println("Which directed hypergraph neural network") + +println("architecture performs best for") + +println("reaction-level regression?") + +println() + +println("Architecture A:") +println("Species -> Reaction") + +println("RMSE: ", rmse_A) +println("MAE : ", mae_A) +println("R² : ", r2_A) + +println() + +println("Architecture B:") +println("Species -> Reaction -> Species -> Reaction") + +println("RMSE: ", rmse_B) +println("MAE : ", mae_B) +println("R² : ", r2_B) + +println() + +println("Best Architecture:") + +println(best_architecture) + +println() + +println("Conclusion:") + +println("The architecture with the lowest RMSE") + +println("and MAE together with the highest R²") + +println("is selected as the preferred") + +println("directed hypergraph message-passing") + +println("architecture for reaction-level") + +println("regression on this chemical") + +println("reaction network.") + +println() + +println("This file provides the experimental") + +println("evidence required to answer the") + +println("primary thesis research question.") + +println() + diff --git a/examples/shonali_prototypes/17 lux directed hypergraph layer.jl b/examples/shonali_prototypes/17 lux directed hypergraph layer.jl new file mode 100644 index 0000000..a061c65 --- /dev/null +++ b/examples/shonali_prototypes/17 lux directed hypergraph layer.jl @@ -0,0 +1,259 @@ +#17 lux directed hypergraph layer + +#purpose: +#first Lux-compatible directed hypergraph message-passing layer. +#this version uses Lux v1.31.4 syntax. + +using Lux +using Random +using LinearAlgebra +using Statistics + +rng = Random.default_rng() +Random.seed!(rng, 42) + +println("Lux Directed Hypergraph Layer") + + +#1.toy species features and incidence matrices + +X_species = Float32[ + 1.0 0.0 2.0; + 0.0 1.0 1.0; + 1.0 1.0 0.0; + 2.0 0.0 1.0 +] + +source_matrix = Float32[ + 1 0 1; + 1 0 0; + 0 1 0; + 0 0 0 +] + +target_matrix = Float32[ + 0 0 0; + 0 1 0; + 1 0 0; + 0 0 1 +] + +println("\nInput species feature size:") +println(size(X_species)) + +println("\nSource matrix size:") +println(size(source_matrix)) + +println("\nTarget matrix size:") +println(size(target_matrix)) + + +#2.normalisation helpers + +function safe_column_normalise(M::AbstractMatrix) + col_sums = sum(M, dims = 1) + safe_sums = similar(col_sums) + + for i in eachindex(col_sums) + safe_sums[i] = col_sums[i] == 0 ? 1 : col_sums[i] + end + + return Float32.(M ./ safe_sums) +end + + +function safe_row_normalise(M::AbstractMatrix) + row_sums = sum(M, dims = 2) + safe_sums = similar(row_sums) + + for i in eachindex(row_sums) + safe_sums[i] = row_sums[i] == 0 ? 1 : row_sums[i] + end + + return Float32.(M ./ safe_sums) +end + + +#3.custom Lux layer + +struct DirectedHypergraphLayer <: Lux.AbstractLuxLayer + species_in_dim::Int + hidden_dim::Int + activation +end + + +#4.parameter initialisation + +function Lux.initialparameters(rng::AbstractRNG, layer::DirectedHypergraphLayer) + return ( + W_species = randn(rng, Float32, layer.species_in_dim, layer.hidden_dim) .* 0.1f0, + b_species = zeros(Float32, 1, layer.hidden_dim), + + W_reaction = randn(rng, Float32, 2 * layer.hidden_dim, layer.hidden_dim) .* 0.1f0, + b_reaction = zeros(Float32, 1, layer.hidden_dim), + + W_update = randn(rng, Float32, 2 * layer.hidden_dim, layer.hidden_dim) .* 0.1f0, + b_update = zeros(Float32, 1, layer.hidden_dim) + ) +end + + +function Lux.initialstates(rng::AbstractRNG, layer::DirectedHypergraphLayer) + return NamedTuple() +end + + +#5.forward pass + +function (layer::DirectedHypergraphLayer)(input, ps, st) + + X_species, source_matrix, target_matrix = input + + membership_matrix = source_matrix .+ target_matrix + + source_norm = safe_column_normalise(source_matrix) + target_norm = safe_column_normalise(target_matrix) + membership_norm = safe_row_normalise(membership_matrix) + + #learnable species transformation + + H_species = + layer.activation.( + X_species * ps.W_species .+ ps.b_species + ) + + #species -> reaction aggregation + + reactant_messages = + transpose(source_norm) * H_species + + product_messages = + transpose(target_norm) * H_species + + directed_reaction_input = + hcat( + reactant_messages, + product_messages + ) + + #learnable reaction transformation + + H_reaction = + layer.activation.( + directed_reaction_input * ps.W_reaction .+ ps.b_reaction + ) + + #reaction -> species aggregation + + species_messages = + membership_norm * H_reaction + + #species update + + species_update_input = + hcat( + H_species, + species_messages + ) + + updated_species = + layer.activation.( + species_update_input * ps.W_update .+ ps.b_update + ) + + output = ( + updated_species = updated_species, + reaction_embeddings = H_reaction + ) + + return output, st + +end + + +#6.create layer and initialise with Lux.setup + +layer = DirectedHypergraphLayer( + size(X_species, 2), + 8, + tanh +) + +ps, st = Lux.setup(rng, layer) + + +#7.forward pass through directed hypergraph layer + +output, st = layer( + ( + X_species, + source_matrix, + target_matrix + ), + ps, + st +) + +println("\nUpdated species embeddings:") +println(output.updated_species) + +println("\nReaction embeddings:") +println(output.reaction_embeddings) + +println("\nUpdated species embedding size:") +println(size(output.updated_species)) + +println("\nReaction embedding size:") +println(size(output.reaction_embeddings)) + + +#8.simple regression head + +regression_head = Lux.Dense(8 => 1) + +ps_head, st_head = Lux.setup(rng, regression_head) + +#Lux Dense expects features × batch +#reaction_embeddings is reactions × features +#so transpose is used + +predictions_matrix, st_head = + regression_head( + transpose(output.reaction_embeddings), + ps_head, + st_head + ) + +predictions = vec(predictions_matrix) + +println("\nReaction-level predictions:") +println(predictions) + + +#9.toy target and loss + +toy_energy_barriers = Float32[ + 12.5, + 18.2, + 7.4 +] + +function mse_loss(y_pred, y_true) + return mean((y_pred .- y_true).^2) +end + +loss = mse_loss(predictions, toy_energy_barriers) + +println("\nToy MSE loss:") +println(loss) + + +#10.summary + +println("\nSummary") +println("This file implements a Lux-compatible directed hypergraph layer.") +println("It contains learnable species, reaction and update transformations.") +println("It performs species-to-reaction and reaction-to-species message passing.") +println("It returns updated species embeddings and reaction embeddings.") +println("The reaction embeddings are passed into a Lux Dense regression head.") \ No newline at end of file From a3362fed5ed7dbde728320363b2e1de386fa1a74 Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Tue, 14 Jul 2026 13:12:36 +0100 Subject: [PATCH 2/7] Add Lux directed hypergraph layer Create layers directory, add reusable DirectedHypergraphLayer, regression model, and initial tests. --- src/HyperGraphNeuralNetworks.jl | 6 + src/layers/DirectedHypergraphLayer.jl | 126 +++++++++++++++++++++ src/layers/DirectedHypergraphRegression.jl | 98 ++++++++++++++++ test/layers/DirectedHypergraphLayer.jl | 47 ++++++++ test/runtests.jl | 1 + 5 files changed, 278 insertions(+) create mode 100644 src/layers/DirectedHypergraphLayer.jl create mode 100644 src/layers/DirectedHypergraphRegression.jl create mode 100644 test/layers/DirectedHypergraphLayer.jl diff --git a/src/HyperGraphNeuralNetworks.jl b/src/HyperGraphNeuralNetworks.jl index 309a889..cdec10f 100644 --- a/src/HyperGraphNeuralNetworks.jl +++ b/src/HyperGraphNeuralNetworks.jl @@ -12,13 +12,18 @@ using MLUtils using SimpleHypergraphs using SimpleDirectedHypergraphs + include("core/abstracttypes.jl") include("core/hypergraphs.jl") +include("layers/DirectedHypergraphLayer.jl") +include("layers/DirectedHypergraphRegression.jl") export AbstractHGNNHypergraph, AbstractHGNNDiHypergraph export HGNNHypergraph, HGNNDiHypergraph export add_vertex, add_vertices, remove_vertex, remove_vertices export add_hyperedge, add_hyperedges, remove_hyperedge, remove_hyperedges +export DirectedHypergraphLayer +export DirectedHypergraphRegression include("core/generate.jl") @@ -50,4 +55,5 @@ include("core/utils.jl") export check_num_vertices, check_num_hyperedges export normalize_graphdata + end \ No newline at end of file diff --git a/src/layers/DirectedHypergraphLayer.jl b/src/layers/DirectedHypergraphLayer.jl new file mode 100644 index 0000000..8fb121d --- /dev/null +++ b/src/layers/DirectedHypergraphLayer.jl @@ -0,0 +1,126 @@ +using Lux +using Random + + +function safe_column_normalise(M::AbstractMatrix) + col_sums = sum(M, dims = 1) + safe_sums = similar(col_sums) + + for i in eachindex(col_sums) + safe_sums[i] = col_sums[i] == 0 ? one(eltype(col_sums)) : col_sums[i] + end + + return M ./ safe_sums +end + +function safe_row_normalise(M::AbstractMatrix) + row_sums = sum(M, dims = 2) + safe_sums = similar(row_sums) + + for i in eachindex(row_sums) + safe_sums[i] = row_sums[i] == 0 ? one(eltype(row_sums)) : row_sums[i] + end + + return M ./ safe_sums +end + +struct DirectedHypergraphLayer{F} <: Lux.AbstractLuxLayer + species_in_dim::Int + hidden_dim::Int + activation::F +end + +function Lux.initialparameters( + rng::AbstractRNG, + layer::DirectedHypergraphLayer +) + return ( + W_species = randn( + rng, + Float32, + layer.species_in_dim, + layer.hidden_dim + ) .* 0.1f0, + + b_species = zeros( + Float32, + 1, + layer.hidden_dim + ), + + W_reaction = randn( + rng, + Float32, + 2 * layer.hidden_dim, + layer.hidden_dim + ) .* 0.1f0, + + b_reaction = zeros( + Float32, + 1, + layer.hidden_dim + ), + + W_update = randn( + rng, + Float32, + 2 * layer.hidden_dim, + layer.hidden_dim + ) .* 0.1f0, + + b_update = zeros( + Float32, + 1, + layer.hidden_dim + ) + ) +end + +Lux.initialstates( + ::AbstractRNG, + ::DirectedHypergraphLayer +) = NamedTuple() + +function (layer::DirectedHypergraphLayer)(input, ps, st) + X_species, source_matrix, target_matrix = input + + membership_matrix = source_matrix .+ target_matrix + + source_norm = safe_column_normalise(source_matrix) + target_norm = safe_column_normalise(target_matrix) + membership_norm = safe_row_normalise(membership_matrix) + + H_species = layer.activation.( + X_species * ps.W_species .+ ps.b_species + ) + + reactant_messages = + transpose(source_norm) * H_species + + product_messages = + transpose(target_norm) * H_species + + directed_reaction_input = + hcat(reactant_messages, product_messages) + + H_reaction = layer.activation.( + directed_reaction_input * ps.W_reaction .+ ps.b_reaction + ) + + species_messages = + membership_norm * H_reaction + + species_update_input = + hcat(H_species, species_messages) + + updated_species = layer.activation.( + species_update_input * ps.W_update .+ ps.b_update + ) + + output = ( + updated_species = updated_species, + reaction_embeddings = H_reaction + ) + + return output, st +end \ No newline at end of file diff --git a/src/layers/DirectedHypergraphRegression.jl b/src/layers/DirectedHypergraphRegression.jl new file mode 100644 index 0000000..9da0bf5 --- /dev/null +++ b/src/layers/DirectedHypergraphRegression.jl @@ -0,0 +1,98 @@ +# directed hypergraph regression model + +#combines a DirectedHypergraphLayer with a Lux Dense regression head. +#the model produces one continuous prediction for every reaction/hyperedge. + +struct DirectedHypergraphRegression{H, R} <: + Lux.AbstractLuxContainerLayer{(:hypergraph_layer, :regression_head)} + + hypergraph_layer::H + regression_head::R +end + + +""" + DirectedHypergraphRegression( + species_in_dim, + hidden_dim; + activation = tanh + ) + + +species features +→ directed hypergraph message passing +→ reaction embeddings +→ dense regression head +→ one scalar prediction per reaction +""" +function DirectedHypergraphRegression( + species_in_dim::Int, + hidden_dim::Int; + activation = tanh +) + hypergraph_layer = DirectedHypergraphLayer( + species_in_dim, + hidden_dim, + activation + ) + + regression_head = Lux.Dense(hidden_dim => 1) + + return DirectedHypergraphRegression( + hypergraph_layer, + regression_head + ) +end + + +# forward pass + +function (model::DirectedHypergraphRegression)(input, ps, st) + X_species, source_matrix, target_matrix = input + + # directed hypergraph message passing + + hypergraph_output, hypergraph_state = model.hypergraph_layer( + ( + X_species, + source_matrix, + target_matrix + ), + ps.hypergraph_layer, + st.hypergraph_layer + ) + + reaction_embeddings = + hypergraph_output.reaction_embeddings + + # Lux Dense expects features × batch. + # The reaction embeddings currently have shape: + # reactions × hidden features. + + reaction_embeddings_for_dense = + transpose(reaction_embeddings) + + prediction_matrix, regression_state = model.regression_head( + reaction_embeddings_for_dense, + ps.regression_head, + st.regression_head + ) + + # Convert the 1 × number_of_reactions matrix + # into a vector containing one prediction per reaction. + + predictions = vec(prediction_matrix) + + output = ( + predictions = predictions, + reaction_embeddings = reaction_embeddings, + updated_species = hypergraph_output.updated_species + ) + + new_state = ( + hypergraph_layer = hypergraph_state, + regression_head = regression_state + ) + + return output, new_state +end \ No newline at end of file diff --git a/test/layers/DirectedHypergraphLayer.jl b/test/layers/DirectedHypergraphLayer.jl new file mode 100644 index 0000000..503f6b3 --- /dev/null +++ b/test/layers/DirectedHypergraphLayer.jl @@ -0,0 +1,47 @@ +using Test +using Lux +using Random +using HyperGraphNeuralNetworks + +@testset "DirectedHypergraphLayer" begin + + rng = Random.default_rng() + + layer = DirectedHypergraphLayer(3, 8, tanh) + + ps, st = Lux.setup(rng, layer) + + X_species = Float32[ + 1.0 0.0 2.0; + 0.0 1.0 1.0; + 1.0 1.0 0.0; + 2.0 0.0 1.0 + ] + + source_matrix = Float32[ + 1 0 1; + 1 0 0; + 0 1 0; + 0 0 0 + ] + + target_matrix = Float32[ + 0 0 0; + 0 1 0; + 1 0 0; + 0 0 1 + ] + + output, st = layer( + (X_species, source_matrix, target_matrix), + ps, + st + ) + + @test size(output.updated_species) == (4, 8) + @test size(output.reaction_embeddings) == (3, 8) + + @test all(isfinite, output.updated_species) + @test all(isfinite, output.reaction_embeddings) + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index c39ad27..1a594b4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ using MLUtils using SimpleHypergraphs using SimpleDirectedHypergraphs using HyperGraphNeuralNetworks +include("layers/DirectedHypergraphLayer.jl") # Necessary for MLDatasets ENV["DATADEPS_ALWAYS_ACCEPT"] = true From 16fef4c55eecbc0e59b408b620af7fe6f7a4b6a3 Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Fri, 17 Jul 2026 16:22:40 +0100 Subject: [PATCH 3/7] Add documentation for directed hypergraph layers --- src/layers/DirectedHypergraphLayer.jl | 63 ++++++++++++++++++++++ src/layers/DirectedHypergraphRegression.jl | 36 +++++++++++++ 2 files changed, 99 insertions(+) diff --git a/src/layers/DirectedHypergraphLayer.jl b/src/layers/DirectedHypergraphLayer.jl index 8fb121d..6aa6152 100644 --- a/src/layers/DirectedHypergraphLayer.jl +++ b/src/layers/DirectedHypergraphLayer.jl @@ -1,6 +1,14 @@ using Lux using Random +""" + safe_column_normalise(M) + +Normalise each column of `M` by its column sum. + +Columns with a sum of zero use a denominator of one, preventing division by +zero while leaving those columns unchanged. +""" function safe_column_normalise(M::AbstractMatrix) col_sums = sum(M, dims = 1) @@ -12,6 +20,15 @@ function safe_column_normalise(M::AbstractMatrix) return M ./ safe_sums end +""" + safe_row_normalise(M) + +Normalise each row of `M` by its row sum. + +Rows with a sum of zero use a denominator of one, preventing division by +zero while leaving those rows unchanged. +""" + function safe_row_normalise(M::AbstractMatrix) row_sums = sum(M, dims = 2) @@ -24,6 +41,42 @@ function safe_row_normalise(M::AbstractMatrix) return M ./ safe_sums end +""" + DirectedHypergraphLayer(species_in_dim, hidden_dim, activation) + +A Lux-compatible message-passing layer for directed hypergraphs. + +The layer accepts a species-feature matrix together with source and target +incidence matrices. It performs: + +1. A learnable transformation of species features. +2. Separate aggregation of source and target species into reaction embeddings. +3. A learnable transformation of reaction embeddings. +4. Propagation of reaction messages back to participating species. +5. A learnable update of the species embeddings. + +# Arguments + +- `species_in_dim`: Number of input features associated with each species. +- `hidden_dim`: Size of the hidden species and reaction embeddings. +- `activation`: Element-wise activation function. + +# Input + +A tuple `(X_species, source_matrix, target_matrix)` where: + +- `X_species` has shape `number_of_species × species_in_dim`. +- `source_matrix` has shape `number_of_species × number_of_reactions`. +- `target_matrix` has shape `number_of_species × number_of_reactions`. + +# Output + +A named tuple containing: + +- `updated_species`: Updated species embeddings. +- `reaction_embeddings`: Learned reaction embeddings. +""" + struct DirectedHypergraphLayer{F} <: Lux.AbstractLuxLayer species_in_dim::Int hidden_dim::Int @@ -81,6 +134,16 @@ Lux.initialstates( ::DirectedHypergraphLayer ) = NamedTuple() +""" + (layer::DirectedHypergraphLayer)(input, ps, st) + +Apply one directed-hypergraph message-passing step. + +`ps` contains the learnable Lux parameters and `st` contains the layer state. +The returned state is unchanged because this layer currently has no mutable +state. +""" + function (layer::DirectedHypergraphLayer)(input, ps, st) X_species, source_matrix, target_matrix = input diff --git a/src/layers/DirectedHypergraphRegression.jl b/src/layers/DirectedHypergraphRegression.jl index 9da0bf5..e24d684 100644 --- a/src/layers/DirectedHypergraphRegression.jl +++ b/src/layers/DirectedHypergraphRegression.jl @@ -3,6 +3,7 @@ #combines a DirectedHypergraphLayer with a Lux Dense regression head. #the model produces one continuous prediction for every reaction/hyperedge. + struct DirectedHypergraphRegression{H, R} <: Lux.AbstractLuxContainerLayer{(:hypergraph_layer, :regression_head)} @@ -18,12 +19,22 @@ end activation = tanh ) +Construct a regression model that combines a +`DirectedHypergraphLayer` with a Lux `Dense` regression head. + +The model performs: species features → directed hypergraph message passing → reaction embeddings → dense regression head → one scalar prediction per reaction + +# Arguments + +- `species_in_dim`: Number of input features for each species. +- `hidden_dim`: Size of the hidden embeddings. +- `activation`: Activation function used in the hypergraph layer. """ function DirectedHypergraphRegression( species_in_dim::Int, @@ -46,6 +57,31 @@ end # forward pass +""" + (model::DirectedHypergraphRegression)(input, ps, st) + +Run a forward pass through the regression model. + +The model first computes reaction embeddings using the +`DirectedHypergraphLayer` and then applies a Dense regression head to +produce one scalar prediction for each reaction. + +# Arguments + +- `input`: Tuple `(X_species, source_matrix, target_matrix)` +- `ps`: Lux parameters +- `st`: Lux state + +# Returns + +A tuple containing: + +- `output`, a named tuple with: + - `predictions` + - `reaction_embeddings` + - `updated_species` +- `new_state`, the updated Lux state. +""" function (model::DirectedHypergraphRegression)(input, ps, st) X_species, source_matrix, target_matrix = input From 9c4ee79f66675a4ca707f0b1b4c6944a6b5f2c5b Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Wed, 22 Jul 2026 22:07:11 +0100 Subject: [PATCH 4/7] Generalise directed hypergraph layer --- src/HyperGraphNeuralNetworks.jl | 4 +- src/layers/DirectedHypergraphLayer.jl | 447 ++++++++++++++++----- src/layers/DirectedHypergraphRegression.jl | 134 ------ test/layers/DirectedHypergraphLayer.jl | 30 +- 4 files changed, 361 insertions(+), 254 deletions(-) delete mode 100644 src/layers/DirectedHypergraphRegression.jl diff --git a/src/HyperGraphNeuralNetworks.jl b/src/HyperGraphNeuralNetworks.jl index cdec10f..1574f3e 100644 --- a/src/HyperGraphNeuralNetworks.jl +++ b/src/HyperGraphNeuralNetworks.jl @@ -16,14 +16,14 @@ using SimpleDirectedHypergraphs include("core/abstracttypes.jl") include("core/hypergraphs.jl") include("layers/DirectedHypergraphLayer.jl") -include("layers/DirectedHypergraphRegression.jl") + export AbstractHGNNHypergraph, AbstractHGNNDiHypergraph export HGNNHypergraph, HGNNDiHypergraph export add_vertex, add_vertices, remove_vertex, remove_vertices export add_hyperedge, add_hyperedges, remove_hyperedge, remove_hyperedges export DirectedHypergraphLayer -export DirectedHypergraphRegression + include("core/generate.jl") diff --git a/src/layers/DirectedHypergraphLayer.jl b/src/layers/DirectedHypergraphLayer.jl index 6aa6152..fca5fa0 100644 --- a/src/layers/DirectedHypergraphLayer.jl +++ b/src/layers/DirectedHypergraphLayer.jl @@ -2,187 +2,418 @@ using Lux using Random """ - safe_column_normalise(M) + safe_normalise(M; dims) -Normalise each column of `M` by its column sum. +Normalise `M` along dimension `dims`. -Columns with a sum of zero use a denominator of one, preventing division by -zero while leaving those columns unchanged. +Use `dims = 1` for column-wise normalisation and `dims = 2` for row-wise +normalisation. Zero-sum rows or columns use a denominator of one to avoid +division by zero. """ - -function safe_column_normalise(M::AbstractMatrix) - col_sums = sum(M, dims = 1) - safe_sums = similar(col_sums) - - for i in eachindex(col_sums) - safe_sums[i] = col_sums[i] == 0 ? one(eltype(col_sums)) : col_sums[i] - end +function safe_normalise(M::AbstractMatrix; dims::Int) + dims in (1, 2) || + throw(ArgumentError("`dims` must be either 1 or 2.")) + + matrix_sums = sum(M; dims = dims) + safe_sums = ifelse.( + iszero.(matrix_sums), + one(eltype(matrix_sums)), + matrix_sums, + ) return M ./ safe_sums end -""" - safe_row_normalise(M) -Normalise each row of `M` by its row sum. -Rows with a sum of zero use a denominator of one, preventing division by -zero while leaving those rows unchanged. """ + DirectedHypergraphLayer( + vertex_in_dim, + hyperedge_in_dim, + hidden_dim; + activation = tanh, + normalize = true, + init_weight = Lux.glorot_uniform, + init_bias = Lux.zeros32, + ) +A general-purpose Lux-compatible message-passing layer for directed +hypergraphs. -function safe_row_normalise(M::AbstractMatrix) - row_sums = sum(M, dims = 2) - safe_sums = similar(row_sums) - - for i in eachindex(row_sums) - safe_sums[i] = row_sums[i] == 0 ? one(eltype(row_sums)) : row_sums[i] - end - - return M ./ safe_sums -end +The layer accepts vertex features, optional hyperedge features, and separate +source and target incidence matrices. Source-side and target-side vertex +representations are aggregated separately to preserve hyperedge direction. -""" - DirectedHypergraphLayer(species_in_dim, hidden_dim, activation) +# Arguments -A Lux-compatible message-passing layer for directed hypergraphs. +- `vertex_in_dim`: Number of input features for each vertex. +- `hyperedge_in_dim`: Number of input features for each hyperedge. Use `0` + when no initial hyperedge features are available. +- `hidden_dim`: Size of the updated vertex and hyperedge representations. +- `activation`: Element-wise activation function. +- `normalize`: Whether incidence matrices are normalised before aggregation. +- `init_weight`: Initialiser used for weight matrices. +- `init_bias`: Initialiser used for bias parameters. -The layer accepts a species-feature matrix together with source and target -incidence matrices. It performs: +# Input -1. A learnable transformation of species features. -2. Separate aggregation of source and target species into reaction embeddings. -3. A learnable transformation of reaction embeddings. -4. Propagation of reaction messages back to participating species. -5. A learnable update of the species embeddings. +When `hyperedge_in_dim == 0`: -# Arguments + (X_vertex, source_matrix, target_matrix) -- `species_in_dim`: Number of input features associated with each species. -- `hidden_dim`: Size of the hidden species and reaction embeddings. -- `activation`: Element-wise activation function. +When `hyperedge_in_dim > 0`: -# Input + (X_vertex, X_hyperedge, source_matrix, target_matrix) -A tuple `(X_species, source_matrix, target_matrix)` where: +Expected shapes: -- `X_species` has shape `number_of_species × species_in_dim`. -- `source_matrix` has shape `number_of_species × number_of_reactions`. -- `target_matrix` has shape `number_of_species × number_of_reactions`. +- `X_vertex`: `number_of_vertices × vertex_in_dim` +- `X_hyperedge`: `number_of_hyperedges × hyperedge_in_dim` +- `source_matrix`: `number_of_vertices × number_of_hyperedges` +- `target_matrix`: `number_of_vertices × number_of_hyperedges` # Output A named tuple containing: -- `updated_species`: Updated species embeddings. -- `reaction_embeddings`: Learned reaction embeddings. +- `updated_vertices` +- `updated_hyperedges` """ - -struct DirectedHypergraphLayer{F} <: Lux.AbstractLuxLayer - species_in_dim::Int +struct DirectedHypergraphLayer{F, IW, IB} <: Lux.AbstractLuxLayer + vertex_in_dim::Int + hyperedge_in_dim::Int hidden_dim::Int activation::F + normalize::Bool + init_weight::IW + init_bias::IB +end + + +function DirectedHypergraphLayer( + vertex_in_dim::Int, + hyperedge_in_dim::Int, + hidden_dim::Int; + activation = tanh, + normalize::Bool = true, + init_weight = Lux.glorot_uniform, + init_bias = Lux.zeros32, +) + vertex_in_dim > 0 || + throw(ArgumentError("`vertex_in_dim` must be positive.")) + + hyperedge_in_dim >= 0 || + throw(ArgumentError("`hyperedge_in_dim` cannot be negative.")) + + hidden_dim > 0 || + throw(ArgumentError("`hidden_dim` must be positive.")) + + return DirectedHypergraphLayer( + vertex_in_dim, + hyperedge_in_dim, + hidden_dim, + activation, + normalize, + init_weight, + init_bias, + ) +end + + +""" +Initialise a weight matrix for the row-major feature convention used by this +layer. + +Lux initialisers produce matrices with shape +`output_dimension × input_dimension`. This layer stores features as +`entities × features`, so the result is transposed to +`input_dimension × output_dimension`. +""" +function _initialise_weight( + initializer, + rng::AbstractRNG, + input_dimension::Int, + output_dimension::Int, +) + return permutedims( + initializer(rng, output_dimension, input_dimension), + ) +end + + +""" +Initialise a bias with shape `1 × output_dimension`. +""" +function _initialise_bias( + initializer, + rng::AbstractRNG, + output_dimension::Int, +) + return permutedims( + initializer(rng, output_dimension, 1), + ) end + function Lux.initialparameters( rng::AbstractRNG, - layer::DirectedHypergraphLayer + layer::DirectedHypergraphLayer, ) + hyperedge_update_in_dim = + 2 * layer.hidden_dim + layer.hyperedge_in_dim + return ( - W_species = randn( + W_vertex = _initialise_weight( + layer.init_weight, rng, - Float32, - layer.species_in_dim, - layer.hidden_dim - ) .* 0.1f0, - - b_species = zeros( - Float32, - 1, - layer.hidden_dim + layer.vertex_in_dim, + layer.hidden_dim, ), - W_reaction = randn( + b_vertex = _initialise_bias( + layer.init_bias, rng, - Float32, - 2 * layer.hidden_dim, - layer.hidden_dim - ) .* 0.1f0, + layer.hidden_dim, + ), - b_reaction = zeros( - Float32, - 1, - layer.hidden_dim + W_hyperedge = _initialise_weight( + layer.init_weight, + rng, + hyperedge_update_in_dim, + layer.hidden_dim, ), - W_update = randn( + b_hyperedge = _initialise_bias( + layer.init_bias, + rng, + layer.hidden_dim, + ), + + W_vertex_update = _initialise_weight( + layer.init_weight, rng, - Float32, 2 * layer.hidden_dim, - layer.hidden_dim - ) .* 0.1f0, + layer.hidden_dim, + ), - b_update = zeros( - Float32, - 1, - layer.hidden_dim - ) + b_vertex_update = _initialise_bias( + layer.init_bias, + rng, + layer.hidden_dim, + ), ) end + +function Lux.parameterlength(layer::DirectedHypergraphLayer) + hyperedge_update_in_dim = + 2 * layer.hidden_dim + layer.hyperedge_in_dim + + vertex_parameters = + layer.vertex_in_dim * layer.hidden_dim + + layer.hidden_dim + + hyperedge_parameters = + hyperedge_update_in_dim * layer.hidden_dim + + layer.hidden_dim + + vertex_update_parameters = + 2 * layer.hidden_dim * layer.hidden_dim + + layer.hidden_dim + + return ( + vertex_parameters + + hyperedge_parameters + + vertex_update_parameters + ) +end + + +# The layer has no running statistics or other mutable non-trainable values. Lux.initialstates( ::AbstractRNG, - ::DirectedHypergraphLayer + ::DirectedHypergraphLayer, ) = NamedTuple() +Lux.statelength(::DirectedHypergraphLayer) = 0 + + +function _unpack_input( + layer::DirectedHypergraphLayer, + input::Tuple{Any, Any, Any}, +) + layer.hyperedge_in_dim == 0 || + throw( + ArgumentError( + "Hyperedge features are required because " * + "`hyperedge_in_dim` is $(layer.hyperedge_in_dim).", + ), + ) + + X_vertex, source_matrix, target_matrix = input + + X_hyperedge = similar( + X_vertex, + size(source_matrix, 2), + 0, + ) + + return ( + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ) +end + + +function _unpack_input( + ::DirectedHypergraphLayer, + input::Tuple{Any, Any, Any, Any}, +) + return input +end + + +function _validate_inputs( + layer::DirectedHypergraphLayer, + X_vertex::AbstractMatrix, + X_hyperedge::AbstractMatrix, + source_matrix::AbstractMatrix, + target_matrix::AbstractMatrix, +) + size(source_matrix) == size(target_matrix) || + throw( + DimensionMismatch( + "`source_matrix` and `target_matrix` must have the same shape.", + ), + ) + + number_of_vertices, number_of_hyperedges = + size(source_matrix) + + size(X_vertex, 1) == number_of_vertices || + throw( + DimensionMismatch( + "The number of rows in `X_vertex` must match the number " * + "of vertices in the incidence matrices.", + ), + ) + + size(X_vertex, 2) == layer.vertex_in_dim || + throw( + DimensionMismatch( + "`X_vertex` has $(size(X_vertex, 2)) features, but the " * + "layer expects $(layer.vertex_in_dim).", + ), + ) + + size(X_hyperedge, 1) == number_of_hyperedges || + throw( + DimensionMismatch( + "The number of rows in `X_hyperedge` must match the number " * + "of hyperedges in the incidence matrices.", + ), + ) + + size(X_hyperedge, 2) == layer.hyperedge_in_dim || + throw( + DimensionMismatch( + "`X_hyperedge` has $(size(X_hyperedge, 2)) features, but " * + "the layer expects $(layer.hyperedge_in_dim).", + ), + ) + + return nothing +end + + """ (layer::DirectedHypergraphLayer)(input, ps, st) Apply one directed-hypergraph message-passing step. -`ps` contains the learnable Lux parameters and `st` contains the layer state. -The returned state is unchanged because this layer currently has no mutable -state. -""" +The layer first updates hyperedge representations using separate source-side +and target-side vertex aggregations. It then propagates the updated hyperedge +representations back to the vertices. +The state is returned unchanged because the layer contains no stateful +operations. +""" function (layer::DirectedHypergraphLayer)(input, ps, st) - X_species, source_matrix, target_matrix = input + ( + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ) = _unpack_input(layer, input) + + _validate_inputs( + layer, + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ) + + membership_matrix = + source_matrix .+ target_matrix + + if layer.normalize + source_used = + safe_normalise(source_matrix; dims = 1) - membership_matrix = source_matrix .+ target_matrix + target_used = + safe_normalise(target_matrix; dims = 1) - source_norm = safe_column_normalise(source_matrix) - target_norm = safe_column_normalise(target_matrix) - membership_norm = safe_row_normalise(membership_matrix) + membership_used = + safe_normalise(membership_matrix; dims = 2) + else + source_used = source_matrix + target_used = target_matrix + membership_used = membership_matrix + end - H_species = layer.activation.( - X_species * ps.W_species .+ ps.b_species + # Transform input vertex features. + H_vertex = layer.activation.( + X_vertex * ps.W_vertex .+ ps.b_vertex ) - reactant_messages = - transpose(source_norm) * H_species + # Aggregate source-side and target-side vertex information separately. + source_messages = + transpose(source_used) * H_vertex - product_messages = - transpose(target_norm) * H_species + target_messages = + transpose(target_used) * H_vertex - directed_reaction_input = - hcat(reactant_messages, product_messages) + # Combine directional messages with initial hyperedge features. + hyperedge_update_input = hcat( + source_messages, + target_messages, + X_hyperedge, + ) - H_reaction = layer.activation.( - directed_reaction_input * ps.W_reaction .+ ps.b_reaction + updated_hyperedges = layer.activation.( + hyperedge_update_input * ps.W_hyperedge .+ + ps.b_hyperedge ) - species_messages = - membership_norm * H_reaction + # Propagate hyperedge information back to participating vertices. + vertex_messages = + membership_used * updated_hyperedges - species_update_input = - hcat(H_species, species_messages) + vertex_update_input = hcat( + H_vertex, + vertex_messages, + ) - updated_species = layer.activation.( - species_update_input * ps.W_update .+ ps.b_update + updated_vertices = layer.activation.( + vertex_update_input * ps.W_vertex_update .+ + ps.b_vertex_update ) output = ( - updated_species = updated_species, - reaction_embeddings = H_reaction + updated_vertices = updated_vertices, + updated_hyperedges = updated_hyperedges, ) return output, st diff --git a/src/layers/DirectedHypergraphRegression.jl b/src/layers/DirectedHypergraphRegression.jl deleted file mode 100644 index e24d684..0000000 --- a/src/layers/DirectedHypergraphRegression.jl +++ /dev/null @@ -1,134 +0,0 @@ -# directed hypergraph regression model - -#combines a DirectedHypergraphLayer with a Lux Dense regression head. -#the model produces one continuous prediction for every reaction/hyperedge. - - -struct DirectedHypergraphRegression{H, R} <: - Lux.AbstractLuxContainerLayer{(:hypergraph_layer, :regression_head)} - - hypergraph_layer::H - regression_head::R -end - - -""" - DirectedHypergraphRegression( - species_in_dim, - hidden_dim; - activation = tanh - ) - -Construct a regression model that combines a -`DirectedHypergraphLayer` with a Lux `Dense` regression head. - -The model performs: - -species features -→ directed hypergraph message passing -→ reaction embeddings -→ dense regression head -→ one scalar prediction per reaction - -# Arguments - -- `species_in_dim`: Number of input features for each species. -- `hidden_dim`: Size of the hidden embeddings. -- `activation`: Activation function used in the hypergraph layer. -""" -function DirectedHypergraphRegression( - species_in_dim::Int, - hidden_dim::Int; - activation = tanh -) - hypergraph_layer = DirectedHypergraphLayer( - species_in_dim, - hidden_dim, - activation - ) - - regression_head = Lux.Dense(hidden_dim => 1) - - return DirectedHypergraphRegression( - hypergraph_layer, - regression_head - ) -end - - -# forward pass -""" - (model::DirectedHypergraphRegression)(input, ps, st) - -Run a forward pass through the regression model. - -The model first computes reaction embeddings using the -`DirectedHypergraphLayer` and then applies a Dense regression head to -produce one scalar prediction for each reaction. - -# Arguments - -- `input`: Tuple `(X_species, source_matrix, target_matrix)` -- `ps`: Lux parameters -- `st`: Lux state - -# Returns - -A tuple containing: - -- `output`, a named tuple with: - - `predictions` - - `reaction_embeddings` - - `updated_species` -- `new_state`, the updated Lux state. -""" - -function (model::DirectedHypergraphRegression)(input, ps, st) - X_species, source_matrix, target_matrix = input - - # directed hypergraph message passing - - hypergraph_output, hypergraph_state = model.hypergraph_layer( - ( - X_species, - source_matrix, - target_matrix - ), - ps.hypergraph_layer, - st.hypergraph_layer - ) - - reaction_embeddings = - hypergraph_output.reaction_embeddings - - # Lux Dense expects features × batch. - # The reaction embeddings currently have shape: - # reactions × hidden features. - - reaction_embeddings_for_dense = - transpose(reaction_embeddings) - - prediction_matrix, regression_state = model.regression_head( - reaction_embeddings_for_dense, - ps.regression_head, - st.regression_head - ) - - # Convert the 1 × number_of_reactions matrix - # into a vector containing one prediction per reaction. - - predictions = vec(prediction_matrix) - - output = ( - predictions = predictions, - reaction_embeddings = reaction_embeddings, - updated_species = hypergraph_output.updated_species - ) - - new_state = ( - hypergraph_layer = hypergraph_state, - regression_head = regression_state - ) - - return output, new_state -end \ No newline at end of file diff --git a/test/layers/DirectedHypergraphLayer.jl b/test/layers/DirectedHypergraphLayer.jl index 503f6b3..8112784 100644 --- a/test/layers/DirectedHypergraphLayer.jl +++ b/test/layers/DirectedHypergraphLayer.jl @@ -4,14 +4,19 @@ using Random using HyperGraphNeuralNetworks @testset "DirectedHypergraphLayer" begin - rng = Random.default_rng() - layer = DirectedHypergraphLayer(3, 8, tanh) + layer = DirectedHypergraphLayer( + 3, + 0, + 8; + activation = tanh, + normalize = true, + ) ps, st = Lux.setup(rng, layer) - X_species = Float32[ + X_vertex = Float32[ 1.0 0.0 2.0; 0.0 1.0 1.0; 1.0 1.0 0.0; @@ -32,16 +37,21 @@ using HyperGraphNeuralNetworks 0 0 1 ] - output, st = layer( - (X_species, source_matrix, target_matrix), + output, new_st = layer( + ( + X_vertex, + source_matrix, + target_matrix, + ), ps, - st + st, ) - @test size(output.updated_species) == (4, 8) - @test size(output.reaction_embeddings) == (3, 8) + @test size(output.updated_vertices) == (4, 8) + @test size(output.updated_hyperedges) == (3, 8) - @test all(isfinite, output.updated_species) - @test all(isfinite, output.reaction_embeddings) + @test all(isfinite, output.updated_vertices) + @test all(isfinite, output.updated_hyperedges) + @test new_st == st end \ No newline at end of file From bdd521423500e82fb32be148969df0bb6933f279 Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Wed, 22 Jul 2026 22:53:40 +0100 Subject: [PATCH 5/7] Expand DirectedHypergraphLayer tests --- test/layers/DirectedHypergraphLayer.jl | 378 +++++++++++++++++++++---- 1 file changed, 328 insertions(+), 50 deletions(-) diff --git a/test/layers/DirectedHypergraphLayer.jl b/test/layers/DirectedHypergraphLayer.jl index 8112784..deba805 100644 --- a/test/layers/DirectedHypergraphLayer.jl +++ b/test/layers/DirectedHypergraphLayer.jl @@ -3,55 +3,333 @@ using Lux using Random using HyperGraphNeuralNetworks + +const X_VERTEX = Float32[ + 1.0 0.0 2.0; + 0.0 1.0 1.0; + 1.0 1.0 0.0; + 2.0 0.0 1.0 +] + +const SOURCE_MATRIX = Float32[ + 1 0 1; + 1 0 0; + 0 1 0; + 0 0 0 +] + +const TARGET_MATRIX = Float32[ + 0 0 0; + 0 1 0; + 1 0 0; + 0 0 1 +] + + @testset "DirectedHypergraphLayer" begin - rng = Random.default_rng() - - layer = DirectedHypergraphLayer( - 3, - 0, - 8; - activation = tanh, - normalize = true, - ) - - ps, st = Lux.setup(rng, layer) - - X_vertex = Float32[ - 1.0 0.0 2.0; - 0.0 1.0 1.0; - 1.0 1.0 0.0; - 2.0 0.0 1.0 - ] - - source_matrix = Float32[ - 1 0 1; - 1 0 0; - 0 1 0; - 0 0 0 - ] - - target_matrix = Float32[ - 0 0 0; - 0 1 0; - 1 0 0; - 0 0 1 - ] - - output, new_st = layer( - ( - X_vertex, - source_matrix, - target_matrix, - ), - ps, - st, - ) - - @test size(output.updated_vertices) == (4, 8) - @test size(output.updated_hyperedges) == (3, 8) - - @test all(isfinite, output.updated_vertices) - @test all(isfinite, output.updated_hyperedges) - - @test new_st == st + + @testset "Basic forward pass" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 0, + 8; + activation = tanh, + normalize = true, + ) + + ps, st = Lux.setup(rng, layer) + + output, new_st = layer( + ( + X_VERTEX, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + + @test size(output.updated_vertices) == (4, 8) + @test size(output.updated_hyperedges) == (3, 8) + + @test all(isfinite, output.updated_vertices) + @test all(isfinite, output.updated_hyperedges) + + @test new_st == st + end + + + @testset "Normalisation enabled and disabled" begin + rng = Random.default_rng() + + layer_normalised = DirectedHypergraphLayer( + 3, + 0, + 8; + activation = tanh, + normalize = true, + ) + + layer_unnormalised = DirectedHypergraphLayer( + 3, + 0, + 8; + activation = tanh, + normalize = false, + ) + + ps_normalised, st_normalised = + Lux.setup(rng, layer_normalised) + + ps_unnormalised, st_unnormalised = + Lux.setup(rng, layer_unnormalised) + + output_normalised, _ = layer_normalised( + ( + X_VERTEX, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps_normalised, + st_normalised, + ) + + output_unnormalised, _ = layer_unnormalised( + ( + X_VERTEX, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps_unnormalised, + st_unnormalised, + ) + + @test size(output_normalised.updated_vertices) == (4, 8) + @test size(output_unnormalised.updated_vertices) == (4, 8) + + @test size(output_normalised.updated_hyperedges) == (3, 8) + @test size(output_unnormalised.updated_hyperedges) == (3, 8) + + @test all(isfinite, output_normalised.updated_vertices) + @test all(isfinite, output_unnormalised.updated_vertices) + end + + + @testset "Hyperedge input features" begin + rng = Random.default_rng() + + X_hyperedge = Float32[ + 1.0 0.0; + 0.0 1.0; + 1.0 1.0 + ] + + layer = DirectedHypergraphLayer( + 3, + 2, + 8; + activation = tanh, + normalize = true, + ) + + ps, st = Lux.setup(rng, layer) + + output, new_st = layer( + ( + X_VERTEX, + X_hyperedge, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + + @test size(output.updated_vertices) == (4, 8) + @test size(output.updated_hyperedges) == (3, 8) + + @test all(isfinite, output.updated_vertices) + @test all(isfinite, output.updated_hyperedges) + + @test new_st == st + end + + + @testset "Zero-sum normalisation" begin + M = Float32[ + 0 1 0; + 0 2 0; + 0 3 0 + ] + + column_normalised = + HyperGraphNeuralNetworks.safe_normalise( + M; + dims = 1, + ) + + row_normalised = + HyperGraphNeuralNetworks.safe_normalise( + M; + dims = 2, + ) + + @test all(isfinite, column_normalised) + @test all(isfinite, row_normalised) + + @test column_normalised[:, 1] == zeros(Float32, 3) + @test column_normalised[:, 3] == zeros(Float32, 3) + + @test row_normalised[1, :] == Float32[0, 1, 0] + @test row_normalised[2, :] == Float32[0, 1, 0] + @test row_normalised[3, :] == Float32[0, 1, 0] + end + + + @testset "Invalid normalisation dimension" begin + M = ones(Float32, 3, 3) + + @test_throws ArgumentError begin + HyperGraphNeuralNetworks.safe_normalise( + M; + dims = 3, + ) + end + end + + + @testset "Mismatched incidence matrices" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 0, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + wrong_target_matrix = zeros(Float32, 4, 2) + + @test_throws DimensionMismatch begin + layer( + ( + X_VERTEX, + SOURCE_MATRIX, + wrong_target_matrix, + ), + ps, + st, + ) + end + end + + + @testset "Incorrect number of vertex rows" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 0, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + wrong_X_vertex = rand(Float32, 5, 3) + + @test_throws DimensionMismatch begin + layer( + ( + wrong_X_vertex, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + end + end + + + @testset "Incorrect vertex feature dimension" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 0, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + wrong_X_vertex = rand(Float32, 4, 2) + + @test_throws DimensionMismatch begin + layer( + ( + wrong_X_vertex, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + end + end + + + @testset "Missing required hyperedge features" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 2, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + @test_throws ArgumentError begin + layer( + ( + X_VERTEX, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + end + end + + + @testset "Incorrect hyperedge feature dimensions" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 2, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + wrong_X_hyperedge = rand(Float32, 2, 2) + + @test_throws DimensionMismatch begin + layer( + ( + X_VERTEX, + wrong_X_hyperedge, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + end + end end \ No newline at end of file From 5190e6ed04813b0f1935d1ea236579ef44315f9c Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Thu, 30 Jul 2026 11:59:03 +0100 Subject: [PATCH 6/7] Improve DirectedHypergraphLayer test coverage Added additional unit tests for DirectedHypergraphLayer, including constructor validation, parameter/state initialization, forward-pass behaviour, and invalid input handling. All project tests pass successfully. --- .../DirectedHypergraphAttentionLayer.jl | 572 ++++++++++++++++++ test/layers/DirectedHypergraphLayer.jl | 112 +++- 2 files changed, 674 insertions(+), 10 deletions(-) create mode 100644 src/layers/DirectedHypergraphAttentionLayer.jl diff --git a/src/layers/DirectedHypergraphAttentionLayer.jl b/src/layers/DirectedHypergraphAttentionLayer.jl new file mode 100644 index 0000000..cd373aa --- /dev/null +++ b/src/layers/DirectedHypergraphAttentionLayer.jl @@ -0,0 +1,572 @@ +using Lux +using Random +using NNlib: leakyrelu + +""" + DirectedHypergraphAttentionLayer( + vertex_in_dim, + hyperedge_in_dim, + hidden_dim; + activation = tanh, + attention_activation = leakyrelu, + init_weight = Lux.glorot_uniform, + init_bias = Lux.zeros32, + return_attention = false, + ) + +A single-head, Lux-compatible attention layer for directed hypergraphs. + +The layer learns separate source-side and target-side attention coefficients +when aggregating vertex information into directed hyperedges. Updated +hyperedge representations are then propagated back to participating vertices. + +# Arguments + +- `vertex_in_dim`: Number of input features associated with each vertex. +- `hyperedge_in_dim`: Number of input features associated with each hyperedge. + Set this to `0` if no initial hyperedge features are available. +- `hidden_dim`: Size of the learned vertex and hyperedge representations. +- `activation`: Activation applied to updated representations. +- `attention_activation`: Activation applied to raw attention scores. +- `init_weight`: Initialiser used for weight parameters. +- `init_bias`: Initialiser used for bias parameters. +- `return_attention`: Whether the output should include the learned source and + target attention matrices. + +# Input + +When `hyperedge_in_dim == 0`, the layer accepts: + + (X_vertex, source_matrix, target_matrix) + +When `hyperedge_in_dim > 0`, the layer accepts: + + (X_vertex, X_hyperedge, source_matrix, target_matrix) + +Expected shapes: + +- `X_vertex`: `number_of_vertices × vertex_in_dim` +- `X_hyperedge`: `number_of_hyperedges × hyperedge_in_dim` +- `source_matrix`: `number_of_vertices × number_of_hyperedges` +- `target_matrix`: `number_of_vertices × number_of_hyperedges` + +# Output + +The layer returns a named tuple containing: + +- `updated_vertices` +- `updated_hyperedges` + +When `return_attention = true`, it additionally returns: + +- `source_attention` +- `target_attention` +""" +struct DirectedHypergraphAttentionLayer{F, A, IW, IB} <: + Lux.AbstractLuxLayer + vertex_in_dim::Int + hyperedge_in_dim::Int + hidden_dim::Int + activation::F + attention_activation::A + init_weight::IW + init_bias::IB + return_attention::Bool +end + + +function DirectedHypergraphAttentionLayer( + vertex_in_dim::Int, + hyperedge_in_dim::Int, + hidden_dim::Int; + activation = tanh, + attention_activation = leakyrelu, + init_weight = Lux.glorot_uniform, + init_bias = Lux.zeros32, + return_attention::Bool = false, +) + vertex_in_dim > 0 || + throw( + ArgumentError( + "`vertex_in_dim` must be positive.", + ), + ) + + hyperedge_in_dim >= 0 || + throw( + ArgumentError( + "`hyperedge_in_dim` cannot be negative.", + ), + ) + + hidden_dim > 0 || + throw( + ArgumentError( + "`hidden_dim` must be positive.", + ), + ) + + return DirectedHypergraphAttentionLayer( + vertex_in_dim, + hyperedge_in_dim, + hidden_dim, + activation, + attention_activation, + init_weight, + init_bias, + return_attention, + ) +end + + +function _attention_row_major_weight( + initializer, + rng::AbstractRNG, + input_dimension::Int, + output_dimension::Int, +) + return permutedims( + initializer( + rng, + output_dimension, + input_dimension, + ), + ) +end + + +function _attention_row_major_bias( + initializer, + rng::AbstractRNG, + output_dimension::Int, +) + return permutedims( + initializer( + rng, + output_dimension, + 1, + ), + ) +end + + +function Lux.initialparameters( + rng::AbstractRNG, + layer::DirectedHypergraphAttentionLayer, +) + hyperedge_update_in_dim = + 2 * layer.hidden_dim + layer.hyperedge_in_dim + + return ( + # Transform input vertex features. + W_vertex = _attention_row_major_weight( + layer.init_weight, + rng, + layer.vertex_in_dim, + layer.hidden_dim, + ), + + b_vertex = _attention_row_major_bias( + layer.init_bias, + rng, + layer.hidden_dim, + ), + + # Compute separate source-side and target-side attention scores. + a_source = _attention_row_major_weight( + layer.init_weight, + rng, + layer.hidden_dim, + 1, + ), + + a_target = _attention_row_major_weight( + layer.init_weight, + rng, + layer.hidden_dim, + 1, + ), + + # Update hyperedges from source messages, target messages, + # and optional hyperedge features. + W_hyperedge = _attention_row_major_weight( + layer.init_weight, + rng, + hyperedge_update_in_dim, + layer.hidden_dim, + ), + + b_hyperedge = _attention_row_major_bias( + layer.init_bias, + rng, + layer.hidden_dim, + ), + + # Update vertices using their existing hidden representation + # and messages received from hyperedges. + W_vertex_update = _attention_row_major_weight( + layer.init_weight, + rng, + 2 * layer.hidden_dim, + layer.hidden_dim, + ), + + b_vertex_update = _attention_row_major_bias( + layer.init_bias, + rng, + layer.hidden_dim, + ), + ) +end + + +Lux.initialstates( + ::AbstractRNG, + ::DirectedHypergraphAttentionLayer, +) = NamedTuple() + + +""" + masked_incidence_softmax(raw_scores, incidence_matrix) + +Convert one raw score per vertex into attention coefficients for each +vertex--hyperedge incidence. + +For every hyperedge column: + +- vertices outside the hyperedge receive attention weight zero; +- participating vertices receive softmax-normalised weights; +- the weights of participating vertices sum to one; +- an empty hyperedge column remains zero. +""" +function masked_incidence_softmax( + raw_scores::AbstractVector, + incidence_matrix::AbstractMatrix, +) + number_of_vertices, number_of_hyperedges = + size(incidence_matrix) + + length(raw_scores) == number_of_vertices || + throw( + DimensionMismatch( + "The number of raw attention scores must match " * + "the number of vertices in the incidence matrix.", + ), + ) + + score_matrix = + reshape(raw_scores, number_of_vertices, 1) .+ + zeros( + eltype(raw_scores), + number_of_vertices, + number_of_hyperedges, + ) + + membership_mask = + .!iszero.(incidence_matrix) + + negative_infinity = + convert(eltype(score_matrix), -Inf) + + masked_scores = + ifelse.( + membership_mask, + score_matrix, + negative_infinity, + ) + + column_has_members = + any(membership_mask, dims = 1) + + column_maximum = + maximum(masked_scores, dims = 1) + + safe_column_maximum = + ifelse.( + column_has_members, + column_maximum, + zero(eltype(column_maximum)), + ) + + exponentials = + ifelse.( + membership_mask, + exp.(score_matrix .- safe_column_maximum), + zero(eltype(score_matrix)), + ) + + denominators = + sum(exponentials, dims = 1) + + safe_denominators = + ifelse.( + iszero.(denominators), + one(eltype(denominators)), + denominators, + ) + + return exponentials ./ safe_denominators +end + + +""" + _safe_attention_row_normalise(M) + +Safely normalise each row of `M`. + +Rows with a sum of zero use a denominator of one, preventing division by zero +and leaving those rows unchanged. +""" +function _safe_attention_row_normalise( + M::AbstractMatrix, +) + row_sums = + sum(M, dims = 2) + + safe_row_sums = + ifelse.( + iszero.(row_sums), + one(eltype(row_sums)), + row_sums, + ) + + return M ./ safe_row_sums +end + + +function _unpack_attention_input( + layer::DirectedHypergraphAttentionLayer, + input::Tuple{Any, Any, Any}, +) + layer.hyperedge_in_dim == 0 || + throw( + ArgumentError( + "Hyperedge features are required because " * + "`hyperedge_in_dim` is $(layer.hyperedge_in_dim).", + ), + ) + + X_vertex, source_matrix, target_matrix = + input + + number_of_hyperedges = + size(source_matrix, 2) + + X_hyperedge = + similar( + X_vertex, + number_of_hyperedges, + 0, + ) + + return ( + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ) +end + + +function _unpack_attention_input( + ::DirectedHypergraphAttentionLayer, + input::Tuple{Any, Any, Any, Any}, +) + return input +end + + +function _validate_attention_inputs( + layer::DirectedHypergraphAttentionLayer, + X_vertex::AbstractMatrix, + X_hyperedge::AbstractMatrix, + source_matrix::AbstractMatrix, + target_matrix::AbstractMatrix, +) + size(source_matrix) == size(target_matrix) || + throw( + DimensionMismatch( + "`source_matrix` and `target_matrix` must have " * + "the same shape.", + ), + ) + + number_of_vertices, number_of_hyperedges = + size(source_matrix) + + size(X_vertex, 1) == number_of_vertices || + throw( + DimensionMismatch( + "The number of rows in `X_vertex` must match " * + "the number of vertices in the incidence matrices.", + ), + ) + + size(X_vertex, 2) == layer.vertex_in_dim || + throw( + DimensionMismatch( + "`X_vertex` has $(size(X_vertex, 2)) features, " * + "but the layer expects $(layer.vertex_in_dim).", + ), + ) + + size(X_hyperedge, 1) == number_of_hyperedges || + throw( + DimensionMismatch( + "The number of rows in `X_hyperedge` must match " * + "the number of hyperedges in the incidence matrices.", + ), + ) + + size(X_hyperedge, 2) == layer.hyperedge_in_dim || + throw( + DimensionMismatch( + "`X_hyperedge` has $(size(X_hyperedge, 2)) features, " * + "but the layer expects $(layer.hyperedge_in_dim).", + ), + ) + + return nothing +end + + +""" + (layer::DirectedHypergraphAttentionLayer)(input, ps, st) + +Apply one single-head attention-based directed-hypergraph message-passing +step. + +Source-side and target-side attention coefficients are calculated separately. +The updated hyperedge representations are then propagated back to the +participating vertices. +""" +function ( + layer::DirectedHypergraphAttentionLayer +)( + input, + ps, + st, +) + ( + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ) = _unpack_attention_input( + layer, + input, + ) + + _validate_attention_inputs( + layer, + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ) + + # Transform the original vertex features. + H_vertex = + layer.activation.( + X_vertex * ps.W_vertex .+ + ps.b_vertex + ) + + # Compute one source-side and one target-side score per vertex. + raw_source_scores = + vec( + layer.attention_activation.( + H_vertex * ps.a_source + ), + ) + + raw_target_scores = + vec( + layer.attention_activation.( + H_vertex * ps.a_target + ), + ) + + # Normalise scores only across vertices participating in each hyperedge. + source_attention = + masked_incidence_softmax( + raw_source_scores, + source_matrix, + ) + + target_attention = + masked_incidence_softmax( + raw_target_scores, + target_matrix, + ) + + # Attention-weighted vertex-to-hyperedge aggregation. + source_messages = + transpose(source_attention) * + H_vertex + + target_messages = + transpose(target_attention) * + H_vertex + + hyperedge_update_input = + hcat( + source_messages, + target_messages, + X_hyperedge, + ) + + updated_hyperedges = + layer.activation.( + hyperedge_update_input * + ps.W_hyperedge .+ + ps.b_hyperedge + ) + + # Propagate updated hyperedge information back to vertices. + membership_matrix = + source_matrix .+ + target_matrix + + membership_weights = + _safe_attention_row_normalise( + membership_matrix, + ) + + vertex_messages = + membership_weights * + updated_hyperedges + + vertex_update_input = + hcat( + H_vertex, + vertex_messages, + ) + + updated_vertices = + layer.activation.( + vertex_update_input * + ps.W_vertex_update .+ + ps.b_vertex_update + ) + + basic_output = ( + updated_vertices = updated_vertices, + updated_hyperedges = updated_hyperedges, + ) + + output = + if layer.return_attention + merge( + basic_output, + ( + source_attention = source_attention, + target_attention = target_attention, + ), + ) + else + basic_output + end + + return output, st +end \ No newline at end of file diff --git a/test/layers/DirectedHypergraphLayer.jl b/test/layers/DirectedHypergraphLayer.jl index deba805..95407d9 100644 --- a/test/layers/DirectedHypergraphLayer.jl +++ b/test/layers/DirectedHypergraphLayer.jl @@ -1,6 +1,7 @@ using Test using Lux using Random +using Enzyme using HyperGraphNeuralNetworks @@ -28,6 +29,27 @@ const TARGET_MATRIX = Float32[ @testset "DirectedHypergraphLayer" begin + @testset "Constructor validation" begin + @test_throws ArgumentError DirectedHypergraphLayer( + 0, + 0, + 8, + ) + + @test_throws ArgumentError DirectedHypergraphLayer( + 3, + -1, + 8, + ) + + @test_throws ArgumentError DirectedHypergraphLayer( + 3, + 0, + 0, + ) + end + + @testset "Basic forward pass" begin rng = Random.default_rng() @@ -61,6 +83,32 @@ const TARGET_MATRIX = Float32[ end + @testset "Parameter and state initialisation" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 2, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + @test size(ps.W_vertex) == (3, 8) + @test size(ps.b_vertex) == (1, 8) + + @test size(ps.W_hyperedge) == (18, 8) + @test size(ps.b_hyperedge) == (1, 8) + + @test size(ps.W_vertex_update) == (16, 8) + @test size(ps.b_vertex_update) == (1, 8) + + @test isempty(st) + @test Lux.statelength(layer) == 0 + @test Lux.parameterlength(layer) == 320 + end + + @testset "Normalisation enabled and disabled" begin rng = Random.default_rng() @@ -114,6 +162,9 @@ const TARGET_MATRIX = Float32[ @test all(isfinite, output_normalised.updated_vertices) @test all(isfinite, output_unnormalised.updated_vertices) + + @test all(isfinite, output_normalised.updated_hyperedges) + @test all(isfinite, output_unnormalised.updated_hyperedges) end @@ -179,12 +230,20 @@ const TARGET_MATRIX = Float32[ @test all(isfinite, column_normalised) @test all(isfinite, row_normalised) - @test column_normalised[:, 1] == zeros(Float32, 3) - @test column_normalised[:, 3] == zeros(Float32, 3) + @test column_normalised[:, 1] == + zeros(Float32, 3) + + @test column_normalised[:, 3] == + zeros(Float32, 3) + + @test row_normalised[1, :] == + Float32[0, 1, 0] - @test row_normalised[1, :] == Float32[0, 1, 0] - @test row_normalised[2, :] == Float32[0, 1, 0] - @test row_normalised[3, :] == Float32[0, 1, 0] + @test row_normalised[2, :] == + Float32[0, 1, 0] + + @test row_normalised[3, :] == + Float32[0, 1, 0] end @@ -211,7 +270,8 @@ const TARGET_MATRIX = Float32[ ps, st = Lux.setup(rng, layer) - wrong_target_matrix = zeros(Float32, 4, 2) + wrong_target_matrix = + zeros(Float32, 4, 2) @test_throws DimensionMismatch begin layer( @@ -238,7 +298,8 @@ const TARGET_MATRIX = Float32[ ps, st = Lux.setup(rng, layer) - wrong_X_vertex = rand(Float32, 5, 3) + wrong_X_vertex = + rand(Float32, 5, 3) @test_throws DimensionMismatch begin layer( @@ -265,7 +326,8 @@ const TARGET_MATRIX = Float32[ ps, st = Lux.setup(rng, layer) - wrong_X_vertex = rand(Float32, 4, 2) + wrong_X_vertex = + rand(Float32, 4, 2) @test_throws DimensionMismatch begin layer( @@ -306,7 +368,36 @@ const TARGET_MATRIX = Float32[ end - @testset "Incorrect hyperedge feature dimensions" begin + @testset "Incorrect number of hyperedge rows" begin + rng = Random.default_rng() + + layer = DirectedHypergraphLayer( + 3, + 2, + 8, + ) + + ps, st = Lux.setup(rng, layer) + + wrong_X_hyperedge = + rand(Float32, 2, 2) + + @test_throws DimensionMismatch begin + layer( + ( + X_VERTEX, + wrong_X_hyperedge, + SOURCE_MATRIX, + TARGET_MATRIX, + ), + ps, + st, + ) + end + end + + + @testset "Incorrect hyperedge feature dimension" begin rng = Random.default_rng() layer = DirectedHypergraphLayer( @@ -317,7 +408,8 @@ const TARGET_MATRIX = Float32[ ps, st = Lux.setup(rng, layer) - wrong_X_hyperedge = rand(Float32, 2, 2) + wrong_X_hyperedge = + rand(Float32, 3, 3) @test_throws DimensionMismatch begin layer( From 9ef5ef3ecec0d3b3b5ad54e11b39bcb622a52dcc Mon Sep 17 00:00:00 2001 From: shonalidixit Date: Fri, 31 Jul 2026 02:11:45 +0100 Subject: [PATCH 7/7] Implement DirectedHypergraphAttentionLayer for Lux ## Summary This PR introduces a new `DirectedHypergraphAttentionLayer` for Lux-compatible directed hypergraph neural networks. ### Features - Added a directed hypergraph attention layer supporting: - source and target attention mechanisms - optional hyperedge features - optional return of attention weights - Added helper functions for masked attention computation and safe normalization. - Added comprehensive unit tests covering: - constructor validation - parameter initialization - forward passes with and without hyperedge features - attention weight behaviour - deterministic zero-parameter behaviour - input validation and error handling All tests pass successfully. --- src/HyperGraphNeuralNetworks.jl | 2 + .../DirectedHypergraphAttentionLayer.jl | 78 ++- .../DirectedHypergraphAttentionLayer.jl | 550 ++++++++++++++++++ test/runtests.jl | 1 + 4 files changed, 617 insertions(+), 14 deletions(-) create mode 100644 test/layers/DirectedHypergraphAttentionLayer.jl diff --git a/src/HyperGraphNeuralNetworks.jl b/src/HyperGraphNeuralNetworks.jl index 1574f3e..bdfefec 100644 --- a/src/HyperGraphNeuralNetworks.jl +++ b/src/HyperGraphNeuralNetworks.jl @@ -16,6 +16,7 @@ using SimpleDirectedHypergraphs include("core/abstracttypes.jl") include("core/hypergraphs.jl") include("layers/DirectedHypergraphLayer.jl") +include("layers/DirectedHypergraphAttentionLayer.jl") export AbstractHGNNHypergraph, AbstractHGNNDiHypergraph @@ -23,6 +24,7 @@ export HGNNHypergraph, HGNNDiHypergraph export add_vertex, add_vertices, remove_vertex, remove_vertices export add_hyperedge, add_hyperedges, remove_hyperedge, remove_hyperedges export DirectedHypergraphLayer +export DirectedHypergraphAttentionLayer include("core/generate.jl") diff --git a/src/layers/DirectedHypergraphAttentionLayer.jl b/src/layers/DirectedHypergraphAttentionLayer.jl index cd373aa..9bb563d 100644 --- a/src/layers/DirectedHypergraphAttentionLayer.jl +++ b/src/layers/DirectedHypergraphAttentionLayer.jl @@ -45,10 +45,10 @@ When `hyperedge_in_dim > 0`, the layer accepts: Expected shapes: -- `X_vertex`: `number_of_vertices × vertex_in_dim` -- `X_hyperedge`: `number_of_hyperedges × hyperedge_in_dim` -- `source_matrix`: `number_of_vertices × number_of_hyperedges` -- `target_matrix`: `number_of_vertices × number_of_hyperedges` +- `X_vertex`: `number_of_vertices √ó vertex_in_dim` +- `X_hyperedge`: `number_of_hyperedges √ó hyperedge_in_dim` +- `source_matrix`: `number_of_vertices √ó number_of_hyperedges` +- `target_matrix`: `number_of_vertices √ó number_of_hyperedges` # Output @@ -220,11 +220,44 @@ function Lux.initialparameters( end +function Lux.parameterlength( + layer::DirectedHypergraphAttentionLayer, +) + hyperedge_update_in_dim = + 2 * layer.hidden_dim + layer.hyperedge_in_dim + + vertex_transform_parameters = + layer.vertex_in_dim * layer.hidden_dim + + layer.hidden_dim + + attention_parameters = + 2 * layer.hidden_dim + + hyperedge_update_parameters = + hyperedge_update_in_dim * layer.hidden_dim + + layer.hidden_dim + + vertex_update_parameters = + 2 * layer.hidden_dim * layer.hidden_dim + + layer.hidden_dim + + return ( + vertex_transform_parameters + + attention_parameters + + hyperedge_update_parameters + + vertex_update_parameters + ) +end + + +# The layer has no mutable non-trainable state. Lux.initialstates( ::AbstractRNG, ::DirectedHypergraphAttentionLayer, ) = NamedTuple() +Lux.statelength(::DirectedHypergraphAttentionLayer) = 0 + """ masked_incidence_softmax(raw_scores, incidence_matrix) @@ -376,6 +409,20 @@ function _unpack_attention_input( end +function _unpack_attention_input( + ::DirectedHypergraphAttentionLayer, + input::Tuple, +) + throw( + ArgumentError( + "The layer expects either a 3-tuple " * + "(X_vertex, source_matrix, target_matrix) or a 4-tuple " * + "(X_vertex, X_hyperedge, source_matrix, target_matrix).", + ), + ) +end + + function _validate_attention_inputs( layer::DirectedHypergraphAttentionLayer, X_vertex::AbstractMatrix, @@ -440,13 +487,7 @@ Source-side and target-side attention coefficients are calculated separately. The updated hyperedge representations are then propagated back to the participating vertices. """ -function ( - layer::DirectedHypergraphAttentionLayer -)( - input, - ps, - st, -) +function (layer::DirectedHypergraphAttentionLayer)(input, ps, st) ( X_vertex, X_hyperedge, @@ -524,9 +565,18 @@ function ( ) # Propagate updated hyperedge information back to vertices. + # A vertex participates in a hyperedge if it occurs on either side. + # Using a logical union avoids double-counting a vertex that is present + # in both the source and target incidence matrices. membership_matrix = - source_matrix .+ - target_matrix + convert.( + promote_type( + eltype(source_matrix), + eltype(target_matrix), + ), + (.!iszero.(source_matrix)) .| + (.!iszero.(target_matrix)), + ) membership_weights = _safe_attention_row_normalise( @@ -569,4 +619,4 @@ function ( end return output, st -end \ No newline at end of file +end diff --git a/test/layers/DirectedHypergraphAttentionLayer.jl b/test/layers/DirectedHypergraphAttentionLayer.jl new file mode 100644 index 0000000..6b05853 --- /dev/null +++ b/test/layers/DirectedHypergraphAttentionLayer.jl @@ -0,0 +1,550 @@ +using Test +using Random +using Lux +using HyperGraphNeuralNetworks + +const HGNN = HyperGraphNeuralNetworks + +@testset "DirectedHypergraphAttentionLayer" begin + + @testset "Constructor validation" begin + layer = HGNN.DirectedHypergraphAttentionLayer(3, 0, 4) + + @test layer isa HGNN.DirectedHypergraphAttentionLayer + @test layer.vertex_in_dim == 3 + @test layer.hyperedge_in_dim == 0 + @test layer.hidden_dim == 4 + @test layer.return_attention == false + + attention_layer = HGNN.DirectedHypergraphAttentionLayer( + 3, + 2, + 4; + activation = identity, + attention_activation = identity, + return_attention = true, + ) + + @test attention_layer.activation === identity + @test attention_layer.attention_activation === identity + @test attention_layer.return_attention + + @test_throws ArgumentError HGNN.DirectedHypergraphAttentionLayer(0, 0, 4) + @test_throws ArgumentError HGNN.DirectedHypergraphAttentionLayer(-1, 0, 4) + @test_throws ArgumentError HGNN.DirectedHypergraphAttentionLayer(3, -1, 4) + @test_throws ArgumentError HGNN.DirectedHypergraphAttentionLayer(3, 0, 0) + @test_throws ArgumentError HGNN.DirectedHypergraphAttentionLayer(3, 0, -1) + end + + @testset "Parameter and state initialisation" begin + rng = Random.Xoshiro(1234) + + layer = HGNN.DirectedHypergraphAttentionLayer(3, 2, 4) + ps, st = Lux.setup(rng, layer) + + @test size(ps.W_vertex) == (3, 4) + @test size(ps.b_vertex) == (1, 4) + @test size(ps.a_source) == (4, 1) + @test size(ps.a_target) == (4, 1) + @test size(ps.W_hyperedge) == (10, 4) + @test size(ps.b_hyperedge) == (1, 4) + @test size(ps.W_vertex_update) == (8, 4) + @test size(ps.b_vertex_update) == (1, 4) + + @test st == NamedTuple() + @test Lux.statelength(layer) == 0 + @test Lux.parameterlength(layer) == 104 + @test Lux.parameterlength(layer) == Lux.parameterlength(ps) + + layer2 = HGNN.DirectedHypergraphAttentionLayer(3, 0, 4) + ps2, _ = Lux.setup(Random.Xoshiro(1234), layer2) + + @test size(ps2.W_hyperedge) == (8, 4) + @test Lux.parameterlength(layer2) == 96 + end + + @testset "masked_incidence_softmax" begin + + raw_scores = Float32[0,0,0] + + incidence_matrix = Float32[ + 1 0 0 + 1 1 0 + 0 1 0 + ] + + attention = HGNN.masked_incidence_softmax( + raw_scores, + incidence_matrix, + ) + + expected = Float32[ + 0.5 0.0 0.0 + 0.5 0.5 0.0 + 0.0 0.5 0.0 + ] + + @test size(attention) == size(incidence_matrix) + @test isapprox(attention, expected) + + @test all(attention[incidence_matrix .== 0] .== 0) + + @test isapprox(sum(attention[:,1]),1.0f0) + @test isapprox(sum(attention[:,2]),1.0f0) + @test sum(attention[:,3]) == 0.0f0 + + @test all(isfinite, attention) + + unequal_scores = Float32[0,1,2] + + unequal_attention = HGNN.masked_incidence_softmax( + unequal_scores, + incidence_matrix, + ) + + @test unequal_attention[2,1] > unequal_attention[1,1] + @test unequal_attention[3,2] > unequal_attention[2,2] + + @test isapprox(sum(unequal_attention[:,1]),1.0f0) + @test isapprox(sum(unequal_attention[:,2]),1.0f0) + + weighted_incidence = Float32[ + 2 0 + 0 4 + 3 5 + ] + + weighted_attention = HGNN.masked_incidence_softmax( + Float32[0,0,0], + weighted_incidence, + ) + + @test isapprox( + weighted_attention, + Float32[ + 0.5 0.0 + 0.0 0.5 + 0.5 0.5 + ] + ) + + @test_throws DimensionMismatch HGNN.masked_incidence_softmax( + Float32[1,2], + incidence_matrix, + ) + end + + @testset "Safe row normalisation" begin + + matrix = Float32[ + 1 1 + 0 0 + 1 3 + ] + + normalised = HGNN._safe_attention_row_normalise(matrix) + + @test isapprox( + normalised, + Float32[ + 0.5 0.5 + 0.0 0.0 + 0.25 0.75 + ] + ) + + @test isapprox(sum(normalised[1,:]),1.0f0) + @test sum(normalised[2,:]) == 0.0f0 + @test isapprox(sum(normalised[3,:]),1.0f0) + + @test all(isfinite, normalised) + end + @testset "Forward pass without hyperedge features" begin + rng = Random.Xoshiro(2026) + + layer = HGNN.DirectedHypergraphAttentionLayer(3, 0, 4) + ps, st = Lux.setup(rng, layer) + + X_vertex = Float32[ + 1 0 2 + 0 1 1 + 2 1 0 + 1 1 1 + ] + + source_matrix = Float32[ + 1 0 + 1 1 + 0 1 + 0 0 + ] + + target_matrix = Float32[ + 0 1 + 0 0 + 1 0 + 0 0 + ] + + output, st_out = layer( + (X_vertex, source_matrix, target_matrix), + ps, + st, + ) + + @test haskey(output, :updated_vertices) + @test haskey(output, :updated_hyperedges) + @test !haskey(output, :source_attention) + @test !haskey(output, :target_attention) + + @test size(output.updated_vertices) == (4, 4) + @test size(output.updated_hyperedges) == (2, 4) + + @test eltype(output.updated_vertices) <: AbstractFloat + @test eltype(output.updated_hyperedges) <: AbstractFloat + + @test all(isfinite, output.updated_vertices) + @test all(isfinite, output.updated_hyperedges) + + @test st_out == st + + X_hyperedge = similar(X_vertex, 2, 0) + + explicit_output, explicit_st = layer( + (X_vertex, X_hyperedge, source_matrix, target_matrix), + ps, + st, + ) + + @test isapprox( + explicit_output.updated_vertices, + output.updated_vertices, + ) + + @test isapprox( + explicit_output.updated_hyperedges, + output.updated_hyperedges, + ) + + @test explicit_st == st + end + + @testset "Forward pass with hyperedge features" begin + rng = Random.Xoshiro(77) + + layer = HGNN.DirectedHypergraphAttentionLayer(3, 2, 5) + ps, st = Lux.setup(rng, layer) + + X_vertex = Float32[ + 1 0 2 + 0 1 1 + 2 1 0 + ] + + X_hyperedge = Float32[ + 1 0 + 0 1 + ] + + source_matrix = Float32[ + 1 0 + 1 1 + 0 1 + ] + + target_matrix = Float32[ + 0 1 + 0 0 + 1 0 + ] + + output, st_out = layer( + ( + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ), + ps, + st, + ) + + @test size(output.updated_vertices) == (3, 5) + @test size(output.updated_hyperedges) == (2, 5) + + @test all(isfinite, output.updated_vertices) + @test all(isfinite, output.updated_hyperedges) + + @test st_out == st + end + + @testset "Returned attention properties" begin + rng = Random.Xoshiro(9) + + layer = HGNN.DirectedHypergraphAttentionLayer( + 2, + 0, + 3; + return_attention = true, + ) + + ps, st = Lux.setup(rng, layer) + + X_vertex = Float32[ + 1 0 + 0 1 + 1 1 + 2 1 + ] + + source_matrix = Float32[ + 1 0 0 + 1 1 0 + 0 1 0 + 0 0 0 + ] + + target_matrix = Float32[ + 0 1 0 + 0 0 0 + 1 0 0 + 0 0 0 + ] + + output, _ = layer( + (X_vertex, source_matrix, target_matrix), + ps, + st, + ) + + @test haskey(output, :source_attention) + @test haskey(output, :target_attention) + + @test size(output.source_attention) == size(source_matrix) + @test size(output.target_attention) == size(target_matrix) + + @test all(output.source_attention[source_matrix .== 0] .== 0) + @test all(output.target_attention[target_matrix .== 0] .== 0) + + @test isapprox(sum(output.source_attention[:, 1]), 1.0f0) + @test isapprox(sum(output.source_attention[:, 2]), 1.0f0) + @test sum(output.source_attention[:, 3]) == 0.0f0 + + @test isapprox(sum(output.target_attention[:, 1]), 1.0f0) + @test isapprox(sum(output.target_attention[:, 2]), 1.0f0) + @test sum(output.target_attention[:, 3]) == 0.0f0 + + @test all(isfinite, output.source_attention) + @test all(isfinite, output.target_attention) + end + @testset "Deterministic zero-parameter behaviour" begin + layer = HGNN.DirectedHypergraphAttentionLayer( + 2, + 0, + 3; + return_attention = true, + ) + + ps, st = Lux.setup(Random.Xoshiro(1), layer) + + zero_ps = map(x -> zero.(x), ps) + + X_vertex = Float32[ + 1 2 + 3 4 + 5 6 + ] + + source_matrix = reshape( + Float32[ + 1, + 1, + 0, + ], + 3, + 1, + ) + + target_matrix = reshape( + Float32[ + 0, + 0, + 1, + ], + 3, + 1, + ) + + output, _ = layer( + (X_vertex, source_matrix, target_matrix), + zero_ps, + st, + ) + + @test output.updated_vertices == zeros(Float32, 3, 3) + @test output.updated_hyperedges == zeros(Float32, 1, 3) + + @test isapprox( + output.source_attention, + reshape(Float32[0.5, 0.5, 0.0], 3, 1), + ) + + @test isapprox( + output.target_attention, + reshape(Float32[0.0, 0.0, 1.0], 3, 1), + ) + end + + @testset "Input validation" begin + layer_without_hyperedge_features = + HGNN.DirectedHypergraphAttentionLayer(3, 0, 4) + + ps0, st0 = Lux.setup( + Random.Xoshiro(4), + layer_without_hyperedge_features, + ) + + X_vertex = ones(Float32, 3, 3) + + source_matrix = Float32[ + 1 0 + 1 1 + 0 1 + ] + + target_matrix = Float32[ + 0 1 + 0 0 + 1 0 + ] + + @test_throws DimensionMismatch layer_without_hyperedge_features( + ( + ones(Float32, 4, 3), + source_matrix, + target_matrix, + ), + ps0, + st0, + ) + + @test_throws DimensionMismatch layer_without_hyperedge_features( + ( + ones(Float32, 3, 2), + source_matrix, + target_matrix, + ), + ps0, + st0, + ) + + @test_throws DimensionMismatch layer_without_hyperedge_features( + ( + X_vertex, + source_matrix, + ones(Float32, 3, 3), + ), + ps0, + st0, + ) + + @test_throws DimensionMismatch layer_without_hyperedge_features( + ( + X_vertex, + ones(Float32, 4, 2), + ones(Float32, 4, 2), + ), + ps0, + st0, + ) + + @test_throws DimensionMismatch layer_without_hyperedge_features( + ( + X_vertex, + ones(Float32, 3, 1), + source_matrix, + target_matrix, + ), + ps0, + st0, + ) + + @test_throws ArgumentError layer_without_hyperedge_features( + ( + X_vertex, + source_matrix, + ), + ps0, + st0, + ) + + @test_throws ArgumentError layer_without_hyperedge_features( + ( + X_vertex, + zeros(Float32, 2, 0), + source_matrix, + target_matrix, + :extra, + ), + ps0, + st0, + ) + + layer_with_hyperedge_features = + HGNN.DirectedHypergraphAttentionLayer(3, 2, 4) + + ps2, st2 = Lux.setup( + Random.Xoshiro(5), + layer_with_hyperedge_features, + ) + + X_hyperedge = ones(Float32, 2, 2) + + @test_throws ArgumentError layer_with_hyperedge_features( + ( + X_vertex, + source_matrix, + target_matrix, + ), + ps2, + st2, + ) + + @test_throws DimensionMismatch layer_with_hyperedge_features( + ( + X_vertex, + ones(Float32, 3, 2), + source_matrix, + target_matrix, + ), + ps2, + st2, + ) + + @test_throws DimensionMismatch layer_with_hyperedge_features( + ( + X_vertex, + ones(Float32, 2, 1), + source_matrix, + target_matrix, + ), + ps2, + st2, + ) + + valid_output, valid_state = layer_with_hyperedge_features( + ( + X_vertex, + X_hyperedge, + source_matrix, + target_matrix, + ), + ps2, + st2, + ) + + @test size(valid_output.updated_vertices) == (3, 4) + @test size(valid_output.updated_hyperedges) == (2, 4) + @test valid_state == st2 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 1a594b4..b17b80a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,7 @@ using SimpleHypergraphs using SimpleDirectedHypergraphs using HyperGraphNeuralNetworks include("layers/DirectedHypergraphLayer.jl") +include("layers/DirectedHypergraphAttentionLayer.jl") # Necessary for MLDatasets ENV["DATADEPS_ALWAYS_ACCEPT"] = true