Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add fourier methods and tests #2

Merged
merged 7 commits into from
Jul 25, 2022
Merged

Add fourier methods and tests #2

merged 7 commits into from
Jul 25, 2022

Conversation

Aman-Pandey-afk
Copy link
Collaborator

Added the functions from fourier.py, along with required functions from gti.py and util.py (Remaining will be implemented according to use). Tests are thoroughly covered. DataFrames and HDF5 support for large data also implemented.

@codecov
Copy link

codecov bot commented Jun 23, 2022

Codecov Report

Merging #2 (53d8ffe) into main (70cafd1) will increase coverage by 94.55%.
The diff coverage is 94.55%.

@@            Coverage Diff            @@
##           main       #2       +/-   ##
=========================================
+ Coverage      0   94.55%   +94.55%     
=========================================
  Files         0        3        +3     
  Lines         0      386      +386     
=========================================
+ Hits          0      365      +365     
- Misses        0       21       +21     
Impacted Files Coverage Δ
src/gti.jl 81.63% <81.63%> (ø)
src/fourier.jl 96.39% <96.39%> (ø)
src/utils.jl 100.00% <100.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 70cafd1...53d8ffe. Read the comment docs.

Copy link
Collaborator

@giordano giordano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The main problem is that the code feels like written in Python, not in Julia. It's very unidiomatic. Function signatures are completely non-annotated, lots of keyword arguments with nothing as default value, that's really not a thing in Julia.

Additionally, the code is completely undocumented, I presume that'd come later?

Also, I haven't actually run the code, but I'm certain performance isn't optimal, there are plenty of type-instabilities, did you use @code_warntype?

Project.toml Show resolved Hide resolved
src/utils.jl Outdated
Comment on lines 1 to 6
function sum_if_not_none_or_initialize(A,B)
if isnothing(A)
return deepcopy(copy((B)))
end
return A + B
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you please explain with words (maybe a comment) what this function is supposed to do? It looks very inefficient and I'm having a hard time thinking about a case where this could be useful.

Also, if anything the julian way to define it would be:

sum_if_not_none_or_initialize(::Nothing,B) = B
sum_if_not_none_or_initialize(A,B) = A + B

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I have used copy() in vain, but the function is executed in loop, so I will need a deepcopy of B when A is nothing (otherwise B points to A).

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so I will need a deepcopy of B when A is nothing (otherwise B points to A).

What does this mean?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Basically it is a recursive copy of objects, but here it is just to ensure that we don't pass a reference of A to B. Shallow copy could be used, but I though there must be some motive in using deepcopy in python too (like maybe variance could be an array in some case so shallow copy doesn't guarantee that A and B are independent).

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Specifically I'm missing what

(otherwise B points to A)

means. Maybe a concrete reproducible example would help.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The specific use-case is here: https://github.com/StingraySoftware/stingray/blob/1bd574a284dfe60177e3313d50cc5700f4d1ec42/stingray/fourier.py#L1005-L1007
As you can see if I just use sum_if_not_none_or_initialize(::Nothing,B) = B , it will point to the reference of B. Ex:
image

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As a general remark, never show screenshots, whoever wants to help you can't copy from that, moreover if someone wants to come back to this discussion whatever is written in the screenshot doesn't show up in search. They aren't really much useful.

That said, I'm not sure that's a problem in practice here.

Also, I think this discussion is a bit out of the point now: I believe you should first understand how to avoid the overuse of nothing and write more idiomatic Julia code here, at which point a new solution may arise.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @giordano, @Aman-Pandey-afk :
This is the use case: in another function, I am calculating (e.g.) an average of many periodograms produced over a number of short time series. In each loop, I produce a periodogram that has to be summed to the total periodogram. The summed periodogram is initialized as None. In the first loop, this function will find that the total periodogram is None, and initialize it to a deep copy of the first periodogram. From that point on, all subsequent arrays will be summed to the total periodogram in place.
The reason for the deep copy is to avoid side effects on the first periodogram, in case it is useful later.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, but the deepcopy is a red herring here, my main point is that we should move away from overusing nothing. There are definitely good use cases for nothing in Julia, but I feel like all the uses here come from the fact it was used (for good reasons!) in Python, not that it's necessary in Julia.

