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

Finalize use cases for submission. #119

Merged
merged 4 commits into from
Apr 20, 2023
Merged

Finalize use cases for submission. #119

merged 4 commits into from
Apr 20, 2023

Conversation

alaindanet
Copy link
Contributor

I put my usecase in the dedicated folder following @ismael-lajaaiti usecases!

Happy to have your feedback on it!

@iago-lito iago-lito self-requested a review March 13, 2023 08:03
@iago-lito iago-lito added the documentation Improvements or additions to documentation label Mar 13, 2023
@iago-lito iago-lito linked an issue Mar 13, 2023 that may be closed by this pull request
5 tasks
extinction_threshold = 1e-6 # Set biomass threshold to consider a species extinct
# Standardize total carrying capacity for the number of producers and interspecific
# competition among producers
function stdK(fw; K = 1.0, αij = 1.0)
Copy link
Collaborator

@ismael-lajaaiti ismael-lajaaiti Mar 14, 2023

Choose a reason for hiding this comment

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

Could be renamed standardize_K?
Also the documentation of the function could be under the form of a docstring?
Also the keyword arguments could become mandatory arguments, as if I'm not wrong you always specify them below (i.e. you don't use the defaults).
Lastly, a naive question out of curiosity, why do you need this standardize the carrying capacity?

n_foodweb = 60 # number of replicates of trophic networks
foodweb_backbone = []
for i in 1:n_foodweb
for C in connectance_values
Copy link
Collaborator

@ismael-lajaaiti ismael-lajaaiti Mar 14, 2023

Choose a reason for hiding this comment

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

You use connectance_values before defining it.

Copy link
Collaborator

Choose a reason for hiding this comment

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

.. uh, does it imply that the script fails?

Copy link
Collaborator

Choose a reason for hiding this comment

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

In its current state yes, but it can be fixed just by moving this loop below the definition of connectance_values. Maybe @alaindanet you run the script chunk by chunk and not as a whole, which could explain why you didn't spot this?

foodweb_backbone = []
for i in 1:n_foodweb
for C in connectance_values
fw = FoodWeb(nichemodel, S, C = C, tol_C = 0.01, Z = Z)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is the default tol_C not good? If so, you could remove it, it would simplify the call of this functions that has already lot of arguments.
Also I think that you didn't format your code because the kwargs are not separated by a ';'.
Lastly note that FoodWeb(nichemodel, S; C = C, Z = Z) can be rewritten FoodWeb(nichemodel, S; C, Z). 🚀

Copy link
Collaborator

Choose a reason for hiding this comment

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

Yupe, that is a cool Julia feature: the awkward ; C = C, Z = Z can be shortened to ; C, Z provided it happens after the ;.

for i in 1:n_foodweb
for C in connectance_values
fw = FoodWeb(nichemodel, S, C = C, tol_C = 0.01, Z = Z)
push!(foodweb_backbone, (rep = i, C = C, fw = fw))
Copy link
Collaborator

Choose a reason for hiding this comment

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

rep could be replaced with the more explicit replicat?

Copy link
Collaborator

@iago-lito iago-lito Mar 14, 2023

Choose a reason for hiding this comment

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

I think replicat is a french word. Either rep or replicate or id then?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Sorry I meant replicatebut maybe rep is just fine, I let you choose @alaindanet.

push!(foodweb_backbone, (rep = i, C = C, fw = fw))
end
end
inter_prod_competition_range = .8:0.05:1.2 # interspecific competition among producers
Copy link
Collaborator

Choose a reason for hiding this comment

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

Write 0.8 instead of .8 for consistency. Again, I assume that the code formatter should correct this automatically.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Exceptionally I think that you could shorten the variable name, e.g. αij_range.
Also be consistent with the use of a or α, I would maybe go for α has it is the variable name use in the paper. But both are good, as long as it commented clearly.

end
end
inter_prod_competition_range = .8:0.05:1.2 # interspecific competition among producers
connectance_values = [0.05, 0.1, 0.2] # Connectance values
Copy link
Collaborator

Choose a reason for hiding this comment

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

No caps.

for a in inter_prod_competition_range
fw = foodweb_backbone[i].fw
connectance = foodweb_backbone[i].C
rep = foodweb_backbone[i].rep # foodweb replicate id
Copy link
Collaborator

Choose a reason for hiding this comment

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

If you rename rep with replicat this comment is not necessary anymore IMO.

# Generate a food-web
foodweb = FoodWeb(nichemodel, 20, C = C, tol = 0.01, Z = 100)
Threads.@threads for i in 1:length(foodweb_backbone) # parallelize computation when possible
for a in inter_prod_competition_range
Copy link
Collaborator

Choose a reason for hiding this comment

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

Replace a by αij?

nprod = length(producers(foodweb))
K_alpha_corrected = 1.0 * (1 + (a * (nprod - 1))) / nprod
env = Environment(foodweb, K = K_alpha_corrected)
prod_comp = ProducerCompetition(fw, αii = αii, αij = a)
Copy link
Collaborator

@ismael-lajaaiti ismael-lajaaiti Mar 14, 2023

