sampling/wang_landau/

wang_landau.rs

use {
    crate::{traits::*, *},
    num_traits::{identities::*, ops::wrapping::*, Bounded},
    rand::Rng,
    std::{io::Write, marker::PhantomData, num::*},
};

#[cfg(feature = "serde_support")]
use serde::{Deserialize, Serialize};

/// # The 1/t Wang Landau approach comes from this paper
/// > R. E. Belardinelli and V. D. Pereyra,
/// > “Fast algorithm to calculate density of states,”
/// > Phys. Rev. E **75**: 046701 (2007), DOI [10.1103/PhysRevE.75.046701](https://doi.org/10.1103/PhysRevE.75.046701)
///
/// * The original Wang Landau algorithm comes from this paper
/// > F. Wang and D. P. Landau,
/// > “Efficient, multiple-range random walk algorithm to calculate the density of states,”
/// > Phys. Rev. Lett. **86**, 2050–2053 (2001), DOI [10.1103/PhysRevLett.86.2050](https://doi.org/10.1103/PhysRevLett.86.2050)
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde_support", derive(Serialize, Deserialize))]
pub struct WangLandau1T<Hist, Rng, Ensemble, S, Res, Energy> {
    pub(crate) ensemble: Ensemble,
    pub(crate) rng: Rng,
    pub(crate) marker_res: PhantomData<Res>,
    pub(crate) steps: Vec<S>,
    mode: WangLandauMode,
    pub(crate) log_density: Vec<f64>,
    pub(crate) log_f: f64,
    pub(crate) log_f_threshold: f64,
    pub(crate) step_size: usize,
    step_count: usize,
    accepted_steps_total: usize,
    recected_steps_total: usize,
    accepted_steps_current: usize,
    recected_steps_current: usize,
    pub(crate) old_bin: usize,
    pub(crate) hist: Hist,
    pub(crate) old_energy: Option<Energy>,
    check_refine_every: usize,
}

impl<Hist, R, E, S, Res, Energy> GlueAble<Hist> for WangLandau1T<Hist, R, E, S, Res, Energy>
where
    Hist: Clone,
{
    fn push_glue_entry_ignoring(&self, job: &mut GlueJob<Hist>, ignore_idx: &[usize]) {
        if !ignore_idx.contains(&0) {
            let sim_progress = SimProgress::LogF(self.log_f);
            let rejected = self.total_steps_rejected() as u64;
            let accepted = self.total_steps_accepted() as u64;

            let stats = IntervalSimStats {
                sim_progress,
                interval_sim_type: SimulationType::WangLandau1T,
                rejected_steps: rejected,
                accepted_steps: accepted,
                replica_exchanges: None,
                proposed_replica_exchanges: None,
                merged_over_walkers: NonZeroUsize::new(1).unwrap(),
            };

            let glue_entry = GlueEntry {
                hist: self.hist.clone(),
                prob: self.log_density.clone(),
                log_base: LogBase::BaseE,
                interval_stats: stats,
            };
            job.collection.push(glue_entry);
        }
    }
}

impl<Hist, Rng, Ensemble, S, Res, Energy> WangLandau1T<Hist, Rng, Ensemble, S, Res, Energy> {
    /// Returns internal ensemble, histogram and Rng
    pub fn into_inner(self) -> (Ensemble, Hist, Rng) {
        (self.ensemble, self.hist, self.rng)
    }
}

impl<Hist, R, E, S, Res, Energy> WangLandau for WangLandau1T<Hist, R, E, S, Res, Energy> {
    #[inline(always)]
    fn log_f(&self) -> f64 {
        self.log_f
    }

    #[inline(always)]
    fn log_f_threshold(&self) -> f64 {
        self.log_f_threshold
    }

    fn set_log_f_threshold(&mut self, log_f_threshold: f64) -> Result<f64, WangLandauErrors> {
        if !log_f_threshold.is_finite() || log_f_threshold.is_sign_negative() {
            return Err(WangLandauErrors::InvalidLogFThreshold);
        }
        let old_threshold = self.log_f_threshold;
        self.log_f_threshold = log_f_threshold;
        Ok(old_threshold)
    }

