diff --git a/Project.toml b/Project.toml index a8edb55..83e195d 100644 --- a/Project.toml +++ b/Project.toml @@ -4,13 +4,16 @@ authors = ["Aman Pandey"] version = "0.1.0" [deps] +CFITSIO = "3b1b4be9-1499-4b22-8d78-7db3344d1961" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568" @@ -19,13 +22,16 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] +CFITSIO = "1.7.1" DataFrames = "1.3" Distributions = "0.25" +DocStringExtensions = "0.9.5" FFTW = "1.4" FITSIO = "0.16" HDF5 = "0.16" Intervals = "1.8" JuliaFormatter = "1.0.62" +LinearAlgebra = "1.11.0" Logging = "1.11.0" NaNMath = "0.3, 1" ProgressBars = "1.4" diff --git a/src/Stingray.jl b/src/Stingray.jl index 6d8e74f..b1ef437 100644 --- a/src/Stingray.jl +++ b/src/Stingray.jl @@ -3,6 +3,8 @@ module Stingray using ResumableFunctions, StatsBase, Statistics, DataFrames using FFTW, NaNMath, FITSIO, Intervals using ProgressBars: tqdm as show_progress +using DocStringExtensions +using LinearAlgebra include("fourier.jl") export positive_fft_bins @@ -33,6 +35,27 @@ export bin_intervals_from_gtis include("utils.jl") include("events.jl") -export readevents, EventList, DictMetadata +export FITSMetadata, + EventList, + times, + energies, + has_energies, + filter_time!, + filter_energy!, + filter_time, + filter_energy, + colnames, + read_energy_column, + readevents, + summary, + filter_on! +include("missionSupport.jl") +export MissionSupport, + get_mission_support, + apply_calibration, + patch_mission_info, + SIMPLE_CALIBRATION_FUNCS, + interpret_fits_data!, + AbstractMissionSupport end diff --git a/src/events.jl b/src/events.jl index 70c54c7..2630b7a 100644 --- a/src/events.jl +++ b/src/events.jl @@ -1,85 +1,745 @@ """ - DictMetadata +$(TYPEDEF) -A structure containing metadata from FITS file headers. +Metadata associated with a FITS or events file. -## Fields +$(TYPEDFIELDS) -- `headers::Vector{Dict{String,Any}}`: A vector of dictionaries containing header information from each HDU. +# Examples +```julia +# Metadata is typically created automatically when reading events +ev = readevents("data.fits") +println(ev.meta.filepath) # Shows the file path +println(ev.meta.energy_units) # Shows "PI", "ENERGY", or "PHA" +``` """ -struct DictMetadata - headers::Vector{Dict{String,Any}} +struct FITSMetadata{H} + "Path to the FITS file" + filepath::String + "HDU index that the metadata was read from" + hdu::Int + "Units of energy (column name: ENERGY, PI, or PHA)" + energy_units::Union{Nothing,String} + "Extra columns that were requested during read" + extra_columns::Dict{String,Vector} + "FITS headers from the selected HDU" + headers::H end +function Base.show(io::IO, ::MIME"text/plain", m::FITSMetadata) + println( + io, + "FITSMetadata for $(basename(m.filepath))[$(m.hdu)] with $(length(m.extra_columns)) extra column(s)", + ) +end + +""" +$(TYPEDEF) + +Container for an events list storing times, energies, and associated metadata. + +$(TYPEDFIELDS) + +# Constructors +```julia +# Read from FITS file (recommended) +ev = readevents("events.fits") + +# Create directly for testing (simplified constructor) +ev = EventList([1.0, 2.0, 3.0], [0.5, 1.2, 2.1]) # times and energies +ev = EventList([1.0, 2.0, 3.0]) # times only +``` + +# Interface +- `length(ev)`: Number of events +- `times(ev)`: Access times vector +- `energies(ev)`: Access energies vector (may be `nothing`) +- `has_energies(ev)`: Check if energies are present + +Generally should not be directly constructed, but read from file using [`readevents`](@ref). + +See also: [`filter_time!`](@ref), [`filter_energy!`](@ref) for filtering operations. +""" +struct EventList{TimeType<:AbstractVector,MetaType<:FITSMetadata} + "Vector with recorded times" + times::TimeType + "Vector with recorded energies (else `nothing`)" + energies::Union{Nothing,TimeType} + "Metadata from FITS file" + meta::MetaType +end + +""" + EventList(times::Vector{T}, energies::Union{Nothing,Vector{T}}=nothing) where T + +Simple constructor for testing without FITS files. Creates an EventList with dummy metadata. + +# Arguments +- `times::Vector{T}`: Vector of event times +- `energies::Union{Nothing,Vector{T}}`: Optional vector of event energies + +# Examples +```julia +# Times only +ev = EventList([1.0, 2.0, 3.0]) + +# Times and energies +ev = EventList([1.0, 2.0, 3.0], [0.5, 1.2, 2.1]) +``` +""" +function EventList(times::Vector{T}, energies::Union{Nothing,Vector{T}} = nothing) where {T} + dummy_meta = FITSMetadata( + "[no file]", # filepath + 1, # hdu + nothing, # energy_units + Dict{String,Vector}(), # extra_columns + Dict{String,Any}(), # headers + ) + EventList(times, energies, dummy_meta) +end + +function Base.show(io::IO, ::MIME"text/plain", ev::EventList) + print(io, "EventList with $(length(ev.times)) times") + if !isnothing(ev.energies) + print(io, " and energies") + end + println(io) +end + +# ============================================================================ +# Interface Methods +# ============================================================================ + +""" + length(ev::EventList) + +Return the number of events in the EventList. +""" +Base.length(ev::EventList) = length(ev.times) + +""" + size(ev::EventList) + +Return the size of the EventList as a tuple. +""" +Base.size(ev::EventList) = (length(ev),) + +""" + times(ev::EventList) + +Access the times vector of an EventList. + +# Examples +```julia +ev = readevents("data.fits") +time_data = times(ev) # Get times vector +``` +""" +times(ev::EventList) = ev.times + +""" + energies(ev::EventList) + +Access the energies vector of an EventList. Returns `nothing` if no energies are present. + +# Examples +```julia +ev = readevents("data.fits") +energy_data = energies(ev) # May be nothing +if !isnothing(energy_data) + println("Energy range: \$(extrema(energy_data))") +end +``` +""" +energies(ev::EventList) = ev.energies + +""" + has_energies(ev::EventList) + +Check whether the EventList contains energy information. + +# Examples +```julia +ev = readevents("data.fits") +if has_energies(ev) + println("Energy data available") +end +``` +""" +has_energies(ev::EventList) = !isnothing(ev.energies) + +# ============================================================================ +# Filtering Functions (Composable and In-Place) +# ============================================================================ + """ - EventList{T} +$(TYPEDSIGNATURES) + +Filter all columns of the EventList based on a predicate `f` applied to the times. +Modifies the EventList in-place for efficiency. + +Returns the modified EventList (for chaining operations). -A structure containing event data from a FITS file. +# Filtering Operations +EventList supports composable filtering operations: +```julia +# Filter by time (in-place) +filter_time!(t -> t > 100.0, ev) -## Fields +# Filter times greater than some minimum using function composition +min_time = 100.0 +filter_time!(x -> x > min_time, ev) -- `filename::String`: Path to the source FITS file. -- `times::Vector{T}`: Vector of event times. -- `energies::Vector{T}`: Vector of event energies. -- `metadata::DictMetadata`: Metadata information extracted from the FITS file headers. +# Chaining filters +filter_energy!(x -> x < 10.0, filter_time!(t -> t > 100.0, ev)) + +# Non-mutating version +ev_filtered = filter_time(t -> t > 100.0, ev) +``` + +See also [`filter_energy!`](@ref), [`filter_time`](@ref). """ -struct EventList{T} - filename::String - times::Vector{T} - energies::Vector{T} - metadata::DictMetadata +filter_time!(f, ev::EventList) = filter_on!(f, ev.times, ev) + +""" +$(TYPEDSIGNATURES) + +Filter all columns of the EventList based on a predicate `f` applied to the energies. +Modifies the EventList in-place for efficiency. + +Returns the modified EventList (for chaining operations). + +# Examples +# Filtering Operations +EventList supports composable filtering operations: +```julia +# Filter energies less than 10 keV +filter_energy!(energy_val -> energy_val < 10.0, ev) + +# With function composition +max_energy = 10.0 +filter_energy!(x -> x < max_energy, ev) + +# Chaining with time filter +filter_energy!(x -> x < 10.0, filter_time!(t -> t > 100.0, ev)) + +# Non-mutating version +ev_filtered = filter_energy(energy_val -> energy_val < 10.0, ev) +``` + +# Throws +- `AssertionError`: If the EventList has no energy data + +See also [`filter_time!`](@ref), [`filter_energy`](@ref). +""" +function filter_energy!(f, ev::EventList) + @assert has_energies(ev) "No energies present in the EventList." + filter_on!(f, ev.energies, ev) end """ - readevents(path; T = Float64) + filter_on!(f, src_col::AbstractVector, ev::EventList) -Read event data from a FITS file into an EventList structure. The `path` is a -string that points to the location of the FITS file. `T` is used to specify -which numeric type to convert the data to. +Internal function to filter EventList based on predicate applied to source column. +Uses efficient in-place filtering adapted from Base.filter! implementation. -Returns an [`EventList`](@ref) containing the extracted data. +This function maintains consistency across all columns (times, energies, extra_columns) +by applying the same filtering mask derived from the source column. -## Notes +# Arguments +- `f`: Predicate function applied to elements of `src_col` +- `src_col::AbstractVector`: Source column to generate filtering mask from +- `ev::EventList`: EventList to filter -The function extracts `TIME` and `ENERGY` columns from any TableHDU in the FITS -file. All headers from each HDU are collected into the metadata field. +# Implementation Notes +- Uses `eachindex` for portable iteration over array indices +- Modifies arrays in-place using a two-pointer technique +- Resizes all arrays to final filtered length +- Maintains type stability with `::Bool` annotation on predicate result """ -function readevents(path; T = Float64) - headers = Dict{String,Any}[] - times = T[] - energies = T[] +function filter_on!(f, src_col::AbstractVector, ev::EventList) + @assert size(src_col) == size(ev.times) "Source column size must match times size" + + # Modified from Base.filter! implementation for multiple arrays + # Use two pointers: i for reading, j for writing + j = firstindex(ev.times) + + for i in eachindex(ev.times) + predicate = f(src_col[i])::Bool - FITS(path, "r") do f - for i = 1:length(f) # Iterate over HDUs - hdu = f[i] - header_dict = Dict{String,Any}() - for key in keys(read_header(hdu)) - header_dict[string(key)] = read_header(hdu)[key] + if predicate + # Copy elements to new position + ev.times[j] = ev.times[i] + + if !isnothing(ev.energies) + ev.energies[j] = ev.energies[i] end - push!(headers, header_dict) - - # Check if the HDU is a table - if isa(hdu, TableHDU) - # Get column names using the correct FITSIO method - colnames = FITSIO.colnames(hdu) - - if "TIME" in colnames - times = convert(Vector{T}, read(hdu, "TIME")) - end - if "ENERGY" in colnames - energies = convert(Vector{T}, read(hdu, "ENERGY")) - end + + # Handle extra columns + for (_, col) in ev.meta.extra_columns + col[j] = col[i] end + + j = nextind(ev.times, j) end end - if isempty(times) - @warn "No TIME data found in FITS file $(path). Time series analysis will not be possible." + + # Resize all arrays to new length + if j <= lastindex(ev.times) + new_length = j - 1 + resize!(ev.times, new_length) + + if !isnothing(ev.energies) + resize!(ev.energies, new_length) + end + + for (_, col) in ev.meta.extra_columns + resize!(col, new_length) + end + end + + ev +end + +# ============================================================================ +# Non-mutating Filter Functions +# ============================================================================ + +""" + filter_time(f, ev::EventList) + +Return a new EventList with events filtered by predicate `f` applied to times. +This is the non-mutating version of [`filter_time!`](@ref). + +# Arguments +- `f`: Predicate function that takes a time value and returns a Boolean +- `ev::EventList`: EventList to filter (not modified) + +# Returns +New EventList with filtered events + +# Examples +```julia +# Create filtered copy +ev_filtered = filter_time(t -> t > 100.0, ev) + +# Original EventList is unchanged +println(length(ev)) # Original length +println(length(ev_filtered)) # Filtered length +``` + +See also [`filter_time!`](@ref), [`filter_energy`](@ref). +""" +function filter_time(f, ev::EventList) + new_ev = deepcopy(ev) + filter_time!(f, new_ev) +end + +""" + filter_energy(f, ev::EventList) + +Return a new EventList with events filtered by predicate `f` applied to energies. +This is the non-mutating version of [`filter_energy!`](@ref). + +# Arguments +- `f`: Predicate function that takes an energy value and returns a Boolean +- `ev::EventList`: EventList to filter (not modified) + +# Returns +New EventList with filtered events + +# Examples +```julia +# Create filtered copy +ev_filtered = filter_energy(energy_val -> energy_val < 10.0, ev) + +# Original EventList is unchanged +println(length(ev)) # Original length +println(length(ev_filtered)) # Filtered length +``` + +# Throws +- `AssertionError`: If the EventList has no energy data + +See also [`filter_energy!`](@ref), [`filter_time`](@ref). +""" +function filter_energy(f, ev::EventList) + new_ev = deepcopy(ev) + filter_energy!(f, new_ev) +end + +# ============================================================================ +# File Reading Functions +# ============================================================================ + +""" + colnames(file::AbstractString; hdu = 2) + +Return a vector of all column names of a FITS file, reading from the specified HDU. + +# Arguments +- `file::AbstractString`: Path to FITS file +- `hdu::Int`: HDU index to read from (default: 2, typical for event data) + +# Returns +Vector of column names as strings + +# Examples +```julia +cols = colnames("events.fits") +println(cols) # ["TIME", "PI", "RAWX", "RAWY", ...] + +# Check if energy column exists +if "ENERGY" in colnames("events.fits") + println("Energy data available") +end +``` +""" +function colnames(file::AbstractString; hdu = 2) + FITS(file) do f + selected_hdu = f[hdu] + FITSIO.colnames(selected_hdu) + end +end + +""" + read_energy_column(hdu; energy_alternatives = ["ENERGY", "PI", "PHA"], T = Float64) + +Attempt to read the energy column of an HDU from a list of alternative names. + +This function provides a robust way to read energy data from FITS files, as different +missions and instruments use different column names for energy information. + +# Arguments +- `hdu`: FITS HDU object to read from +- `energy_alternatives::Vector{String}`: List of column names to try (default: ["ENERGY", "PI", "PHA"]) +- `T::Type`: Type to convert energy data to (default: Float64) + +# Returns +`(column_name, data)` tuple where: +- `column_name::Union{Nothing,String}`: Name of the column that was successfully read, or `nothing` +- `data::Union{Nothing,Vector{T}}`: Energy data as Vector{T}, or `nothing` if no column found + +# Examples +```julia +FITS("events.fits") do f + hdu = f[2] + col_name, energy_data = read_energy_column(hdu) + if !isnothing(energy_data) + println("Found energy data in column: \$col_name") + println("Energy range: \$(extrema(energy_data))") + end +end +``` + +# Implementation Notes +- Tries columns in order until one is successfully read +- Uses case-insensitive matching for column names +- Handles read errors gracefully by trying the next column +- Separated from main reading function for testability and clarity +- Type-stable with explicit return type annotation +- Added case_sensitive=false parameter: This tells FITSIO.jl to use the old behavior for backward compatibility +""" +function read_energy_column( + hdu; + energy_alternatives::Vector{String} = ["ENERGY", "PI", "PHA"], + T::Type = Float64, +)::Tuple{Union{Nothing,String},Union{Nothing,Vector{T}}} + + # Get actual column names from the file + all_cols = FITSIO.colnames(hdu) + + for col_name in energy_alternatives + # Find matching column name (case-insensitive) + actual_col = findfirst(col -> uppercase(col) == uppercase(col_name), all_cols) + + if !isnothing(actual_col) + actual_col_name = all_cols[actual_col] + try + # Use the actual column name from the file + data = read(hdu, actual_col_name, case_sensitive = false) + return actual_col_name, convert(Vector{T}, data) + catch + # If this column exists but can't be read, try the next one + continue + end + end + end + + return nothing, nothing +end +""" + readevents(path; kwargs...) + +Read an [`EventList`](@ref) from a FITS file with optional mission-specific support. +Will attempt to read an energy column if one exists, with mission-specific calibration +and interpretation capabilities. + +This is the primary function for loading X-ray event data from FITS files. +It handles the complexities of different file formats, provides mission-specific +energy calibration, and offers a consistent interface for accessing event data. + +# Arguments +- `path::AbstractString`: Path to the FITS file + +# Keyword Arguments +- `mission::Union{String,Nothing} = nothing`: Mission name for mission-specific support (e.g., "nustar", "chandra", "xmm") +- `instrument::Union{String,Nothing} = nothing`: Instrument identifier for mission-specific calibration +- `epoch::Union{Float64,Nothing} = nothing`: Observation epoch in MJD for time-dependent calibrations +- `hdu::Int = 2`: HDU index to read from (typically 2 for event data) +- `T::Type = Float64`: Type to cast the time and energy columns to +- `sort::Bool = false`: Whether to sort by time if not already sorted +- `extra_columns::Vector{String} = []`: Extra columns to read from the same HDU +- `energy_alternatives::Vector{String} = ["ENERGY", "PI", "PHA"]`: Energy column alternatives to try (overridden by mission-specific preferences) + +# Returns +`EventList{Vector{T}, FITSMetadata{FITSIO.FITSHeader}}`: EventList containing the event data + +# Mission Support +When a mission is specified, the function will: +- Use mission-specific energy column preferences (overrides `energy_alternatives`) +- Apply mission-specific calibration functions to convert PI channels to energy +- Apply mission-specific header patches and interpretations +- Handle mission-specific GTI (Good Time Interval) extensions + +Supported missions: +- `"nustar"`: Nuclear Spectroscopic Telescope Array +- `"xmm"`: X-ray Multi-Mirror Mission +- `"nicer"`: Neutron star Interior Composition Explorer +- `"ixpe"`: Imaging X-ray Polarimetry Explorer +- `"chandra"` / `"axaf"`: Chandra X-ray Observatory +- `"xte"` / `"rxte"`: Rossi X-ray Timing Explorer + +# Examples +```julia +# Basic usage +ev = readevents("events.fits") + +# With mission-specific support +ev = readevents("events.fits", mission="nustar", instrument="FPM_A") + +# Mission support with epoch for time-dependent calibration +ev = readevents("events.fits", mission="xte", instrument="PCA", epoch=50000.0) + +# With custom options +ev = readevents("events.fits", hdu=3, sort=true, T=Float32) + +# Reading extra columns +ev = readevents("events.fits", extra_columns=["RAWX", "RAWY", "DETX", "DETY"]) + +# Accessing the data +println("Number of events: \$(length(ev))") +println("Time range: \$(extrema(times(ev)))") +if has_energies(ev) + println("Energy range: \$(extrema(energies(ev)))") + println("Energy column: \$(ev.meta.energy_units)") +end +``` + +# Mission-Specific Examples +```julia +# NuSTAR data with automatic PI to energy conversion +ev = readevents("nustar_events.fits", mission="nustar") +# Uses mission-specific energy alternatives: ["PI", "ENERGY", "PHA"] +# Applies calibration: E(keV) = PI * 0.04 + 1.62 + +# Chandra data with mission-specific handling +ev = readevents("chandra_events.fits", mission="chandra") +# Uses energy alternatives: ["ENERGY", "PI", "PHA"] +# Applies Chandra-specific header interpretations + +# XMM data with mission patches +ev = readevents("xmm_events.fits", mission="xmm") +# Handles XMM-specific GTI extensions: ["GTI", "GTI0", "STDGTI"] +``` + +# Error Handling +- Throws `AssertionError` if time and energy vectors have different sizes +- Throws `AssertionError` if times are not sorted and `sort=false` +- Throws `ArgumentError` if mission name is empty string +- FITS reading errors are propagated from the FITSIO.jl library +- Warns if mission is not recognized (uses identity calibration function) + +# Implementation Notes +- Uses type-stable FITS reading with explicit type conversions +- Handles missing energy data gracefully +- Supports efficient multi-column sorting when `sort=true` +- Creates metadata with all relevant file information +- Validates data consistency before returning +- Mission-specific energy alternatives override the default parameter +- Applies mission-specific calibration to PI channel data automatically +- Uses case-insensitive column matching for robustness +""" +function readevents( + path::AbstractString; + mission::Union{String,Nothing} = nothing, + instrument::Union{String,Nothing} = nothing, + epoch::Union{Float64,Nothing} = nothing, + hdu::Int = 2, + T::Type = Float64, + sort::Bool = false, + extra_columns::Vector{String} = String[], + energy_alternatives::Vector{String} = ["ENERGY", "PI", "PHA"], + kwargs..., +)::EventList{Vector{T},FITSMetadata{FITSIO.FITSHeader}} + + # Get mission support if specified + mission_support = if !isnothing(mission) + ms = get_mission_support(mission, instrument, epoch) + # Use mission-specific energy alternatives if available + energy_alternatives = ms.energy_alternatives + ms + else + nothing + end + + # Read data from FITS file with type-stable operations + time::Vector{T}, + energy::Union{Nothing,Vector{T}}, + energy_col::Union{Nothing,String}, + header::FITSIO.FITSHeader, + extra_data::Dict{String,Vector} = FITS(path, "r") do f + + selected_hdu = f[hdu] + + # Apply mission-specific FITS interpretation if available + if !isnothing(mission_support) && !isnothing(mission_support.interpretation_func) + interpret_fits_data!(f, mission_support) + end + + # Read header (type-stable) + header = read_header(selected_hdu) + + # Get actual column names to find the correct TIME column + all_cols = FITSIO.colnames(selected_hdu) + time = convert(Vector{T}, read(selected_hdu, "TIME", case_sensitive = false)) + + # Read energy column using separated function with mission-specific alternatives + energy_column, energy = read_energy_column( + selected_hdu; + T = T, + energy_alternatives = energy_alternatives, + ) + + # Apply mission-specific calibration if we have PI data and mission support + if !isnothing(energy) && !isnothing(mission_support) && + !isnothing(energy_column) && uppercase(energy_column) == "PI" + energy = convert(Vector{T}, apply_calibration(mission_support, energy)) + # Update the energy column name to reflect that it's now calibrated + energy_column = "ENERGY" + end + + # Read extra columns with case-insensitive option + extra_data = Dict{String,Vector}() + for col_name in extra_columns + # Find actual column name (case-insensitive) + actual_col_idx = + findfirst(col -> uppercase(col) == uppercase(col_name), all_cols) + if !isnothing(actual_col_idx) + actual_col_name = all_cols[actual_col_idx] + extra_data[col_name] = + read(selected_hdu, actual_col_name, case_sensitive = false) + else + @warn "Column '$col_name' not found in FITS file" + end + end + + (time, energy, energy_column, header, extra_data) + end + + # Apply mission-specific header patches if available + if !isnothing(mission_support) + # Convert header to dictionary for patching + header_dict = Dict{String,Any}() + + # Use the proper way to access FITSHeader keys and values + for key in keys(header) + header_dict[key] = header[key] + end + + # Apply mission patches + patched_header_dict = patch_mission_info(header_dict, mission) + + # Note: We keep the original header structure but could extend this + # to update the header with patched information if needed + end + + # Validate energy-time consistency + if !isnothing(energy) + @assert size(time) == size(energy) "Time and energy do not match sizes ($(size(time)) != $(size(energy)))" + end + + # Handle sorting if requested + if !issorted(time) + if sort + # Efficient sorting of multiple arrays + sort_indices = sortperm(time) + time = time[sort_indices] + + if !isnothing(energy) + energy = energy[sort_indices] + end + + # Sort extra columns + for (col_name, col_data) in extra_data + extra_data[col_name] = col_data[sort_indices] + end + else + @assert false "Times are not sorted (pass `sort = true` to force sorting)" + end + end + + # Create metadata - record the column name that was found (possibly updated by calibration) + meta = FITSMetadata(path, hdu, energy_col, extra_data, header) + + # Return type-stable EventList + EventList(time, energy, meta) +end + +# ============================================================================ +# Utility Functions +# ============================================================================ + +""" + summary(ev::EventList) + +Provide a comprehensive summary of the EventList contents. + +# Arguments +- `ev::EventList`: EventList to summarize + +# Returns +String with summary information including: +- Number of events +- Time span +- Energy range (if available) +- Energy units (if available) +- Number of extra columns + +# Examples +```julia +ev = readevents("events.fits") +println(summary(ev)) +# Output: "EventList: 1000 events over 3600.0 time units, energies: 0.5 - 12.0 (PI), 2 extra columns" +``` +""" +function Base.summary(ev::EventList) + n_events = length(ev) + time_span = isempty(ev.times) ? 0.0 : maximum(ev.times) - minimum(ev.times) + + summary_str = "EventList: $n_events events over $(time_span) time units" + + if has_energies(ev) + energy_range = extrema(ev.energies) + summary_str *= ", energies: $(energy_range[1]) - $(energy_range[2])" + if !isnothing(ev.meta.energy_units) + summary_str *= " ($(ev.meta.energy_units))" + end end - if isempty(energies) - @warn "No ENERGY data found in FITS file $(path). Energy spectrum analysis will not be possible." + if !isempty(ev.meta.extra_columns) + summary_str *= ", $(length(ev.meta.extra_columns)) extra columns" end - metadata = DictMetadata(headers) - return EventList{T}(path, times, energies, metadata) + return summary_str end diff --git a/src/missionSupport.jl b/src/missionSupport.jl new file mode 100644 index 0000000..9ac0089 --- /dev/null +++ b/src/missionSupport.jl @@ -0,0 +1,226 @@ +""" +Dictionary of simple conversion functions for different missions. + +This dictionary provides PI (Pulse Invariant) to energy conversion functions +for various X-ray astronomy missions. Each function takes a PI channel value +and returns the corresponding energy in keV. + +Supported missions: +- NuSTAR: Nuclear Spectroscopic Telescope Array +- XMM: X-ray Multi-Mirror Mission +- NICER: Neutron star Interior Composition Explorer +- IXPE: Imaging X-ray Polarimetry Explorer +- AXAF/Chandra: Advanced X-ray Astrophysics Facility +- XTE/RXTE: Rossi X-ray Timing Explorer +""" +const SIMPLE_CALIBRATION_FUNCS = Dict{String, Function}( + "nustar" => (pi) -> pi * 0.04 + 1.62, + "xmm" => (pi) -> pi * 0.001, + "nicer" => (pi) -> pi * 0.01, + "ixpe" => (pi) -> pi / 375 * 15, + "axaf" => (pi) -> (pi - 1) * 14.6e-3, # Chandra/AXAF + "chandra" => (pi) -> (pi - 1) * 14.6e-3, # Explicit chandra entry + "xte" => (pi) -> pi * 0.025 # RXTE/XTE +) + +""" +Abstract type for mission-specific calibration and interpretation. + +This serves as the base type for all mission support implementations, +allowing for extensibility and type safety in mission-specific operations. +""" +abstract type AbstractMissionSupport end + +""" + MissionSupport{T} <: AbstractMissionSupport + +Structure containing mission-specific calibration and interpretation information. + +This structure encapsulates all the necessary information for handling +data from a specific X-ray astronomy mission, including calibration +functions, energy column alternatives, and GTI extension preferences. + +# Fields +- `name::String`: Mission name (normalized to lowercase) +- `instrument::Union{String, Nothing}`: Instrument identifier +- `epoch::Union{T, Nothing}`: Observation epoch in MJD (for time-dependent calibrations) +- `calibration_func::Function`: PI to energy conversion function +- `interpretation_func::Union{Function, Nothing}`: Mission-specific FITS interpretation function +- `energy_alternatives::Vector{String}`: Preferred energy column names in order of preference +- `gti_extensions::Vector{String}`: GTI extension names in order of preference + +# Type Parameters +- `T`: Type of the epoch parameter (typically Float64) +""" +struct MissionSupport{T} <: AbstractMissionSupport + name::String + instrument::Union{String, Nothing} + epoch::Union{T, Nothing} + calibration_func::Function + interpretation_func::Union{Function, Nothing} + energy_alternatives::Vector{String} + gti_extensions::Vector{String} +end + +""" + get_mission_support(mission::String, instrument=nothing, epoch=nothing) -> MissionSupport + +Create mission support object with mission-specific parameters. + +This function creates a MissionSupport object containing all the necessary +information for processing data from a specified X-ray astronomy mission. +It handles mission aliases (e.g., Chandra/AXAF) and provides appropriate +defaults for each mission. + +# Arguments +- `mission::String`: Mission name (case-insensitive) +- `instrument::Union{String, Nothing}=nothing`: Instrument identifier +- `epoch::Union{Float64, Nothing}=nothing`: Observation epoch in MJD + +# Returns +- `MissionSupport{Float64}`: Mission support object + +# Throws +- `ArgumentError`: If mission name is empty + +# Examples +```julia +# Basic usage +ms = get_mission_support("nustar") + +# With instrument specification +ms = get_mission_support("nustar", "FPM_A") + +# With epoch for time-dependent calibrations +ms = get_mission_support("xte", "PCA", 50000.0) +``` +""" +function get_mission_support(mission::String, + instrument::Union{String, Nothing}=nothing, + epoch::Union{Float64, Nothing}=nothing) + + # Check for empty mission string + if isempty(mission) + throw(ArgumentError("Mission name cannot be empty")) + end + + mission_lower = lowercase(mission) + + if mission_lower in ["chandra", "axaf"] + mission_lower = "chandra" + end + + calib_func = if haskey(SIMPLE_CALIBRATION_FUNCS, mission_lower) + SIMPLE_CALIBRATION_FUNCS[mission_lower] + else + @warn "Mission $mission not recognized, using identity function" + identity + end + + energy_alts = if mission_lower in ["chandra", "axaf"] + ["ENERGY", "PI", "PHA"] + elseif mission_lower == "xte" + ["PHA", "PI", "ENERGY"] + elseif mission_lower == "nustar" + ["PI", "ENERGY", "PHA"] + else + ["ENERGY", "PI", "PHA"] + end + + # Mission-specific GTI extensions + gti_exts = if mission_lower == "xmm" + ["GTI", "GTI0", "STDGTI"] + elseif mission_lower in ["chandra", "axaf"] + ["GTI", "GTI0", "GTI1", "GTI2", "GTI3"] + else + ["GTI", "STDGTI"] + end + + MissionSupport{Float64}(mission_lower, instrument, epoch, calib_func, nothing, energy_alts, gti_exts) +end + +""" + apply_calibration(mission_support::MissionSupport, pi_channels::AbstractArray) -> Vector{Float64} + +Apply calibration function to PI channels. + +Converts PI (Pulse Invariant) channel values to energies in keV using +the mission-specific calibration function stored in the MissionSupport object. + +# Arguments +- `mission_support::MissionSupport`: Mission support object containing calibration function +- `pi_channels::AbstractArray{T}`: Array of PI channel values + +# Returns +- `Vector{Float64}`: Array of energy values in keV + +# Examples +```julia +ms = get_mission_support("nustar") +pi_values = [100, 500, 1000] +energies = apply_calibration(ms, pi_values) +``` +""" +function apply_calibration(mission_support::MissionSupport, pi_channels::AbstractArray{T}) where T + if isempty(pi_channels) + return similar(pi_channels, Float64) + end + return mission_support.calibration_func.(pi_channels) +end + +""" + patch_mission_info(info::Dict{String,Any}, mission=nothing) -> Dict{String,Any} + +Apply mission-specific patches to header information. + +This function applies mission-specific modifications to FITS header information +to handle mission-specific quirks and conventions. It's based on the Python +implementation in Stingray's mission interpretation module. + +# Arguments +- `info::Dict{String,Any}`: Dictionary containing header information +- `mission::Union{String,Nothing}=nothing`: Mission name + +# Returns +- `Dict{String,Any}`: Patched header information dictionary + +# Examples +```julia +info = Dict("gti" => "STDGTI", "ecol" => "PHA") +patched = patch_mission_info(info, "xmm") # Adds GTI0 to gti field +``` +""" +function patch_mission_info(info::Dict{String,Any}, mission::Union{String,Nothing}=nothing) + if isnothing(mission) + return info + end + + mission_lower = lowercase(mission) + patched_info = copy(info) + + # Normalize chandra/axaf + if mission_lower in ["chandra", "axaf"] + mission_lower = "chandra" + end + + if mission_lower == "xmm" && haskey(patched_info, "gti") + patched_info["gti"] = string(patched_info["gti"], ",GTI0") + elseif mission_lower == "xte" && haskey(patched_info, "ecol") + patched_info["ecol"] = "PHA" + patched_info["ccol"] = "PCUID" + elseif mission_lower == "chandra" + # Chandra-specific patches + if haskey(patched_info, "DETNAM") + patched_info["detector"] = patched_info["DETNAM"] + end + if haskey(patched_info, "TIMESYS") + patched_info["time_system"] = patched_info["TIMESYS"] + end + end + + return patched_info +end +function interpret_fits_data!(f::FITS, mission_support::MissionSupport) + # Placeholder for mission-specific interpretation + return nothing +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 4275ebf..16dd9d2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,13 @@ using Stingray using Test using FFTW, Distributions, Statistics, StatsBase, HDF5, FITSIO -using Logging +using Logging ,LinearAlgebra +using CFITSIO include("test_fourier.jl") include("test_gti.jl") -include("test_events.jl") +@testset verbose=true "Eventlist" begin + include("test_events.jl") +end +@testset verbose=true "Synthetic Events Tests" begin + include("test_missionSupport.jl") +end \ No newline at end of file diff --git a/test/test_events.jl b/test/test_events.jl index ff1c91b..b014915 100644 --- a/test/test_events.jl +++ b/test/test_events.jl @@ -1,149 +1,439 @@ -@testset "EventList Tests" begin - # Test 1: Create a sample FITS file for testing - @testset "Sample FITS file creation" begin - test_dir = mktempdir() - sample_file = joinpath(test_dir, "sample.fits") - f = FITS(sample_file, "w") - write(f, Int[]) # Empty primary array - # Create a binary table HDU with TIME and ENERGY columns - times = Float64[1.0, 2.0, 3.0, 4.0, 5.0] - energies = Float64[10.0, 20.0, 15.0, 25.0, 30.0] - # Add a binary table extension - table = Dict{String,Array}() - table["TIME"] = times - table["ENERGY"] = energies - write(f, table) - close(f) - @test isfile(sample_file) - # Test reading the sample file - data = readevents(sample_file) - @test data.filename == sample_file - @test length(data.times) == 5 - @test length(data.energies) == 5 - @test eltype(data.times) == Float64 - @test eltype(data.energies) == Float64 - end +# Test data path (if using test data directory) +# TEST_DATA_PATH = joinpath(@__DIR__, "_data", "testdata") - # Test 2: Test with different data types - @testset "Different data types" begin - test_dir = mktempdir() - sample_file = joinpath(test_dir, "sample_float32.fits") - f = FITS(sample_file, "w") - write(f, Int[]) - # Create data - times = Float64[1.0, 2.0, 3.0] - energies = Float64[10.0, 20.0, 30.0] - table = Dict{String,Array}() - table["TIME"] = times - table["ENERGY"] = energies - write(f, table) - close(f) - # Test with Float32 - data_f32 = readevents(sample_file, T = Float32) - @test eltype(data_f32.times) == Float32 - @test eltype(data_f32.energies) == Float32 - # Test with Int64 - data_i64 = readevents(sample_file, T = Int64) - @test eltype(data_i64.times) == Int64 - @test eltype(data_i64.energies) == Int64 - end +# Helper function to generate mock data +function mock_data(times, energies; energy_column = "ENERGY") + test_dir = mktempdir() + sample_file = joinpath(test_dir, "sample.fits") + + # Create a sample FITS file + FITS(sample_file, "w") do f + # Create primary HDU with a small array instead of empty + # Use a single element array instead of empty + write(f, [0]) - # Test 3: Test with missing columns - @testset "Missing columns" begin - test_dir = mktempdir() - sample_file = joinpath(test_dir, "sample_no_energy.fits") - # Create a sample FITS file with only TIME column - f = FITS(sample_file, "w") - write(f, Int[]) - times = Float64[1.0, 2.0, 3.0] + # Create event table in HDU 2 table = Dict{String,Array}() table["TIME"] = times + table[energy_column] = energies + table["INDEX"] = collect(1:length(times)) write(f, table) - close(f) - local data - @test_logs (:warn, "No ENERGY data found in FITS file $(sample_file). Energy spectrum analysis will not be possible.") begin - data = readevents(sample_file) - end - @test length(data.times) == 3 - @test length(data.energies) == 0 - #create a file with only ENERGY column - sample_file2 = joinpath(test_dir, "sample_no_time.fits") - f = FITS(sample_file2, "w") - write(f, Int[]) # Empty primary array - energies = Float64[10.0, 20.0, 30.0] - table = Dict{String,Array}() - table["ENERGY"] = energies - write(f, table) - close(f) - local data2 - @test_logs (:warn, "No TIME data found in FITS file $(sample_file2). Time series analysis will not be possible.") begin - data2 = readevents(sample_file2) - end - @test length(data2.times) == 0 # No TIME column - @test length(data2.energies) == 3 end + sample_file +end - # Test 4: Test with multiple HDUs - @testset "Multiple HDUs" begin - test_dir = mktempdir() - sample_file = joinpath(test_dir, "sample_multi_hdu.fits") - # Create a sample FITS file with multiple HDUs - f = FITS(sample_file, "w") - write(f, Int[]) - times1 = Float64[1.0, 2.0, 3.0] - energies1 = Float64[10.0, 20.0, 30.0] - table1 = Dict{String,Array}() - table1["TIME"] = times1 - table1["ENERGY"] = energies1 - write(f, table1) - # Second table HDU (with OTHER column) - other_data = Float64[100.0, 200.0, 300.0] - table2 = Dict{String,Array}() - table2["OTHER"] = other_data - write(f, table2) - # Third table HDU (with TIME only) - times3 = Float64[4.0, 5.0, 6.0] - table3 = Dict{String,Array}() - table3["TIME"] = times3 - write(f, table3) - close(f) - data = readevents(sample_file) - @test length(data.metadata.headers) == 4 # Primary + 3 table HDUs - # Should read the first HDU with both TIME and ENERGY - @test length(data.times) == 3 - @test length(data.energies) == 3 - end +# Test basic EventList creation and validation +let + # Test valid construction with simplified constructor + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 15.0, 25.0, 30.0] - # Test 5: Test with monol_testA.evt - @testset "test monol_testA.evt" begin - test_filepath = joinpath("data", "monol_testA.evt") - if isfile(test_filepath) - @testset "monol_testA.evt" begin - old_logger = global_logger(ConsoleLogger(stderr, Logging.Error)) - try - data = readevents(test_filepath) - @test data.filename == test_filepath - @test length(data.metadata.headers) > 0 - finally - global_logger(old_logger) - end - end - else - @info "Test file '$(test_filepath)' not found. Skipping this test." - end - end + ev = EventList(times, energies) + @test ev.times == times + @test ev.energies == energies + @test length(ev) == 5 + @test has_energies(ev) + + # Test with no energies + ev_no_energy = EventList(times) + @test ev_no_energy.times == times + @test isnothing(ev_no_energy.energies) + @test !has_energies(ev_no_energy) + +end + +# Test accessor functions +let + times_vec = [1.0, 2.0, 3.0] + energies_vec = [10.0, 20.0, 30.0] + + ev = EventList(times_vec, energies_vec) + + # Test times() accessor + @test times(ev) === ev.times + @test times(ev) == times_vec + + # Test energies() accessor + @test energies(ev) === ev.energies + @test energies(ev) == energies_vec + + # Test with nothing energies + ev_no_energy = EventList(times_vec) + @test isnothing(energies(ev_no_energy)) + +end + +# Test Base interface methods +let + times_vec = [1.0, 2.0, 3.0, 4.0] + + ev = EventList(times_vec) + + # Test length + @test length(ev) == 4 + @test length(ev) == length(times_vec) + + # Test size + @test size(ev) == (4,) + @test size(ev) == (length(times_vec),) + +end + +# Test filter_time! function (in-place filtering by time) +let + # Test basic time filtering - keep times >= 3.0 + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + + ev = EventList(times, energies) + filter_time!(t -> t >= 3.0, ev) + + @test ev.times == [3.0, 4.0, 5.0] + @test ev.energies == [30.0, 40.0, 50.0] + @test length(ev) == 3 + + # Test filtering that removes all elements + ev_empty = EventList([1.0, 2.0], [10.0, 20.0]) + filter_time!(t -> t > 10.0, ev_empty) + @test length(ev_empty) == 0 + @test ev_empty.times == Float64[] + @test ev_empty.energies == Float64[] + + # Test filtering with no energies + ev_no_energy = EventList([1.0, 2.0, 3.0, 4.0]) + filter_time!(t -> t > 2.0, ev_no_energy) + @test ev_no_energy.times == [3.0, 4.0] + @test isnothing(ev_no_energy.energies) + @test length(ev_no_energy) == 2 + + # Test filtering with extra columns + times_extra = [1.0, 2.0, 3.0, 4.0] + energies_extra = [10.0, 20.0, 30.0, 40.0] + dummy_meta = FITSMetadata{Dict{String,Any}}( + "", + 1, + nothing, + Dict("INDEX" => [1, 2, 3, 4]), + Dict{String,Any}(), + ) + ev_extra = EventList(times_extra, energies_extra, dummy_meta) + + filter_time!(t -> t >= 2.5, ev_extra) + @test ev_extra.times == [3.0, 4.0] + @test ev_extra.energies == [30.0, 40.0] + @test ev_extra.meta.extra_columns["INDEX"] == [3, 4] + +end + +# Test filter_energy! function (in-place filtering by energy) +let + # Test basic energy filtering - keep energies >= 25.0 + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + + ev = EventList(times, energies) + filter_energy!(e -> e >= 25.0, ev) + + @test ev.times == [3.0, 4.0, 5.0] + @test ev.energies == [30.0, 40.0, 50.0] + @test length(ev) == 3 + + # Test filtering that removes all elements + ev_all_removed = EventList([1.0, 2.0], [10.0, 20.0]) + filter_energy!(e -> e > 100.0, ev_all_removed) + @test length(ev_all_removed) == 0 + @test ev_all_removed.times == Float64[] + @test ev_all_removed.energies == Float64[] + + # Test error when no energies present + ev_no_energy = EventList([1.0, 2.0, 3.0]) + @test_throws AssertionError filter_energy!(e -> e > 10.0, ev_no_energy) + + # Test filtering with extra columns + times_extra = [1.0, 2.0, 3.0, 4.0] + energies_extra = [15.0, 25.0, 35.0, 45.0] + dummy_meta = FITSMetadata{Dict{String,Any}}( + "", + 1, + nothing, + Dict("DETX" => [0.1, 0.2, 0.3, 0.4]), + Dict{String,Any}(), + ) + ev_extra = EventList(times_extra, energies_extra, dummy_meta) + + filter_energy!(e -> e >= 30.0, ev_extra) + @test ev_extra.times == [3.0, 4.0] + @test ev_extra.energies == [35.0, 45.0] + @test ev_extra.meta.extra_columns["DETX"] == [0.3, 0.4] + +end + +# Test filter_on! function (generic in-place filtering) +let + # Test filtering on times using filter_on! + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + + ev = EventList(times, energies) + filter_on!(t -> t % 2 == 0, ev.times, ev) # Keep even time values (2.0, 4.0) + + @test ev.times == [2.0, 4.0] + @test ev.energies == [20.0, 40.0] + @test length(ev) == 2 + + # Test filtering on energies using filter_on! + times2 = [1.0, 2.0, 3.0, 4.0] + energies2 = [15.0, 25.0, 35.0, 45.0] + ev2 = EventList(times2, energies2) + filter_on!(e -> e > 30.0, ev2.energies, ev2) # Keep energies > 30 + + @test ev2.times == [3.0, 4.0] + @test ev2.energies == [35.0, 45.0] + + # Test assertion error for mismatched sizes + times3 = [1.0, 2.0, 3.0] + energies3 = [10.0, 20.0, 30.0] + ev3 = EventList(times3, energies3) + wrong_size_col = [1.0, 2.0] # Wrong size + @test_throws AssertionError filter_on!(x -> x > 1.0, wrong_size_col, ev3) + + # Test with extra columns + times_extra = [1.0, 2.0, 3.0, 4.0, 5.0] + energies_extra = [10.0, 20.0, 30.0, 40.0, 50.0] + dummy_meta = FITSMetadata{Dict{String,Any}}( + "", + 1, + nothing, + Dict("FLAG" => [1, 0, 1, 0, 1]), + Dict{String,Any}(), + ) + ev_extra = EventList(times_extra, energies_extra, dummy_meta) - @testset "Error handling" begin - # Test with non-existent file - using a more generic approach - @test_throws Exception readevents("non_existent_file.fits") + # Filter based on FLAG column (keep where FLAG == 1) + filter_on!(flag -> flag == 1, ev_extra.meta.extra_columns["FLAG"], ev_extra) + @test ev_extra.times == [1.0, 3.0, 5.0] + @test ev_extra.energies == [10.0, 30.0, 50.0] + @test ev_extra.meta.extra_columns["FLAG"] == [1, 1, 1] - # Test with invalid FITS file - invalid_file = tempname() - open(invalid_file, "w") do io - write(io, "This is not a FITS file") - end - @test_throws Exception readevents(invalid_file) +end + +# Test non-mutating filter functions (filter_time and filter_energy) +let + # Test filter_time (non-mutating) + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + + ev_original = EventList(times, energies) + ev_filtered = filter_time(t -> t >= 3.0, ev_original) + + # Original should be unchanged + @test ev_original.times == [1.0, 2.0, 3.0, 4.0, 5.0] + @test ev_original.energies == [10.0, 20.0, 30.0, 40.0, 50.0] + @test length(ev_original) == 5 + + # Filtered should have new values + @test ev_filtered.times == [3.0, 4.0, 5.0] + @test ev_filtered.energies == [30.0, 40.0, 50.0] + @test length(ev_filtered) == 3 + + # Test filter_energy (non-mutating) + ev_original2 = EventList(times, energies) + ev_filtered2 = filter_energy(e -> e <= 30.0, ev_original2) + + # Original should be unchanged + @test ev_original2.times == times + @test ev_original2.energies == energies + + # Filtered should have new values + @test ev_filtered2.times == [1.0, 2.0, 3.0] + @test ev_filtered2.energies == [10.0, 20.0, 30.0] + @test length(ev_filtered2) == 3 + + # Test with no energies + ev_no_energy = EventList([1.0, 2.0, 3.0, 4.0]) + ev_filtered_no_energy = filter_time(t -> t > 2.0, ev_no_energy) + + @test ev_no_energy.times == [1.0, 2.0, 3.0, 4.0] # Original unchanged + @test ev_filtered_no_energy.times == [3.0, 4.0] + @test isnothing(ev_filtered_no_energy.energies) + + # Test filter_energy error with no energies + @test_throws AssertionError filter_energy(e -> e > 10.0, ev_no_energy) + +end + +# Test complex filtering scenarios +let + # Test multiple sequential filters + times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] + energies = [10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0] + + ev = EventList(times, energies) + + # First filter by time (keep t >= 3.0) + filter_time!(t -> t >= 3.0, ev) + @test length(ev) == 6 + @test ev.times == [3.0, 4.0, 5.0, 6.0, 7.0, 8.0] + @test ev.energies == [20.0, 25.0, 30.0, 35.0, 40.0, 45.0] + + # Then filter by energy (keep e <= 35.0) + filter_energy!(e -> e <= 35.0, ev) + @test length(ev) == 4 + @test ev.times == [3.0, 4.0, 5.0, 6.0] + @test ev.energies == [20.0, 25.0, 30.0, 35.0] + + # Test edge case: empty result after filtering + ev_edge = EventList([1.0, 2.0], [10.0, 20.0]) + filter_time!(t -> t > 5.0, ev_edge) + @test length(ev_edge) == 0 + @test ev_edge.times == Float64[] + @test ev_edge.energies == Float64[] + + # Test edge case: no filtering needed (all pass) + ev_all_pass = EventList([1.0, 2.0, 3.0], [10.0, 20.0, 30.0]) + original_times = copy(ev_all_pass.times) + original_energies = copy(ev_all_pass.energies) + filter_time!(t -> t > 0.0, ev_all_pass) # All should pass + @test ev_all_pass.times == original_times + @test ev_all_pass.energies == original_energies + @test length(ev_all_pass) == 3 + +end + +# Test readevents basic functionality with mock data +let + mock_times = Float64[1.0, 2.0, 3.0, 4.0, 5.0] + mock_energies = Float64[10.0, 20.0, 15.0, 25.0, 30.0] + sample = mock_data(mock_times, mock_energies) + + # Try reading mock data + data = readevents(sample) + @test data.times == mock_times + @test data.energies == mock_energies + @test length(data.meta.headers) >= 1 # Should have at least primary header + @test length(data.meta.extra_columns) == 0 + + # Test reading with extra columns + data = readevents(sample; extra_columns = ["INDEX"]) + @test length(data.meta.extra_columns) == 1 + @test haskey(data.meta.extra_columns, "INDEX") + @test data.meta.extra_columns["INDEX"] == collect(1:5) + +end + +# Test readevents HDU handling +let + # Test with events in HDU 3 instead of default HDU 2 + test_dir = mktempdir() + sample_file = joinpath(test_dir, "hdu3_sample.fits") + f = FITS(sample_file, "w") + write(f, [0]) # Primary HDU with non-empty array + + # Empty table in HDU 2 + empty_table = Dict{String,Array}() + empty_table["OTHER"] = Float64[1.0, 2.0] + write(f, empty_table) + + # Event data in HDU 3 + times = Float64[1.0, 2.0, 3.0] + energies = Float64[10.0, 20.0, 30.0] + event_table = Dict{String,Array}() + event_table["TIME"] = times + event_table["ENERGY"] = energies + write(f, event_table) + close(f) + + # Test specifying specific HDU + data_hdu3 = readevents(sample_file, hdu = 3) + @test data_hdu3.times == times + @test data_hdu3.energies == energies + +end + +# Test readevents alternative energy columns +let + # Test with PI column + times = Float64[1.0, 2.0, 3.0] + pi_values = Float64[100.0, 200.0, 300.0] + + pi_file = mock_data(times, pi_values; energy_column = "PI") + + data = readevents(pi_file) + @test data.times == times + @test data.energies == pi_values + @test data.meta.energy_units == "PI" + + # Test with PHA column + pha_file = mock_data(times, pi_values; energy_column = "PHA") + + data_pha = readevents(pha_file) + @test data_pha.times == times + @test data_pha.energies == pi_values + @test data_pha.meta.energy_units == "PHA" + +end + +# Test readevents missing columns +let + # File with only TIME column + test_dir = mktempdir() + time_only_file = joinpath(test_dir, "time_only.fits") + f = FITS(time_only_file, "w") + write(f, [0]) # Non-empty primary HDU + + times = Float64[1.0, 2.0, 3.0] + table = Dict{String,Array}() + table["TIME"] = times + write(f, table) + close(f) + + data = readevents(time_only_file) + @test data.times == times + @test isnothing(data.energies) + @test isnothing(data.meta.energy_units) + +end + +# Test error handling +let + # Test non-existent file + @test_throws Exception readevents("non_existent_file.fits") + + # Test invalid FITS file + test_dir = mktempdir() + invalid_file = joinpath(test_dir, "invalid.fits") + open(invalid_file, "w") do io + write(io, "This is not a FITS file") end + @test_throws Exception readevents(invalid_file) + +end + +# Test case insensitive column names +let + test_dir = mktempdir() + sample_file = joinpath(test_dir, "case_test.fits") + + f = FITS(sample_file, "w") + + # Create primary HDU with valid data + primary_data = reshape([1.0], 1, 1) + write(f, primary_data) + + # Use lowercase column names + times = Float64[1.0, 2.0, 3.0] + energies = Float64[10.0, 20.0, 30.0] + + table = Dict{String,Array}() + table["time"] = times # lowercase + table["energy"] = energies # lowercase + write(f, table) + close(f) + + data = readevents(sample_file) + @test data.times == times + @test data.energies == energies + end diff --git a/test/test_missionSupport.jl b/test/test_missionSupport.jl new file mode 100644 index 0000000..66a060e --- /dev/null +++ b/test/test_missionSupport.jl @@ -0,0 +1,572 @@ +""" + create_synthetic_events(n_events::Int=1000; mission::String="nustar", + time_range::Tuple{Float64,Float64}=(0.0, 1000.0), + pi_range::Tuple{Int,Int}=(1, 4096), + T::Type=Float64) -> EventList{Vector{T}, FITSMetadata} + +Create synthetic X-ray event data for testing purposes. + +This function generates synthetic X-ray event data with proper time ordering, +energy calibration, and mission-specific metadata. The generated events include: +- Random event times within the specified range +- Random PI (Pulse Invariant) channels +- Energy values calculated using mission-specific calibration +- Detector ID assignments +- Realistic FITS-style metadata headers + +# Arguments +- `n_events::Int=1000`: Number of events to generate +- `mission::String="nustar"`: Mission name for calibration and metadata +- `time_range::Tuple{Float64,Float64}=(0.0, 1000.0)`: Time range for events (start, stop) +- `pi_range::Tuple{Int,Int}=(1, 4096)`: PI channel range +- `T::Type=Float64`: Numeric type for the data + +# Returns +- `EventList{Vector{T}, FITSMetadata}`: Synthetic event list with times, energies, and metadata + +# Examples +```julia +# Default NuSTAR events +events = create_synthetic_events() + +# Custom XMM events +events = create_synthetic_events(5000; mission="xmm", time_range=(100.0, 200.0)) + +# NICER events with custom PI range +events = create_synthetic_events(2000; mission="nicer", pi_range=(50, 1000)) +``` +""" +function create_synthetic_events(n_events::Int=1000; + mission::String="nustar", + time_range::Tuple{Float64,Float64}=(0.0, 1000.0), + pi_range::Tuple{Int,Int}=(1, 4096), + T::Type=Float64) + + # Input validation + if n_events <= 0 + throw(ArgumentError("Number of events must be positive")) + end + + if time_range[1] >= time_range[2] + throw(ArgumentError("Time range start must be less than end")) + end + + if pi_range[1] > pi_range[2] + throw(ArgumentError("PI range start must be less than or equal to end")) + end + + # Get mission support for calibration + mission_support = get_mission_support(mission) + + # Generate random event times and sort them + times = sort(rand(n_events) * (time_range[2] - time_range[1]) .+ time_range[1]) + times = convert(Vector{T}, times) + + # Generate random PI channels + pi_channels = rand(pi_range[1]:pi_range[2], n_events) + + # Apply mission-specific calibration to get energies + energies = convert(Vector{T}, apply_calibration(mission_support, pi_channels)) + + # Generate detector IDs (assume 4-detector system like NuSTAR) + detector_ids = rand(0:3, n_events) + + # Create extra columns dictionary + extra_columns = Dict{String, Vector}( + "PI" => pi_channels, + "DETID" => detector_ids + ) + + # Create synthetic filename + filepath = "synthetic_$(mission).fits" + + # Create FITS headers using FITSIO.FITSHeader structure + # Note: This is a simplified header for testing - in real usage, + # headers would come from actual FITS files + header_dict = Dict{String, Any}( + "TELESCOP" => uppercase(mission), + "INSTRUME" => "TEST_INSTRUMENT", + "OBJECT" => "TEST_SOURCE", + "RA_NOM" => 180.0, + "DEC_NOM" => 0.0, + "EQUINOX" => 2000.0, + "RADECSYS" => "FK5", + "TSTART" => time_range[1], + "TSTOP" => time_range[2], + "EXPOSURE" => time_range[2] - time_range[1], + "ONTIME" => time_range[2] - time_range[1], + "LIVETIME" => time_range[2] - time_range[1], + "NAXIS2" => n_events, + "TIMESYS" => "TT", + "TIMEREF" => "LOCAL", + "TIMEUNIT" => "s", + "MJDREFI" => 55197, # Standard MJD reference for many missions + "MJDREFF" => 0.00076601852, + "CLOCKCORR" => "T", + "DATE-OBS" => "2020-01-01T00:00:00", + "DATE-END" => "2020-01-01T01:00:00" + ) + + # Add mission-specific header information + if lowercase(mission) == "nustar" + header_dict["DETNAM"] = "TEST_DET" + header_dict["PIFLTCOR"] = "T" + elseif lowercase(mission) == "xmm" + header_dict["FILTER"] = "Medium" + header_dict["SUBMODE"] = "PrimeFullWindow" + elseif lowercase(mission) == "nicer" + header_dict["DETNAM"] = "TEST_NICER" + header_dict["FILTFILE"] = "NONE" + elseif lowercase(mission) in ["chandra", "axaf"] + header_dict["DETNAM"] = "ACIS-S" + header_dict["GRATING"] = "NONE" + elseif lowercase(mission) == "xte" + header_dict["DETNAM"] = "PCA" + header_dict["LAYERS"] = "ALL" + elseif lowercase(mission) == "ixpe" + header_dict["DETNAM"] = "GPD" + header_dict["POLMODE"] = "ON" + end + + # Create a simple header structure (for testing purposes) + # In real usage, this would be a proper FITSIO.FITSHeader + headers = header_dict + + # Create FITSMetadata + metadata = FITSMetadata( + filepath, # filepath + 2, # hdu (typical for event data) + "ENERGY", # energy_units (after calibration) + extra_columns, # extra_columns + headers # headers + ) + + # Create and return EventList + return EventList(times, energies, metadata) +end +# Test: Basic Synthetic Event Creation - Default Parameters +let + events = create_synthetic_events() + + # Test basic structure + @test isa(events, EventList) + @test length(events) == 1000 # Default n_events + @test length(times(events)) == 1000 + @test length(energies(events)) == 1000 + @test events.meta.filepath == "synthetic_nustar.fits" + + # Test times are sorted and in range + @test issorted(times(events)) + @test all(0.0 .<= times(events) .<= 1000.0) + + # Test energies are positive (after calibration) + @test all(energies(events) .> 0) + + # Test extra columns + @test haskey(events.meta.extra_columns, "PI") + @test haskey(events.meta.extra_columns, "DETID") + @test length(events.meta.extra_columns["PI"]) == 1000 + @test length(events.meta.extra_columns["DETID"]) == 1000 + + # Test PI channels are in expected range + @test all(1 .<= events.meta.extra_columns["PI"] .<= 4096) + + # Test detector IDs are in expected range + @test all(0 .<= events.meta.extra_columns["DETID"] .<= 3) + + # Test metadata structure + @test isa(events.meta, FITSMetadata) + @test events.meta.hdu == 2 + @test events.meta.energy_units == "ENERGY" + + # Test metadata content + @test events.meta.headers["TELESCOP"] == "NUSTAR" + @test events.meta.headers["INSTRUME"] == "TEST_INSTRUMENT" + @test events.meta.headers["NAXIS2"] == 1000 + @test events.meta.headers["TSTART"] == 0.0 + @test events.meta.headers["TSTOP"] == 1000.0 +end + +# Test: Basic Synthetic Event Creation - Custom Parameters +let + n_events = 500 + mission = "xmm" + time_range = (100.0, 200.0) + pi_range = (50, 1000) + + events = create_synthetic_events(n_events; + mission=mission, + time_range=time_range, + pi_range=pi_range) + + # Test custom parameters are respected + @test length(events) == n_events + @test length(times(events)) == n_events + @test length(energies(events)) == n_events + @test events.meta.filepath == "synthetic_xmm.fits" + + # Test time range + @test all(time_range[1] .<= times(events) .<= time_range[2]) + @test events.meta.headers["TSTART"] == time_range[1] + @test events.meta.headers["TSTOP"] == time_range[2] + + # Test PI range + @test all(pi_range[1] .<= events.meta.extra_columns["PI"] .<= pi_range[2]) + + # Test mission-specific metadata + @test events.meta.headers["TELESCOP"] == "XMM" +end + +# Test: Different Missions +let + missions = ["nustar", "xmm", "nicer", "ixpe", "axaf", "chandra", "xte"] + + for mission in missions + events = create_synthetic_events(100; mission=mission) + + @test events.meta.filepath == "synthetic_$(mission).fits" + @test events.meta.headers["TELESCOP"] == uppercase(mission) + @test length(events) == 100 + @test length(times(events)) == 100 + @test length(energies(events)) == 100 + + # Test that calibration was applied correctly + ms = get_mission_support(mission) + expected_energies = apply_calibration(ms, events.meta.extra_columns["PI"]) + @test energies(events) ≈ expected_energies + end +end + +# Test: Time Ordering +let + for _ in 1:10 # Test multiple times due to randomness + events = create_synthetic_events(100) + @test issorted(times(events)) + + # Test no duplicate times (very unlikely but possible) + @test length(unique(times(events))) >= 95 # Allow some duplicates due to floating point + end +end + +# Test: Energy Calibration Consistency +let + missions = ["nustar", "xmm", "nicer", "ixpe"] + + for mission in missions + events = create_synthetic_events(200; mission=mission) + + # Manually verify calibration + ms = get_mission_support(mission) + expected_energies = apply_calibration(ms, events.meta.extra_columns["PI"]) + @test energies(events) ≈ expected_energies + + # Test energy ranges are reasonable for each mission + if mission == "nustar" + # NuSTAR: pi * 0.04 + 1.62, PI range 1-4096 + @test minimum(energies(events)) >= 1.66 # 1*0.04 + 1.62 + @test maximum(energies(events)) <= 165.46 # 4096*0.04 + 1.62 + elseif mission == "xmm" + # XMM: pi * 0.001, PI range 1-4096 + @test minimum(energies(events)) >= 0.001 + @test maximum(energies(events)) <= 4.096 + elseif mission == "nicer" + # NICER: pi * 0.01, PI range 1-4096 + @test minimum(energies(events)) >= 0.01 + @test maximum(energies(events)) <= 40.96 + elseif mission == "ixpe" + # IXPE: pi / 375 * 15, PI range 1-4096 + @test minimum(energies(events)) >= 15.0/375 # ≈ 0.04 + @test maximum(energies(events)) <= 4096*15.0/375 # ≈ 163.84 + end + end +end + +# Test: Statistical Properties +let + events = create_synthetic_events(10000) # Large sample for statistics + + # Test time distribution (should be roughly uniform) + time_hist = fit(Histogram, times(events), 0:100:1000) + counts = time_hist.weights + # Expect roughly equal counts in each bin (within statistical fluctuation) + expected_count = 10000 / 10 # 10 bins + @test all(abs.(counts .- expected_count) .< 3 * sqrt(expected_count)) # 3-sigma test + + # Test PI distribution (should be roughly uniform over discrete range) + pi_hist = fit(Histogram, events.meta.extra_columns["PI"], 1:100:4096) + pi_counts = pi_hist.weights + # More lenient test due to discrete uniform distribution + @test std(pi_counts) / mean(pi_counts) < 0.2 # Coefficient of variation < 20% + + # Test detector ID distribution + detid_counts = [count(==(i), events.meta.extra_columns["DETID"]) for i in 0:3] + @test all(abs.(detid_counts .- 2500) .< 3 * sqrt(2500)) # Each detector ~2500 events +end + +# Test: Small Event Counts +let + for n in [1, 2, 5, 10] + events = create_synthetic_events(n) + @test length(events) == n + @test length(times(events)) == n + @test length(energies(events)) == n + @test length(events.meta.extra_columns["PI"]) == n + @test length(events.meta.extra_columns["DETID"]) == n + @test issorted(times(events)) + end +end + +# Test: Large Event Counts +let + events = create_synthetic_events(100000) + @test length(events) == 100000 + @test length(times(events)) == 100000 + @test length(energies(events)) == 100000 + @test issorted(times(events)) + @test events.meta.headers["NAXIS2"] == 100000 +end + +# Test: Extreme Time Ranges +let + # Very short time range + events = create_synthetic_events(100; time_range=(0.0, 0.1)) + @test all(0.0 .<= times(events) .<= 0.1) + @test events.meta.headers["TSTART"] == 0.0 + @test events.meta.headers["TSTOP"] == 0.1 + + # Very long time range + events = create_synthetic_events(100; time_range=(0.0, 1e6)) + @test all(0.0 .<= times(events) .<= 1e6) + @test events.meta.headers["TSTART"] == 0.0 + @test events.meta.headers["TSTOP"] == 1e6 + + # Negative time range + events = create_synthetic_events(100; time_range=(-1000.0, -500.0)) + @test all(-1000.0 .<= times(events) .<= -500.0) + @test issorted(times(events)) +end + +# Test: Extreme PI Ranges +let + # Small PI range + events = create_synthetic_events(100; pi_range=(100, 110)) + @test all(100 .<= events.meta.extra_columns["PI"] .<= 110) + + # Single PI value + events = create_synthetic_events(100; pi_range=(500, 500)) + @test all(events.meta.extra_columns["PI"] .== 500) + @test all(energies(events) .== energies(events)[1]) # All same energy + + # Large PI range + events = create_synthetic_events(100; pi_range=(1, 10000)) + @test all(1 .<= events.meta.extra_columns["PI"] .<= 10000) +end + +# Test: Unknown Mission +let + # Should still work but with warning + @test_logs (:warn, r"Mission unknown_mission not recognized") begin + events = create_synthetic_events(100; mission="unknown_mission") + @test events.meta.filepath == "synthetic_unknown_mission.fits" + @test events.meta.headers["TELESCOP"] == "UNKNOWN_MISSION" + @test length(events) == 100 + # Should use identity calibration + @test energies(events) == Float64.(events.meta.extra_columns["PI"]) + end +end + +# Test: No Missing Data +let + events = create_synthetic_events(1000) + + # Check no NaN or missing values + @test all(isfinite.(times(events))) + @test all(isfinite.(energies(events))) + @test all(isfinite.(events.meta.extra_columns["PI"])) + @test all(isfinite.(events.meta.extra_columns["DETID"])) + + # Check no negative energies (after calibration) + @test all(energies(events) .>= 0) +end + +# Test: Correct Data Types +let + events = create_synthetic_events(100) + + @test eltype(times(events)) == Float64 + @test eltype(energies(events)) == Float64 + @test eltype(events.meta.extra_columns["PI"]) <: Integer + @test eltype(events.meta.extra_columns["DETID"]) <: Integer + @test isa(events.meta, FITSMetadata) + @test isa(events.meta.filepath, String) +end + +# Test: Array Length Consistency +let + for n in [10, 100, 1000, 5000] + events = create_synthetic_events(n) + + @test length(events) == n + @test length(times(events)) == n + @test length(energies(events)) == n + @test length(events.meta.extra_columns["PI"]) == n + @test length(events.meta.extra_columns["DETID"]) == n + @test events.meta.headers["NAXIS2"] == n + end +end + +# Test: Mission Name Handling +let + # Test case sensitivity + missions = ["NUSTAR", "nustar", "NuSTAR", "NuStar"] + for mission in missions + events = create_synthetic_events(100; mission=mission) + @test events.meta.headers["TELESCOP"] == "NUSTAR" + @test events.meta.filepath == "synthetic_$(mission).fits" + end +end + +# Test: Calibration Differences +let + # Same PI values should give different energies for different missions + test_pi_range = (1000, 1000) # Fixed PI value + + results = Dict{String, Float64}() + for mission in ["nustar", "xmm", "nicer", "ixpe"] + events = create_synthetic_events(100; mission=mission, pi_range=test_pi_range) + results[mission] = energies(events)[1] # All should be same since PI is fixed + end + + # Different missions should give different energies + missions = collect(keys(results)) + for i in 1:length(missions) + for j in (i+1):length(missions) + @test results[missions[i]] ≠ results[missions[j]] + end + end +end + +# Test: Metadata Consistency +let + missions = ["nustar", "xmm", "nicer", "ixpe", "axaf", "chandra"] + + for mission in missions + events = create_synthetic_events(100; mission=mission) + + @test events.meta.headers["TELESCOP"] == uppercase(mission) + @test events.meta.headers["INSTRUME"] == "TEST_INSTRUMENT" + @test haskey(events.meta.headers, "NAXIS2") + @test haskey(events.meta.headers, "TSTART") + @test haskey(events.meta.headers, "TSTOP") + end +end + +# Test: EventList Interface Methods +let + events = create_synthetic_events(1000) + + # Test length methods + @test length(events) == 1000 + @test size(events) == (1000,) + + # Test accessor methods + @test times(events) === events.times + @test energies(events) === events.energies + @test has_energies(events) == true + + # Test show methods (should not error) + io = IOBuffer() + show(io, MIME"text/plain"(), events) + @test length(String(take!(io))) > 0 + + show(io, MIME"text/plain"(), events.meta) + @test length(String(take!(io))) > 0 +end + +# Test: Large Dataset Creation +let + # Test that large datasets can be created without issues + @time events = create_synthetic_events(50000) + @test length(events) == 50000 + @test sizeof(times(events)) + sizeof(energies(events)) < 1e6 # Less than 1MB for 50k events +end + +# Test: Memory Efficiency +let + # Test that no unnecessary copies are made + events1 = create_synthetic_events(1000) + events2 = create_synthetic_events(1000) + + # Each should have independent data + @test times(events1) !== times(events2) + @test energies(events1) !== energies(events2) + @test events1.meta.extra_columns["PI"] !== events2.meta.extra_columns["PI"] +end + +# Test: Random Seed Behavior +let + # Test that different calls produce different results (random behavior) + events1 = create_synthetic_events(1000) + events2 = create_synthetic_events(1000) + + # Should be different due to randomness + @test times(events1) != times(events2) + @test events1.meta.extra_columns["PI"] != events2.meta.extra_columns["PI"] + @test events1.meta.extra_columns["DETID"] != events2.meta.extra_columns["DETID"] + + # But same structure + @test length(events1) == length(events2) + @test typeof(events1) == typeof(events2) +end + +# Test: Input Validation - Error Cases +let + # Test negative event count + @test_throws ArgumentError create_synthetic_events(-1) + @test_throws ArgumentError create_synthetic_events(0) + + # Test invalid time range + @test_throws ArgumentError create_synthetic_events(100; time_range=(100.0, 50.0)) + @test_throws ArgumentError create_synthetic_events(100; time_range=(100.0, 100.0)) + + # Test invalid PI range + @test_throws ArgumentError create_synthetic_events(100; pi_range=(1000, 500)) +end + +# Test: FITSMetadata Structure +let + events = create_synthetic_events(100) + meta = events.meta + + # Test all required fields are present + @test isa(meta.filepath, String) + @test isa(meta.hdu, Int) + @test isa(meta.energy_units, Union{String, Nothing}) + @test isa(meta.extra_columns, Dict{String, Vector}) + @test !isnothing(meta.headers) + + # Test metadata consistency + @test meta.hdu == 2 + @test meta.energy_units == "ENERGY" + @test length(meta.extra_columns) == 2 # PI and DETID +end + +# Test: Summary Function +let + events = create_synthetic_events(1000) + summary_str = summary(events) + + @test isa(summary_str, String) + @test occursin("1000 events", summary_str) + # The time span should be close to but less than the full range (1000.0) + actual_time_span = maximum(times(events)) - minimum(times(events)) + expected_pattern = "$(actual_time_span) time units" + @test occursin(expected_pattern, summary_str) + + # Alternative: Use regex to match any floating point number + @test occursin(r"\d+\.\d+ time units", summary_str) + + @test occursin("energies:", summary_str) + @test occursin("(ENERGY)", summary_str) + @test occursin("2 extra columns", summary_str) +end \ No newline at end of file