Choose a reason for hiding this comment

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

Could be simplify to ProducerCompetition(fw, αii, αij) if you rename a by αij

Copy link
Collaborator

@ismael-lajaaiti ismael-lajaaiti Mar 14, 2023

Choose a reason for hiding this comment

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

Also you could rename prod_comp to producer_competition? To simplify in this same way the function call below.

# Measure species persistence, the number of species that have
# a biomass above 0 at the last timestep
pers = species_persistence(sol, threshold = 0, last = 1)
push!(output, (rep = rep, C = connectance, αij = a, persistence = pers))
Copy link
Collaborator

Choose a reason for hiding this comment

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

If you choose variable names such that they match the tuple names you could do:
push!(output, (; replicat, C, αij, persistence))
Because

a, b = 1, 2
(; a, b) # give (a = 1, b = 2) 



df = DataFrame(output)
Copy link
Collaborator

@ismael-lajaaiti ismael-lajaaiti Mar 14, 2023

Choose a reason for hiding this comment

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

You could also directly create a dataframe, e.g. replace the line output = []by df = DataFrame(:replicat = Integer[], :persistence = AbstractFloat[], :C = AbstractFloat[], :αij = AbstractFloat[]) and fill it row by row with push!, as you were already doing.

What I like about this alternative is that it is more direct. You directly store the information in a structure that you want in the end, i.e. a DataFrame and not a vector of named tuple that you will later convert to a DataFrame.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I need to point out that @alaindanet's approach has other benefits here:

  • It is (maybe) more efficient than repeated calls to push!(::DataFrame, ...) (although this depends on the internals of DataFrames.jl which I am unfamiliar with).
  • It is more robust to changes to DataFrames.jl, or to the possible drop of this particular dependency.

So, just pick your poison maybe :) I would argue that having a diversity of approaches can be useful for users to appreciate the package flexibility. On the other hand, using different approaches make the use_cases/ folder less consistent. I have no strong opinion about that.

xlabel!(p1, "Interspecific producer competition (αᵢⱼ)")
ylabel!(p1, "Species persistence")
xlabel!("Interspecific producer competition (αᵢⱼ)")
ylabel!("Species persistence")
save_plot("persistence_c_alpha")
Copy link
Collaborator

Choose a reason for hiding this comment

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

Should be removed.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Agreed, we can leave this to users.

@ismael-lajaaiti
Copy link
Collaborator

Here is my review @alaindanet. Maybe before diving into it wait for @iago-lito's round that will add a layer of comments on top of mine.

Before summarizing my comments, this is really good and I am mainly nitpicking, thank you!

Below are the main points I raised in my comments:

  • remove code that doesn't belong to the scope of the use case (that maybe need still to be refined, see this comment), i.e. saving plots and data.
  • variable names could be still be improve (see my suggestions) by being more explicit or more consistent, and if variable names correspond to function kwargs you can use the following trick: f(x; y = y, z = z)-> f(x; y, z).
  • use directly a DataFrame to store your data
  • I will work on the plotting to use library that works well with DataFrame with the hope to simplify the code section for creating the plot.

@ismael-lajaaiti
Copy link
Collaborator

ismael-lajaaiti commented Mar 20, 2023

Here are my first modifications. I integrated the above suggestions and while re-writing stuff I found that the two loops (loop 1: generation of trophic networks, loop 2: add producer competition and run simulations) could be merged. I also now use a DataFrame to store the output and I removed the unnecessary information that was saved in the output.
Also instead of using println to show simulations progress I now use the @info macro.

Coming next: use Makie and maybe AlgebraOfGraphics for the plot.

EDIT: I am adding below the to do list of the remaining things I aim to address this PR

  • simplify Hana's use case (no need to use biomass extrema)
  • update 'paradox of enrichment' use case to new style (Makie to handle the plots, DataFrames to handle the data).

@ismael-lajaaiti
Copy link
Collaborator

Coming next: use Makie and maybe AlgebraOfGraphics for the plot.

Here it is 🚀! I really like the Makie experience, the code for the plot is a bit longer than before, but I find it more intuitive and more customizable. Also I used the theme set_aog_theme!() from AlgebraOfGraphics, so we can ensure a consistent style between the different figures with one function (and this theme is quite nice).

@iago-lito you can go for your review 🤓

@ismael-lajaaiti ismael-lajaaiti force-pushed the usecase_competition3 branch 2 times, most recently from 309aaa7 to 8287404 Compare March 21, 2023 13:59
@ismael-lajaaiti
Copy link
Collaborator

I've reviewed and added Hana's use case. I've also added the part of code to plot the figure with the theme used for the other use case figures. I still need to add some comments to clarify some parameter choices, I'll update it as soon as I have answers from Hana.

@ismael-lajaaiti ismael-lajaaiti force-pushed the usecase_competition3 branch 2 times, most recently from 0598318 to 4ff4f85 Compare March 22, 2023 10:30
@ismael-lajaaiti
Copy link
Collaborator

ismael-lajaaiti commented Mar 22, 2023

I am quite happy with the result, use cases are now all consistent in terms of code style and of the figures they produce. I think I've done all I could done for the moment. I am ready for your review @iago-lito.