    #[inline(always)]
    fn log_density(&self) -> &Vec<f64> {
        &self.log_density
    }

    fn write_log<W: Write>(&self, mut writer: W) -> Result<(), std::io::Error> {
        writeln!(writer,
            "#Acceptance prob_total: {}\n#Acceptance prob current: {}\n#total_steps: {}\n#log_f: {:e}\n#Current_mode {:?}",
            self.fraction_accepted_total(),
            self.fraction_accepted_current(),
            self.step_counter(),
            self.log_f(),
            self.mode
        )?;
        writeln!(
            writer,
            "#total_steps_accepted: {}\n#total_steps_rejected: {}\n#current_accepted_steps: {}\n#current_rejected_steps: {}",
            self.accepted_steps_total,
            self.recected_steps_total,
            self.accepted_steps_current,
            self.recected_steps_current
        )
    }

    #[inline(always)]
    fn mode(&self) -> WangLandauMode {
        self.mode
    }

    #[inline(always)]
    fn step_counter(&self) -> usize {
        self.step_count
    }

    #[inline(always)]
    fn total_steps_rejected(&self) -> usize {
        self.recected_steps_total
    }

    #[inline(always)]
    fn total_steps_accepted(&self) -> usize {
        self.accepted_steps_total
    }
}

impl<Hist, R, E, S, Res, Energy> WangLandau1T<Hist, R, E, S, Res, Energy>
where
    Hist: Histogram + HistogramVal<Energy>,
{
    /// # Check if `self` is initialized
    /// * if this returns true, you can begin the WangLandau simulation
    /// * otherwise call one of the `self.init*` methods
    pub fn is_initialized(&self) -> bool {
        match &self.old_energy {
            None => false,
            Some(e) => self.hist.is_inside(e),
        }
    }
}

/// Possible errors when setting initial guess
#[derive(Clone, Copy, Debug)]
pub enum SetInitialError {
    /// # Dimensions do not match!
    /// The length of the initial guess and the amount of bins have to be the same
    DimensionError,
    /// All values inside the initial guess have to be finite
    NonFiniteEncountered,
    /// log_f has to fullfill 0.0 < log_f < 10.0
    InvalidLogF,
}

impl<Hist, R, E, S, Res, Energy> WangLandauEnsemble<E>
    for WangLandau1T<Hist, R, E, S, Res, Energy>
{
    #[inline(always)]
    fn ensemble(&self) -> &E {
        &self.ensemble
    }

    #[inline(always)]
    unsafe fn ensemble_mut(&mut self) -> &mut E {
        &mut self.ensemble
    }
}

impl<Hist, R, E, S, Res, Energy> WangLandauEnergy<Energy>
    for WangLandau1T<Hist, R, E, S, Res, Energy>
{
    #[inline(always)]
    fn energy(&self) -> Option<&Energy> {
        self.old_energy.as_ref()
    }
}

impl<Hist, R, E, S, Res, Energy> WangLandauHist<Hist> for WangLandau1T<Hist, R, E, S, Res, Energy> {
    #[inline(always)]
    fn hist(&self) -> &Hist {
        &self.hist
    }
}

impl<Hist, R, E, S, Res, Energy> WangLandau1T<Hist, R, E, S, Res, Energy> {
    /// # Acceptance rate
    /// Fraction of performed wang landau steps, that were accepted
    fn fraction_accepted_total(&self) -> f64 {
        let sum = self.accepted_steps_total + self.recected_steps_total;
        self.accepted_steps_total as f64 / sum as f64
    }

    /// # Acceptance rate since last Refinement
    /// Fraction of performed wang landau steps since
    /// the last time, the factor f was refined, that were accepted
    fn fraction_accepted_current(&self) -> f64 {
        let total = self.accepted_steps_current + self.recected_steps_current;
        if total == 0 {
            f64::NAN
        } else {
            self.accepted_steps_current as f64 / total as f64
        }
    }

