JointFitter and two examples.

[KriSun95]

Adds the first version of the new JointFitter and two examples for the example gallery.

This is still a work in progress. My local docs building is not working.

To add:

• changelog
• get local docs build working
• make tests
• run pre-commit

Thoughts

This is workable but there are still somethings to do or clean-up. I'm not sure if the following should be this PR or another one.

• We have to decide how to choose and optimiser and statistic. At the minute, JointFitter needs to be initialised with one of each since I'm trying to follow their set-up. The issue is, not all of their other fitters follow this and end up just bypassing this and hardcoding in their optimiser and statistic choices anyway. How do we want to do this?
• Unit support! Right now, JointFitter just handles units briefly and crudely to get stuff to work. What I've written isn't complicated and we should look into altering the Astropy decorator that handles units.
• The thick target function has its issues with input+output units and also returning NaNs for some very high precision input values. I've not included the thick target in any example here for that reason. This should be another PR but I wasn't sure if that should be one on its own or with one of Jake's PRs. I just wanted to mention it here.

@jajmitchell , @samaloney , let me know what you both think about this. I think the JointFitter doc-strings could use a little work but other than that, what would be crucial to get this initial stuff merged in your eyes?

[settwi]

Nice. I want to add a sensible “what if” to think about.

What if a user wants to use a prior distribution that is not a uniform prior? Bounding the parameters is equivalent to a uniform prior. Say I want to use a truncated normal distribution on a parameter, or a Jeffreys prior. How hard would that be to implement as this code is written, or as Astropy modeling is implemented?

[KriSun95]

@settwi , that's a good point and worth thinking about.

The permanent fix would be to PR Astropy

Personally, I think that the place to allow this would be in the module level fitter_to_model_params_array function but in Astropy main. This function is the one responsible for building up the parameter list for a model from the partial list of fitted parameters and is the one that enforces the bounds in the first place.

If we could just allow the bounds attribute for the parameter object to accept a prior function of the user's choice (either in place of a simple length 2 tuple or as well as) then I think the code-block that enforces the given bounds at the minute (here) could just run the value through the parameter's "bound function" to get the corresponding value.

This would make the bounds on any parameter far more powerful for any fitting.

The temporary fix (while we wait for Astropy to accept)

I can think of two ways:

1. We can edit our fitter_to_model_params_array to do this.
2. Add a prior attribute to our model parameters (like temperature, etc.) that we can check for once the parameters are returned from fitter_to_model_params_array in our own fitter objects.

[samaloney]

This is really nice!

I wonder would it make sense to extract an ABC which defines the API then make this a concrete implementation?

The tests cover:

• unitless models,
• forward-direction tying
  but misses:
• fixed params
• bounds enforcement check,
• weights,
• Quantity units.

The units-handling branch and the weights/fkwarg paths used in are untested. I do think tests for any feature expose by the API would be needed - could be done in a later PR but would be better here.

It would be good to try and find some specific examples of the thick target function issues so we can pin it down and hopefully fix.

On the optimiser and statistic I think the way it set up now makes sense if need to define our own API and then map the the underling fitters - other subclasses can override the methods and raise error is not supported

[samaloney]

As partially noted in the PR description the __init__ validates and stores self._opt_method/self._stat_method here always calls scipy.optimize.minimize(..., method="Nelder-Mead") with a hardcoded sum-of-squares.

Both gallery examples call JointFitter(optimizer=SLSQP, statistic=leastsquare) which doesn't actually use the optim/stats functions.

[KriSun95]

Absolutely, and I don't like this either. I wanted to get the code up here to discuss first but I definitely thought this is something to fix. I'm hoping my pseudo code gets at this and I'm happy to build on this.

[samaloney]

mod_params is only assigned when parameter_units is not None. It's unreachable via the public call path (_get_param_units never returns None for non-empty model lists)

[samaloney]

Should see if we can eventually PR this to astropy

[KriSun95]

Absolutely. I Think that pretty much all the JointFitter code and fitters.py code as it stands could live in Astropy eventually. The only thing I've essentially done (after a lot of work to get there) was to:

• Get the fitter to accept (model, x, y) groups instead of just one
  • This mainly meant I had to add in some loops to loop through the groups (+minimal other things)
• Make those two module level functions (only because one of these functions uses the other) accept this model_list parameter

I don't see why the other fitting classes in Astropy couldn't just be made to do this. Then we only have to have our own fitters like the NelderMeadLeastSquaresFitter fitter with our own optimiser and statistics modules to let us do whatever we want.

[KriSun95]

This is really nice!

I wonder would it make sense to extract an ABC which defines the API then make this a concrete implementation?

Thanks! I've been happy with how it's turned out so far.

I would love to do that! It's something that the API in Astropy isn't very solid in and I don't like exactly how I've done it here.

The way Astropy does it is by defining the base Fitter class which any other fitter needs to inherit from (I also inherit from Fitter for our JointFitter). These other classes then redefine the same functions as I did, like objective_function, _run_fitter, and __call__. The thing that annoys me with the Astropy fitters is that a lot of them inherit from Fitter but never actually use it (i.e., they never call super in their init) then just hard code everything.

