From fc25de254d08da69bffcfd325fa1208fb393a091 Mon Sep 17 00:00:00 2001 From: Aman-Pandey-afk <77021852+Aman-Pandey-afk@users.noreply.github.com> Date: Mon, 29 Aug 2022 13:17:28 +0530 Subject: [PATCH 1/3] Implement EventList --- Project.toml | 2 + src/Stingray.jl | 8 ++- src/events.jl | 136 ++++++++++++++++++++++++++++++++++++++++++++ src/gti.jl | 10 +++- src/io.jl | 50 ++++++++++++++++ src/lightcurve.jl | 53 +++++++++++++++++ test/runtests.jl | 5 +- test/test_events.jl | 101 ++++++++++++++++++++++++++++++++ 8 files changed, 359 insertions(+), 6 deletions(-) create mode 100644 src/events.jl create mode 100644 src/io.jl create mode 100644 src/lightcurve.jl create mode 100644 test/test_events.jl diff --git a/Project.toml b/Project.toml index cd910a8..de20406 100644 --- a/Project.toml +++ b/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" diff --git a/src/Stingray.jl b/src/Stingray.jl index 740e849..c5177e3 100644 --- a/src/Stingray.jl +++ b/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,12 @@ export get_btis export time_intervals_from_gtis export bin_intervals_from_gtis +include("events.jl") +export EventList + +include("io.jl") +export load_events_from_fits + include("utils.jl") end diff --git a/src/events.jl b/src/events.jl new file mode 100644 index 0000000..5ef7137 --- /dev/null +++ b/src/events.jl @@ -0,0 +1,136 @@ +@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{Int64}=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(filename::String, format::String) + if format=="fits" + return load_events_from_fits(filename) + else + throw(ArgumentError("File 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 not supported still.")) + end +end + +function to_lc(ev::EventList, dt) + tstart = ev.gti[1][1] + tseg = ev.gti[end][end] + return to_lc(ev, dt, tstart, tseg) +end + +to_lc(ev::EventList, dt::Real, + tstart::Real, tseg::Real) = make_lightcurves(ev.time, dt, tstart=tstart, + gti=ev.gti, 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([],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 \ No newline at end of file diff --git a/src/gti.jl b/src/gti.jl index 1d82065..6ffc3d3 100644 --- a/src/gti.jl +++ b/src/gti.jl @@ -21,8 +21,8 @@ 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)]) @@ -152,7 +152,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}) diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 0000000..30b8598 --- /dev/null +++ b/src/io.jl @@ -0,0 +1,50 @@ +function load_events_from_fits(filename::String) + a = 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 "MJDREFI" in header.keys + ev.mjdref += header["MJDREFI"] + end + if "MJDREFF" in header.keys + ev.mjdref += header["MJDREFF"] + end + if "INSTRUME" in header.keys + ev.instr = header["INSTRUME"] + end + if "TIMEREF" in header.keys + ev.timeref = header["TIMEREF"] + end + if "TIMESYS" in header.keys + ev.timesys = header["TIMESYS"] + end + if "PLEPHEM" in header.keys + ev.ephem = header["PLEPHEM"] + end + ev + end + return a +end + +# TODO +function write_events_to_fits(filename::String, ev::EventList) + # FITS(filename,"w") do event + # FITSIO.write(event, ev) + # end +end \ No newline at end of file diff --git a/src/lightcurve.jl b/src/lightcurve.jl new file mode 100644 index 0000000..5d1d8ed --- /dev/null +++ b/src/lightcurve.jl @@ -0,0 +1,53 @@ +@with_kw mutable struct LightCurve{T1<:AbstractVector,T2<:AbstractMatrix} + time::T1 + counts::Int=0 + err::Vector{Float64}=Float64[] + input_counts::Bool=true + gti::T2=reshape(Float64[],0,2) + err_dist::String="" + mjdref::Float64=0 + skip_checks::Bool=false + low_memory::Bool=false + mission::String="" + instr::String="" + header::String="" +end + +function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0.0, + tstart::Real=nothing, gti::AbstractMatrix{<:Real}=nothing, + mjdref::Real=0, use_hist::Bool=false) + sort(toa) + if isnothing(tstart) + # if tstart is not set, assume light curve starts with first photon + tstart = toa[0] + end + + # compute the number of bins in the light curve + # for cases where tseg/dt is not integer. + # TODO: check that this is always consistent and that we + # are not throwing away good events. + + if isnothing(tseg) + tseg = toa[-1] - tstart + end + + timebin = tseg ÷ dt + # If we are missing the next bin by just 1%, let's round up: + if tseg / dt - timebin >= 0.99 + timebin += 1 + end + + tend = tstart + timebin * dt + good = (tstart <= toa) & (toa < tend) + if !use_hist + binned_toas = (toa[good] - tstart) ÷ dt + counts = bincount(binned_toas, minlength=timebin) + time = tstart + range(0.5, 0.5 + len(counts)) * dt + else + histbins = range(tstart, tend + dt, dt) + counts, histbins = histogram(toa[good], bins=histbins) + time = histbins[begin:end] + 0.5 * dt + + return Lightcurve(time, counts; gti=gti, mjdref=mjdref, dt=dt, + skip_checks=True, err_dist="poisson") +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index a4460f4..880e8da 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,5 +2,6 @@ using Stingray using Test using FFTW, Distributions, Statistics, StatsBase, Metadata, HDF5 -include("test_fourier.jl") -include("test_gti.jl") +# include("test_fourier.jl") +# include("test_gti.jl") +include("test_events.jl") diff --git a/test/test_events.jl b/test/test_events.jl new file mode 100644 index 0000000..8a564df --- /dev/null +++ b/test/test_events.jl @@ -0,0 +1,101 @@ +@testset "test_events" begin + times = [0.5, 1.5, 2.5, 3.5] + counts = [3000, 2000, 2200, 3600] + counts_flat = [3000, 3000, 3000, 3000] + spectrum = [[1, 2, 3, 4, 5, 6], [1000, 2040, 1000, 3000, 4020, 2070]] + gti = [0, 4] + + @testset "test_initialize_eventList" begin + times = sort(rand(Uniform(1e8,1e8+100),101)) + ev = EventList(time=times, mjdref=54600) + @test ev.time ≈ times atol = 1e-15 + @test ev.mjdref == 54600 + end + + @testset "test_join_events" begin + @testset "test_join_different_dt" begin + ev = EventList(time=[10, 20, 30], dt=1) + ev_other = EventList(time=[40, 50, 60], dt=3) + ev_new = Stingray.join(ev,ev_other) + @test ev_new.dt == 3 + end + @testset "test_join_different_instr" begin + ev = EventList(time=[10, 20, 30], instr="fpma") + ev_other = EventList(time=[40, 50, 60], instr="fpmb") + ev_new = Stingray.join(ev,ev_other) + @test ev_new.instr == "fpma,fpmb" + end + @testset "test_join_without_energy" begin + ev = EventList(time=[1, 2, 3], energy=[3, 3, 3]) + ev_other = EventList(time=[4, 5]) + ev_new = Stingray.join(ev, ev_other) + @test ev_new.energy == [3, 3, 3, 0, 0] + end + @testset "test_join_without_pi" begin + ev = EventList(time=[1, 2, 3], PI=[3, 3, 3]) + ev_other = EventList(time=[4, 5]) + ev_new = Stingray.join(ev, ev_other) + @test ev_new.PI == [3, 3, 3, 0, 0] + end + @testset "test_join_with_gti_none" begin + ev = EventList(time=[1, 2, 3]) + ev_other = EventList(time=[4, 5], gti=[3.5 5.5]) + ev_new = Stingray.join(ev, ev_other) + @test ev_new.gti == [[1 3]; [3.5 5.5];] + + ev = EventList(time=[1, 2, 3], gti=[0.5 3.5]) + ev_other = EventList(time=[4, 5]) + ev_new = Stingray.join(ev,ev_other) + @test ev_new.gti == [[0.5 3.5]; [4 5];] + + ev = EventList(time=[1, 2, 3]) + ev_other = EventList(time=[4, 5]) + ev_new = Stingray.join(ev, ev_other) + @test isempty(ev_new.gti) + end + @testset "test_non_overlapping_join" begin + ev = EventList(time=[1, 1, 2, 3, 4], + energy=[3, 4, 7, 4, 3], gti=[[1 2]; [3 4];]) + ev_other = EventList(time=[5, 6, 6, 7, 10], + energy=[4, 3, 8, 1, 2], gti=[6 7]) + ev_new = Stingray.join(ev, ev_other) + + @test ev_new.time == [1, 1, 2, 3, 4, 5, 6, 6, 7, 10] + @test ev_new.energy == [3, 4, 7, 4, 3, 4, 3, 8, 1, 2] + @test ev_new.gti == [[1 2]; [3 4]; [6 7];] + end + @testset "test_overlapping_join" begin + ev = EventList(time=[1, 1, 10, 6, 5], + energy=[10, 6, 3, 11, 2], gti=[[1 3]; [5 6];]) + ev_other = EventList(time=[5, 7, 6, 6, 10], + energy=[2, 3, 8, 1, 2], gti=[[5 7]; [8 10];]) + ev_new = Stingray.join(ev, ev_other) + + @test ev_new.time == [1, 1, 5, 5, 6, 6, 6, 7, 10, 10] + @test ev_new.energy == [10, 6, 2, 2, 11, 8, 1, 3, 3, 2] + @test ev_new.gti == [5 6] + end + end + @testset "test_sort" begin + ev = EventList(time=[2, 4, 1, 4, 3], energy=[3, 4, 7, 4, 3]) + @testset "test_sort_not_inplace" begin + new_ev = Stingray.sort(ev); + @test new_ev.time == [1, 2, 3, 4, 4] + @test new_ev.energy == [7, 3, 3, 4, 4] + @test ev.time == [2, 4, 1, 4, 3] + end + @testset "test_sort_inplace" begin + Stingray.sort!(ev); + @test ev.time == [1, 2, 3, 4, 4] + @test ev.energy == [7, 3, 3, 4, 4] + end + end + @testset "Input/Output" begin + @testset "test_fits_with_standard_file" begin + fname = joinpath(@__DIR__ ,"data","monol_testA.evt") + ev = Stingray.read(fname, "fits") + @test ev.mjdref == 55197.00076601852 + end + end +end + From 6daf452716b30238a355c0ea6ac38f0714f81f8b Mon Sep 17 00:00:00 2001 From: Aman-Pandey-afk <77021852+Aman-Pandey-afk@users.noreply.github.com> Date: Fri, 2 Sep 2022 12:00:29 +0530 Subject: [PATCH 2/3] Implement File writing --- src/gti.jl | 29 ++++++++++++++++++++++++----- src/io.jl | 37 +++++++++++++++++++++++++++---------- src/lightcurve.jl | 13 +++++++------ test/runtests.jl | 4 ++-- test/test_events.jl | 10 +++++++++- 5 files changed, 69 insertions(+), 24 deletions(-) diff --git a/src/gti.jl b/src/gti.jl index 6ffc3d3..6b380bd 100644 --- a/src/gti.jl +++ b/src/gti.jl @@ -24,8 +24,13 @@ function get_gti_from_hdu(gtihdu::TableHDU) 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}}, @@ -189,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; diff --git a/src/io.jl b/src/io.jl index 30b8598..50b35a7 100644 --- a/src/io.jl +++ b/src/io.jl @@ -19,22 +19,22 @@ function load_events_from_fits(filename::String) end ev.gti = load_gtis(filename) header = FITSIO.read_header(eventHDU) - if "MJDREFI" in header.keys + if haskey(header,"MJDREFI") ev.mjdref += header["MJDREFI"] end - if "MJDREFF" in header.keys + if haskey(header,"MJDREFF") ev.mjdref += header["MJDREFF"] end - if "INSTRUME" in header.keys + if haskey(header,"INSTRUME") ev.instr = header["INSTRUME"] end - if "TIMEREF" in header.keys + if haskey(header,"TIMEREF") ev.timeref = header["TIMEREF"] end - if "TIMESYS" in header.keys + if haskey(header,"TIMESYS") ev.timesys = header["TIMESYS"] end - if "PLEPHEM" in header.keys + if haskey(header,"PLEPHEM") ev.ephem = header["PLEPHEM"] end ev @@ -42,9 +42,26 @@ function load_events_from_fits(filename::String) return a end -# TODO function write_events_to_fits(filename::String, ev::EventList) - # FITS(filename,"w") do event - # FITSIO.write(event, ev) - # end + 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 \ No newline at end of file diff --git a/src/lightcurve.jl b/src/lightcurve.jl index 5d1d8ed..64051c9 100644 --- a/src/lightcurve.jl +++ b/src/lightcurve.jl @@ -19,7 +19,7 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0 sort(toa) if isnothing(tstart) # if tstart is not set, assume light curve starts with first photon - tstart = toa[0] + tstart = toa[begin] end # compute the number of bins in the light curve @@ -28,7 +28,7 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0 # are not throwing away good events. if isnothing(tseg) - tseg = toa[-1] - tstart + tseg = toa[end-1] - tstart end timebin = tseg ÷ dt @@ -38,15 +38,16 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0 end tend = tstart + timebin * dt - good = (tstart <= toa) & (toa < tend) + good = tstart <= toa < tend if !use_hist - binned_toas = (toa[good] - tstart) ÷ dt + binned_toas = (keepat!(toa,good) .- tstart) .÷ dt counts = bincount(binned_toas, minlength=timebin) time = tstart + range(0.5, 0.5 + len(counts)) * dt else histbins = range(tstart, tend + dt, dt) - counts, histbins = histogram(toa[good], bins=histbins) - time = histbins[begin:end] + 0.5 * dt + counts, histbins = histogram(keepat!(toa,good), bins=histbins) + time = @view(histbins[begin:end-1]) + 0.5 * dt + end return Lightcurve(time, counts; gti=gti, mjdref=mjdref, dt=dt, skip_checks=True, err_dist="poisson") diff --git a/test/runtests.jl b/test/runtests.jl index 880e8da..9b9f4b6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,6 @@ using Stingray using Test using FFTW, Distributions, Statistics, StatsBase, Metadata, HDF5 -# include("test_fourier.jl") -# include("test_gti.jl") +include("test_fourier.jl") +include("test_gti.jl") include("test_events.jl") diff --git a/test/test_events.jl b/test/test_events.jl index 8a564df..5390879 100644 --- a/test/test_events.jl +++ b/test/test_events.jl @@ -1,5 +1,5 @@ @testset "test_events" begin - times = [0.5, 1.5, 2.5, 3.5] + time = [0.5, 1.5, 2.5, 3.5] counts = [3000, 2000, 2200, 3600] counts_flat = [3000, 3000, 3000, 3000] spectrum = [[1, 2, 3, 4, 5, 6], [1000, 2040, 1000, 3000, 4020, 2070]] @@ -96,6 +96,14 @@ ev = Stingray.read(fname, "fits") @test ev.mjdref == 55197.00076601852 end + + @testset "test_io_with_fits" begin + ev = EventList(time=time, mjdref=54000) + Stingray.write(ev, "ev.fits", "fits") + new_ev = Stingray.read("ev.fits", "fits") + @test new_ev.time == time + rm("ev.fits") + end end end From cdeaf8ea0714365a63d606f0dca0e9a89a5bca8a Mon Sep 17 00:00:00 2001 From: Aman-Pandey-afk <77021852+Aman-Pandey-afk@users.noreply.github.com> Date: Wed, 7 Sep 2022 01:40:46 +0530 Subject: [PATCH 3/3] Implement LightCurve Rebin --- src/Stingray.jl | 6 +++ src/events.jl | 40 +++++++++++++------ src/io.jl | 5 +-- src/lightcurve.jl | 82 +++++++++++++++++++++++++++++--------- src/utils.jl | 87 +++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/test_events.jl | 23 +++++++++-- test/test_lightcurve.jl | 49 +++++++++++++++++++++++ 8 files changed, 255 insertions(+), 38 deletions(-) create mode 100644 test/test_lightcurve.jl diff --git a/src/Stingray.jl b/src/Stingray.jl index c5177e3..414a981 100644 --- a/src/Stingray.jl +++ b/src/Stingray.jl @@ -30,8 +30,14 @@ 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 diff --git a/src/events.jl b/src/events.jl index 5ef7137..6624a2a 100644 --- a/src/events.jl +++ b/src/events.jl @@ -1,3 +1,5 @@ +import Base: read, write, join, sort, sort! + @with_kw mutable struct EventList{T1<:AbstractVector,T2<:AbstractVector,T3<:AbstractMatrix} time::T1 energy::T2=Float64[] @@ -6,7 +8,7 @@ dt::Float64=0 notes::String="" gti::T3=reshape(Float64[],0,2) - PI::Vector{Int64}=Int[] + PI::Vector{Int}=Int[] high_precision::Bool=false mission::String="" instr::String="" @@ -17,11 +19,11 @@ timesys::String="" end -function read(filename::String, format::String) +function read(::Type{EventList},filename::String, format::String) if format=="fits" return load_events_from_fits(filename) else - throw(ArgumentError("File not supported still.")) + throw(ArgumentError("File format $(format) not supported still.")) end end @@ -29,19 +31,29 @@ function write(ev::EventList, filename::String, format::String) if format=="fits" write_events_to_fits(filename, ev) else - throw(ArgumentError("File not supported still.")) + 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) - tstart = ev.gti[1][1] - tseg = ev.gti[end][end] + 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, tstart=tstart, - gti=ev.gti, tseg=tseg, + tstart::Real, tseg::Real) = make_lightcurves(ev.time, dt, ev.gti; + tstart=tstart, tseg=tseg, mjdref=ev.mjdref) function join(ev1::EventList, ev2::EventList) @@ -93,15 +105,17 @@ function join(ev1::EventList, ev2::EventList) 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" + @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([],0,2) + new_ev.gti = reshape(Float64[],0,2) end - for attr in [:mission, :instr] + for attr in (:mission, :instr) if getfield(ev1, attr) != getfield(ev2, attr) setfield!(new_ev, attr, string(getfield(ev1, attr),',',getfield(ev2, attr))) else @@ -133,4 +147,4 @@ function apply_mask(ev::EventList, mask::AbstractVector{T};) where {T<:Union{Boo setfield!(ev, field, deepcopy(fval)[mask]) end end -end \ No newline at end of file +end diff --git a/src/io.jl b/src/io.jl index 50b35a7..57e2e57 100644 --- a/src/io.jl +++ b/src/io.jl @@ -1,5 +1,5 @@ function load_events_from_fits(filename::String) - a = FITS(filename) do hduList + FITS(filename) do hduList eventHDU = hduList["EVENTS"] cols = FITSIO.colnames(eventHDU) df = DataFrame(eventHDU) @@ -39,7 +39,6 @@ function load_events_from_fits(filename::String) end ev end - return a end function write_events_to_fits(filename::String, ev::EventList) @@ -64,4 +63,4 @@ function write_events_to_fits(filename::String, ev::EventList) gtiData = Dict("START"=>ev.gti[:,1], "STOP"=>ev.gti[:,2]) FITSIO.write(hduList, gtiData, name = "GTI") end -end \ No newline at end of file +end diff --git a/src/lightcurve.jl b/src/lightcurve.jl index 64051c9..cc57cb1 100644 --- a/src/lightcurve.jl +++ b/src/lightcurve.jl @@ -1,11 +1,12 @@ @with_kw mutable struct LightCurve{T1<:AbstractVector,T2<:AbstractMatrix} time::T1 - counts::Int=0 + counts::Vector{Int}=Int[] err::Vector{Float64}=Float64[] input_counts::Bool=true gti::T2=reshape(Float64[],0,2) err_dist::String="" mjdref::Float64=0 + dt::Float64=0 skip_checks::Bool=false low_memory::Bool=false mission::String="" @@ -13,11 +14,11 @@ header::String="" end -function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0.0, - tstart::Real=nothing, gti::AbstractMatrix{<:Real}=nothing, - mjdref::Real=0, use_hist::Bool=false) - sort(toa) - if isnothing(tstart) +function make_lightcurves(toa::AbstractVector{<:Real}, dt::Real, + gti::AbstractMatrix{<:Real}; tseg::Real=-1, + tstart::Real=-1, mjdref::Real=0) + sort!(toa) + if tstart == -1 # if tstart is not set, assume light curve starts with first photon tstart = toa[begin] end @@ -27,7 +28,7 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0 # TODO: check that this is always consistent and that we # are not throwing away good events. - if isnothing(tseg) + if tseg == -1 tseg = toa[end-1] - tstart end @@ -38,17 +39,62 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0 end tend = tstart + timebin * dt - good = tstart <= toa < tend - if !use_hist - binned_toas = (keepat!(toa,good) .- tstart) .÷ dt - counts = bincount(binned_toas, minlength=timebin) - time = tstart + range(0.5, 0.5 + len(counts)) * dt + good = tstart .<= toa .< tend + + histbins = range(tstart, tend, step = dt) + counts= fit(Histogram, keepat!(toa,good), histbins).weights + time = @view(histbins[begin:end-1]) .+ 0.5 .* dt + + return LightCurve(time=time, counts=counts, gti=gti, mjdref=mjdref, dt=dt, + skip_checks=true, err_dist="poisson") +end + +function rebin(lc::LightCurve, dt_new::Real; method::String = "sum") + if dt_new 0 + pop!(ybin) + pop!(ybinerr) + pop!(step_size) + end + + dx_var = Statistics.var(dx_old) / Statistics.mean(dx_old) + + if size(dx_old) == 1 || dx_var < 1e-6 + new_step = step_size[1] + end + + new_x0 = (x[1] - (0.5 * dx_old[1])) + (0.5 * dx_new) + xbin = 1:length(ybin) * dx_new .+ new_x0 + + return xbin, ybin, ybinerr, new_step +end diff --git a/test/runtests.jl b/test/runtests.jl index 9b9f4b6..a864cae 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,3 +5,4 @@ using FFTW, Distributions, Statistics, StatsBase, Metadata, HDF5 include("test_fourier.jl") include("test_gti.jl") include("test_events.jl") +include("test_lightcurve.jl") diff --git a/test/test_events.jl b/test/test_events.jl index 5390879..a86cd0f 100644 --- a/test/test_events.jl +++ b/test/test_events.jl @@ -3,7 +3,7 @@ counts = [3000, 2000, 2200, 3600] counts_flat = [3000, 3000, 3000, 3000] spectrum = [[1, 2, 3, 4, 5, 6], [1000, 2040, 1000, 3000, 4020, 2070]] - gti = [0, 4] + gti = [0 4] @testset "test_initialize_eventList" begin times = sort(rand(Uniform(1e8,1e8+100),101)) @@ -12,6 +12,21 @@ @test ev.mjdref == 54600 end + @testset "test_to_and_from_lightcurve" begin + @testset "test_from_lc" begin + lc = LightCurve(time=[0.5, 1.5, 2.5], counts=[2, 1, 2]) + ev = from_lc(lc) + @test ev.time == [0.5, 0.5, 1.5, 2.5, 2.5] + end + + @testset "test_to_lc" begin + ev = EventList(time=time, gti=gti) + lc = to_lc(ev, 1) + @test lc.time == [0.5, 1.5, 2.5, 3.5] + @test lc.gti == gti + end + end + @testset "test_join_events" begin @testset "test_join_different_dt" begin ev = EventList(time=[10, 20, 30], dt=1) @@ -93,14 +108,14 @@ @testset "Input/Output" begin @testset "test_fits_with_standard_file" begin fname = joinpath(@__DIR__ ,"data","monol_testA.evt") - ev = Stingray.read(fname, "fits") + ev = read(EventList, fname, "fits") @test ev.mjdref == 55197.00076601852 end @testset "test_io_with_fits" begin ev = EventList(time=time, mjdref=54000) - Stingray.write(ev, "ev.fits", "fits") - new_ev = Stingray.read("ev.fits", "fits") + write(ev, "ev.fits", "fits") + new_ev = read(EventList, "ev.fits", "fits") @test new_ev.time == time rm("ev.fits") end diff --git a/test/test_lightcurve.jl b/test/test_lightcurve.jl new file mode 100644 index 0000000..f9a0018 --- /dev/null +++ b/test/test_lightcurve.jl @@ -0,0 +1,49 @@ +@testset "test_lightcurve" begin + @testset "rebin_lightcurve" begin + dt = 0.0001220703125 + n = 1384132 + mean_counts = 2.0 + times = range(dt/2, dt/2 + n*dt, step = dt) + counts = zero(times) .+ mean_counts + lc = LightCurve(time=times, counts=counts) + + @testset "test_rebin_even" begin + dt_new = 2.0 + lc_binned = rebin(lc, dt_new) + @test lc_binned.dt == dt_new + counts_test = zero(lc_binned.time) .+ lc.counts[1]*dt_new/lc.dt + @test lc_binned.counts == counts_test + end + + @testset "test_rebin_odd" begin + dt_new = 1.5 + lc_binned = rebin(lc, dt_new) + @test lc_binned.dt == dt_new + counts_test = zero(lc_binned.time) .+ lc.counts[1]*dt_new/lc.dt + @test lc_binned.counts == counts_test + end + + @testset "test_rebin_several" for dt in [2,3,pi,5] + lc_binned = rebin(lc, dt) + @test length(lc_binned.time) == length(lc_binned.counts) + end + + @testset "test_rebin_with_gtis" begin + times = collect(range(0, 99.9, step=0.1)) + + counts = round.(rand(Normal(100, 0.1), length(times))) + gti = [[0 40]; [60 100];] + + good, gti = create_gti_mask(times, gti) + + keepat!(times, good) + keepat!(counts, good) + + lc = LightCurve(time=times,counts=counts, gti=gti, skip_checks=true, dt=0.1) + + lc_rebin = rebin(lc, 1.0) + + @test (lc_rebin.time[40] - lc_rebin.time[39]) >= 1.0 + end + end +end