    /// # Set the initial guess for the non-normalized probability estimate
    /// * `new_guess` your new guess for the probability estimate. Its length has to equal the number of bins of the internal histogram
    ///   which is the same as the length of the old estimate which you can get by calling [log_density](Self::log_density). All contained values have
    ///   to be finite
    /// * `new_log_f`: Which log_f to start at? 0.0 < log_f <= 10.0 has to be true.
    ///   If you don't know what's best I recommend starting with log_f=1.0, the better your probability estimate is, the smaller this value can be
    /// # Note
    /// This will reset the calculation. Meaning you will have to call one of the initializing functions like `init_greedy_heuristic`again
    /// and all internal counters are reset to 0
    pub fn set_initial_probability_guess(
        mut self,
        new_guess: Vec<f64>,
        new_log_f: f64,
    ) -> Result<Self, SetInitialError>
    where
        Hist: Histogram,
    {
        if 0.0 >= new_log_f || new_log_f > 10.0 {
            Err(SetInitialError::InvalidLogF)
        } else if new_guess.len() != self.log_density.len() {
            Err(SetInitialError::DimensionError)
        } else if new_guess.iter().any(|val| !val.is_finite()) {
            Err(SetInitialError::NonFiniteEncountered)
        } else {
            self.log_density = new_guess;
            self.log_f = new_log_f;
            self.step_count = 0;
            self.accepted_steps_current = 0;
            self.accepted_steps_total = 0;
            self.recected_steps_current = 0;
            self.recected_steps_total = 0;
            self.mode = WangLandauMode::RefineOriginal;
            self.hist.reset();
            self.old_energy = None;
            self.old_bin = usize::MAX;
            Ok(self)
        }
    }
}

