[MRG] updated visualisation methods and unit handling #264

cjayb · 2021-01-27T14:25:30Z

NB 81ff283 removes testing of MPIBackend due to CircleCI shitting kittens (see #215)

This PR updates methods for visualisation of Dipole

improve time series smoothing using convolution (and add Savitzky-Golay and better docstrings)
power spectrum estimation and plotting (simple periodogram, no Welch)
handling of units in plots (alternative to a scale-operation) (closes Units missing in dipole plots #243 )
update alpha-example to demonstrate
implement test for copy method
add new test_viz.py for testing visualisations
closes plot_tfr_morlet does not work for list of Dipoles #267 by allowing passing list of dipoles to plot_tfr_morlet (and plot_psd)

Closes #265 : baseline normalisation left in to have parity with HNN GUI, this PR addresses units-issues

To follow up on a discussion in #249 (@jasmainak @rythorpe): I think this example might benefit from smoothing the waveforms? I tried with smooth(1000), which looks quite illustrative for both cases:

We could leave it for 0.2 release, though. There's a lot we need to figure out regarding visualisation.

examples/plot_simulate_alpha.py

ntolley · 2021-01-27T14:46:29Z

Thanks @cjayb, I'll go ahead and merge once the CI's go green.

I actually really like the smoothing for illustrative purposes. Might require a quick statement in the docstring like "MEG/EEG signals undergo low-pass filtering when passing through the skull. We can apply a smoothing to the simulated current dipole to better reflect experimentally recorded signals."

Totally fine to wait to merge if you want to go ahead and add the smoothing.

jasmainak · 2021-01-27T15:02:58Z

MEG/EEG signals undergo low-pass filtering when passing through the skull.

yikes, Matti will disown me if he sees this sentence in hnn-core :) I believe the low-pass hypothesis has never been proven experimentally ... but I am happy to be proven wrong

jasmainak · 2021-01-27T15:05:06Z

I also saw this sentence in the example:

The dpl[trial_idx].scale() call relates to the amount of cortical tissue necessary to observe the electric current dipole outside the head with M/EEG.

but no call to "scale". Remove the sentence?

jasmainak · 2021-01-27T15:15:59Z

I tried with smooth(1000), which looks quite illustrative for both cases

I wonder if the smooth function should take input in seconds instead of samples

ntolley · 2021-01-27T15:20:21Z

yikes, Matti will disown me if he sees this sentence in hnn-core :) I believe the low-pass hypothesis has never been proven experimentally ... but I am happy to be proven wrong

Ah good to know, quick reading suggests that it may reflect how synchronous high frequency activity occurs with smaller volumes of cortical tissue.

Point stills stands, if we want to include smoothing do you have a suggested explanation that won't offend too many people? Perhaps just "Several factors contribute to the appearance of low-pass filtered signal when recording neural signals outside the head. We can apply a smoothing to the simulated current dipole to better reflect experimentally recorded signals."

jasmainak · 2021-01-27T15:26:47Z

So -- low-pass filtering the data will create a corner frequency whereas I have mostly seen the spectrum decay (what's usually called a 1/f but that's a bit of a misnomer too). I am also not exactly sure what the smoothing with Hamming window does. It low-pass filters but what's the corner frequency? Okay I have said enough ... floor is yours @cjayb to pitch opinions :)

codecov-io · 2021-01-27T15:57:33Z

Codecov Report

Merging #264 (db0555c) into master (a4a8ac0) will decrease coverage by 0.09%.
The diff coverage is 37.35%.

@@            Coverage Diff             @@
##           master     #264      +/-   ##
==========================================
- Coverage   52.92%   52.83%   -0.10%     
==========================================
  Files          23       24       +1     
  Lines        2749     2843      +94     
==========================================
+ Hits         1455     1502      +47     
- Misses       1294     1341      +47