src/Stingray.jl Outdated
export time_intervals_from_gtis

include("utils.jl")
export sum_if_not_none_or_initialize
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are so many exports here. I think it's good to export what's useful to end users, otherwise you just clutter the scope needlessly. sum_if_not_none_or_initialize looks a rather internal function (and see my comment about the need for this function in the first place).

src/fourier.jl Outdated
@@ -0,0 +1,755 @@
function positive_fft_bins(n_bin; include_zero = false)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please, put type annotation to the signatures of the functions. It's really hard to guess what the code is supposed to do if I have no idea what the input should be. I appreciate the original code probably doesn't have type annotations because it was written in Python, but this is a different language, and the goal is to write something idiomatic in Julia.

include_zero would be a Bool, so use include_zero::Bool, n_bin is (my guess) an Integer?

src/fourier.jl Outdated
minbin = 1
end

if n_bin%2 == 0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please use the semantically relevant functions. Julia has iseven

Suggested change
if n_bin%2 == 0
if iseven(n_bin)

src/fourier.jl Show resolved Hide resolved
src/fourier.jl Outdated
n_bin = nothing

for flux in local_show_progress(flux_iterable)
if isnothing(flux) || all(flux .== 0)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if isnothing(flux) || all(flux .== 0)
if isnothing(flux) || all(iszero, flux)

has two advantages: you don't allocate the array flux .== 0, and also all is able to return early at the first non-zero element of flux.

src/fourier.jl Outdated
sum_of_photons1 = sum_of_photons2 = 0

for (flux1, flux2) in zip(flux_iterable1, flux_iterable2)
if isnothing(flux1) || isnothing(flux2) || all(flux1 .== 0) || all(flux2 .== 0)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if isnothing(flux1) || isnothing(flux2) || all(flux1 .== 0) || all(flux2 .== 0)
if isnothing(flux1) || isnothing(flux2) || all(iszero, flux1) || all(iszero, flux2)

src/fourier.jl Outdated

for (flux1, flux2) in local_show_progress(zip(flux_iterable1,
flux_iterable2))
if isnothing(flux1) || isnothing(flux2) || all(flux1 == 0) || all(flux2 == 0)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if isnothing(flux1) || isnothing(flux2) || all(flux1 == 0) || all(flux2 == 0)
if isnothing(flux1) || isnothing(flux2) || all(iszero, flux1) || all(iszero, flux2)

src/fourier.jl Outdated
Comment on lines 722 to 723
is_events = all([isnothing(val) for val in (fluxes1, fluxes2, errors1,
errors2)])
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
is_events = all([isnothing(val) for val in (fluxes1, fluxes2, errors1,
errors2)])
is_events = all(isnothing, (fluxes1, fluxes2, errors1,errors2))

@giordano
Copy link
Collaborator

