fixed calculation of energy

[kpannek]

Hi @jdtournier,
As discussed, these are the fixes for the calculation of the energy.
dirmerge, for some reason, used 1.2/r. I've changed it to 1.0/r. Please have a look if this is correct as well.
Thanks

[jdtournier]

Thanks @kpannek, most welcome. Just to give a bit of context here, Kerstin pointed out that the energy term for electrostatic repulsion (F ∝ 1/r²) should actually scale as E ∝ 1/r, whereas I had the energy scaling down as 1/r²... So this fixes this in the places where it's used. I'll merge in the equivalent changes to dirgen, which do lead to slightly different results - I'll post these here in a minute.

[jdtournier]

dirmerge, for some reason, used 1.2/r. I've changed it to 1.0/r.

The reason for that is given in the preceding comment - although that might not be clear. The 1.2 factor means that there is a 0.2 contribution from the unipolar model, which is there to try to promote uniformity over the sphere, in an effort to avoid any biases due to eddy-currents. I'll make that a more explicit option...

[jdtournier]

OK, I think I've done all I wanted to do. Still need to double-check that all these commands are still behaving sensibly, but I doubt that would be a huge issue. Simplest will be to run gen_scheme and see what comes out - that script exercises most of these commands.

To give you an idea of what impact these changes have on dirgen, after some of the changes above (including multiple restarts, as originally suggested by @thijsdhollander way back when), command is:

 $ dirgen 60 dir60_new.txt

Using previous version:

$ dirstat dir60_old.txt
dir60_old.txt [ 60 directions ]

  Bipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 18.8232, range [ 18.302 - 19.5847 ]
    energy: total = 3222.42, mean = 107.414, range [ 107.235 - 107.7 ]

  Unipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 19.1198, range [ 18.302 - 23.0698 ]
    energy: total = 1614.56, mean = 53.8188, range [ 46.8849 - 60.2786 ]

  Spherical Harmonic fit:
    condition numbers for lmax = 2 -> 8: 1.00612 1.01637 1.03319 1.15016

and using new version:

$ dirstat dir60_new.txt
dir60_new.txt [ 60 directions ]

  Bipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 18.7958, range [ 18.2769 - 19.5421 ]
    energy: total = 3222.41, mean = 107.414, range [ 107.269 - 107.654 ]

  Unipolar electrostatic repulsion model:
    nearest-neighbour angles: mean = 19.3551, range [ 18.2769 - 31.3355 ]
    energy: total = 1609.95, mean = 53.665, range [ 44.9445 - 59.834 ]

  Spherical Harmonic fit:
    condition numbers for lmax = 2 -> 8: 1.00423 1.01308 1.02902 1.15597

So the main difference in a slight improvement in the SH fit for the lower orders, slightly worse for the highest one. Not too sure what's best here... But there's nothing stopping anyone from using a higher power if they want to.

[jdtournier]

OK, just run a couple of instances of gen_scheme with this branch versus master, there are minor differences as expected, but they are relatively small, bordering on negligible. I think this is safe to merge, certainly to the dev branch with a view to it becoming available in the next tagged release.

@kpannek, have a look through to make sure you're happy with the changes, and assuming everyone's OK with this, I'll merge to dev tomorrow.

[kpannek]

Thanks @jdtournier . Looks all good to me!

[thijsdhollander]

Looks all good to me! Very happy with the multiple restarts (it's something I do in practice every time already anyway; and would always recommend). The change / "fix" to the repulsion law makes sense. I actually remember wondering about this myself when I started using this functionality in MRtrix a long time ago (I used to just use my own Matlab script with the same kind of functionality for this). Is there a good (or any) reference for including such a power parameter to optimise this problem, rather than going with default Coulomb for the Thomson problem? Interesting impact on the condition numbers; especially it being different (in direction) across lmax'es. The difference (and the direction of it for power = 1 vs 2) may very well depend on the number of directions you're asking for too.

There's a few very particularly challenging numbers for this problem by the way. After some very thorough investigations a year (or something) ago, I would for instance advise (and have advised a practical case) against going for N=45. It's a particularly annoying case. There's a very high risk of ending up with a particularly annoying local configuration somewhere on the sphere for that case. It's these things that may impact condition numbers in pretty surprising (but then again: not really surprising, if you think about it) ways. But then again, it's also these things that aren't always reflected very well by the condition numbers. N=44 and N=46 are much more well behaved. As with many things, a combination of quantitative as well as qualitative checks are the best recipe to end up with an optimal outcome (set of directions).

It does look like the test data for dirgen will have to be regenerated for the test to pass.

With the current pull request to go to master being tag_3.0_RC2 rather than dev, it may be good to rebase this pull request to that branch. As far as I'm concerned, I'd say these changes are entirely ready to end up on master as soon as possible.

[jdtournier]

OK, test data updated, and this pull request is now set to merge to tag_3.0_RC2. Not keen to merge to master directly since this does modify results (even if only a little), and isn't actually much of a problem anyway...

No literature that I'm aware of that looks into the effects of the power, but it is used in an early paper by Papadakis et al.. There's not much discussion of the effect of the power though, other than just mentioning that a low power has been shown to be suboptimal, apparently. This is what I'd originally also observed based on the nearest-neighbour angle - but as we now know, that's not the whole story: lower powers typically give up better conditioned SH fits, which is a much better indication of practical performance in my opinion...

Anyway, I reckon this is probably good to go, as soon as the tests have passed.