impl<Hist, R, E, S, Res, Energy> WangLandau1T<Hist, R, E, S, Res, Energy>
where
    R: Rng,
    E: MarkovChain<S, Res>,
    Energy: Clone,
    Hist: Histogram + HistogramVal<Energy>,
{
    /// # Create a new WangLandau simulation
    /// **IMPORTANT** You have to call one of the `init*` functions,
    /// to create a valid state, before you can start the simulation
    /// ## Parameter
    /// * `log_f_threshold`: how small should the ln(f) (see paper) become
    ///   until the simulation is finished?
    /// * `ensemble`: The ensemble to explore.
    ///   Current state of ensemble will be used as inital condition for the `init*` functions
    /// * `step_size`: The markov steps will be performed with this step size, e.g.,
    ///   `ensemble.m_steps(step_size)`
    /// * `histogram`: Provides the binning. You can either use one of the already implemented
    ///   histograms, like `HistU32Fast`, `HistU32`, `HistF64` etc. or implement your own by
    ///   implementing the traits `Histogram + HistogramVal<Energy>` yourself
    /// * `check_refine_every`: how often to check, if every bin in the histogram was hit.
    ///   Needs to be at least 1. Good values depend on the problem at hand, but if you are
    ///   unsure, you can start with a value like 1000
    pub fn new(
        log_f_threshold: f64,
        ensemble: E,
        rng: R,
        step_size: usize,
        histogram: Hist,
        check_refine_every: usize,
    ) -> Result<Self, WangLandauErrors> {
        if !log_f_threshold.is_finite() || log_f_threshold.is_sign_negative() {
            return Err(WangLandauErrors::InvalidLogFThreshold);
        } else if check_refine_every == 0 {
            return Err(WangLandauErrors::CheckRefineEvery0);
        }
        let log_density = vec![0.0; histogram.bin_count()];
        let steps = Vec::with_capacity(step_size);

        Ok(Self {
            ensemble,
            step_count: 0,
            step_size,
            hist: histogram,
            rng,
            marker_res: PhantomData::<Res>,
            log_f: 1.0,
            log_density,
            log_f_threshold,
            mode: WangLandauMode::RefineOriginal,
            recected_steps_current: 0,
            recected_steps_total: 0,
            accepted_steps_current: 0,
            accepted_steps_total: 0,
            old_bin: usize::MAX,
            old_energy: None,
            check_refine_every,
            steps,
        })
    }

    fn init<F>(&mut self, energy_fn: F, step_limit: Option<u64>) -> Result<(), WangLandauErrors>
    where
        F: Fn(&mut E) -> Option<Energy>,
    {
        self.old_energy = energy_fn(&mut self.ensemble);
        if self.old_energy.is_some() {
            return Ok(());
        }

        match step_limit {
            None => loop {
                self.ensemble.m_steps_quiet(self.step_size);
                self.old_energy = energy_fn(&mut self.ensemble);

                if self.old_energy.is_some() {
                    self.count_accepted();
                    return Ok(());
                }
                self.count_rejected();
            },
            Some(limit) => {
                for _ in 0..limit {
                    self.ensemble.m_steps_quiet(self.step_size);
                    self.old_energy = energy_fn(&mut self.ensemble);

                    if self.old_energy.is_some() {
                        self.count_accepted();
                        return Ok(());
                    }
                    self.count_rejected();
                }
                Err(WangLandauErrors::InitFailed)
            }
        }
    }

    fn greedy_helper<F, H, J>(&mut self, old_distance: &mut J, energy_fn: F, distance_fn: H)
    where
        F: Fn(&mut E) -> Option<Energy> + Copy,
        H: Fn(&Hist, &Energy) -> J,
        J: PartialOrd,
    {
        self.ensemble.m_steps(self.step_size, &mut self.steps);

        if let Some(energy) = energy_fn(&mut self.ensemble) {
            let distance = distance_fn(&self.hist, &energy);
            if distance <= *old_distance {
                self.old_energy = Some(energy);
                *old_distance = distance;
                self.count_accepted();

                return;
            }
        }

        self.count_rejected();
        self.ensemble.undo_steps_quiet(&self.steps);
    }

    /// # Find a valid starting Point
    /// * if the ensemble is already at a valid starting point,
    ///   the ensemble is left unchanged (as long as your energy calculation does not change the ensemble)
    /// * Uses a greedy heuristic. Performs markov steps. If that brought us closer to the target interval,
    ///   the step is accepted. Otherwise it is rejected
    /// # Parameter
    /// * `step_limit`: Some(val) -> val is max number of steps tried, if no valid state is found, it will return an Error. None -> will loop until either
    ///   a valid state is found or forever
    /// * `energy_fn` function calculating `Some(energy)` of the system
    ///   or rather the Parameter of which you wish to obtain the probability distribution.
    ///   Has to be the same function as used for the wang landau simulation later.
    ///   If there are any states, for which the calculation is invalid, `None` should be returned
    /// * steps resulting in ensembles for which `energy_fn(&mut ensemble)` is `None`
    ///   will always be rejected
    pub fn init_greedy_heuristic<F>(
        &mut self,
        energy_fn: F,
        step_limit: Option<u64>,
    ) -> Result<(), WangLandauErrors>
    where
        F: Fn(&mut E) -> Option<Energy>,
    {
        self.init(&energy_fn, step_limit)?;
        let mut old_distance = self.hist.distance(self.old_energy_ref());
        let mut step_count = 0;
        while old_distance != 0.0 {
            self.greedy_helper(&mut old_distance, &energy_fn, |hist, energy| {
                hist.distance(energy)
            });
            if let Some(limit) = step_limit {
                if limit == step_count {
                    return Err(WangLandauErrors::InitFailed);
                }
                step_count += 1;
            }
        }
        self.end_init();
        Ok(())
    }

    /// # Find a valid starting Point
    /// * if the ensemble is already at a valid starting point,
    ///   the ensemble is left unchanged (as long as your energy calculation does not change the ensemble)
    /// * Uses overlapping intervals. Accepts a step, if the resulting ensemble is in the same interval as before,
    ///   or it is in an interval closer to the target interval
    /// * Take a look at the [`HistogramIntervalDistance` trait](`crate::HistogramIntervalDistance`)
    /// # Parameter
    /// * `step_limit`: Some(val) -> val is max number of steps tried, if no valid state is found, it will return an Error. None -> will loop until either
    ///   a valid state is found or forever
    /// * `energy_fn` function calculating `Some(energy)` of the system
    ///   or rather the Parameter of which you wish to obtain the probability distribution.
    ///   Has to be the same function as used for the wang landau simulation later.
    ///   If there are any states, for which the calculation is invalid, `None` should be returned
    /// * steps resulting in ensembles for which `energy_fn(&mut ensemble)` is `None`
    ///   will always be rejected
    pub fn init_interval_heuristik<F>(
        &mut self,
        overlap: NonZeroUsize,
        energy_fn: F,
        step_limit: Option<u64>,
    ) -> Result<(), WangLandauErrors>
    where
        F: Fn(&mut E) -> Option<Energy>,
        Hist: HistogramIntervalDistance<Energy>,
    {
        self.init(&energy_fn, step_limit)?;
        let mut old_dist = self
            .hist
            .interval_distance_overlap(self.old_energy_ref(), overlap);

        let dist = |h: &Hist, val: &Energy| h.interval_distance_overlap(val, overlap);
        let mut step_count = 0;
        while old_dist != 0 {
            self.greedy_helper(&mut old_dist, &energy_fn, dist);
            if let Some(limit) = step_limit {
                if limit == step_count {
                    return Err(WangLandauErrors::InitFailed);
                }
                step_count += 1;
            }
        }
        self.end_init();
        Ok(())
    }

    /// # Find a valid starting Point
    /// * if the ensemble is already at a valid starting point,
    ///   the ensemble is left unchanged (as long as your energy calculation does not change the ensemble)
    /// * `overlap` - see [`HistogramIntervalDistance` trait](`crate::HistogramIntervalDistance`)
    ///   Should be greater than 0 and smaller than the number of bins in your histogram. E.g. `overlap = 3` if you have 200 bins
    /// * `mid` - should be something like `128u8`, `0i8` or `0i16`. It is very unlikely that using a type with more than 16 bit makes sense for mid
    /// * `step_limit`: Some(val) -> val is max number of steps tried, if no valid state is found, it will return an Error. None -> will loop until either
    ///   a valid state is found or forever
    /// * alternates between greedy and interval heuristic every time a wrapping counter passes `mid` or `U::min_value()`
    /// * I recommend using this heuristic, if you do not know which one to use
    /// # Parameter
    /// * `energy_fn` function calculating `Some(energy)` of the system
    ///   or rather the Parameter of which you wish to obtain the probability distribution.
    ///   Has to be the same function as used for the wang landau simulation later.
    ///   If there are any states, for which the calculation is invalid, `None` should be returned
    /// * steps resulting in ensembles for which `energy_fn(&mut ensemble)` is `None`
    ///   will always be rejected
    pub fn init_mixed_heuristik<F, U>(
        &mut self,
        overlap: NonZeroUsize,
        mid: U,
        energy_fn: F,
        step_limit: Option<u64>,
    ) -> Result<(), WangLandauErrors>
    where
        F: Fn(&mut E) -> Option<Energy>,
        Hist: HistogramIntervalDistance<Energy>,
        U: One + Bounded + WrappingAdd + Eq + PartialOrd,
    {
        self.init(&energy_fn, step_limit)?;
        if self.hist.is_inside(self.old_energy_ref()) {
            self.end_init();
            return Ok(());
        }

        let mut old_dist = f64::INFINITY;
        let mut old_dist_interval = usize::MAX;
        let mut counter: U = U::min_value();
        let min_val = U::min_value();
        let one = U::one();
        let dist_interval = |h: &Hist, val: &Energy| h.interval_distance_overlap(val, overlap);
        let mut step_count = 0;
        loop {
            if counter == min_val {
                let current_energy = self.old_energy_ref();
                old_dist = self.hist.distance(current_energy);
            } else if counter == mid {
                let current_energy = self.old_energy_ref();
                old_dist_interval = dist_interval(&self.hist, current_energy);
            }
            if counter < mid {
                self.greedy_helper(&mut old_dist, &energy_fn, |hist, val| hist.distance(val));
                if old_dist == 0.0 {
                    break;
                }
            } else {
                self.greedy_helper(&mut old_dist_interval, &energy_fn, dist_interval);
                if old_dist_interval == 0 {
                    break;
                }
            }
            counter = counter.wrapping_add(&one);
            if let Some(limit) = step_limit {
                if limit == step_count {
                    return Err(WangLandauErrors::InitFailed);
                }
                step_count += 1;
            }
        }
        self.end_init();
        Ok(())
    }

    fn end_init(&mut self) {
        self.old_bin = self
            .hist
            .get_bin_index(self.old_energy_ref())
            .expect("Error in heuristic - old bin invalid");
    }

    fn old_energy_clone(&self) -> Energy {
        self.old_energy_ref().clone()
    }

    fn old_energy_ref(&self) -> &Energy {
        self.old_energy.as_ref().unwrap()
    }

    fn count_accepted(&mut self) {
        self.ensemble.steps_accepted(&self.steps);
        self.accepted_steps_current += 1;
        self.accepted_steps_total += 1;
    }

    fn count_rejected(&mut self) {
        self.ensemble.steps_rejected(&self.steps);
        self.recected_steps_current += 1;
        self.recected_steps_total += 1;
    }

    fn check_refine(&mut self) {
        match self.mode {
            WangLandauMode::Refine1T => {
                self.log_f = self.log_f_1_t();
            }
            WangLandauMode::RefineOriginal => {
                if self.step_count % self.check_refine_every == 0 && !self.hist.any_bin_zero() {
                    self.recected_steps_current = 0;
                    self.accepted_steps_current = 0;
                    let ref_1_t = self.log_f_1_t();
                    self.log_f *= 0.5;
                    if self.log_f < ref_1_t {
                        self.log_f = ref_1_t;
                        self.mode = WangLandauMode::Refine1T;
                    }
                    self.hist.reset();
                }
            }
        }
    }

    fn wl_step_helper(&mut self, energy: Option<Energy>) {
        let current_energy = match energy {
            Some(energy) => energy,
            None => {
                self.count_rejected();
                self.hist.increment_index(self.old_bin).unwrap();
                self.log_density[self.old_bin] += self.log_f;
                self.ensemble.undo_steps_quiet(&self.steps);
                return;
            }
        };

        match self.hist.get_bin_index(&current_energy) {
            Ok(current_bin) => {
                let accept_prob = self.metropolis_acception_prob(current_bin);

                if self.rng.random::<f64>() > accept_prob {
                    // reject step
                    self.count_rejected();
                    self.ensemble.undo_steps_quiet(&self.steps);
                } else {
                    // accept step
                    self.count_accepted();

                    self.old_energy = Some(current_energy);
                    self.old_bin = current_bin;
                }
            }
            _ => {
                // invalid step -> reject
                self.count_rejected();
                self.ensemble.undo_steps_quiet(&self.steps);
            }
        };

        self.hist.increment_index(self.old_bin).unwrap();
        self.log_density[self.old_bin] += self.log_f;
    }

    /// # Wang Landau Step
    /// * performs a single Wang Landau step
    /// # Parameter
    /// * `energy_fn` function calculating `Some(energy)` of the system
    ///   or rather the Parameter of which you wish to obtain the probability distribution.
    ///   If there are any states, for which the calculation is invalid, `None` should be returned
    /// * steps resulting in ensembles for which `energy_fn(&mut ensemble)` is `None`
    ///   will always be rejected
    /// # Important
    /// * You have to call one of the `self.init*` functions before calling this one -
    ///   **will panic otherwise**
    pub fn wang_landau_step<F>(&mut self, energy_fn: F)
    where
        F: Fn(&E) -> Option<Energy>,
    {
        unsafe { self.wang_landau_step_unsafe(|e| energy_fn(e)) }
    }

    /// # Wang Landau Step
    /// * if possible, use `self.wang_landau_step()` instead - it is safer
    /// * performs a single Wang Landau step
    /// # Parameter
    /// * `energy_fn` function calculating `Some(energy)` of the system
    ///   or rather the Parameter of which you wish to obtain the probability distribution.
    ///   If there are any states, for which the calculation is invalid, `None` should be returned
    /// * steps resulting in ensembles for which `energy_fn(&mut ensemble)` is `None`
    ///   will always be rejected
    /// # Safety
    /// * You have to call one of the `self.init*` functions before calling this one -
    ///   **will panic otherwise**
    /// * unsafe, because you have to make sure, that the `energy_fn` function
    ///   does not change the state of the ensemble in such a way, that the result of
    ///   `energy_fn` changes when called again. Maybe do cleanup at the beginning of the energy function?
    pub unsafe fn wang_landau_step_unsafe<F>(&mut self, mut energy_fn: F)
    where
        F: FnMut(&mut E) -> Option<Energy>,
    {
        debug_assert!(
            self.old_energy.is_some(),
            "Error - self.old_energy invalid - Did you forget to call one of the `self.init*` members for initialization?"
        );

        self.step_count += 1;

        self.ensemble.m_steps(self.step_size, &mut self.steps);

        self.check_refine();
        let current_energy = energy_fn(&mut self.ensemble);

        self.wl_step_helper(current_energy);
    }

    /// # Wang Landau Step
    /// * performs a single Wang Landau step
    /// # Parameter
    /// * `energy_fn` function calculating the energy of the system **on the fly**
    /// * **steps resulting in invalid ensembles are not allowed!**
    /// # Important
    /// * You have to call one of the `self.init*` functions before calling this one -
    ///   **will panic otherwise**
    pub fn wang_landau_step_acc<F>(&mut self, energy_fn: F)
    where
        F: FnMut(&E, &S, &mut Energy),
    {
        debug_assert!(
            self.old_energy.is_some(),
            "Error - self.old_energy invalid - Did you forget to call one of the `self.init*` members for initialization?"
        );

        self.step_count += 1;

        let mut new_energy = self.old_energy_clone();

        self.ensemble
            .m_steps_acc(self.step_size, &mut self.steps, &mut new_energy, energy_fn);

        self.check_refine();

        self.wl_step_helper(Some(new_energy));
    }

    /// # Wang Landau
    /// * perform Wang Landau simulation
    /// * calls `self.wang_landau_step(energy_fn, valid_ensemble)` until `self.is_finished()`
    pub fn wang_landau_convergence<F>(&mut self, energy_fn: F)
    where
        F: Fn(&E) -> Option<Energy>,
    {
        while !self.is_finished() {
            self.wang_landau_step(&energy_fn);
        }
    }

    /// # Wang Landau - efficient energy calculation
    /// * perform Wang Landau simulation
    /// * calls `self.wang_landau_step_acc(energy_fn, valid_ensemble)` until `self.is_finished()`
    pub fn wang_landau_convergence_acc<F>(&mut self, mut energy_fn: F)
    where
        F: FnMut(&E, &S, &mut Energy),
    {
        while !self.is_finished() {
            self.wang_landau_step_acc(&mut energy_fn);
        }
    }

    /// # Wang Landau
    /// * if possible, use `self.wang_landau_convergence()` instead - it is safer
    /// * perform Wang Landau simulation
    /// * calls `self.wang_landau_step_unsafe(energy_fn, valid_ensemble)` until `self.is_finished()`
    /// # Safety
    /// * You have mutable access to your ensemble, which is why this function is unsafe.
    ///   If you do anything, which changes the future outcome of the energy function, the results will be wrong!
    ///   I use the unsafe keyword here to force the user to acknowledge that.
    pub unsafe fn wang_landau_convergence_unsafe<F>(&mut self, mut energy_fn: F)
    where
        F: FnMut(&mut E) -> Option<Energy>,
    {
        while !self.is_finished() {
            self.wang_landau_step_unsafe(&mut energy_fn);
        }
    }

    /// # Wang Landau
    /// * perform Wang Landau simulation
    /// * calls `self.wang_landau_step(energy_fn)` until `self.is_finished()`
    ///   or `condition(&self)` is false
    pub fn wang_landau_while<F, W>(&mut self, energy_fn: F, mut condition: W)
    where
        F: Fn(&E) -> Option<Energy>,
        W: FnMut(&Self) -> bool,
    {
        while !self.is_finished() && condition(self) {
            self.wang_landau_step(&energy_fn);
        }
    }

    /// # Wang Landau
    /// * perform Wang Landau simulation
    /// * calls `self.wang_landau_step(energy_fn)` until `self.is_finished()`
    ///   or `condition(&self)` is false
    pub fn wang_landau_while_acc<F, W>(&mut self, mut energy_fn: F, mut condition: W)
    where
        F: FnMut(&E, &S, &mut Energy),
        W: FnMut(&Self) -> bool,
    {
        while !self.is_finished() && condition(self) {
            self.wang_landau_step_acc(&mut energy_fn);
        }
    }

    /// # Wang Landau
    /// * if possible, use `self.wang_landau_while()` instead - it is safer
    /// * perform Wang Landau simulation
    /// * calls `self.wang_landau_step(energy_fn)` until `self.is_finished()`
    ///   or `condition(&self)` is false
    /// # Safety
    /// * You have mutable access to your ensemble, which is why this function is unsafe.
    ///   If you do anything, which changes the future outcome of the energy function, the results will be wrong!
    ///   I use the unsafe keyword here to force the user to acknowledge that
    pub unsafe fn wang_landau_while_unsafe<F, W>(&mut self, mut energy_fn: F, mut condition: W)
    where
        F: FnMut(&mut E) -> Option<Energy>,
        W: FnMut(&Self) -> bool,
    {
        while !self.is_finished() && condition(self) {
            self.wang_landau_step_unsafe(&mut energy_fn);
        }
    }

    /// **panics** if index is invalid
    #[inline(always)]
    fn metropolis_acception_prob(&self, new_bin: usize) -> f64 {
        (self.log_density[self.old_bin] - self.log_density[new_bin]).exp()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::examples::coin_flips::*;
    use rand::SeedableRng;
    use rand_pcg::Pcg64Mcg;
    #[test]
    #[cfg_attr(miri, ignore)]
    fn wl_simulations_equal() {
        let mut rng = Pcg64Mcg::seed_from_u64(2239790);
        let ensemble = CoinFlipSequence::new(100, Pcg64Mcg::from_rng(&mut rng));
        let histogram = HistogramFast::new_inclusive(0, 100).unwrap();
        let mut wl = WangLandau1T::new(0.0075, ensemble, rng, 1, histogram, 30).unwrap();

        wl.init_mixed_heuristik(
            NonZeroUsize::new(3).unwrap(),
            6400i16,
            |e| Some(e.head_count()),
            None,
        )
        .unwrap();

        let mut wl_backup = wl.clone();
        let start_wl = std::time::Instant::now();
        wl.wang_landau_convergence(|e| Some(e.head_count()));
        let dur_1 = start_wl.elapsed();
        let start_wl_acc = std::time::Instant::now();
        wl_backup.wang_landau_convergence_acc(CoinFlipSequence::update_head_count);
        let dur_2 = start_wl_acc.elapsed();
        println!(
            "WL: {}, WL_ACC: {}, difference: {}",
            dur_1.as_nanos(),
            dur_2.as_nanos(),
            dur_1.as_nanos() - dur_2.as_nanos()
        );

        // assert, that the different methods lead to the same result
        for (&log_value, &log_value_acc) in
            wl.log_density().iter().zip(wl_backup.log_density().iter())
        {
            assert_eq!(log_value, log_value_acc);
        }
    }
}