More statistics functions #115
Conversation
src/Maths.jl
Outdated
| spl = CSpline(x, y) | ||
| lb = xmax - 2*(xmax-left) | ||
| ub = xmax + 2*(right-xmax) | ||
| f(x) = spl(x) - val |
There was a problem hiding this comment.
if you make the spline of y - val we can take the analytic roots of the spline, possibly faster. We haven;t implemented this yet, but it is easy to get from the roots of each piecewise polynomial.
There was a problem hiding this comment.
should be much faster, yes. I tried to look up how to do that but didn't get very far--do you have a reference where it's explained?
There was a problem hiding this comment.
There is an old fun fortran 77 implementation here: http://www.netlib.org/dierckx/sproot.f
The comments explain the methodology.
Scipy has an implementation. Their spline is at: https://github.com/scipy/scipy/blob/ab2c51ca7d9eaf2ee93f89a9dd00ed67546999f4/scipy/interpolate/_cubic.py#L449
This inherits from PPoly which is implemented here: https://github.com/scipy/scipy/blob/ab2c51ca7d9eaf2ee93f89a9dd00ed67546999f4/scipy/interpolate/interpolate.py#L1165
| function energy_λ(grid, energyfun_ω, λlims; label=nothing, winwidth=0) | ||
| λlims = collect(λlims) | ||
| ωmin, ωmax = extrema(wlfreq.(λlims)) | ||
| window = Maths.planck_taper(grid.ω, ωmin-winwidth, ωmin, ωmax, ωmax+winwidth) |
There was a problem hiding this comment.
a luxury item would be to be able to specify the window type/shape
There was a problem hiding this comment.
there is now a version which you can pass your own window function to, I figured that's the most flexible way
| return Eta, analytic! | ||
| end | ||
|
|
||
| function plan_analytic(grid::Grid.RealGrid, Eω) |
There was a problem hiding this comment.
can you use your new plan_hilbert here?
There was a problem hiding this comment.
you could, but it's two extra transforms--you'd first need to get the field from Eω and then do the Hilbert transform. getting the analytic field is the main cost of the statistics functions so I wanted it to be as fast as possible
There was a problem hiding this comment.
I thought about interleaving this with the windowing function (defined in Luna.run), so you already have the time-domain field, but that ties together very separate bits of the code
src/Stats.jl
Outdated
|
|
||
| function electrondensity(Presp, dfun) | ||
| function addstat!(d, Eω, Et, z, dz) | ||
| d["electrondensity"] = maximum(Presp.fraction)*dfun(z) |
There was a problem hiding this comment.
You're relying on internal stats of the Plasma polarisation here. This is not ideal. Along with the fact that you get just the last step. I would propose that we supply the nonlinear responses to the stats functions, along with a way to get the peak field. That way they can interrogate the nonlinear response at that step. In fnfep I just take the on axis electric field and calculate the electron density from that.
There was a problem hiding this comment.
by "peak" you mean in space, right? seeing as you need the whole time axis to calculate the electron density. then, if you want to be precise (i.e. have exactly the value used in the propagation), you'd also need to oversample the field onto the finer time-grid first.
There was a problem hiding this comment.
I've tested it for the mode-averaged case, and it turns out the way I had it here doesn't work at all--compared to calculating it from the field, the density plot is shifted and the magnitudes are wrong.
The oversampling doesn't seem to matter too much, though this may depend on the degree of ionisation:
julia> Maths.fwhm(grid.t, It, :nearest, minmax=:min)*1e15
30.000177980675797
julia> Maths.fwhm(grid.t, It, :spline, minmax=:min)*1e15
30.00000000548264
|
src/Maths.jl
Outdated
| #spline method | ||
| try | ||
| spl = BSpline(x, y .- val) | ||
| r = roots(spl) |
There was a problem hiding this comment.
For horrible functions the number of roots could be very high, we probably need to increase the argument to the Dierckx roots function. See here:
https://github.com/kbarbary/Dierckx.jl/blob/master/src/Dierckx.jl#L378
There was a problem hiding this comment.
any feeling for how many might be enough?
There was a problem hiding this comment.
Well for a "noise burst" it could be thousands. That might come with a performance hit though. To be clear, passing a high maxn to Dierckx doesn't cost anything if there are not many roots. But if there are thousands then it will find all of them.
This isn't just academic, when looking for pulse convergence in say modelocking (or a FOPO!), a key metric is if the min and max FWHM match. We could be doing a lot of that quite soon. I suppose in that case we'd specify a different method here, and only use the spline for more precise requirements.
It might be worth adding a "noise burst" test though.
There was a problem hiding this comment.
With this pulse:
it finds around 600 roots, FWHM max still works, and it still only takes 1.5 ms with 2^14 points
There was a problem hiding this comment.
funnily, the warning in Dierckx.jl when more roots than maxn were found is now an exception, because warn no longer exists. Medium to long term it would be good to no longer depend on that package.
There was a problem hiding this comment.
Agreed, but Dierckx (the fortran part) is by far the highest quality and most tested spline package I've found so far. And the wrapper code is rather minimal, so we can help maintain it. If I had sufficient free time I'd translate it, but alas I can;t see that happening very soon.
|
Do you plan to implement a peak intensity function too? I use that quite a lot, but it can be added later. |
I was until I forgot. Just adding it now. I take it you won't want to assume that the peak is at the centre of the waveguide? |
In fnfep I do assume this, but I had an optimistic TODO to fix this. If we can not make this assumption it is better, but it quickly becomes a 3D global optimisation problem to solve this in general, so I'm not sure how easy it is. |
For now I'm going to assume the same thing. I'll open an issue to keep track |
hmmmm good point. I guess I could just save both the on-axis and mode-averaged intensity, since one of them is just |
|
Honestly I'd just choose one convention and stick to it everywhere for consistency. |
|
This is basically done, but there are a few things we're not saving explicitly, e.g. the name of the gas species. Of course you could add that manually, but I think the only way to really cover everything is to move to a higher-level interface like #129 , i.e. something more declarative with a fixed set of input parameters. Most of the time, that is what we will use anyway. Then, running scripts like we do right now is the "advanced" interface with a big "no warranty provided" sign, where the user has to make sure they save everything they need. Does that make sense? |
|
Agreed. As I understand it we can already save arbitrary extra info that we choose to anyway. Note that to convert types to dicts two really useful functions are
The first easily produces a text representation which could be saved as a reference. The second could be useful for iterating over the fields of a struct. We could have a quite simple function which iterates over the fields of an arbitrary object and saves them if they have a reasonable type. Though that would be close to a lightweight JLD. |
|
Anyway, go ahead and merge and we can add more features later. |
Apologies, I see that this is exactly what you use! So much for my reviewing! |
|
Using |
|
Tests still pass |
This adds statistics functions for the following:
zand stepsizedzWhat's missing at the moment is anything to do with space, i.e. peak intensity, beam size, and anything for free-space propagation