Impacted Files	Coverage Δ
hnn_core/parallel_backends.py	`15.25% <0.00%> (-0.27%)`	⬇️
hnn_core/tests/test_network.py	`100.00% <ø> (ø)`
hnn_core/viz.py	`4.68% <2.89%> (-0.59%)`	⬇️
hnn_core/dipole.py	`24.13% <31.11%> (-0.34%)`	⬇️
hnn_core/tests/test_dipole.py	`100.00% <100.00%> (ø)`
hnn_core/tests/test_viz.py	`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 a4a8ac0...2bc577f. Read the comment docs.

cjayb · 2021-01-27T19:35:06Z

I think Matti and @jasmainak are vindicated: http://dx.doi.org/10.1523/ENEURO.0291-16.2016 :)

Any real evoked or even rhythmic response is riding on top of that 1/f-like sea of, you know, that ninety-something% of "brain activity" we rarely focus on. So in essence we have hundreds of thousands of 'other' HNN networks active, creating a stochastic smoothing. Furthermore, one of our 200 pyramidal neuron-strong networks is probably only on the order of a percent of the total number of microcolumns creating the macroscopic signals we measure. And each of these is also doing other computations, involving other dynamics that interact and interfere.

I see smoothing here as a way to 'simulate' the combined effects of 1) non-zero background and 2) hundreds or thousands of more n_trials. I don't think we can 'justify it' with reference to any biophysical properties of neurons or tissue. The former surely generate HH-dynamics at a very high rate, and the latter is essentially a linear conductor with a flat impedance spectrum.

I'll take a look at the smoothing code (agree it's currently a little ad hoc), and see if I can write some prose about the background (without drowning the example in too many details). @ntolley had a delightfully 'Swiss' proposal above, maybe best to dodge the issue like that :)

jasmainak · 2021-01-27T20:38:19Z

yes, agree with @cjayb 's explanation. Another point I can't resist making is that intracranial EEG signals look pretty smooth as well ... at least what I've seen so far from recent explorations. So the destructive/constructive interference between multiple HH networks seems to make more sense to me. I'm okay with either @ntolley 's swiss explanation or a more detailed explanation that is less wrong than the first explanation we had :)

jasmainak · 2021-01-28T06:04:53Z

I like the returning of self :)

what's the benefit of savgol_filter over smooth?

EDIT: okay I see you can specify h_freq :)

jasmainak · 2021-01-28T06:09:52Z

OK to use periodogram, or should it be Welch with a winlen-parameter?

I guess it depends what you're trying to find. Welch will dampen any local effects but periodogram may show even artifacts

jasmainak · 2021-01-28T06:10:37Z

should the method be called plot_spectrum, for simplicity?

sure. Okay I am deferring the review of the rest of the PR to @ntolley . Need to go underground for grant writing again :)

jasmainak · 2021-01-28T06:11:45Z

hnn_core/dipole.py

@@ -318,3 +329,38 @@ def write(self, fname):
                   self.data['L5']]].T
        np.savetxt(fname, X, fmt=['%3.3f', '%5.4f', '%5.4f', '%5.4f'],
                   delimiter='\t')
+
+    # Savitzky-Golay filtering, lifted and adapted from mne-python (0.22)
+    def savgol_filter(self, h_freq):


you have to lift the test as well :) getting lazy? :P

I like the idea of a Savitzky–Golay filter, but I think this code should go in the smooth() method. I'm late to the party here so if this has already been discussed and resolved, my suggestion isn't a big deal.

cjayb · 2021-01-28T21:34:26Z

@rythorpe @ntolley could you weigh in on these? I'd like to get you immediate reactions before moving forward (yes: even with tests ;)

I guess it depends what you're trying to find. Welch will dampen any local effects but periodogram may show even artifacts

Our windows are so short that there's hardly a case where we'd benefit from Welch, or is there? The Hamming window applied in the periodogram should take care of our simulation edge-artefacts, too! Though it might be prudent to let the tmin-parameter actually crop the time used to calculate the spectrum.

ntolley · 2021-01-29T02:24:54Z

should the method be called plot_spectrum, for simplicity?

Agree that this sounds reasonable as well.

rather than reading from Dipole.units, might be better to have an explicit units-argument? Dipole.scale doesn't update the units, so we can easily end up in a mismatch. But we don't currently expose the scale-method, so I guess we can defer the decision for later?

Unit argument sounds good since a large amount of fitting waveforms to real data involves appropriate scaling. Perhaps we could rename Dipole.units, or at least emphasize in the docs that it reflects the original units from a single cortical column before scaling.

I'll need to refresh my memory on these spectral methods but will comment on the proposed functions soon.

ntolley · 2021-01-29T16:27:31Z

Back from some light reading. Definitely supportive of the savgol filter. Specifying the frequency seems quite useful, but I am particularly a fan of the shape preserving properties.

Tough decision with welch vs. periodogram. While this isn't the ideal example, the periodogram may prove more useful in practice. PSD's from this sort of experiment would typically be averaged over multiple trials. I believe this would have a similar benefit to welch's method if applied to a single longer time series?

To stay consistent with the lab's literature, I'll plan to put together an example for beta events that shows where averaged spectral methods can misleading (after the 0.1 release).

rythorpe · 2021-01-29T20:04:38Z

hnn_core/viz.py

-def plot_spectrogram(dpl, *, fmin, fmax, winlen=None, tmin=None, tmax=None,
-                     layer='agg', ax=None, show=True):
-    """Plot Welch spectrogram (power spectrum) of dipole time course
+def plot_periodogram(dpl, *, fmin, fmax, tmin=None, tmax=None, layer='agg',


What about plot_psd() to be explicit that we're plotting spectral density in units of power rather than e.g. spectral mass in units of energy?

rythorpe · 2021-01-29T20:14:11Z

hnn_core/viz.py

@@ -110,7 +110,7 @@ def plot_dipole(dpl, tmin=None, tmax=None, ax=None, layer='agg', decim=None,

    ax.ticklabel_format(axis='both', scilimits=(-2, 3))
    ax.set_xlabel('Time (ms)')
-    ax.set_ylabel('Dipole moment')
+    ax.set_ylabel(f'Dipole moment ({dpl[0].units})')


Suggested change

ax.set_ylabel(f'Dipole moment ({dpl[0].units})')

ax.set_ylabel(f'Dipole moment (scaled {dpl[0].units})')

Regarding the discussion about how to handle units, what if we add 'scaled' in the y-label and then add a label+legend to the figure that specifies the scaling factor. Here's an example snippet of how I've handled it in some of my figures:

jasmainak · 2021-01-29T20:18:05Z

might be a radical solution, but can we get rid of units attribute altogether? And we store always in the SI unit, Am like in MNE. For plotting, you can have a param called scalings which is 1e9 by default indicating that it's in nAm

cjayb · 2021-02-01T14:26:57Z

5dee224 re-implements smooth: behaviour depends on argument passed (post_proc needed a minor change, otherwise simple)
ab26e26 implements Welch PSD: I'm more convinced than ever that it's not what we want! Note that nfft is much longer than the trial duration (16k vs 11k). For meaningful smoothing at low frequencies, we'd need a lot longer trials
5d31f4b shows some how we might demonstrate the effects of smoothing, and the two PSDs.

Are we OK with smooth having this dual role? I may be a slight antipattern, but I feel it's semantically pretty clear to a user who just wants to 'smooth' their data for visualisation. I could improve the docstring, if that would help?

How about dropping Welch and re-implementing plot_periodogram to take in a list of Dipoles? That would be equivalent to a Welch windowing with no overlap and nfft==nperseg. WDYT?

Next up: units scaling.

jasmainak · 2021-02-01T19:43:51Z

Before we go into the rabbit hole of implementing more sophisticated periodograms, let's ask ourselves ... do we really need this or is a TFR plot good enough for everything? Because that's what I have seen in HNN so far. I feel the efforts might be better spent in trying to copy the TFR code from MNE. It's also easier to interpret in HNN contexts where we can have transient oscillations. The periodogram that I had in the alpha oscillation was a lazy effort to prevent an MNE dependency while showing alpha for cheap with scipy

jasmainak · 2021-02-01T19:46:47Z

Are we OK with smooth having this dual role?

Need to sleep on it a bit. You could consider the following alternatives:

dipole.smooth_savgol(30)
dipole.smooth_hamming(500)

or

dipole.smooth('savgol', 30)
dipole.smooth('hamming', 500)

also not entirely sure about forcing named arguments in functions with only a few arguments. Might be annoying for a new user who doesn't know Python that much.

cjayb · 2021-02-01T20:19:35Z

why does it return an instance if it works in-place?

Enables chaining, e.g.: `dpl.smooth(...).scale(...).plot(...)1