I can try and make JointFitter a sort of base class that needs to be initialised with an optimiser and statistic. Then inherit from it in a new class, like NelderMeadLeastSquaresFitter, that we can use in the same way. This might mean we can just be a bit more rigid the Astropy in how new fitters should be created.

E.g.,
In fitters.py:

from scipy.optimize import minimize
from some_module import least_squares_func

class JointFitter(Fitter):
    def __init__(self, optimizer, statistic):
        super().__init__(optimizer=optimizer, statistic=statistic)

    def joint_model_to_fit_params(self, models):
        ...

    def _update_model_params(self, models, fps, jfit_param_indices):
        ...

    def _verify_input(self, args):
        ...

    @staticmethod
    def _get_param_units(models):
        ...

    def _run_fitter(self, models, farg, fkwarg=None):
        raise NotImplementedError("Subclasses should implement this method.")

    def objective_function(self, models, farg, fkwarg=None):
        raise NotImplementedError("Subclasses should implement this method.")

    def __call__(self, *args, fkwarg=None, inplace=False):
        # will check inputs then return output of `self._run_fitter`
        ...
        return self._run_fitter(...)


class NelderMeadLeastSquaresFitter(JointFitter):
    def __init__(self):
        super().__init__(optimizer=minimize, statistic=least_squares_func)

    def objective_function(self, fps, *args, weights=None, jfit_param_indices=None, parameter_units=None):
        # gathers parameters and evaluates model to run`self._stat_method` which is `least_squares_func` here
        ...
        self._stat_method()

    def _run_fitter(self, models, farg, fkwarg=None):
        # sets up and runs `self._opt_method` which is `minimize` here
        ...
        self._opt_method(...)
        ...

The user could then just:

fitter = NelderMeadLeastSquaresFitter()
result = fitter(model0, x0, y0, ..., modelN, xN, yN)

The JointFitter class could realistically replace the base Fitter class from Astropy eventually...I think.

The tests cover:

• unitless models,
• forward-direction tying
  but misses:
• fixed params
• bounds enforcement check,
• weights,
• Quantity units.

The units-handling branch and the weights/fkwarg paths used in are untested. I do think tests for any feature expose by the API would be needed - could be done in a later PR but would be better here.

Agreed! I'll work through this list and update the tests to account for them.

It would be good to try and find some specific examples of the thick target function issues so we can pin it down and hopefully fix.

Agreed again and I'm working on another branch to do this.

On the optimiser and statistic I think the way it set up now makes sense if need to define our own API and then map the the underling fitters - other subclasses can override the methods and raise error is not supported

Have I answered this with the pseudo code above?

[KriSun95]

New tests to add:

• unitless models,
• forward-direction tying
• fixed params
• bounds enforcement check,
• weights,
• Quantity units

[KriSun95]

Don't believe the failed ubuntu-latest test. I've seen it fail then pass with no change to the code here and in PR #289

[KriSun95]

Minimum bar to get this merged:

• Update joint thermal fitting example to use the Constant model to account for differently observed photon spectra rather than just using different emission measure.
  • @jajmitchell and @samaloney , the new example is live using the constant. I fixed a small bug in the model too while I was at it.
• Update the API to the pseudo code above. Might as well start out closer to something more useable.
  • [~] I want to include another example (maybe even a docs page) showing how the new API can be used to easily work with different fit statistics, optimisation methods, and beyond. Never mind, let's just wait until later, not need to overcomplicate things just now.
• @samaloney will look at the current tests and changes may be requested.

[KriSun95]

I know what's broken, I moved some stuff while making our starting API clearer and it broke some obvious things I forgot were there in the examples. I will fix this soon.

[samaloney]

Few small comments, fixes required, really like the example now but main topic is the tests.

Only test_JointFitter_initial_values isolates a single unit (joint_model_to_fit_params).

The other three each run a full scipy.optimize.minimize fit and assert on final parameter values, in particular test_ScipyMinimizeJointFitter_units_bounds_fixing_weights bundles four features at once (units + bounds + fixing + weights). If it fails, you can't tell which one broke without debugging. There's also no direct test of _update_model_params, _evaluate_models, _assign_param_units/_get_param_units, or objective_function in isolation from the optimizer loop.

I'd suggest splitting the combined test into one narrow test per feature (units-only, bounds-only, fixed-only, weights-only) plus keeping one "everything together" integration test, and adding a couple of direct unit tests for _update_model_params/_evaluate_models/fitter_to_model_params_array that check parameter reconstruction without invoking scipy at all.

Obviviosy this is more test code to maintain, but you get much faster failure localisation and don't need a full optimiser convergence to verify parameter bookkeeping.

This could be handled in a separate issue (that would have to be done before a release) to not block progress.

[samaloney]

hum should this be deleted as other wise we have two chi_squared functions the new tuple based (below) one and the old array in this renamed file?

[samaloney]

Missing from the __all__ above

[samaloney]

Suggested change

	raise ValueError("Inputs `data_ys_tuple` and `model_ys_tuple` must be tuples of of same length.")
	raise ValueError("Inputs `data_ys_tuple` and `model_ys_tuple` must be tuples of same length.")

[samaloney]

Hum I wonder what the performance hit is for this, maybe open an issue to profile later.