Implement EventList

Project.toml

@@ -16,6 +16,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 
 [compat]
 julia = "1.7"
@@ -30,6 +31,7 @@ HDF5 = "0.16"
 FITSIO = "0.16"
 Intervals = "1.8"
 NaNMath = "0.3, 1"
+Parameters = "0.12"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/Stingray.jl

@@ -1,7 +1,7 @@
 module Stingray
 
 using ResumableFunctions, StatsBase, Statistics, DataFrames
-using FFTW, Metadata, NaNMath, FITSIO, Intervals
+using FFTW, Metadata, NaNMath, FITSIO, Intervals, Parameters
 using ProgressBars: tqdm as show_progress 
 
 include("fourier.jl")
@@ -30,6 +30,18 @@ export get_btis
 export time_intervals_from_gtis
 export bin_intervals_from_gtis
 
+include("lightcurve.jl")
+export LightCurve
+export rebin
+
+include("events.jl")
+export EventList
+export to_lc
+export from_lc
+
+include("io.jl")
+export load_events_from_fits
+
 include("utils.jl")
 
 end 

src/events.jl

@@ -0,0 +1,150 @@
+import Base: read, write, join, sort, sort!
+
+@with_kw mutable struct EventList{T1<:AbstractVector,T2<:AbstractVector,T3<:AbstractMatrix}
+    time::T1
+    energy::T2=Float64[]
+    ncounts::Int=0
+    mjdref::Float64=0
+    dt::Float64=0
+    notes::String=""
+    gti::T3=reshape(Float64[],0,2)
+    PI::Vector{Int}=Int[]
+    high_precision::Bool=false
+    mission::String=""
+    instr::String=""
+    header::String=""
+    detector_id::Vector{String}=String[]
+    ephem::String=""
+    timeref::String=""
+    timesys::String=""
+end
+
+function read(::Type{EventList},filename::String, format::String)
+    if format=="fits"
+        return load_events_from_fits(filename)
+    else
+        throw(ArgumentError("File format $(format) not supported still."))
+    end
+end
+
+function write(ev::EventList, filename::String, format::String)
+    if format=="fits"
+        write_events_to_fits(filename, ev)
+    else
+        throw(ArgumentError("File format $(format) not supported still."))
+    end
+end
+
+function from_lc(lc::LightCurve)
+    times = [lc.time[i] for i in 1:length(lc.time) for _ in 1:lc.counts[i]]
+    return EventList(time=times,gti=lc.gti)
+end
+
+function to_lc(ev::EventList, dt)
+    if isempty(ev.gti)
+        tstart=0
+        tseg=0
+    else
+        tstart = ev.gti[begin][begin]
+        tseg = ev.gti[end][end]-tstart
+    end
+    return to_lc(ev, dt, tstart, tseg)
+end
+
+to_lc(ev::EventList, dt::Real,
+      tstart::Real, tseg::Real) = make_lightcurves(ev.time, dt, ev.gti;
+                                                   tstart=tstart, tseg=tseg,
+                                                   mjdref=ev.mjdref)
+
+function join(ev1::EventList, ev2::EventList)
+    new_ev = EventList(time=vcat(ev1.time, ev2.time))
+    if ev1.dt!=ev2.dt
+        new_ev.dt = max(ev1.dt,ev2.dt);
+    end
+
+    # Tolerance for MJDREF:1 microsecond
+    if !isapprox(ev1.mjdref, ev2.mjdref, atol=1e-6 / 86400)
+        ev2.mjdref = ev1.mjdref
+    end
+
+    order = sortperm(new_ev.time)
+    new_ev.time = new_ev.time[order]
+
+    if !isempty(ev1.PI) || !isempty(ev2.PI)
+        if isempty(ev1.PI)
+            ev1.PI = zero(ev1.time)
+        end
+        if isempty(ev2.PI)
+            ev2.PI = zero(ev2.time)
+        end
+        new_ev.PI = vcat(ev1.PI, ev2.PI)
+        new_ev.PI = new_ev.PI[order]
+    else
+        new_ev.PI = Float64[]
+    end
+
+    if !isempty(ev1.energy) || !isempty(ev2.energy)
+        if isempty(ev1.energy)
+            ev1.energy = zero(ev1.time)
+        end
+        if isempty(ev2.energy)
+            ev2.energy = zero(ev2.time)
+        end
+        new_ev.energy = vcat(ev1.energy, ev2.energy)
+        new_ev.energy = new_ev.energy[order]
+    else
+        new_ev.energy = Float64[]
+    end
+
+    if !isempty(ev1.gti) || !isempty(ev2.gti)
+        if isempty(ev1.gti)
+            ev1.gti = [ev1.time[begin] - ev1.dt/2 ev1.time[end] + ev1.dt/2]
+        end
+        if isempty(ev2.gti)
+            ev2.gti = [ev2.time[begin] - ev2.dt/2 ev2.time[end] + ev2.dt/2]
+        end
+        new_ev.gti = operations_on_gtis([ev1.gti,ev2.gti], intersect)
+        if isempty(new_ev.gti)
+            @warn """
+            GTIs in these two event lists do not overlap at all.
+            Merging instead of returning an overlap.
+            """
+            new_ev.gti = operations_on_gtis([ev1.gti,ev2.gti], union)
+        end
+    else
+        new_ev.gti = reshape(Float64[],0,2)
+    end
+
+    for attr in (:mission, :instr)
+        if getfield(ev1, attr) != getfield(ev2, attr)
+            setfield!(new_ev, attr, string(getfield(ev1, attr),',',getfield(ev2, attr)))
+        else
+            setfield!(new_ev, attr, getfield(ev1, attr))
+        end
+    end
+
+    new_ev.mjdref = ev1.mjdref
+
+    return new_ev
+end
+
+function sort(ev::EventList)
+    order = sortperm(ev.time)
+    new_ev = deepcopy(ev)
+    apply_mask(new_ev, order)
+    return new_ev
+end
+
+function sort!(ev::EventList)
+    order = sortperm(ev.time)
+    apply_mask(ev, order)
+end
+
+function apply_mask(ev::EventList, mask::AbstractVector{T};) where {T<:Union{Bool,Int}}
+    for field in fieldnames(EventList)
+        fval = getfield(ev, field)
+        if fval isa AbstractVector && !isempty(fval) 
+            setfield!(ev, field, deepcopy(fval)[mask])
+        end
+    end
+end