do we really need this or is a TFR plot good enough for everything? Because that's what I have seen in HNN so far.

Fair enough, but I do believe HNN has PSD plots too? The fundamental problem of too short data segments applies to either approach. With 300 ms of data it makes little sense even pretending to estimate power <10 Hz. The simulations we run are sampled at 40 kHz, but we know our aggregate dipoles are effectively silent above 200 Hz. I think it would be great to have a simple PSD routine available, such as the periodogram that could be used to average over trials to obtain a little of that sweet Welch feeling.

It's also easier to interpret in HNN contexts where we can have transient oscillations. The periodogram that I had in the alpha oscillation was a lazy effort to prevent an MNE dependency while showing alpha for cheap with scipy

Agreed, but I wouldn't sell the simple spectrum quite that short: many non-engineers balk at TFRs, but can easily get behind a power spectrum.

Need to sleep on it a bit. You could consider the following alternatives:

I guess I see what you're trying to avoid (e.g. forced named arguments), but having to memorise different smoothing methods is also an overhead. I still like my proposal best, but can easily be swayed :)

(P.S. Travis will have his pep8 with the next commit)

jasmainak · 2021-02-02T04:21:10Z

Enables chaining, e.g.: `dpl.smooth(...).scale(...).plot(...)1

Take a look at what happens in MNE. I don't think chaining is supposed to work with in-place operations, e.g., raw.filter

EDIT: I take it back. Looks like it's not how I was thinking :)

jasmainak · 2021-02-02T04:25:44Z

Travis will have his pep8 with the next commit

We should separate style checks into a separate CI build so it fails fast ...

hnn_core/dipole.py

hnn_core/viz.py

jasmainak · 2021-02-18T04:37:59Z

hnn_core/parallel_backends.py

@@ -371,6 +376,7 @@ def simulate(self, net, n_trials, postproc=True):
        """

        # just use the joblib backend for a single core
+        # XXX this does not work: if n_procs==1, empty dipole list returned


something feels fishy here. I'd rather fix it here than in the example if it's easy enough to do?

When you say "I'd" do you really mean "you'd"? ;) I'm a bit stumped about that behaviour. I don't know enough about the context manager to understand how you can nest a JoblibBackend().simulate inside a with MPIBackend...?

okay let me give this a shot after lunch

see db0555c

what would you do without my reviews? ;-)

oops I forgot to remove the XXX comment

jasmainak · 2021-02-18T04:38:28Z

anything to update in whats_new.rst?

hnn_core/viz.py

cjayb · 2021-02-18T12:27:20Z

Did some work on making the visualisations more robust to different label widths and subplot combos (this is a great resource).

I made all the implicit post-processing stuff private, so we can easily remove them in the next round. Behaviour remains GUI-like (smoothing and scaling read from param-file by default).

Beginning to feel pretty done here...

Edit: I'm not impressed by CodeCov, it doesn't seem to be updating and I can't figure out how to trigger it!

ntolley · 2021-02-18T14:57:47Z

examples/plot_firing_pattern.py

@@ -100,7 +100,7 @@
 ###############################################################################
 # We can plot the firing pattern of individual cells by indexing with the gid
 gid = 170
-plt.figure(figsize=(4, 4))
+plt.figure(figsize=(4, 4), constrained_layout=True)


Maybe use plt.tightlayout() instead? It seems that this setting may break our figures in a matplotlib update:

unfortunately this doesn't work for more complicated plots (like the one in plot_simulate_gamma). Without some extra effort, many of our plots get the ylabel cut out (see e.g. the firing rate demo in our master docs). The challenge is to find something that respects sharex=True, and allows variable-width labels (1-2 lines), and understands colorbars... Any suggestions welcome! (tight_layout doesn't cut it)

I wonder if we need to worry about that warning though? I think hnn-core will be developing a lot faster than matplotlib, so I'm sure we'll come you with new plots faster than they can remove functionality from us ;)

can't we just get rid of this argument? Does it look too bad? The aim of the examples is not to make beautiful figures. People can always add plt.tight_layout to their scripts.

But the warnings don't look too bad. I am fine with letting it stay :)

hnn_core/dipole.py

jasmainak · 2021-02-18T17:14:16Z

Other than MPI issue, PR looks good to me. I'll give it a shot myself after lunch. Thanks @cjayb for the awesome PR as always!

jasmainak · 2021-02-24T13:45:21Z

.circleci/config.yml

@@ -26,7 +26,7 @@ jobs:
              command: |
                  export PATH=~/miniconda/bin:$PATH
                  conda update --yes --quiet conda
-                  conda create -n testenv --yes pip python=3.6
+                  conda create -n testenv --yes pip python=3.7


Are you certain this has nothing to do with this PR? One way to verify is to make a "blank PR" with a tiny doc change and see if you have the same CircleCI problem. If so, I can live with this. It's not good to create such tight bleeding edge Python dependency. Because Neuron doesn't work on 3.8, it means we have only one version of Python for which MPI backend would work

Switching back to 3.6, works on master (#272)

jasmainak · 2021-02-24T13:48:48Z

hnn_core/parallel_backends.py

+        dpls[-1]._baseline_renormalize(N_pyr_x, N_pyr_y)  # XXX cf. #270
+        dpls[-1]._convert_fAm_to_nAm()  # always applied, cf. #264


this is different from the definition of a "raw dipole" in the GUI? I am fine with this but perhaps GUI users should weigh in.

@kohl-carmen @rythorpe any opinions?

I'm not sure what a "raw dipole" is?

jasmainak · 2021-02-24T13:50:20Z

hnn_core/dipole.py

        for key in self.data.keys():
            self.data[key] = _hammfilt(self.data[key], winsz)

+        return self
+
+    def savgol_filter(self, h_freq):


no test for this function?

Only implicitly (test_viz.py:L36). Eric implemented some actual tests on this in mne, but I decided against copying those, since (again) we don't have that much data to make FFTs of...

jasmainak · 2021-02-24T13:51:40Z

@cjayb I left a few more comments. But let me know if you're too saturated with this PR, we can move those discussions elsewhere

jasmainak · 2021-02-25T12:52:00Z

examples/plot_simulate_mpi_backend.py

@@ -63,7 +63,8 @@
 # ``openmpi``, which must be installed on the system
 from hnn_core import MPIBackend

-with MPIBackend(n_procs=2, mpi_cmd='mpiexec'):
+n_procs = 1


think we should revert this to n_procs=2 as it was before? not sure if that was causing the CircleCI issue

Will do. In fact, #272 passed just fine with python=3.6, so I'll revert it here too. I'm hoping it's just been a fluke. I really don't see how this PR could have broken MPI...

jasmainak · 2021-02-25T13:28:57Z

All good merging. Thanks a lot @cjayb !! 🥳

cjayb · 2021-02-25T13:30:29Z

Sweet!

cjayb mentioned this pull request Jan 27, 2021

[MRG] DOC: Update alpha example #249

Merged

ntolley reviewed Jan 27, 2021

View reviewed changes

examples/plot_simulate_alpha.py Show resolved Hide resolved

cjayb changed the title ~~[MRG] minor alpha COSMITs (after #249)~~ [WIP] minor alpha COSMITs (after #249) Jan 27, 2021

jasmainak reviewed Jan 28, 2021

View reviewed changes

rythorpe reviewed Jan 29, 2021

View reviewed changes

cjayb changed the title ~~[WIP] minor alpha COSMITs (after #249)~~ [WIP] updated visualisation methods and unit handling; alpha update Feb 1, 2021

MAINT small cosmits after jonescompneurolab#249

f049bff

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/dipole.py Show resolved Hide resolved

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/dipole.py Outdated Show resolved Hide resolved

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/dipole.py Outdated Show resolved Hide resolved

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/dipole.py Outdated Show resolved Hide resolved

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/viz.py Outdated Show resolved Hide resolved

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/viz.py Show resolved Hide resolved

cjayb added 3 commits February 18, 2021 12:20

ENH visualisations and add whats new

dbbcf75

missed a plot (constrained_layout)

5e6f9a4

MAINT make post proc methods private and adjust tests

97d9d47

cjayb changed the title ~~[MRG] updated visualisation methods and unit handling; alpha update~~ [MRG] updated visualisation methods and unit handling Feb 18, 2021

ntolley reviewed Feb 18, 2021

View reviewed changes

jasmainak reviewed Feb 18, 2021

View reviewed changes

hnn_core/dipole.py Outdated Show resolved Hide resolved

Cosmits

db0555c

jasmainak force-pushed the minor_alpha_nitpicks branch from e1e0be3 to db0555c Compare February 24, 2021 13:36

jasmainak reviewed Feb 24, 2021

View reviewed changes

cjayb mentioned this pull request Feb 25, 2021

trigger CI run on master #272

Closed

jasmainak reviewed Feb 25, 2021

View reviewed changes

cjayb and others added 2 commits February 25, 2021 14:05

MPI n_procs back to 2; CicleCI python=3.6

2bc577f

DOC: remove superfluous comment

49c9e33

jasmainak merged commit 7f050b9 into jonescompneurolab:master Feb 25, 2021

cjayb mentioned this pull request Mar 16, 2021

[MRG] Add ability to modify cell-cell connectivity #276

Merged

	ax.set_ylabel(f'Dipole moment ({dpl[0].units})')
	ax.set_ylabel(f'Dipole moment (scaled {dpl[0].units})')

		dpls[-1]._baseline_renormalize(N_pyr_x, N_pyr_y) # XXX cf. #270
		dpls[-1]._convert_fAm_to_nAm() # always applied, cf. #264

[MRG] updated visualisation methods and unit handling #264

[MRG] updated visualisation methods and unit handling #264

Conversation

cjayb commented Jan 27, 2021 • edited Loading

ntolley commented Jan 27, 2021 • edited Loading

jasmainak commented Jan 27, 2021 • edited Loading

jasmainak commented Jan 27, 2021

jasmainak commented Jan 27, 2021

ntolley commented Jan 27, 2021

jasmainak commented Jan 27, 2021

codecov-io commented Jan 27, 2021 • edited Loading

Codecov Report

cjayb commented Jan 27, 2021

jasmainak commented Jan 27, 2021

jasmainak commented Jan 28, 2021 • edited Loading

jasmainak commented Jan 28, 2021

jasmainak commented Jan 28, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cjayb commented Jan 28, 2021 • edited Loading

ntolley commented Jan 29, 2021

ntolley commented Jan 29, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jasmainak commented Jan 29, 2021

cjayb commented Feb 1, 2021

jasmainak commented Feb 1, 2021

jasmainak commented Feb 1, 2021

cjayb commented Feb 1, 2021

jasmainak commented Feb 2, 2021 • edited Loading

jasmainak commented Feb 2, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jasmainak commented Feb 18, 2021

cjayb commented Feb 18, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jasmainak commented Feb 18, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jasmainak commented Feb 24, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jasmainak commented Feb 25, 2021

cjayb commented Feb 25, 2021

cjayb commented Jan 27, 2021 •

edited

Loading

ntolley commented Jan 27, 2021 •

edited

Loading

jasmainak commented Jan 27, 2021 •

edited

Loading

codecov-io commented Jan 27, 2021 •

edited

Loading

jasmainak commented Jan 28, 2021 •

edited

Loading

cjayb commented Jan 28, 2021 •

edited

Loading

jasmainak commented Feb 2, 2021 •

edited

Loading

cjayb commented Feb 18, 2021 •

edited

Loading