diff --git a/.gitignore b/.gitignore index d1458d2..bff5b40 100644 --- a/.gitignore +++ b/.gitignore @@ -36,3 +36,6 @@ benchmark/* /profile/* /statprof/* /debug/* + +# Compiled packages +JetReconstructionCompiled diff --git a/compile/Project.toml b/compile/Project.toml new file mode 100644 index 0000000..e93077e --- /dev/null +++ b/compile/Project.toml @@ -0,0 +1,11 @@ +[deps] +JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196" +PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" + +[sources] +JetReconstruction = {path = ".."} + +[compat] +JetReconstruction = "0.4" +PackageCompiler = "2" +julia = "1.11" diff --git a/compile/README.md b/compile/README.md new file mode 100644 index 0000000..838e3f1 --- /dev/null +++ b/compile/README.md @@ -0,0 +1,65 @@ +# Compiling JetReconstruction.jl to a C-library + +Minimal C bindings for JetReconstruction.jl + +- [C-header](JetReconstruction.h) +- shared library compiled with [PackageCompiler.jl](https://github.com/JuliaLang/PackageCompiler.jl) + +## Building library + +To build the library, run the following command from the package root directory: + +```sh +julia --project=compile compile/build.jl +``` + +## Usage example + +### Example source file + +Example usage of C bindings in an application: + +```C +#include "JetReconstruction.h" +#include "julia_init.h" /*Should be automatically generated by PackageCompiler.jl and distributed together with the "JetReconstruction.h" header file*/ + +int main(int argc, char *argv[]) { + init_julia(argc, argv); /*initialization of julia runtime*/ + + /*Prepare array of pseudo jets*/ + size_t particles_len; + jetreconstruction_PseudoJet* particles; + /*Initialize with desired values*/ + + /*Call jet reconstruction*/ + jetreconstruction_JetAlgorithm algorithm = JETRECONSTRUCTION_JETALGORITHM_CA; + double R = 3.0; + jetreconstruction_RecoStrategy strategy = JETRECONSTRUCTION_RECOSTRATEGY_BEST; + + jetreconstruction_ClusterSequence cluster_seq; + jetreconstruction_jet_reconstruct(particles, len, algorithm, R, strategy, + &cluster_seq); + + /*Use the cluster sequence in your application + then free memory allocations done by library*/ + jetreconstruction_ClusterSequence_free_members(&cluster_seq); + shutdown_julia(0); /*teardown of julia runtime*/ + return 0; +} + +``` + +### Example compilation + +To build an example application run the following command: + +```shell +cc -o jetreconstruction_test compile/test/jetreconstruction_test.c -IJetReconstructionCompiled/include -LJetReconstructionCompiled/lib -ljetreconstruction -ljulia +``` + +In case the compiled library resides in non-standard location, add its location to `LD_LIBRARY_PATH` when running example application: + +```shell +LD_LIBRARY_PATH=JetReconstructionCompiled/lib/:${LD_LIBRARY_PATH} ./jetreconstruction_test +``` + diff --git a/compile/build.jl b/compile/build.jl new file mode 100644 index 0000000..4508cf7 --- /dev/null +++ b/compile/build.jl @@ -0,0 +1,16 @@ +using PackageCompiler + +output_directory = joinpath(splitdir(@__DIR__) |> first, "JetReconstructionCompiled") + +@info "Creating library in $output_directory" +PackageCompiler.create_library(".", output_directory; + lib_name = "jetreconstruction", + header_files = [joinpath(@__DIR__, "include", + "JetReconstruction.h")], + # precompile_execution_file = [joinpath(@__DIR__, + # "precompile_execution.jl")], + # precompile_statements_file = [jointpath(@__DIR__, + # "precompile_statements.jl")], + incremental = false, + filter_stdlibs = true, + force = true) diff --git a/compile/include/JetReconstruction.h b/compile/include/JetReconstruction.h new file mode 100644 index 0000000..c9c79c7 --- /dev/null +++ b/compile/include/JetReconstruction.h @@ -0,0 +1,145 @@ +/* + * SPDX-FileCopyrightText: 2022-2024 JetReconstruction.jl authors, CERN + * SPDX-License-Identifier: MIT + */ + +#include + +/* +Enumeration representing different jet algorithms used in +the JetReconstruction module. +*/ +typedef enum { + JETRECONSTRUCTION_JETALGORITHM_ANTIKT = 0, /* The Anti-Kt algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_CA = 1, /* The Cambridge/Aachen algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_KT = 2, /* The Inclusive-Kt algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_GENKT = + 3, /* The Generalised Kt algorithm (with arbitrary power). */ + JETRECONSTRUCTION_JETALGORITHM_EEKT = + 4, /* The Generalised e+e- kt algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_DURHAM = + 5 /* The e+e- kt algorithm, aka Durham. */ +} jetreconstruction_JetAlgorithm; + +/* +Scoped enumeration (using EnumX) representing the different strategies for jet +reconstruction. +*/ +typedef enum { + JETRECONSTRUCTION_RECOSTRATEGY_BEST = 0, /* The best strategy. */ + JETRECONSTRUCTION_RECOSTRATEGY_N2PLAIN = 1, /* The plain N2 strategy. */ + JETRECONSTRUCTION_RECOSTRATEGY_N2TILTED = 2 /* The tiled N2 strategy. */ +} jetreconstruction_RecoStrategy; + +/* The `PseudoJet` struct represents a pseudojet, a four-momentum object used in +jet reconstruction algorithms. Additional information for the link back into the +history of the clustering is stored in the `_cluster_hist_index` field. There is +caching of the more expensive calculations for rapidity and azimuthal angle.*/ +typedef struct { + double px; /* The x-component of the momentum. */ + double py; /* The y-component of the momentum. */ + double pz; /* The z-component of the momentum. */ + double E; /* The energy component of the momentum. */ + long _cluster_hist_index; /* The index of the cluster history. */ + double _pt2; /* The squared transverse momentum. */ + double _inv_pt2; /* The inverse squared transverse momentum. */ + double _rap; /* The rapidity. */ + double _phi; /* The azimuthal angle. */ +} jetreconstruction_PseudoJet; + +int jetreconstruction_PseudoJet_init(jetreconstruction_PseudoJet *ptr, + double px, double py, double pz, double E); +/* +A struct holding a record of jet mergers and finalisations +*/ +typedef struct { + long parent1; /* Index in history where first parent of this jet was + created (NonexistentParent if this jet is an original + particle) */ + long parent2; /* Index in history where second parent of this jet was + created (NonexistentParent if this jet is an original + particle); BeamJet if this history entry just labels the + fact that the jet has recombined with the beam */ + long child; /* Index in history where the current jet is recombined with + another jet to form its child. It is Invalid + if this jet does not further recombine. */ + long jetp_index; /* Index in the jets vector where we will find the + PseudoJet object corresponding to this jet (i.e. the + jet created at this entry of the history). NB: if this + element of the history corresponds to a beam + recombination, then `jetp_index=Invalid`. */ + double dij; /* The distance corresponding to the recombination at this + stage of the clustering. */ + double max_dij_so_far; /* The largest recombination distance seen so far + in the clustering history */ +} jetreconstruction_HistoryElement; + +/* +A struct holding the full history of a jet clustering sequence, including the +final jets. +*/ +typedef struct { + jetreconstruction_JetAlgorithm + algorithm; /* The algorithm used for clustering. */ + double power; /* The power value used for the clustering algorithm (note that + this value is always stored as a Float64 to be + type stable) */ + double R; /* The R parameter used for the clustering algorithm. */ + jetreconstruction_RecoStrategy + strategy; /* The strategy used for clustering. */ + jetreconstruction_PseudoJet + *jets; /* The actual jets in the cluster sequence, which are of type `T + <: FourMomentum`. */ + size_t jets_length; /* Length of jets. */ + long n_initial_jets; /* The initial number of particles used for exclusive + jets. */ + jetreconstruction_HistoryElement + *history; /* The branching history of the cluster sequence. Each + stage in the history indicates where to look in the + jets vector to get the physical PseudoJet. */ + size_t history_length; /* Length of history. */ + double Qtot; /* The total energy of the event. */ +} jetreconstruction_ClusterSequence; + +void jetreconstruction_ClusterSequence_free_members_( + jetreconstruction_ClusterSequence *ptr); +static inline void jetreconstruction_ClusterSequence_free_members( + jetreconstruction_ClusterSequence *ptr) { + jetreconstruction_ClusterSequence_free_members_(ptr); + ptr->jets = NULL; + ptr->jets_length = 0; + ptr->history = NULL; + ptr->history_length = 0; +} + +int jetreconstruction_jet_reconstruct( + const jetreconstruction_PseudoJet *particles, size_t particles_length, + jetreconstruction_JetAlgorithm algorithm, double R, + jetreconstruction_RecoStrategy strategy, + jetreconstruction_ClusterSequence *result); + +typedef struct { + jetreconstruction_PseudoJet *data; + size_t length; +} jetreconstruction_JetsResult; + +void jetreconstruction_JetsResult_free_members_( + jetreconstruction_JetsResult *ptr); +static inline void +jetreconstruction_JetsResult_free_members(jetreconstruction_JetsResult *ptr) { + jetreconstruction_JetsResult_free_members_(ptr); + ptr->data = NULL; + ptr->length = 0; +} + +int jetreconstruction_exclusive_jets_dcut( + const jetreconstruction_ClusterSequence *clustersequence, double dcut, + jetreconstruction_JetsResult *result); + +int jetreconstruction_exclusive_jets_njets( + const jetreconstruction_ClusterSequence *clustersequence, size_t njets, + jetreconstruction_JetsResult *result); + +int jetreconstruction_inclusive_jets( + const jetreconstruction_ClusterSequence *clustersequence, double ptmin, + jetreconstruction_JetsResult *result); diff --git a/compile/precompile_execution.jl b/compile/precompile_execution.jl new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/compile/precompile_execution.jl @@ -0,0 +1 @@ + diff --git a/compile/precompile_statements.jl b/compile/precompile_statements.jl new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/compile/precompile_statements.jl @@ -0,0 +1 @@ + diff --git a/compile/test/jetreconstruction_test.c b/compile/test/jetreconstruction_test.c new file mode 100644 index 0000000..0a04917 --- /dev/null +++ b/compile/test/jetreconstruction_test.c @@ -0,0 +1,67 @@ +#include "JetReconstruction.h" +#include "julia_init.h" +#include +#include +#include + +void printPseudoJet(const jetreconstruction_PseudoJet *jet) { + assert(jet != NULL); + printf("PseudoJet(%f %f %f %f %ld %f %f %f %f)\n", + jet->px, jet->py, jet->pz, jet->E, jet->_cluster_hist_index, jet->_pt2, + jet->_inv_pt2, jet->_rap, jet->_phi); +} + +void printHistoryElement(const jetreconstruction_HistoryElement *history) { + assert(history != NULL); + printf("HistoryElement(%ld %ld %ld %ld %lf %lf)\n", + history->parent1, history->parent2, history->child, + history->jetp_index, history->dij, history->max_dij_so_far); +} + +void printClusterSequence(const jetreconstruction_ClusterSequence *sequence) { + printf("Cluster Sequence Information:\n" + "Algorithm: %d\n" + "Power: %f\n" + "R parameter: %f\n" + "Strategy: %d\n" + "Initial number of jets: %ld\n" + "Total event energy (Qtot): %f\n", + sequence->algorithm, sequence->power, sequence->R, sequence->strategy, + sequence->n_initial_jets, sequence->Qtot); + printf("Number of jets: %zu\n", sequence->jets_length); + if (sequence->jets != NULL) { + for (size_t i = 0; i < sequence->jets_length; i++) { + printPseudoJet(sequence->jets + i); + } + } + printf("History length: %zu\n", sequence->history_length); + if (sequence->history != NULL) { + for (size_t i = 0; i < sequence->history_length; i++) { + printHistoryElement(sequence->history + i); + } + } +} + +int main(int argc, char *argv[]) { + init_julia(argc, argv); + size_t len = 2; + jetreconstruction_PseudoJet particles[2]; + jetreconstruction_PseudoJet_init(&particles[0], 0.0, 1.0, 2.0, 3.0); + jetreconstruction_PseudoJet_init(&particles[1], 1.0, 2.0, 3.0, 4.0); + + jetreconstruction_JetAlgorithm algorithm = JETRECONSTRUCTION_JETALGORITHM_CA; + double R = 3.0; + jetreconstruction_RecoStrategy strategy = JETRECONSTRUCTION_RECOSTRATEGY_BEST; + + jetreconstruction_ClusterSequence cluster_seq; + jetreconstruction_jet_reconstruct(particles, len, algorithm, R, strategy, + &cluster_seq); + + printClusterSequence(&cluster_seq); + jetreconstruction_JetsResult result; + jetreconstruction_exclusive_jets_njets(&cluster_seq, 2, &result); + jetreconstruction_JetsResult_free_members(&result); + jetreconstruction_ClusterSequence_free_members(&cluster_seq); + shutdown_julia(0); + return 0; +} diff --git a/src/C_JetReconstruction/C_JetReconstruction.jl b/src/C_JetReconstruction/C_JetReconstruction.jl new file mode 100644 index 0000000..4c12334 --- /dev/null +++ b/src/C_JetReconstruction/C_JetReconstruction.jl @@ -0,0 +1,191 @@ +""" +Minimal C bindings for JetReconstruction.jl +""" +module C_JetReconstruction + +using ..JetReconstruction: JetAlgorithm, RecoStrategy, PseudoJet, ClusterSequence, + HistoryElement, jet_reconstruct, + exclusive_jets, inclusive_jets + +""" + unsafe_wrap_c_array(ptr::Ptr{T}, array_length::Csize_t) where {T} + +Wraps a C array into a Julia `Vector` for both bits and non-bits types. + +# Arguments +- `ptr::Ptr{T}`: A pointer to the C array. +- `array_length::Csize_t`: The length of the C array. + +# Returns +- A Julia `Vector{T}` containing the elements of the C array. + +# Safety +This function use 'unsafe' methods and has undefined behaviour +if pointer isn't valid or length isn't correct. +""" +function unsafe_wrap_c_array(ptr::Ptr{T}, array_length::Csize_t) where {T} + if isbitstype(T) + return unsafe_wrap(Vector{T}, ptr, array_length) + end + + vec = Vector{T}(undef, array_length) + for i in eachindex(vec) + @inbounds vec[i] = unsafe_load(ptr, i) + end + return vec +end + +""" + make_c_array(v::Vector{T}) where {T} + +Helper function for converting a Julia vector to a C-style array. +A C-style array is dynamically allocated and contents of input vector copied to it. + +# Arguments +- `v::Vector{T}`: A Julia vector of type `T`. + +# Returns +- `ptr::Ptr{T}`: A pointer to the allocated C-style array. +- `len::Int`: The length of the vector. + +# Throws +- `OutOfMemoryError`: If memory allocation fails. + +# Notes +- The caller is responsible for freeing the allocated memory using `Libc.free(ptr)`. +""" +function make_c_array(v::Vector{T}) where {T} + len = length(v) + ptr = Ptr{T}(Libc.malloc(len * sizeof(T))) + if ptr == C_NULL + throw(OutOfMemoryError("Libc.malloc failed to allocate memory")) + end + for i in 1:len + @inbounds unsafe_store!(ptr, v[i], i) + end + return ptr, len +end + +Base.@ccallable function jetreconstruction_PseudoJet_init(ptr::Ptr{PseudoJet}, px::Cdouble, + py::Cdouble, pz::Cdouble, + E::Cdouble)::Cint + pseudojet = PseudoJet(px, py, pz, E) + unsafe_store!(ptr, pseudojet) + return 0 +end + +struct C_ClusterSequence{T} + algorithm::JetAlgorithm.Algorithm + power::Cdouble + R::Cdouble + strategy::RecoStrategy.Strategy + jets::Ptr{T} + jets_length::Csize_t + n_initial_jets::Clong + history::Ptr{HistoryElement} + history_length::Csize_t + Qtot::Cdouble +end + +function free_members(ptr::Ptr{C_ClusterSequence{T}}) where {T} + if ptr != C_NULL + clusterseq = unsafe_load(ptr) + Libc.free(clusterseq.jets) + Libc.free(clusterseq.history) + # Struct is immutable so pointers can't be assigned NULL and lengths updated to zero (without making a copy) + end +end +Base.@ccallable function jetreconstruction_ClusterSequence_free_members_(ptr::Ptr{C_ClusterSequence{PseudoJet}})::Cvoid + free_members(ptr) + return nothing +end + +function ClusterSequence{T}(c::C_ClusterSequence{T}) where {T} + jets_v = unsafe_wrap_c_array(c.jets, c.jets_length) + history_v = unsafe_wrap_c_array(c.history, c.history_length) + return ClusterSequence{T}(c.algorithm, c.power, c.R, c.strategy, jets_v, + c.n_initial_jets, + history_v, c.Qtot) +end + +function C_ClusterSequence{T}(clustersequence::ClusterSequence{T}) where {T} + jets_ptr, jets_length = make_c_array(clustersequence.jets) + history_ptr, history_length = make_c_array(clustersequence.history) + return C_ClusterSequence{T}(clustersequence.algorithm, clustersequence.power, + clustersequence.R, clustersequence.strategy, jets_ptr, + jets_length, clustersequence.n_initial_jets, history_ptr, + history_length, clustersequence.Qtot) +end + +function c_jet_reconstruct(particles::Ptr{T}, + particles_length::Csize_t, + algorithm::JetAlgorithm.Algorithm, + R::Cdouble, + strategy::RecoStrategy.Strategy, + result::Ptr{C_ClusterSequence{U}}) where {T, U} + particles_v = unsafe_wrap_c_array(particles, particles_length) + clusterseq = jet_reconstruct(particles_v; p = nothing, algorithm = algorithm, R = R, + strategy = strategy) + c_clusterseq = C_ClusterSequence{U}(clusterseq) + unsafe_store!(result, c_clusterseq) + return 0 +end + +Base.@ccallable function jetreconstruction_jet_reconstruct(particles::Ptr{PseudoJet}, + particles_length::Csize_t, + algorithm::JetAlgorithm.Algorithm, + R::Cdouble, + strategy::RecoStrategy.Strategy, + result::Ptr{C_ClusterSequence{PseudoJet}})::Cint + c_jet_reconstruct(particles, particles_length, algorithm, R, strategy, result) + return 0 +end + +struct C_JetsResult{T} + data::Ptr{T} + length::Csize_t +end + +function free_members(ptr::Ptr{C_JetsResult{T}}) where {T} + if ptr != C_NULL + result = unsafe_load(ptr) + Libc.free(result.data) + # Struct is immutable so pointer can't be assigned NULL and length updated to zero (without making a copy) + end +end + +Base.@ccallable function jetreconstruction_JetsResult_free_members_(ptr::Ptr{C_JetsResult{PseudoJet}})::Cvoid + free_members(ptr) + return nothing +end + +function jets_selection(selector, clustersequence::Ptr{C_ClusterSequence{T}}, + result::Ptr{C_JetsResult{U}}; kwargs...)::Cint where {T, U} + c_clusterseq = unsafe_load(clustersequence) + clusterseq = ClusterSequence{T}(c_clusterseq) + jets_result = selector(clusterseq; kwargs...) + println(jets_result) + c_results = C_JetsResult{U}(C_NULL, 0) # TODO convert and write to result + unsafe_store!(result, c_results) + return 0 +end + +Base.@ccallable function jetreconstruction_exclusive_jets_dcut(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + dcut::Cdouble, + result::Ptr{C_JetsResult{PseudoJet}})::Cint + return jets_selection(exclusive_jets, clustersequence, result; dcut = dcut) +end + +Base.@ccallable function jetreconstruction_exclusive_jets_njets(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + njets::Csize_t, + result::Ptr{C_JetsResult{PseudoJet}})::Cint + return jets_selection(exclusive_jets, clustersequence, result; njets = njets) +end + +Base.@ccallable function jetreconstruction_inclusive_jets(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + ptmin::Cdouble, + result::Ptr{C_JetsResult{PseudoJet}})::Cint + return jets_selection(inclusive_jets, clustersequence, result; ptmin = ptmin) +end + +end # module C_JetReconstruction diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 9929187..f7ff618 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -89,4 +89,7 @@ export jetsplot, animatereco include("JSONresults.jl") export FinalJet, FinalJets +# C bindings sub-module +include("C_JetReconstruction/C_JetReconstruction.jl") + end