src/gti.jl

@@ -21,11 +21,16 @@ function get_gti_from_hdu(gtihdu::TableHDU)
         stopstr = "Stop"
     end
 
-    gtistart = read(gtihdu,startstr)
-    gtistop = read(gtihdu,stopstr)
+    gtistart = FITSIO.read(gtihdu,startstr)
+    gtistop = FITSIO.read(gtihdu,stopstr)
 
-    return mapreduce(permutedims, vcat, 
-    [[a, b] for (a,b) in zip(gtistart, gtistop)])
+    final_gti = [[a, b] for (a,b) in zip(gtistart, gtistop)]
+
+    if isempty(final_gti)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, final_gti)
+    end
 end
 
 function check_gtis(gti::AbstractMatrix)
@@ -96,7 +101,13 @@ function create_gti_mask(times::AbstractVector{<:Real},gtis::AbstractMatrix{<:Re
         end
     end
 
-    return mask, mapreduce(permutedims, vcat, keepat!(new_gtis,new_gti_mask))
+    keepat!(new_gtis,new_gti_mask)
+
+    if isempty(new_gtis)
+        return mask, reshape(Float64[],0,2)
+    else
+        return mask, mapreduce(permutedims, vcat, new_gtis)
+    end
 end
 
 function create_gti_from_condition(time::AbstractVector{<:Real}, condition::AbstractVector{Bool};
@@ -124,7 +135,11 @@ function create_gti_from_condition(time::AbstractVector{<:Real}, condition::Abst
         end
         push!(gtis,[t0, t1])
     end
-    return mapreduce(permutedims, vcat, gtis)
+    if isempty(gtis)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, gtis)
+    end
 end
 
 function operations_on_gtis(gti_list::AbstractVector{<:AbstractMatrix{T}}, 
@@ -152,7 +167,11 @@ function operations_on_gtis(gti_list::AbstractVector{<:AbstractMatrix{T}},
         push!(final_gti, [first(interval), last(interval)])
     end
 
-    return mapreduce(permutedims, vcat, final_gti)
+    if isempty(final_gti)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, final_gti)
+    end
 end
 
 function get_btis(gtis::AbstractMatrix{<:Real})
@@ -185,7 +204,11 @@ function get_btis(gtis::AbstractMatrix{T}, start_time, stop_time) where {T<:Real
         push!(btis, [first(interval), last(interval)])
     end
 
-    return mapreduce(permutedims, vcat, btis)
+    if isempty(btis)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, btis)
+    end
 end
 
 function time_intervals_from_gtis(gtis::AbstractMatrix{<:Real}, segment_size::Real;

src/io.jl

@@ -0,0 +1,66 @@
+function load_events_from_fits(filename::String)
+    FITS(filename) do hduList
+        eventHDU = hduList["EVENTS"]
+        cols = FITSIO.colnames(eventHDU)
+        df = DataFrame(eventHDU)
+        if "TIME" in cols
+            time = df[!,"TIME"]
+        else
+            throw(KeyError("Time stamp data not provided or column names not appropriate in the file"))
+        end
+        ev = EventList(time=time)
+        if "PI" in cols
+            PI = df[!,"PI"]
+            ev.PI = PI
+        end
+        if "ENERGY" in cols
+            energy = df[!,"ENERGY"]
+            ev.energy = energy
+        end
+        ev.gti = load_gtis(filename)
+        header = FITSIO.read_header(eventHDU)
+        if haskey(header,"MJDREFI") 
+            ev.mjdref += header["MJDREFI"]
+        end
+        if haskey(header,"MJDREFF") 
+            ev.mjdref += header["MJDREFF"]
+        end
+        if haskey(header,"INSTRUME") 
+            ev.instr = header["INSTRUME"]
+        end
+        if haskey(header,"TIMEREF") 
+            ev.timeref = header["TIMEREF"]
+        end
+        if haskey(header,"TIMESYS") 
+            ev.timesys = header["TIMESYS"]
+        end
+        if haskey(header,"PLEPHEM") 
+            ev.ephem = header["PLEPHEM"]
+        end
+        ev
+    end
+end
+
+function write_events_to_fits(filename::String, ev::EventList)
+    FITS(filename,"w") do hduList
+
+        tkeys = String[]
+        values = []
+
+        for field in fieldnames(EventList)
+            fval = getfield(ev, field)
+            if !(fval isa AbstractVecOrMat)
+                push!(tkeys, String(field))
+                push!(values, fval)
+            end
+        end
+
+        header = FITSHeader(tkeys, values,fill("",length(tkeys)))
+
+        eventData = Dict("TIME"=>ev.time, "ENERGY"=>ev.energy, "PI"=>ev.PI);
+        FITSIO.write(hduList, eventData, header = header, name = "EVENTS")
+
+        gtiData = Dict("START"=>ev.gti[:,1], "STOP"=>ev.gti[:,2])
+        FITSIO.write(hduList, gtiData, name = "GTI")
+    end
+end