BTW, another thing I just realised is that in the code there are a lot of all(arr .== 0) (which I suggested to replace with all(iszero, arr) because it's faster and it avoids allocation of temporary arrays), but with floating point numbers it isn't generally a good idea to check exact equality with 0, you'd use for example isapprox (with a non-zero absolute tolerance). Is exact equality with zero really what you want to check?

@Aman-Pandey-afk
Copy link
Collaborator Author

I have used ≈ wherever required, the places you have suggested to use iszero requires the array sum to be not equal to zero (leading to unnecessary calculations) thus, exact equality is needed to be checked. I will use iszero instead of flux1 .== 0.

src/gti.jl Outdated
Comment on lines 23 to 26
st = st[1:end-1]
end
if st[end] + nbin > stopbin
return st[1:end-1]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is equivalent to

Suggested change
st = st[1:end-1]
end
if st[end] + nbin > stopbin
return st[1:end-1]
pop!(st)
end
if st[end] + nbin > stopbin
pop!(st)

right?

Also, I have already mentioned that slicing creates copies, so in st = st[1:end-1] you are creating a copy of st[1:end-1] before reassigning it to the st identifier. pop! simply removes the last element of the collection, without assuming anything about indexing of the array (Julia is 1-based by default, but arbitrary indices are allowed).

src/gti.jl Outdated
Comment on lines 43 to 44
gti_low = (@view gtis[:, 1]) .+ (dt / 2 - epsilon_times_dt)
gti_up = (@view gtis[:, 2]) .- (dt / 2 - epsilon_times_dt)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
gti_low = (@view gtis[:, 1]) .+ (dt / 2 - epsilon_times_dt)
gti_up = (@view gtis[:, 2]) .- (dt / 2 - epsilon_times_dt)
gti_low = @view(gtis[:, 1]) .+ (dt ./ 2 .- epsilon_times_dt)
gti_up = @view(gtis[:, 2]) .- (dt ./ 2 .- epsilon_times_dt)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both dt and epsilon_times_dt are float, so I guess no need to dot broadcast

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as #2 (comment)

src/gti.jl Outdated

function calculate_segment_bin_start(startbin::Integer, stopbin::Integer,
nbin::Integer; fraction_step::Real=1)
st = floor.(range(startbin, stopbin, step=Int64(nbin * fraction_step)))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless you have a specific reason for using 64-bit integers, it's usually more natural to use Int, which is the integer type with the same number of bits as the current architecture (32-bit on 32-bit systems, 64-bit on 64-bit integers). Int is also the type of integer literals in Julia

Suggested change
st = floor.(range(startbin, stopbin, step=Int64(nbin * fraction_step)))
st = floor.(range(startbin, stopbin, step=Int(nbin * fraction_step)))

src/gti.jl Outdated
end

epsilon_times_dt = epsilon * dt
nbin = Int64(round(segment_size / dt))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
nbin = Int64(round(segment_size / dt))
nbin = round(Int, segment_size / dt)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But actually, integer division (already suggested in #2 (comment)) is much idiomatic

Suggested change
nbin = Int64(round(segment_size / dt))
nbin = segment_size ÷ dt

src/gti.jl Outdated
epsilon_times_dt = epsilon * dt
nbin = Int64(round(segment_size / dt))

spectrum_start_bins = Int64[]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
spectrum_start_bins = Int64[]
spectrum_start_bins = Int[]

src/fourier.jl Outdated
Comment on lines 102 to 108
raw_coherence(cross_power::T1, power1::T2, power2::T2, power1_noise::Real,
power2_noise::Real, n_ave::Integer;
intrinsic_coherence::Real=1.0) where
{T1<:Union{AbstractVector{<:Number},Number},T2<:Union{AbstractVector{<:Real},Real}} =
_raw_coherence_single.(cross_power, power1, power2,
power1_noise, power2_noise, n_ave;
intrinsic_coherence)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, it sounds like this should be defined on scalars, and use broadcasting on the call site

src/fourier.jl Outdated
Comment on lines 131 to 136
estimate_intrinsic_coherence(cross_power:: AbstractVector{<:Complex},
power1:: AbstractVector{<:Real}, power2:: AbstractVector{<:Real},
power1_noise::Real, power2_noise::Real, n_ave::Integer)=
_estimate_intrinsic_coherence_single.(
cross_power, power1, power2,
power1_noise, power2_noise, n_ave)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this estimate_intrinsic_coherence only calling _estimate_intrinsic_coherence_single broadcasted?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I think I should define a single function only and call it broadcasted everywhere, the thing was API is meant to be used with vector inputs and not just broadcasted. We can mention in the docstrings though.

src/fourier.jl Outdated
Comment on lines 172 to 177
function cross_to_covariance(cross_power:: AbstractVector{<:Complex},
ref_power:: AbstractVector{<:Real},
ref_power_noise::Real,
delta_nu::T) where {T<:Union{AbstractVector{<:Real},Real}}
return cross_power .* sqrt.(delta_nu ./ (ref_power .- ref_power_noise))
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Once again, the pattern in Julia is to define a function f(x::Real) = ... only on scalars, and then call it f.(V) where V::AbstractVector{<:Real}. There is no need to define the function on scalars and vectors, it's a useless and confusing duplicate.

src/fourier.jl Outdated
if isnothing(segment_size)
segment_size = max(gti) - min(gti)
end
n_bin = Int64(round(segment_size / dt))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
n_bin = Int64(round(segment_size / dt))
n_bin = segment_size ÷ dt

Copy link
Collaborator Author

@Aman-Pandey-afk Aman-Pandey-afk Jul 9, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well Integer Division always rounds to smaller integer right? (EX: round(0.7)=1 is the result needed)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, then the first suggestion in #2 (comment) stands.

src/fourier.jl Outdated
if isnothing(segment_size)
segment_size = max(gti) - min(gti)
end
n_bin = Int64(round(segment_size / dt))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
n_bin = Int64(round(segment_size / dt))
n_bin = segment_size ÷ dt

src/utils.jl Outdated
return A + B
end

sqrt(x::Real) = x < 0.0 ? NaN : Base.sqrt(x)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I previously wrote a comment (in the wrong page, so I deleted it) that this was overriding a Base method. This wasn't true, you're instead shadowing it (less bad, but still a bit). What's the use case for redefining sqrt in this way? If the argument is negative, an error is usually the right thing to do, instead of poising everything with NaN, which is hard to debug.

@Aman-Pandey-afk
Copy link
Collaborator Author

Methods like Keepat and broadcast are surely great in reducing allocations, as I can observe them with tests. I have added all the suggestions plus some changes influenced by this discussion. But, somehow the time taken is more than that before these 3 commits. (As per the ss on slack, the allocations were around 2.3 million, now 2.5 with slight increase in GC%), although opposite should happen. I think a test on other OS will provide good measure. (I'm sure this code is more efficient).

return spectrum_start_bins, spectrum_start_bins.+nbin
end

@resumable function generate_indices_of_segment_boundaries_unbinned(times::AbstractVector{<:Real},
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is showing not to be covered by tests, but it is being called on line 177 in the fourier.jl file which is showing to be covered. I have tried logging some statements in this function during test runs and it is actually being executed (otherwise the tests won't pass also).

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

src/fourier.jl Outdated
intrinsic_coherence::Real=1.0)

if n_ave > 500
return 0.0 * power1
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
return 0.0 * power1
return float(zero(power1))

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually I should just return 0.0 now as the function is meant to be used with dot broadcast for vectors.

src/fourier.jl Show resolved Hide resolved
src/fourier.jl Outdated
Comment on lines 217 to 218
# astype here serves to avoid integer rounding issues in Windows,
# where long is a 32-bit integer.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was just a np.astype comment. I had kept them to understand the codes. Will change/remove them.

src/fourier.jl Outdated
Comment on lines 216 to 221
event_times = times[idx0:idx1-1]
# astype here serves to avoid integer rounding issues in Windows,
# where long is a 32-bit integer.
cts = fit(Histogram,Float64.(event_times .- s);nbins=n_bin).weights
else
cts = Float64.(fluxes[idx0+1:idx1])
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
event_times = times[idx0:idx1-1]
# astype here serves to avoid integer rounding issues in Windows,
# where long is a 32-bit integer.
cts = fit(Histogram,Float64.(event_times .- s);nbins=n_bin).weights
else
cts = Float64.(fluxes[idx0+1:idx1])
event_times = @view times[idx0:idx1-1]
# astype here serves to avoid integer rounding issues in Windows,
# where long is a 32-bit integer.
cts = fit(Histogram, float.(event_times .- s); nbins=n_bin).weights
else
cts = float.(@view fluxes[idx0+1:idx1])

src/fourier.jl Outdated
else
cts = Float64.(fluxes[idx0+1:idx1])
if !isnothing(errors)
cts = cts, errors[idx0+1:idx1]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this correct? cts is changing in a hot loop from a vector to a tuple of vectors?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, this can be both tuple or a simple vector which is element of the flux_iterable. If this is costly we may separate errors into another array and isTuple check then becomes an isnothing check.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed

src/fourier.jl Outdated
end

# No need for the negative frequencies
unnorm_power = keepat!(unnorm_power,fgt0)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Am I missing something or fgt0 here can be also nothing?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cross is initialized with nothing which means fgt0 = positive_fft_bins(n_bin) will always be executed first time.

src/fourier.jl Outdated

# Calculate the unnormalized cross spectrum
unnorm_power = ft1 .* conj.(ft2)
unnorm_pd1 = unnorm_pd2 = 0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Below unnorm_pd1 and unnorm_pd2 can become vectors, no?

Copy link
Collaborator Author

@Aman-Pandey-afk Aman-Pandey-afk Jul 13, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I think it is a initialization problem, these are only used when return_auxil is true and as variables have a scope outside the if block too in Julia, I will remove these initialization.

src/fourier.jl Outdated
Comment on lines 538 to 541
unnorm_power = keepat!(unnorm_power,fgt0)
if return_auxil
unnorm_pd1 = keepat!(unnorm_pd1,fgt0)
unnorm_pd2 = keepat!(unnorm_pd2,fgt0)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

keepat! modifies the vector in-place, no need to reassign the output to itself.

Suggested change
unnorm_power = keepat!(unnorm_power,fgt0)
if return_auxil
unnorm_pd1 = keepat!(unnorm_pd1,fgt0)
unnorm_pd2 = keepat!(unnorm_pd2,fgt0)
keepat!(unnorm_power,fgt0)
if return_auxil
keepat!(unnorm_pd1,fgt0)
keepat!(unnorm_pd2,fgt0)

src/fourier.jl Outdated
Comment on lines 619 to 620
unnorm_pds1 /= n_ave
unnorm_pds2 /= n_ave
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With

Suggested change
unnorm_pds1 /= n_ave
unnorm_pds2 /= n_ave
unnorm_pds1 ./= n_ave
unnorm_pds2 ./= n_ave

you modify unnorm_pds1 and unnorm_pds2 in-place without allocating temporary arrays and assigning them to these variables.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So any place I work with operation on both scalar and vectors, I should use dot broadcast right? (Shouldn't be this then the natural behaviour)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Uhm, no, you can't do x .= ... with x being a scalar. And no, in-place assignment is not going to be the "natural" (if you mean default) behaviour, because that's potentially disruptive.

Copy link
Collaborator

@giordano giordano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think main issues have been addressed. However I believe the code could be much more idiomatic (I gave a concrete example for poisson_level), but that's for another time, although I'd like to see more code like that in the future.

Comment on lines +9 to +21
function poisson_level(norm::String; meanrate = nothing, n_ph = nothing, backrate::Real = 0.0)
if norm == "abs"
return 2.0 * meanrate
elseif norm == "frac"
return 2.0 / (meanrate - backrate)^2 * meanrate
elseif norm == "leahy"
return 2.0
elseif norm == "none"
return float(n_ph)
else
throw(ArgumentError("Unknown value for norm: $norm"))
end
end
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a better way to represent this is to define a "normalisation" type. Something along these lines (also with better names than Abs and Frac, I'm just copying the short names used here):

abstract type Normalization end

struct Abs{T} <: Normalization
    meanrate::T
end
poisson_level(n::Abs) = 2 * n.meanrate

struct Frac{T} <: Normalization
    meanrate::T
    backrate::T
end
poisson_level(n::Frac) = 2 / (n.meanrate - n.backrate) ^ 2 * n.meanrate

struct Leahy <: Normalization end
poisson_level(::Leahy) = 2

struct Unnormalized{T} <: Normalization
    n_ph::T
end
poisson_level(n::Unnormalized) = n_ph

It's idiomatic, it makes natural use of multiple dispatch (instead of dispatching on strings....), it's optimisable (the current implementation of poisson_level can't possibly be optimised by the compiler, since all arguments are basically Any), no need to abuse nothing, etc.

Similar things can be done elsewhere throughout the code, to make it more idiomatic and actually use Julia features, instead of writing Python code with Julia syntax. But let's not focus on this now, this is mainly an idea for the future, to give a concrete suggestion for how to make the code more idiomatic.

@matteobachetti
Copy link
Member

Thanks @Aman-Pandey-afk and @giordano, I learned a lot from this PR and its review.

@matteobachetti matteobachetti merged commit 0a7b39f into main Jul 25, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants