Add imf function #19
Add imf function #19
Conversation
* TODO.md: Remove 'imf' from list. * docs/src/ref.md: Add "imf" entry to the manual. * src/utils.jl: Add imf entry to "utils". * src/imf.jl: Contains code of imf function. * test/util-tests.jl: Include tests for imf. Values now match (with the normal difference in least significant digits) with those of IDL AstroLib.
Codecov Report
@@ Coverage Diff @@
## master #19 +/- ##
==========================================
+ Coverage 99.55% 99.56% +<.01%
==========================================
Files 63 64 +1
Lines 1125 1144 +19
==========================================
+ Hits 1120 1139 +19
Misses 5 5
Continue to review full report at Codecov.
|
src/imf.jl
Outdated
| println("Length of array mass_range is not one more than that of expon") | ||
| return zeros(mass) | ||
| end | ||
| integ = ones(T, ne_comp) |
There was a problem hiding this comment.
The initial value of integ doesn't look to be relevant. You can save some nanoseconds by initializing the vector with Vector{T}(ne_comp).
src/imf.jl
Outdated
| end | ||
| integ = ones(T, ne_comp) | ||
| joint = ones(T, ne_comp) | ||
| norm = ones(T, ne_comp) |
There was a problem hiding this comment.
I don't think you need to initialize the vector here, this line can be removed.
src/imf.jl
Outdated
| joint[i] = joint[i-1]*(mass_range[i]^(expon[i-1] - expon[i])) | ||
| end | ||
| end | ||
| norm = (1/sum_kbn(integ.*joint))*joint |
There was a problem hiding this comment.
Is there a particular reason why you used sum_kbn instead of the simpler sum? Did you find inaccurate results with it? Anyway, you can save an operation by using joint as numerator:
norm = joint ./ sum(integ .* joint)There was a problem hiding this comment.
I thought sum_kbn would give a significant boost in accuracy (as I was brought to believe from the docs). But apparently sum_kbn advantage was almost completely nerfed by JuliaLang/julia#4039. So yeah, using sum would be better.
There was a problem hiding this comment.
sum_kbn should be indeed more accurate, for example in the case of elements spanning several orders of magnitude, but if this isn't the case we shouldn't worry too much ;-)
Anyway, I just realized that sum(integ .* joint) is simply the scalar product between the two vectors, so dot(integ, joint) is better. Sorry for not having noticed before.
src/imf.jl
Outdated
|
|
||
| ### Example ### | ||
|
|
||
| ```julia |
There was a problem hiding this comment.
This line should probably go below ;-)
src/imf.jl
Outdated
| mass_range::AbstractVector{T}) | ||
| ne_comp = length(expon) | ||
| if length(mass_range) != ne_comp + 1 | ||
| println("Length of array mass_range is not one more than that of expon") |
There was a problem hiding this comment.
IDL AstroLib throws an error here, instead of returning a vector of zeros. Why the different behavior?
There was a problem hiding this comment.
The option of using exception handling completely slipped my mind. Too much worrying about keeping return-type stable ;-)
Change sum() to dot() as only scalar product is involved
I thought of making mass arguement as scalar and then writing another imf function that'll take the mass vector indices one by one (as in precess_xyz), but just restricting the mass argument to vector made things simpler.