Bayesian extend/resume#257
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## develop #257 +/- ##
===========================================
+ Coverage 82.92% 83.30% +0.37%
===========================================
Files 62 63 +1
Lines 5084 5253 +169
===========================================
+ Hits 4216 4376 +160
- Misses 868 877 +9
Flags with carried forward coverage won't be shown. Click here to find out more.
|
compatibility for resampling API
damskii9992
left a comment
There was a problem hiding this comment.
Just gonna stop my review here as there is already plenty to work on.
As we discussed in-person, I also think we should refactor the bayesian sampling to a seperate "samplers" module. This should help us when we add other Bayesian samplers in the future as well as provide a cleaner user API (the user doesn't create a "Fitter" when doing sampling).
It also more cleanly separates responsibility in the code. A Minimizer minimizes a function to the data, a Sampler samples the parameters to create distributions.
To avoid duplicating code, it might be smart to move shared functionality to new base-classes.
damskii9992
left a comment
There was a problem hiding this comment.
I'm gonna stop my review again here as the changes I have proposed are again substantial and another review would likely be required afterwards.
My biggest issue lies with the class/file structure. I think we should copy the factory/minimizer structure for the samplers, as this is an easily extendable and logically sounds structure and it makes reviews easier as the structure would be the same across our codebase. But we probably need an ADR for this and to have discussion about it, probably in-person.
| "extended_draws = extended_results.draws\n", | ||
| "\n", | ||
| "print(f'Short chain: {short_draws.shape[0]} retained draws')\n", | ||
| "print(f'Extended chain: {extended_draws.shape[0]} retained draws')" |
There was a problem hiding this comment.
The output here is that the extended chain has 7000 draws. Adding 5000 with a thin of 2 should make the chain 2500 longer. Clearly we're applying the burn-in of 500 draws again, but why are we doing that when we extend the chain? Shouldn't that only happen when you initiate the chain?
There was a problem hiding this comment.
Nope. the burn-in is not re-applied. Sampler.extend() passes burn=0, and Bumps.mcmc_sample sets burn=0 whenever resume_state is not None, warning if the caller asked otherwise. There is no path that re-burns a resumed chain.
There was a problem hiding this comment.
I'll update the notebook/docs with proper statements
There was a problem hiding this comment.
2-parameter model, sample(samples=10000, burn=500, thin=2) then extend(additional_samples=5000, thin=2)):
| Config | sample() draws |
extend() draws |
Delta |
|---|---|---|---|
default (trim=True) |
4520 | 6600 | +2080 |
sampler_kwargs={'trim': False} |
5000 | 7500 | +2500 |
With trimming off, the numbers are exactly as you expect. The ring buffer also grows exactly as intended: Nthin * Npop == 7500, and state.draw(portion=1.0, outliers=True) returns all 7500 points. No draws are lost, and none are re-burned.
Running the actual tutorial (4 parameters, not the 2 above) reproduces your number exactly and shows where it comes from:
| chain | results.draws (trimmed) |
full chain |
|---|---|---|
| original | 4400 | 5000 |
| extended | 7000 | 7600 |
So the "missing" draws were never missing: the chain grew from 5000 to 7600.
There was a problem hiding this comment.
Okay, but this still doesn't explain everything.
So the "Full" chain included the burn-in but not the thinned-out, since the first full chain is 5000 steps, but the extended is 2600 larger, which is more than half of the 5000 asked for. You write that DREAM rounds up to nearest whole generation, so I guess that must be 200 steps? Is that the default or is that a parameter we can change? It is weird that it runs longer than you ask it to.
| Construct directly with a configured ``Fitter`` (or ``MultiFitter``) whose | ||
| minimizer has been switched to ``AvailableMinimizers.Bumps``. |
There was a problem hiding this comment.
When the Sampler is constructed from the Fitter it appears like you have to do a fit first, which you don't.
Instead of having the Sampler take a Fitter as input, I would make the Sampler an independent class.
Maybe it can have a constructor Sampler.from_fitter which Constructs a Sampler from a Fitter object, but I wouldn't make it required.
There was a problem hiding this comment.
Also, we should probably copy over the structure of the minimizers. I.e. with a factory, a list of available samplers (for now just DREAM) and separate backend classes.
There was a problem hiding this comment.
Think of the cost: Sampler reaches into Fitter internals (_precompute_reshaping, _fit_function_wrapper, and assigning fitter.fit_function to swap the wrapper.
Making it standalone means extracting that machinery into a shared base, which will affect the basic functionality of the fitter. I woudn't want to do it unless really, really needed.
I will reword the docstring/notebook to say no prior fit is needed as a solution.
There was a problem hiding this comment.
I understand why you're reluctant, this is the easy solution, but it is also the messy solution which will compound our technical debt. By not doing this now, you're just postponing this task and making it harder in the future. It also means that WHEN we do this task in the future, it is gonna be an API-breaking change, whereas by adding it correctly from the get-go isn't API-breaking.
There was a problem hiding this comment.
Following our discussion on Slack, I will reluctantly relent this issue to a (soon) future PR.
There is exactly one sampler (DREAM) and one backend. |
|
|
I added the plan for refactoring in #280 This follows discussion in Slack. I would prefer addressing this in a separate PR. |
damskii9992
left a comment
There was a problem hiding this comment.
Answers to comments, will continue with my review.
| "extended_draws = extended_results.draws\n", | ||
| "\n", | ||
| "print(f'Short chain: {short_draws.shape[0]} retained draws')\n", | ||
| "print(f'Extended chain: {extended_draws.shape[0]} retained draws')" |
There was a problem hiding this comment.
Okay, but this still doesn't explain everything.
So the "Full" chain included the burn-in but not the thinned-out, since the first full chain is 5000 steps, but the extended is 2600 larger, which is more than half of the 5000 asked for. You write that DREAM rounds up to nearest whole generation, so I guess that must be 200 steps? Is that the default or is that a parameter we can change? It is weird that it runs longer than you ask it to.
| except Exception: | ||
| return None |
There was a problem hiding this comment.
You should AT LEAST report a warning then, maybe using your fancy-schmancy new logging module? The user should be notified that the saved state lost its identification fingerprint.
| from .minimizers.minimizer_bumps import save_sampler_state | ||
|
|
||
| save_sampler_state(state, str(path)) |
There was a problem hiding this comment.
So you don't import a save function now, but you still have both a function and a method on the class, why not consolidate them? They're quite small functions/methods both.
| Construct directly with a configured ``Fitter`` (or ``MultiFitter``) whose | ||
| minimizer has been switched to ``AvailableMinimizers.Bumps``. |
There was a problem hiding this comment.
I understand why you're reluctant, this is the easy solution, but it is also the messy solution which will compound our technical debt. By not doing this now, you're just postponing this task and making it harder in the future. It also means that WHEN we do this task in the future, it is gonna be an API-breaking change, whereas by adding it correctly from the get-go isn't API-breaking.
| Construct directly with a configured ``Fitter`` (or ``MultiFitter``) whose | ||
| minimizer has been switched to ``AvailableMinimizers.Bumps``. |
There was a problem hiding this comment.
Following our discussion on Slack, I will reluctantly relent this issue to a (soon) future PR.
|
Issues from the most recent round addressed. |
damskii9992
left a comment
There was a problem hiding this comment.
I did not look closely at the tests, neither did I pay too much attention to the structure since we plan to change that, but here are the last comments I found after going through the PR again.
| from .utils import FitError | ||
| from .utils import FitResults | ||
|
|
||
| __all__ = [MinimizerBase, Bumps, DFO, LMFit, FitError, FitResults] |
There was a problem hiding this comment.
I recently learned that all the __all__ attributed does, is change what gets imported when you do star imports, i.e.:
from easyscience.fitters.minimizers import *The default for that is to import all modules, so simply adding all modules to an __all__ is completely redundant . . .
| burn: int = 2000, | ||
| thin: int = 10, | ||
| population: int | None = None, | ||
| resume_state: Any | None = None, |
There was a problem hiding this comment.
I assume this shouldn't be typehinted Any?
| # Recover the original scale factor from the state's Npop. | ||
| # Try ceil-division first, then fall back to floor. | ||
| recovered = int(math.ceil(resume_state.Npop / n_params)) | ||
| if int(math.ceil(recovered * n_params)) != resume_state.Npop: |
There was a problem hiding this comment.
You KNOW that recovered and n_params are integers, hence you don't need the int and math.ceil calls here.
Also, what are you even trying to do here? You get the floor rounding always, UNLESS resume_state.Npop / n_params is an integer-like float, in which case the ceil function does the same as int?
This section kinda seems redundant.
|
|
||
| return _outer(self) | ||
|
|
||
| def mcmc_sample( |
There was a problem hiding this comment.
This method is getting stupidly long now with the resume functionality, please try to factor out some functionality into helper methods to ease readability.
| # from .multi_fitter import MultiFitter # noqa: F401, E402 | ||
|
|
||
| all = [AvailableMinimizers, Fitter, FitResults] | ||
| all = [AvailableMinimizers, Fitter, FitResults, Sampler, SamplingResults] |
There was a problem hiding this comment.
As per my other comment.
There was a problem hiding this comment.
There is currently a lot of BUMPS-specific code in this file, which makes sense since that is the only sampler we currently have, but bear in mind that WHEN we have another sampler, there should be NO bumps-specific code or syntax in this file, those should all live in the minimizer_bumps.py file.
All syntax here should be sampler-independent.
| Calling ``sample()`` on a sampler that already holds a chain starts a | ||
| **fresh** chain — the previous state and results are replaced. Use | ||
| ``extend()`` to continue an existing chain. |
There was a problem hiding this comment.
Maybe add a warning if a user tries to run sample() when there already is a result, something like:
"Warning, replacing previous chain with fresh chain, if you meant to continue the chain, please use the extend() method"
Or something along those lines.
| warnings.warn( | ||
| 'The data bound to this Sampler does not match the data ' | ||
| 'fingerprint stored with the saved chain. Extending a ' | ||
| 'chain against different data is undefined behaviour.', | ||
| UserWarning, | ||
| stacklevel=2, |
There was a problem hiding this comment.
You should probably use your own logging system here.
| # NOTE: The Bayesian MCMC tests formerly here (TestMultiFitterMcmcSample) were | ||
| # moved to test_sampler.py and adapted to ``Sampler(f, ...)``. |
There was a problem hiding this comment.
We probably don't need this comment lol.
There was a problem hiding this comment.
Most of these tests belong in unit tests. Remember, the unit tests folder should mimic the src folder. So sampler.py should have a mirror test_samlper.py in unit tests.
Integration tests are only for full runs, like running an actual sampling and checking the resulting R value or convergence or the value of some parameter or some such.
It would also be the place to have a full save/load/resume roundtrip test.
Support for resuming and extending Bayesian MCMC chains using the BUMPS DREAM sampler.
Added ability to save and reload sampler state and continue chains across sessions.
Resume and extend MCMC chains:
resume_stateparameter tomcmc_sampleandsamplemethods in bothfitter.pyandminimizer_bumps.py, allowing users to continue sampling from a saved chain state. This includes documentation of the "ring-buffer contract" for extending chains without losing samples.Population/chains argument resolution:
_resolve_population_aliasmethod to handle bothchainsandpopulationarguments, raising an error if both are set with different values, improving API clarity and preventing silent bugs.