final_touches

Author: kashish2210

src/fourier.jl

@@ -224,24 +224,15 @@ end
     end
 end
 
-"""
-    avg_pds_from_iterable(flux_iterable, dt::Real; kwargs...)
-
-Compute averaged power density spectrum from an iterable of flux segments.
-
-## Parameters
-- `flux_iterable`: Iterator of flux segments
-- `dt::Real`: Time resolution
-- `norm::String="frac"`: Normalization type
-- `use_common_mean::Bool=true`: Use common mean across segments
-- `silent::Bool=false`: Suppress progress output
-"""
-function avg_pds_from_iterable(flux_iterable, dt::Real; 
-                             norm::String="frac", 
-                             use_common_mean::Bool=true,
-                             silent::Bool=false)
-    
-    # Initialize accumulators
+function avg_pds_from_iterable(flux_iterable, dt::Real; norm::String="frac", 
+                               use_common_mean::Bool=true,
+                               silent::Bool=false)
+    local_show_progress = show_progress
+    if silent
+        local_show_progress = identity
+    end
+
+    # Initialize stuff
     cross = unnorm_cross = nothing
     n_ave = 0
     freq = nothing
@@ -251,9 +242,8 @@ function avg_pds_from_iterable(flux_iterable, dt::Real;
     fgt0 = nothing
     n_bin = nothing
 
-    for flux in flux_iterable
-        # Validate segment
-        if isnothing(flux) || all(iszero, flux)
+    for flux in local_show_progress(flux_iterable)
+        if isnothing(flux) || all(iszero,flux)
             continue
         end
 
@@ -262,7 +252,7 @@ function avg_pds_from_iterable(flux_iterable, dt::Real;
         variance = nothing
         if flux isa Tuple
             flux, err = flux
-            variance = Statistics.mean(err)^2
+            variance = Statistics.mean(err) ^ 2
         end
 
         # Calculate the FFT
@@ -279,43 +269,40 @@ function avg_pds_from_iterable(flux_iterable, dt::Real;
         sum_of_photons += n_ph
 
         if !isnothing(variance)
-            common_variance = isnothing(common_variance) ? variance : common_variance + variance
+            common_variance = sum_if_not_none_or_initialize(common_variance, variance)
         end
 
         # In the first loop, define the frequency and the freq. interval > 0
-        # Initialize frequency arrays
         if isnothing(cross)
             fgt0 = positive_fft_bins(n_bin)
-            freq = FFTW.fftfreq(n_bin, dt)[fgt0]
+            freq = fftfreq(n_bin, dt)[fgt0]
         end
 
-        keepat!(unnorm_power, fgt0)
+        # No need for the negative frequencies
+        keepat!(unnorm_power,fgt0)
 
         # If the user wants to normalize using the mean of the total
         # lightcurve, normalize it here
         cs_seg = unnorm_power
-        if !use_common_mean
-            mean_flux = n_ph / n_bin
+        if !(use_common_mean)
+            mean = n_ph / n_bin
+
             cs_seg = normalize_periodograms(
-                unnorm_power,
-                dt,
-                n_bin;
-                mean_flux = mean_flux,
-                n_ph = n_ph,
-                norm = norm,
-                power_type = "real"
+                unnorm_power, dt, n_bin; mean_flux = mean, n_ph=n_ph,
+                norm=norm, variance=variance,
             )
         end
 
         # Accumulate the total sum cross spectrum
         cross = sum_if_not_none_or_initialize(cross, cs_seg)
-        unnorm_cross = sum_if_not_none_or_initialize(unnorm_cross, unnorm_power)
+        unnorm_cross = sum_if_not_none_or_initialize(unnorm_cross,
+                                                     unnorm_power)
 
         n_ave += 1
     end
 
-    # Validate results
-    if isnothing(cross) || n_ave == 0
+    # If there were no good intervals, return nothing
+    if isnothing(cross)
         return nothing
     end
 
@@ -337,28 +324,19 @@ function avg_pds_from_iterable(flux_iterable, dt::Real;
     # Final normalization (If not done already!)
     if use_common_mean
         cross = normalize_periodograms(
-            unnorm_cross,
-            dt,
-            n_bin;
-            mean_flux = common_mean,
-            n_ph = n_ph,
-            norm = norm,
-            power_type = "real"
+            unnorm_cross, dt, n_bin; mean_flux=common_mean, n_ph=n_ph,
+            norm=norm, variance=common_variance
         )
     end
 
-    # Calculate power errors
-    power_errors = cross ./ sqrt(n_ave)
+    results = DataFrame()
+    results[!,"freq"] = freq
+    results[!,"power"] = cross
+    results[!,"unnorm_power"] = unnorm_cross
+    results = (; results , n = n_bin, m = n_ave, dt, norm, df = 1 / (dt * n_bin), nphots = n_ph, setment_size = dt * n_bin,  mean = common_mean, variance = common_variance)
+    results = DataFrame(results.results)
+    return results
 
-    # Return AveragedPowerspectrum
-    return AveragedPowerspectrum(
-        freq[2:end],
-        cross[2:end],
-        power_errors[2:end],
-        n_ave,
-        n_bin,
-        dt
-    )
 end
 
 function avg_cs_from_iterables_quick(flux_iterable1 ,flux_iterable2,
@@ -698,24 +676,33 @@ function avg_cs_from_iterables(
     return results
 
 end
-function avg_pds_from_events(times::AbstractVector{<:Real}, gti::AbstractMatrix{<:Real}, 
-                           segment_size::Real, dt::Real; norm::String="frac",
-                           use_common_mean::Bool=true, silent::Bool=false, 
-                           fluxes=nothing, errors=nothing)
+
+function avg_pds_from_events(times:: AbstractVector{<:Real}, gti::AbstractMatrix{<:Real}, 
+                             segment_size::Real, dt::Real; norm::String="frac",
+                             use_common_mean::Bool=true, silent::Bool=false, 
+                             fluxes=nothing, errors=nothing)
     if isnothing(segment_size)
         segment_size = max(gti) - min(gti)
     end
-    n_bin = round(Int, segment_size / dt)
+    n_bin = round(Int,segment_size / dt)
     dt = segment_size / n_bin
 
     flux_iterable = get_flux_iterable_from_segments(times, gti, segment_size;
-                                                  n_bin, fluxes=fluxes,
-                                                  errors=errors)
+                                                    n_bin, fluxes=fluxes,
+                                                    errors=errors)
+    cross = avg_pds_from_iterable(flux_iterable, dt, norm=norm,
+                                  use_common_mean=use_common_mean,
+                                  silent=silent)
+
+    if !isnothing(cross)
+        cross = (; cross, gti=gti)
+        cross = DataFrame(cross.cross)
+    end
+
+    return cross
 
-    return avg_pds_from_iterable(flux_iterable, dt, norm=norm,
-                               use_common_mean=use_common_mean,
-                               silent=silent)
 end
+
 function avg_cs_from_events(times1:: AbstractVector{<:Real}, times2:: AbstractVector{<:Real}, 
                             gti::AbstractMatrix{<:Real}, segment_size::Real, dt::Real; 
                             norm::String="frac", use_common_mean::Bool=true, 
@@ -771,4 +758,4 @@ function avg_cs_from_events(times1:: AbstractVector{<:Real}, times2:: AbstractVe
     end
 
     return results
-end
+end

test/test_fourier.jl

@@ -1,22 +1,3 @@
-function compare_tablespds(ps1::AveragedPowerspectrum, ps2::AveragedPowerspectrum; rtol=0.001, discard=[])
-    test_result = true
-
-    if !isapprox(ps1.freqs, ps2.freqs, rtol=rtol) ||
-       !isapprox(ps1.power, ps2.power, rtol=rtol) ||
-       !isapprox(ps1.power_errors, ps2.power_errors, rtol=rtol) ||
-       ps1.m != ps2.m ||
-       ps1.n != ps2.n ||
-       !isapprox(ps1.dt, ps2.dt, rtol=rtol)
-        test_result = false
-    end
-    
-    @test test_result
-end
-# ----------------------------------------------------------------
-# here there are two function's compare_tables(){"with dataframe"} and compare_tablespds(){"without dataframe"} which are meant to be merged
-# in future while converting cross spectrum with Struct{}
-# ----------------------------------------------------------------
-
 function compare_tables(table1, table2; rtol=0.001, discard = [])
 
     s_discard = Symbol.(discard)
@@ -64,6 +45,7 @@ function compare_tables(table1, table2; rtol=0.001, discard = [])
     end
     @test test_result
 end
+
 @testset "positive_fft_bins" begin
     freq = fftfreq(11)
     goodbins = positive_fft_bins(11)
@@ -181,7 +163,7 @@ end
         @test isnothing(out_ev) 
     end
 
-   @testset "test_avg_pds_use_common_mean_similar_stats" for norm in ["frac", "abs", "none", "leahy"]
+    @testset "test_avg_pds_use_common_mean_similar_stats" for norm in ["frac", "abs", "none", "leahy"]
         out_comm = avg_pds_from_events(
             times,
             gti,
@@ -191,7 +173,7 @@ end
             use_common_mean=true,
             silent=true,
             fluxes=nothing,
-        )
+        ).power
         out = avg_pds_from_events(
             times,
             gti,
@@ -201,8 +183,8 @@ end
             use_common_mean=false,
             silent=true,
             fluxes=nothing,
-        )
-        @test std(out_comm.power) ≈ std(out.power) rtol=0.1
+        ).power
+        @test std(out_comm)≈std(out) rtol=0.1
     end
 
     @testset "test_avg_cs_use_common_mean_similar_stats" for 
@@ -256,7 +238,7 @@ end
             silent=true,
             fluxes=counts,
         )
-        compare_tablespds(out_ev, out_ct)
+        compare_tables(out_ev, out_ct)
     end
 
     @testset "test_avg_pds_cts_and_err_and_events_are_equal" for use_common_mean in [true,false], norm in ["frac", "abs", "none", "leahy"]       
@@ -281,7 +263,12 @@ end
             fluxes=counts,
             errors=errs,
         )
-        compare_tablespds(out_ev, out_ct)
+        # The variance is not _supposed_ to be equal, when we specify errors
+        if use_common_mean
+            compare_tables(out_ev, out_ct, rtol=0.01, discard=["variance"])
+        else
+            compare_tables(out_ev, out_ct, rtol=0.1, discard=["variance"])
+        end
     end
 
     @testset "test_avg_cs_cts_and_events_are_equal" for use_common_mean in [true,false], norm in ["frac", "abs", "none", "leahy"]
@@ -435,4 +422,4 @@ end
         pdsnorm = normalize_periodograms(pds, dt, N; mean_flux=mean, n_ph=nph, norm="none")
         @test Statistics.mean(pdsnorm)≈Statistics.mean(pds) rtol=0.01
     end
-end
+end