The next step would is to incorporate these use case in the documentation, instead of having them as raw Julia scripts in a separate folder. If I have time I may tackle this before next week.

EDIT: I can't integrate them in the documentation at least for the NTI use case, because of #118 😬. I think one trick could be to slightly modify the use case so that it runs faster (e.g. decrease community size). I think it is not an issue if the use case in the online documentation do not perfectly corresponds to the one in the paper. As we can still keep the scripts to produce the figures of the paper in the use_case/, there the code should remain untouched, or rather we should commit ourselves to maintain this code so it always produce the same output (in case we introduce breaking changes). But we can let ourselves some flexibility for the online due to the computational constraints. IMO this is the best compromise.

@ismael-lajaaiti ismael-lajaaiti changed the title Usecase producer competition Finalize use cases for submission. Mar 24, 2023
Copy link
Collaborator

@iago-lito iago-lito left a comment

Choose a reason for hiding this comment

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

Hi @ismael-lajaaiti @alaindanet! Here is my review for the use cases as they are standing now. Unfortunately, we won't be able to merge them as long as #87 has not landed, because the commits here lie downstream some commits of simulation_output_utils branch. This would not be a problem if the use cases used no features from #87 though, do they?

The comments are mostly style nitpicks, supposed to polish the use cases very much because there will be readerrrs reading them ;)

However, I have one big concern that the pattern used here to fill vectors of data frames from parallel julia threads is ⚠ unsound ⚠. You cannot just push! to a shared vector from multiple threads without protecting it with a mutex / a lock, or it causes data races. Details here and here in particular. The parallel scripts here look like they run fine (and indeed they will do most of the time because simulations are rather long so there are very few chances of simultaneous writes), but they are actually exhibiting undefined behaviour. Bad luck could lead users to severe bugs like results loss, inconsistent vector state or julia memory corruption in general.

Fortunately, the fix is rather simple: just protect your shared vector with a mutex, or (as suggested in my comments) pre-fill it with dummy values and make sure that each thread only writes into its own dedicated slot.

See you next round ;)

use_cases/diversity-vs-nti-intensity.jl Show resolved Hide resolved
use_cases/diversity-vs-nti-intensity.jl Outdated Show resolved Hide resolved
use_cases/diversity-vs-nti-intensity.jl Outdated Show resolved Hide resolved
use_cases/diversity-vs-nti-intensity.jl Outdated Show resolved Hide resolved
use_cases/diversity-vs-nti-intensity.jl Outdated Show resolved Hide resolved
use_cases/trophic-classes-vs-mortality.jl Show resolved Hide resolved
use_cases/trophic-classes-vs-mortality.jl Outdated Show resolved Hide resolved
use_cases/trophic-classes-vs-mortality.jl Show resolved Hide resolved
use_cases/trophic-classes-vs-mortality.jl Outdated Show resolved Hide resolved
use_cases/trophic-classes-vs-mortality.jl Outdated Show resolved Hide resolved
@ismael-lajaaiti
Copy link
Collaborator

Thanks again @iago-lito for your review and bringing up the data race issue and all your other suggestions. I really have the feeling to improve a lot thanks to them 🙏. I am going to rerun the use cases to ensure that they all run correctly. Once it is done, I will push my corrections. 🚀

@ismael-lajaaiti
Copy link
Collaborator

Here are my corrections. Also it seems that I forgot to previously push the paradox_enrichment.jl use case, so you haven't seen this script yet. Still I tried to adapt the suggestions you made for the other use cases in it. Let me know, what you think about it. 😉

Copy link
Collaborator

@iago-lito iago-lito left a comment

Choose a reason for hiding this comment

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

Great! Thank you @ismael-lajaaiti for polishing this up :) I have no blocking concern anymore with these use cases. Unfortunately, we cannot land them now because they depend on #87 which is almost ready but not ready yet.

I usually check that the scripts run well on a separate machine without a screen that I can access over ssh, but today I have tried to run it on my machine just because I was curious to see the actual plots.. but the plots did not show up. Is that expected? Do I need to save() or to show() them or whatever so I can have a look at them?

use_cases/paradox-enrichment.jl Outdated Show resolved Hide resolved
use_cases/tritrophic-biomass-vs-temperature.jl Outdated Show resolved Hide resolved
use_cases/paradox-enrichment.jl Show resolved Hide resolved
use_cases/tritrophic-biomass-vs-temperature.jl Outdated Show resolved Hide resolved
@ismael-lajaaiti
Copy link
Collaborator

Thank you :) No worries, I'll wait #87 to be merged. I incorporated your latest suggestions.

To see the plot see my comment above.

If not working from an IDE, CairoMakie won't display plots contrary to the Plots package. You need to save the figure to see it. Try:

save("path/to/save/the/figure", fig; resolution = (450, 350), px_per_unit = 3)

The kwargs are optional, but ensure you that the figure will look very nice. 😆

@iago-lito
Copy link
Collaborator

#87 Has landed! I'll rebase this on top of that and land it asap :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation Improvements or additions to documentation
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add use cases to the documentation
4 participants