System.h

Go to the documentation of this file.

 
#ifndef GOOSEEPM_SYSTEM_H
#define GOOSEEPM_SYSTEM_H
 
#include <algorithm>
#include <chrono>
#include <prrng.h>
#include <xtensor/xset_operation.hpp>
#include <xtensor/xsort.hpp>
 
#include "config.h"
#include "version.h"
 
namespace GooseEPM {
 
namespace detail {
 
inline std::string get_namespace()
{
    std::string ret = "GooseEPM";
#ifdef GOOSEEPM_USE_XTENSOR_PYTHON
    return ret + ".";
#else
    return ret + "::";
#endif
}
 
template <class T>
inline std::string shape_to_string(const T& shape)
{
    std::string ret = "[";
    for (size_t i = 0; i < shape.size(); ++i) {
        ret += std::to_string(shape[i]);
        if (i < shape.size() - 1) {
            ret += ", ";
        }
    }
    ret += "]";
    return ret;
}
 
inline void unravel_index(
    ptrdiff_t idx,
    const std::array<ptrdiff_t, 2>& strides,
    std::array<ptrdiff_t, 2>& indices
)
{
    indices[0] = idx / strides[0];
    indices[1] = idx % strides[0];
}
 
template <class T>
inline bool check_distances(const T& distance)
{
    using value_type = typename T::value_type;
    static_assert(std::numeric_limits<value_type>::is_integer, "Distances must be integer");
    static_assert(std::numeric_limits<value_type>::is_signed, "Distances must be signed");
 
    if (distance.dimension() != 1) {
        return false;
    }
 
    value_type lower = xt::amin(distance)();
    value_type upper = xt::amax(distance)() + 1;
    auto d = xt::arange<value_type>(lower, upper);
 
    if (!xt::all(xt::in1d(distance, d))) {
        return false;
    }
 
    if (static_cast<value_type>(distance.size()) != upper - lower) {
        return false;
    }
 
    return true;
}
 
template <class T>
inline T create_distance_lookup(const T& distance)
{
    using value_type = typename T::value_type;
    GOOSEEPM_REQUIRE(detail::check_distances(distance), std::out_of_range);
 
    value_type lower = xt::amin(distance)();
    value_type upper = xt::amax(distance)() + 1;
 
    value_type N = static_cast<value_type>(distance.size());
    T ret = xt::empty<value_type>({2 * N - 1});
 
    for (value_type i = 0; i < upper; ++i) {
        ret(i) = xt::argmax(xt::equal(distance, i))();
    }
    for (value_type i = upper; i < N; ++i) {
        ret(i) = xt::argmax(xt::equal(distance, i - N))();
    }
    for (value_type i = -1; i >= lower; --i) {
        ret.periodic(i) = xt::argmax(xt::equal(distance, i))();
    }
    for (value_type i = lower; i > -N; --i) {
        ret.periodic(i) = xt::argmax(xt::equal(distance, N + i))();
    }
 
    return ret;
}
 
template <size_t Dim, class T, class D>
inline T reorganise_popagator(const T& propagator, const D& distances)
{
    static_assert(Dim > 0 && Dim < 3, "Only 1d/2d systems are supported");
    GOOSEEPM_REQUIRE(propagator.dimension() == Dim, std::out_of_range);
 
    T ret = xt::empty_like(propagator);
 
    if constexpr (Dim == 1) {
        auto d = xt::argsort(distances[0]);
        for (size_t i = 0; i < propagator.size(); ++i) {
            ret(i) = propagator(d(i));
        }
    }
    else if constexpr (Dim == 2) {
        auto dr = xt::argsort(distances[0]);
        auto dc = xt::argsort(distances[1]);
 
        for (size_t i = 0; i < propagator.shape(0); ++i) {
            for (size_t j = 0; j < propagator.shape(1); ++j) {
                ret(i, j) = propagator(dr(i), dc(j));
            }
        }
    }
 
    return ret;
}
 
template <size_t Dim, class T, class D>
inline T popagator_fix_stress(const T& propagator, const D& distances)
{
    static_assert(Dim > 0 && Dim < 3, "Only 1d/2d systems are supported");
    GOOSEEPM_REQUIRE(propagator.dimension() == Dim, std::out_of_range);
 
    T ret = propagator;
 
    if constexpr (Dim == 1) {
        auto i = xt::argwhere(xt::equal(distances[0], 0));
        GOOSEEPM_REQUIRE(i.size() == 1, std::out_of_range);
        for (size_t k = 0; k < 20; ++k) { // loop to increase precision
            ret -= xt::mean(ret);
            ret(i[0][0]) = -1.0;
        }
    }
    else if constexpr (Dim == 2) {
        auto i = xt::argwhere(xt::equal(distances[0], 0));
        auto j = xt::argwhere(xt::equal(distances[1], 0));
        GOOSEEPM_REQUIRE(i.size() == 1, std::out_of_range);
        GOOSEEPM_REQUIRE(j.size() == 1, std::out_of_range);
        for (size_t k = 0; k < 20; ++k) { // loop to increase precision
            ret -= xt::mean(ret);
            ret(i[0][0], j[0][0]) = -1.0;
        }
    }
 
    return reorganise_popagator<Dim>(ret, distances);
}
 
template <size_t Dim, class T, class D>
inline T popagator_fix_strain(const T& propagator, const D& distances)
{
    static_assert(Dim > 0 && Dim < 3, "Only 1d/2d systems are supported");
    GOOSEEPM_REQUIRE(propagator.dimension() == Dim, std::out_of_range);
 
    T ret = propagator;
 
    if constexpr (Dim == 1) {
        auto i = xt::argwhere(xt::equal(distances[0], 0));
        GOOSEEPM_REQUIRE(i.size() == 1, std::out_of_range);
        for (size_t k = 0; k < 20; ++k) { // loop to increase precision
            ret -= xt::mean(ret) + 1.0 / static_cast<double>(ret.size());
            ret(i[0][0]) = -1.0;
        }
    }
    else if constexpr (Dim == 2) {
        auto i = xt::argwhere(xt::equal(distances[0], 0));
        auto j = xt::argwhere(xt::equal(distances[1], 0));
        GOOSEEPM_REQUIRE(i.size() == 1, std::out_of_range);
        GOOSEEPM_REQUIRE(j.size() == 1, std::out_of_range);
        for (size_t k = 0; k < 20; ++k) { // loop to increase precision
            ret -= xt::mean(ret) + 1.0 / static_cast<double>(ret.size());
            ret(i[0][0], j[0][0]) = -1.0;
        }
    }
 
    return reorganise_popagator<Dim>(ret, distances);
}
 
class Default {
public:
    Default()
    {
    }
};
 
class Depinning {
public:
    Depinning()
    {
    }
};
 
class Finite {
public:
    Finite()
    {
    }
};
 
class FixStrain {
public:
    FixStrain()
    {
    }
};
 
class FixStress {
private:
    double m_sigbar; 
 
public:
    FixStress(double sigbar) : m_sigbar(sigbar)
    {
    }
 
    template <class T>
    void set_average_stress(T& sig)
    {
        sig -= xt::mean(sig)() - m_sigbar;
    }
 
    void set_sigmabar(double sigbar)
    {
        m_sigbar = sigbar;
    }
 
    double sigmabar() const
    {
        return m_sigbar;
    }
};
 
class Spring {
private:
    double m_u; 
    double m_k_div_n; 
    double m_kp1_div_k; 
 
public:
    Spring(double k, double u, size_t size)
        : m_u(u), m_k_div_n(k / static_cast<double>(size)), m_kp1_div_k((k + 1.0) / k)
    {
    }
 
    template <class T>
    void stress_match_frame(T& sig, double depsp_f)
    {
        sig -= depsp_f * m_k_div_n;
    }
 
    void frame_match_stress(double delta_sigma)
    {
        m_u += delta_sigma * m_kp1_div_k;
    }
 
    double u() const
    {
        return m_u;
    }
 
    void set_u(double u)
    {
        m_u = u;
    }
};
 
} // namespace detail
 
template <size_t Dimension, class StressFixer, class Temperature, class FailureRules>
class SystemBase {
public:
    static constexpr size_t Dim = Dimension; 
 
protected:
    prrng::pcg32 m_gen; 
    array_type::tensor<double, Dim> m_propagator; 
    std::array<ptrdiff_t, Dim> m_offset; 
    array_type::tensor<double, Dim> m_sig; 
    array_type::tensor<double, Dim> m_sigy; 
    array_type::tensor<double, Dim> m_sigy_mu; 
    array_type::tensor<double, Dim> m_sigy_std; 
    array_type::tensor<double, Dim> m_epsp; 
    array_type::tensor<size_t, Dim> m_nfails; 
    array_type::tensor<bool, Dim> m_unstable; 
    array_type::tensor<double, Dim> m_tau_thermal; 
    double m_t; 
    double m_alpha; 
    StressFixer* m_stress_fixer; 
    double m_temperature; 
    double m_inv_temp; 
    double m_tau; 
    double m_failure_rate; 
    size_t m_failure_idx; 
    double m_failure_x; 
 
protected:
    void initSystem(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_mean,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate,
        double alpha,
        bool random_stress,
        StressFixer* stress_fixer = nullptr
    )
    {
        static_assert(Dim > 0 && Dim < 3, "Only 1d/2d systems are supported");
 
        m_stress_fixer = stress_fixer;
        m_temperature = 0.0;
        m_inv_temp = std::numeric_limits<double>::infinity();
 
        GOOSEEPM_REQUIRE(propagator.dimension() == Dim, std::out_of_range);
        for (size_t i = 0; i < Dim; ++i) {
            GOOSEEPM_REQUIRE(distances[i].size() == propagator.shape(i), std::out_of_range);
            GOOSEEPM_REQUIRE(detail::check_distances(distances[i]), std::out_of_range);
            m_offset[i] = xt::amin(distances[i])();
        }
        GOOSEEPM_REQUIRE(xt::has_shape(sigmay_mean, sigmay_std.shape()), std::out_of_range);
 
        if constexpr (std::is_same<StressFixer, detail::FixStrain>::value) {
            m_propagator = detail::popagator_fix_strain<Dim>(propagator, distances);
        }
        else {
            m_propagator = detail::popagator_fix_stress<Dim>(propagator, distances);
        }
 
        m_gen.seed(seed);
        m_alpha = alpha;
        m_failure_rate = failure_rate;
        m_tau = 1.0 / m_failure_rate;
        m_sigy_mu = sigmay_mean;
        m_sigy_std = sigmay_std;
        m_sig = xt::zeros_like(m_sigy_mu);
        m_epsp = xt::zeros_like(m_sigy_mu);
        m_nfails = xt::zeros<size_t>(m_epsp.shape());
        m_sigy = xt::empty_like(m_sigy_mu);
        m_unstable = xt::empty<bool>(m_epsp.shape());
        m_tau_thermal = xt::empty_like(m_sigy_mu);
 
        for (size_t i = 0; i < m_sigy.size(); ++i) {
            this->drawYieldStress(i);
        }
 
        if (random_stress) {
            this->initRandomStress(propagator, distances);
        }
        else {
            this->stress_updated();
        }
 
        m_t = 0.0;
        m_epsp.fill(0.0);
        m_nfails.fill(0);
        m_failure_idx = m_epsp.size();
        m_failure_x = 0.0;
    }
 
private:
    template <class S, class T>
    void initRandomStress(const S& propagator, const T& distances)
    {
        if constexpr (std::is_same<StressFixer, detail::FixStrain>::value) {
            m_propagator = detail::popagator_fix_stress<Dim>(propagator, distances);
        }
 
        this->initSigmaPropogator(0.1);
 
        if constexpr (std::is_same<StressFixer, detail::FixStrain>::value) {
            m_propagator = detail::popagator_fix_strain<Dim>(propagator, distances);
        }
 
        for (size_t i = 0; i < 20; ++i) {
            size_t n = this->relaxAthermal();
            if (n == 0) {
                return;
            }
 
            if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
                m_stress_fixer->set_u(0.0);
            }
            m_sig -= xt::mean(m_sig)();
            this->stress_updated();
        }
        GOOSEEPM_WARNING("Initialisation failed to converge");
    }
 
    void stress_updated()
    {
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            xt::noalias(m_unstable) = m_sig >= m_sigy;
            if constexpr (std::is_same<Temperature, detail::Finite>::value) {
                xt::noalias(m_tau_thermal) = xt::where(
                    m_unstable,
                    m_tau,
                    xt::exp(xt::pow(m_sigy - m_sig, m_alpha) * m_inv_temp) * m_tau
                );
            }
        }
        else {
            xt::noalias(m_unstable) = xt::abs(m_sig) >= m_sigy;
            if constexpr (std::is_same<Temperature, detail::Finite>::value) {
                xt::noalias(m_tau_thermal) = xt::where(
                    m_unstable,
                    m_tau,
                    xt::exp(xt::pow(m_sigy - xt::abs(m_sig), m_alpha) * m_inv_temp) * m_tau
                );
            }
        }
    }
 
    void stress_updated(size_t idx)
    {
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            m_unstable.flat(idx) = m_sig.flat(idx) >= m_sigy.flat(idx);
        }
        else {
            m_unstable.flat(idx) = std::abs(m_sig.flat(idx)) >= m_sigy.flat(idx);
        }
 
        if constexpr (std::is_same<Temperature, detail::Finite>::value) {
            if (m_unstable.flat(idx)) {
                m_tau_thermal.flat(idx) = m_tau;
            }
            else {
                double E;
                if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
                    E = std::pow(m_sigy.flat(idx) - m_sig.flat(idx), m_alpha);
                }
                else {
                    E = std::pow(m_sigy.flat(idx) - std::abs(m_sig.flat(idx)), m_alpha);
                }
                m_tau_thermal.flat(idx) = std::exp(E * m_inv_temp) * m_tau;
            }
        }
    }
 
    void drawYieldStress(size_t idx)
    {
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            m_sigy.flat(idx) = std::abs(m_gen.normal(m_sigy_mu.flat(idx), m_sigy_std.flat(idx)));
        }
        else {
            m_sigy.flat(idx) = m_gen.normal(m_sigy_mu.flat(idx), m_sigy_std.flat(idx));
            // avoid negative yield stresses. todo: permanent fix with truncated normal distribution
            while (m_sigy.flat(idx) < 0) {
                m_sigy.flat(idx) = m_gen.normal(m_sigy_mu.flat(idx), m_sigy_std.flat(idx));
            }
        }
    }
 
public:
    std::string repr() const
    {
        std::string ret = detail::get_namespace() + "System";
 
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            ret += "Depinning";
        }
 
        if constexpr (std::is_same<Temperature, detail::Finite>::value) {
            ret += "Thermal";
        }
 
        if constexpr (std::is_same<StressFixer, detail::FixStrain>::value) {
            ret += "StrainControl";
        }
        else if constexpr (std::is_same<StressFixer, detail::FixStress>::value) {
            ret += "StressControl";
        }
        else if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
            ret += "SpringLoading";
        }
 
        return ret + std::to_string(Dim) + " " + detail::shape_to_string(m_sig.shape());
    }
 
    const auto& shape() const
    {
        return m_epsp.shape();
    }
 
    auto size() const
    {
        return m_epsp.size();
    }
 
    const array_type::tensor<double, Dim>& propagator() const
    {
        return m_propagator;
    }
 
    array_type::tensor<ptrdiff_t, 1> distances(size_t axis = 0) const
    {
        GOOSEEPM_REQUIRE(axis < Dim, std::out_of_range);
        return m_offset[axis] + xt::arange<ptrdiff_t>(m_propagator.shape(axis));
    }
 
    auto alpha() const
    {
        return m_alpha;
    }
 
    auto failure_rate() const
    {
        return m_failure_rate;
    }
 
    auto failure_idx() const
    {
        return m_failure_idx;
    }
 
    auto failure_x() const
    {
        return m_failure_x;
    }
 
    void set_t(double t)
    {
        m_t = t;
    }
 
    double t() const
    {
        return m_t;
    }
 
    void set_state(uint64_t state)
    {
        m_gen.restore(state);
    }
 
    uint64_t state() const
    {
        return m_gen.state();
    }
 
    void set_epsp(const array_type::tensor<double, Dim>& epsp)
    {
        m_epsp = epsp;
    }
 
    const array_type::tensor<double, Dim>& epsp() const
    {
        return m_epsp;
    }
 
    void set_nfails(const array_type::tensor<size_t, Dim>& nfails)
    {
        m_nfails = nfails;
    }
 
    const array_type::tensor<size_t, Dim>& nfails() const
    {
        return m_nfails;
    }
 
    void set_sigmay(const array_type::tensor<double, Dim>& sigmay)
    {
        m_sigy = sigmay;
        this->stress_updated();
    }
 
    const array_type::tensor<double, Dim>& sigmay() const
    {
        return m_sigy;
    }
 
    void set_sigmay_mean(const array_type::tensor<double, Dim>& value)
    {
        GOOSEEPM_ASSERT(xt::has_shape(value, m_sigy_mu.shape()), std::out_of_range);
        m_sigy_mu = value;
    }
 
    const array_type::tensor<double, Dim>& sigmay_mean() const
    {
        return m_sigy_mu;
    }
 
    void set_sigmay_std(const array_type::tensor<double, Dim>& value)
    {
        GOOSEEPM_ASSERT(xt::has_shape(value, m_sigy_std.shape()), std::out_of_range);
        m_sigy_std = value;
    }
 
    const array_type::tensor<double, Dim>& sigmay_std() const
    {
        return m_sigy_std;
    }
 
    void set_sigma(const array_type::tensor<double, Dim>& sigma)
    {
        m_sig = sigma;
        this->stress_updated();
    }
 
    const array_type::tensor<double, Dim>& sigma() const
    {
        return m_sig;
    }
 
    void set_sigmabar(double sigmabar)
    {
        if constexpr (std::is_same<StressFixer, detail::FixStress>::value) {
            m_stress_fixer->set_sigmabar(sigmabar);
            m_stress_fixer->set_average_stress(m_sig);
        }
        else if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
            double dsigbar = sigmabar - xt::mean(m_sig)();
            m_sig += dsigbar;
            m_stress_fixer->frame_match_stress(dsigbar);
        }
        else {
            m_sig -= xt::mean(m_sig)() - sigmabar;
        }
        this->stress_updated();
    }
 
    double sigmabar() const
    {
        if constexpr (std::is_same<StressFixer, detail::FixStress>::value) {
            return m_stress_fixer->sigmabar();
        }
        else {
            return xt::mean(m_sig)();
        }
    }
 
    double epsframe() const
    {
        if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
            return m_stress_fixer->u();
        }
        throw std::invalid_argument("Only implemented for Spring");
    }
 
    void set_epsframe(double u)
    {
        if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
            return m_stress_fixer->set_u(u);
        }
        throw std::invalid_argument("Only implemented for Spring");
    }
 
    void initSigmaFast(double sigma_std, size_t delta_r)
    {
        m_sig.fill(0);
        ptrdiff_t d = static_cast<ptrdiff_t>(delta_r);
        auto dsig = m_gen.normal<decltype(m_sig)>(m_sig.shape(), 0, sigma_std);
 
        if constexpr (Dim == 1) {
            for (ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(m_sig.size()); ++i) {
                m_sig(i) += dsig(i);
                m_sig.periodic(i - d) -= 0.5 * dsig(i);
                m_sig.periodic(i + d) -= 0.5 * dsig(i);
            }
        }
        else if constexpr (Dim == 2) {
            for (ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(m_sig.shape(0)); ++i) {
                for (ptrdiff_t j = 0; j < static_cast<ptrdiff_t>(m_sig.shape(1)); ++j) {
                    m_sig(i, j) += dsig(i, j);
                    m_sig.periodic(i - d, j) -= 0.5 * dsig(i, j);
                    m_sig.periodic(i + d, j) -= 0.5 * dsig(i, j);
                    m_sig.periodic(i, j - d) -= 0.5 * dsig(i, j);
                    m_sig.periodic(i, j + d) -= 0.5 * dsig(i, j);
                    m_sig.periodic(i + d, j + d) += 0.25 * dsig(i, j);
                    m_sig.periodic(i - d, j + d) += 0.25 * dsig(i, j);
                    m_sig.periodic(i + d, j - d) += 0.25 * dsig(i, j);
                    m_sig.periodic(i - d, j - d) += 0.25 * dsig(i, j);
                }
            }
            m_sig *= 0.5;
        }
        this->stress_updated();
    }
 
    void initSigmaPropogator(double sigma_std)
    {
        m_sig.fill(0);
        auto dsig = m_gen.normal<decltype(m_sig)>(m_sig.shape(), 0, sigma_std);
 
        if constexpr (Dim == 1) {
            ptrdiff_t ni = static_cast<ptrdiff_t>(m_sig.size());
            ptrdiff_t nj = static_cast<ptrdiff_t>(m_propagator.size());
            ptrdiff_t dj = m_offset[0];
            for (ptrdiff_t i = 0; i < ni; ++i) {
                for (ptrdiff_t j = 0; j < nj; ++j) {
                    m_sig.periodic(i + dj + j) -= dsig(i) * m_propagator(j);
                }
            }
        }
        else if constexpr (Dim == 2) {
            ptrdiff_t ni = static_cast<ptrdiff_t>(m_sig.shape(0));
            ptrdiff_t nj = static_cast<ptrdiff_t>(m_sig.shape(1));
            ptrdiff_t nk = static_cast<ptrdiff_t>(m_propagator.shape(0));
            ptrdiff_t nl = static_cast<ptrdiff_t>(m_propagator.shape(1));
            ptrdiff_t dk = m_offset[0];
            ptrdiff_t dl = m_offset[1];
 
            for (ptrdiff_t i = 0; i < ni; ++i) {
                for (ptrdiff_t j = 0; j < nj; ++j) {
                    for (ptrdiff_t k = 0; k < nk; ++k) {
                        for (ptrdiff_t l = 0; l < nl; ++l) {
                            m_sig.periodic(i + dk + k, j + dl + l) -=
                                dsig(i, j) * m_propagator(k, l);
                        }
                    }
                }
            }
        }
 
        m_sig /= xt::sqrt(xt::sum(xt::pow(m_propagator, 2.0)));
        this->stress_updated();
    }
 
    void makeAthermalFailureStep()
    {
        size_t nfailing = xt::sum(m_unstable)(); // todo: costly, cache?
 
        if (nfailing == 0) {
            m_failure_idx = m_epsp.size();
            return;
        }
 
        // == m_gen.exponential() / (m_failure_rate * static_cast<double>(nfailing));
        m_t -= std::log(1.0 - m_gen.random()) / (m_failure_rate * static_cast<double>(nfailing));
 
        size_t i = 0;
        if (nfailing > 1) {
            i = m_gen.randint(nfailing);
        }
 
        size_t idx;
        for (idx = 0; idx < m_sig.size(); ++idx) {
            if (m_unstable.flat(idx)) {
                if (i == 0) {
                    break;
                }
                --i;
            }
        }
 
        this->spatialParticleFailure(idx);
    }
 
    void makeAthermalFailureSteps(size_t n)
    {
        for (size_t i = 0; i < n; ++i) {
            this->makeAthermalFailureStep();
        }
    }
 
    size_t makeAthermalFailureSteps_chrono(size_t n, double elapsed)
    {
        std::chrono::steady_clock::time_point tic = std::chrono::steady_clock::now();
        auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(tic - tic).count();
        dt = static_cast<decltype(dt)>(elapsed * 1e3);
 
        for (size_t i = 0; i < n; ++i) {
            this->makeAthermalFailureStep();
 
            std::chrono::steady_clock::time_point toc = std::chrono::steady_clock::now();
            if (std::chrono::duration_cast<std::chrono::milliseconds>(toc - tic).count() > dt) {
                return i + 1;
            }
        }
 
        return n;
    }
 
    void makeWeakestFailureStep(bool allow_stable = false)
    {
        size_t idx;
        double x;
 
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            idx = xt::argmin(m_sigy - m_sig)(); // todo: costly, cache?
            x = m_sigy.flat(idx) - m_sig.flat(idx);
        }
        else {
            idx = xt::argmax(xt::abs(m_sig) - m_sigy)(); // todo: costly, cache?
            x = m_sigy.flat(idx) - std::abs(m_sig.flat(idx));
        }
 
        if (x > 0 && !allow_stable) {
            m_failure_idx = m_epsp.size();
            return;
        }
 
        m_t += 1.0;
        this->spatialParticleFailure(idx);
    }
 
    void makeWeakestFailureSteps(size_t n, bool allow_stable = false)
    {
        for (size_t i = 0; i < n; ++i) {
            this->makeWeakestFailureStep(allow_stable);
        }
    }
 
    size_t makeWeakestFailureSteps_chrono(size_t n, double elapsed, bool allow_stable = false)
    {
        std::chrono::steady_clock::time_point tic = std::chrono::steady_clock::now();
        auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(tic - tic).count();
        dt = static_cast<decltype(dt)>(elapsed * 1e3);
 
        for (size_t i = 0; i < n; ++i) {
            this->makeWeakestFailureStep(allow_stable);
 
            std::chrono::steady_clock::time_point toc = std::chrono::steady_clock::now();
            if (std::chrono::duration_cast<std::chrono::milliseconds>(toc - tic).count() > dt) {
                return i + 1;
            }
        }
 
        return n;
    }
 
    double temperature() const
    {
        if constexpr (!std::is_same<Temperature, detail::Finite>::value) {
            throw std::invalid_argument("Only implemented for thermal systems");
        }
        return m_temperature;
    }
 
    void set_temperature(double temperature)
    {
        if constexpr (!std::is_same<Temperature, detail::Finite>::value) {
            throw std::invalid_argument("Only implemented for thermal systems");
        }
        m_temperature = temperature;
        m_inv_temp = 1.0 / m_temperature;
        this->stress_updated();
    }
 
    array_type::tensor<double, Dim> time_to_failure() const
    {
        if constexpr (std::is_same<Temperature, detail::Finite>::value) {
            return m_tau_thermal;
        }
        else {
            return m_tau * xt::ones_like(m_sig);
        }
    }
 
    void makeThermalFailureStep()
    {
        if constexpr (!std::is_same<Temperature, detail::Finite>::value) {
            throw std::invalid_argument("Only implemented for thermal systems");
        }
        // == m_gen.exponential<decltype(m_sig)>(m_sig.shape()) * m_tau_thermal (but faster)
        auto dt =
            xt::eval(-xt::log(1.0 - m_gen.random<decltype(m_sig)>(m_sig.shape())) * m_tau_thermal);
        size_t idx = xt::argmin(dt)();
        m_t += dt.flat(idx);
        this->spatialParticleFailure(idx);
    }
 
    void makeThermalFailureSteps(size_t n)
    {
        for (size_t i = 0; i < n; ++i) {
            this->makeThermalFailureStep();
        }
    }
 
    size_t makeThermalFailureSteps_chrono(size_t n, double elapsed)
    {
        std::chrono::steady_clock::time_point tic = std::chrono::steady_clock::now();
        auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(tic - tic).count();
        dt = static_cast<decltype(dt)>(elapsed * 1e3);
 
        for (size_t i = 0; i < n; ++i) {
            this->makeThermalFailureStep();
 
            std::chrono::steady_clock::time_point toc = std::chrono::steady_clock::now();
            if (std::chrono::duration_cast<std::chrono::milliseconds>(toc - tic).count() > dt) {
                return i + 1;
            }
        }
 
        return n;
    }
 
    const array_type::tensor<double, Dim> x() const
    {
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            return m_sigy - m_sig;
        }
        else {
            return m_sigy - xt::abs(m_sig);
        }
    }
 
    const array_type::tensor<bool, Dim>& unstable() const
    {
        return m_unstable;
    }
 
    void spatialParticleFailure(size_t idx)
    {
        double dsig;
 
        m_failure_idx = idx;
        if constexpr (std::is_same<FailureRules, detail::Depinning>::value) {
            m_failure_x = m_sigy.flat(idx) - m_sig.flat(idx);
            dsig = 1.0;
        }
        else {
            m_failure_x = m_sigy.flat(idx) - std::abs(m_sig.flat(idx));
            dsig = m_sig.flat(idx) + m_gen.normal(0.0, 0.01);
        }
 
        m_nfails.flat(idx) += 1;
        m_epsp.flat(idx) += dsig;
        this->drawYieldStress(idx);
 
        if constexpr (Dim == 1) {
            ptrdiff_t di = static_cast<ptrdiff_t>(idx) + m_offset[0];
            ptrdiff_t n = static_cast<ptrdiff_t>(m_sig.size());
 
            for (ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(m_propagator.size()); ++i) {
                // index corrected for periodicity
                ptrdiff_t ii = (n + di + i) % n;
                // propagate
                m_sig.flat(ii) += m_propagator(i) * dsig;
                if constexpr (!std::is_same<StressFixer, detail::Spring>::value) {
                    this->stress_updated(ii);
                }
            }
        }
        else if constexpr (Dim == 2) {
            auto index = xt::unravel_index(idx, m_sig.shape());
            ptrdiff_t di = static_cast<ptrdiff_t>(index[0]) + m_offset[0];
            ptrdiff_t dj = static_cast<ptrdiff_t>(index[1]) + m_offset[1];
            ptrdiff_t n = static_cast<ptrdiff_t>(m_sig.shape(0));
            ptrdiff_t m = static_cast<ptrdiff_t>(m_sig.shape(1));
 
            for (ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(m_propagator.shape(0)); ++i) {
                for (ptrdiff_t j = 0; j < static_cast<ptrdiff_t>(m_propagator.shape(1)); ++j) {
                    // index corrected for periodicity
                    ptrdiff_t ii = (n + di + i) % n;
                    ptrdiff_t jj = (m + dj + j) % m;
                    ptrdiff_t iidx = ii * m + jj;
                    // propagate
                    m_sig.flat(iidx) += m_propagator(i, j) * dsig;
                    if constexpr (!std::is_same<StressFixer, detail::Spring>::value) {
                        this->stress_updated(iidx);
                    }
                }
            }
        }
 
        if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
            m_stress_fixer->stress_match_frame(m_sig, dsig);
            this->stress_updated();
        }
    }
 
    double shiftImposedShear(int direction = 1)
    {
        double dsig;
 
        if (direction > 0) {
            dsig = xt::amin(m_sigy - m_sig)() + 2.0 * std::numeric_limits<double>::epsilon();
        }
        else {
            dsig = -xt::amin(m_sig + m_sigy)() + 2.0 * std::numeric_limits<double>::epsilon();
        }
 
        m_sig += dsig;
        this->stress_updated();
 
        if constexpr (std::is_same<StressFixer, detail::Spring>::value) {
            m_stress_fixer->frame_match_stress(dsig);
        }
 
        return dsig;
    }
 
    size_t relaxAthermal(size_t max_steps = 1e9, bool max_steps_is_error = true)
    {
        for (size_t i = 0; i < max_steps; ++i) {
            this->makeAthermalFailureStep();
            if (m_failure_idx == m_epsp.size()) {
                return i;
            }
        }
 
        if (max_steps_is_error) {
            throw std::runtime_error("Failed to converge.");
        }
 
        return max_steps;
    }
 
    size_t relaxWeakest(size_t max_steps = 1e9, bool max_steps_is_error = true)
    {
        for (size_t i = 0; i < max_steps; ++i) {
            this->makeWeakestFailureStep();
            if (m_failure_idx == m_epsp.size()) {
                return i;
            }
        }
 
        if (max_steps_is_error) {
            throw std::runtime_error("Failed to converge.");
        }
 
        return max_steps;
    }
};
 
template <size_t Dim>
class SystemStrainControl : public SystemBase<Dim, detail::FixStrain, void, void> {
public:
    SystemStrainControl(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_mean,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
    {
        this->initSystem(
            propagator, distances, sigmay_mean, sigmay_std, seed, failure_rate, alpha, random_stress
        );
    }
};
 
template <size_t Dim>
class SystemStressControl : public SystemBase<Dim, detail::FixStress, void, void> {
private:
    detail::FixStress m_my_stress_fixer;
 
public:
    SystemStressControl(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_mean,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
        : m_my_stress_fixer(0.0)
    {
        this->initSystem(
            propagator,
            distances,
            sigmay_mean,
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress,
            &m_my_stress_fixer
        );
    }
};
 
template <size_t Dim>
class SystemSpringLoading : public SystemBase<Dim, detail::Spring, void, void> {
private:
    detail::Spring m_my_stress_fixer;
 
public:
    SystemSpringLoading(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_mean,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double kframe,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
        : m_my_stress_fixer(kframe, 0.0, sigmay_mean.size())
    {
        this->initSystem(
            propagator,
            distances,
            sigmay_mean,
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress,
            &m_my_stress_fixer
        );
    }
};
 
template <size_t Dim>
class SystemThermalStressControl : public SystemBase<Dim, detail::FixStress, detail::Finite, void> {
private:
    detail::FixStress m_my_stress_fixer;
 
public:
    SystemThermalStressControl(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_mean,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
        : m_my_stress_fixer(0.0)
    {
        this->initSystem(
            propagator,
            distances,
            sigmay_mean,
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress,
            &m_my_stress_fixer
        );
    }
};
 
template <size_t Dim>
class SystemDepinningStrainControl
    : public GooseEPM::
          SystemBase<Dim, GooseEPM::detail::FixStrain, void, GooseEPM::detail::Depinning> {
public:
    SystemDepinningStrainControl(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
    {
        this->initSystem(
            propagator,
            distances,
            xt::zeros_like(sigmay_std),
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress
        );
    }
};
 
template <size_t Dim>
class SystemDepinningStressControl
    : public GooseEPM::
          SystemBase<Dim, GooseEPM::detail::FixStress, void, GooseEPM::detail::Depinning> {
private:
    GooseEPM::detail::FixStress m_my_stress_fixer;
 
public:
    SystemDepinningStressControl(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
        : m_my_stress_fixer(0.0)
    {
        this->initSystem(
            propagator,
            distances,
            xt::zeros_like(sigmay_std),
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress,
            &m_my_stress_fixer
        );
    }
};
 
template <size_t Dim>
class SystemDepinningSpringLoading
    : public GooseEPM::
          SystemBase<Dim, GooseEPM::detail::Spring, void, GooseEPM::detail::Depinning> {
private:
    GooseEPM::detail::Spring m_my_stress_fixer;
 
public:
    SystemDepinningSpringLoading(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double kframe,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
        : m_my_stress_fixer(kframe, 0.0, sigmay_std.size())
    {
        this->initSystem(
            propagator,
            distances,
            xt::zeros_like(sigmay_std),
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress,
            &m_my_stress_fixer
        );
    }
};
 
template <size_t Dim>
class SystemDepinningThermalStressControl : public GooseEPM::SystemBase<
                                                Dim,
                                                GooseEPM::detail::FixStress,
                                                GooseEPM::detail::Finite,
                                                GooseEPM::detail::Depinning> {
private:
    GooseEPM::detail::FixStress m_my_stress_fixer;
 
public:
    SystemDepinningThermalStressControl(
        const array_type::tensor<double, Dim>& propagator,
        const std::vector<array_type::tensor<ptrdiff_t, 1>>& distances,
        const array_type::tensor<double, Dim>& sigmay_std,
        uint64_t seed,
        double failure_rate = 1,
        double alpha = 1.5,
        bool random_stress = true
    )
        : m_my_stress_fixer(0.0)
    {
        this->initSystem(
            propagator,
            distances,
            xt::zeros_like(sigmay_std),
            sigmay_std,
            seed,
            failure_rate,
            alpha,
            random_stress,
            &m_my_stress_fixer
        );
    }
};
 
inline std::tuple<
    array_type::tensor<size_t, 1>,
    array_type::tensor<size_t, 1>,
    array_type::tensor<size_t, 1>>
segment_avalanche(
    const array_type::tensor<bool, 1>& condition,
    array_type::tensor<size_t, 1> idx,
    bool first = true,
    bool last = true
)
{
    GOOSEEPM_REQUIRE(xt::has_shape(condition, idx.shape()), std::out_of_range);
 
    if (condition.size() == 0) {
        return std::make_tuple(
            xt::empty<size_t>({0}), xt::empty<size_t>({0}), xt::empty<size_t>({0})
        );
    }
 
    if (xt::all(condition)) {
        if (!first || !last) {
            return std::make_tuple(
                xt::empty<size_t>({0}), xt::empty<size_t>({0}), xt::empty<size_t>({0})
            );
        }
        std::sort(idx.begin(), idx.end());
        size_t a = std::unique(idx.begin(), idx.end()) - idx.begin();
        array_type::tensor<size_t, 1> S = {condition.size()};
        array_type::tensor<size_t, 1> A = {a};
        array_type::tensor<size_t, 1> index = {0};
        return std::make_tuple(S, A, index);
    }
 
    size_t n = condition.size();
    if (!last) {
        for (size_t i = condition.size() - 1; i > 0; --i) {
            if (!condition(i)) {
                n = i + 1;
                break;
            }
        }
    }
 
    size_t nevent = static_cast<size_t>(condition(0));
    for (size_t i = 0; i < n - 1; ++i) {
        if (!condition(i) && condition(i + 1)) {
            ++nevent;
        }
    }
 
    size_t i = 0;
    if (!first && condition(0)) {
        --nevent;
        for (; i < condition.size() - 1; ++i) {
            if (!condition(i)) {
                break;
            }
        }
    }
 
    array_type::tensor<size_t, 1> S = xt::empty<size_t>({nevent});
    array_type::tensor<size_t, 1> A = xt::empty<size_t>({nevent});
    array_type::tensor<size_t, 1> index = xt::empty<size_t>({nevent});
    size_t event = 0;
    size_t j;
 
    for (; i < n; ++i) {
        if (condition(i)) {
            for (j = i + 1; j < n; ++j) {
                if (!condition(j)) {
                    break;
                }
            }
            if (j < n) {
                auto begin = &idx(i);
                auto end = &idx(j);
                std::sort(begin, end);
                S(event) = j - i;
                A(event) = std::unique(begin, end) - begin;
                index(event) = i;
                ++event;
                i = j;
            }
            else {
                break;
            }
        }
    }
 
    if (last && j == n) {
        if (condition(j - 1)) {
            auto begin = &idx(i);
            auto end = &idx(j - 1) + 1;
            std::sort(begin, end);
            S(event) = j - i;
            A(event) = std::unique(begin, end) - begin;
            index(event) = i;
        }
    }
 
    return std::make_tuple(S, A, index);
}
 
inline array_type::tensor<size_t, 1> cumsum_n_unique(const array_type::tensor<size_t, 1>& idx)
{
    array_type::tensor<size_t, 1> ret = xt::empty_like(idx);
    xt::xtensor<bool, 1> failed = xt::zeros<bool>({xt::amax(idx)() + 1});
 
    if (idx.size() == 0) {
        return ret;
    }
 
    ret(0) = 1;
    failed(idx(0)) = true;
 
    for (size_t i = 1; i < idx.size(); ++i) {
        if (!failed(idx(i))) {
            failed(idx(i)) = true;
            ret(i) = ret(i - 1) + 1;
        }
        else {
            ret(i) = ret(i - 1);
        }
    }
 
    return ret;
}
 
template <size_t Dimension>
class [[deprecated("Use GooseEYE instead")]] AvalancheSegmenter {
public:
    static constexpr size_t Dim = Dimension; 
 
private:
    std::array<size_t, Dim> m_shape; 
    array_type::tensor<size_t, 1> m_idx; 
    array_type::tensor<double, 1> m_time; 
    array_type::tensor<ptrdiff_t, Dim> m_dx; 
    array_type::tensor<size_t, Dim> m_s; 
    array_type::tensor<ptrdiff_t, Dim> m_label; 
    ptrdiff_t m_new_label = 1; 
    size_t m_step = 0; 
 
    std::vector<ptrdiff_t> m_renum;
 
    std::vector<ptrdiff_t> m_next;
 
    std::vector<ptrdiff_t> m_connected; 
 
public:
    AvalancheSegmenter(
        const std::array<size_t, Dim>& shape,
        const array_type::tensor<size_t, 1>& idx,
        const array_type::tensor<double, 1>& t
    )
    {
        GOOSEEPM_WARNING_PYTHON("AvalancheSegmenter is deprecated, use GooseEYE instead.");
        m_shape = shape;
        m_idx = idx;
        m_time = t;
        m_s = xt::zeros<size_t>(m_shape);
        m_label = xt::zeros<ptrdiff_t>(m_shape);
        m_renum.resize(m_s.size() + 1);
        m_next.resize(m_s.size() + 1);
        std::iota(m_renum.begin(), m_renum.end(), 0);
        this->clean_next();
 
        GOOSEEPM_REQUIRE(xt::has_shape(m_idx, t.shape()), std::invalid_argument);
        GOOSEEPM_REQUIRE(std::is_sorted(m_time.begin(), m_time.end()), std::invalid_argument);
        GOOSEEPM_REQUIRE(xt::amax(m_idx)() < m_s.size(), std::invalid_argument);
 
        if constexpr (Dim == 1) {
            // kernel = {1, 1, 1}
            m_dx = {-1, 1};
        }
        else if constexpr (Dim == 2) {
            // kernel = {{0, 1, 0}, {1, 1, 1}, {0, 1, 0}};
            m_dx = {{-1, 0}, {0, -1}, {0, 1}, {1, 0}};
        }
        m_connected.resize(m_dx.shape(0));
    }
 
    size_t nlabels()
    {
        return m_new_label;
    }
 
    void prune()
    {
        ptrdiff_t n = static_cast<ptrdiff_t>(m_new_label);
        m_new_label = 1;
        m_renum[0] = 0;
        for (ptrdiff_t i = 1; i < n; ++i) {
            if (m_renum[i] == i) {
                m_renum[i] = m_new_label;
                ++m_new_label;
            }
        }
        this->private_renumber(m_renum);
        std::iota(m_renum.begin(), m_renum.begin() + n, 0);
    }
 
private:
    void clean_next()
    {
        std::fill(m_next.begin(), m_next.end(), -1);
    }
 
    template <class T>
    void private_renumber(const T& renum)
    {
        for (size_t i = 0; i < m_label.size(); ++i) {
            m_label.flat(i) = renum[m_label.flat(i)];
        }
    }
 
public:
    void change_labels(const array_type::tensor<size_t, 1>& new_label)
    {
        GOOSEEPM_ASSERT(new_label.size() >= this->nlabels(), std::out_of_range);
        GOOSEEPM_ASSERT(
            *std::max_element(new_label.begin(), new_label.end()) < this->size() + 1,
            std::out_of_range
        );
 
        this->private_renumber(new_label);
        m_new_label = xt::amax(m_label)() + 1;
        std::fill(m_renum.begin(), m_renum.begin() + m_new_label, 0);
 
        for (size_t i = 0; i < new_label.size(); ++i) {
            m_renum[new_label[i]] = new_label[i];
        }
    }
 
    void reorder(const array_type::tensor<ptrdiff_t, 1>& labels)
    {
        auto maxlab = *std::max_element(labels.begin(), labels.end());
 
#ifdef GOOSEEPM_ENABLE_ASSERT
        GOOSEEPM_ASSERT(maxlab < m_new_label, std::out_of_range);
 
        std::vector<bool> label_in_use(m_new_label, false);
        for (size_t i = 0; i < m_label.size(); ++i) {
            label_in_use[m_label.flat(i)] = true;
        }
 
        std::vector<bool> label_in_order(m_new_label, false);
        for (size_t i = 0; i < labels.size(); ++i) {
            label_in_order[labels[i]] = true;
        }
 
        for (ptrdiff_t i = 0; i < m_new_label; ++i) {
            if (label_in_use[i]) {
                GOOSEEPM_ASSERT(label_in_order[i], std::out_of_range);
            }
        }
#endif
 
        decltype(m_renum) tmp(m_new_label);
 
        for (size_t i = 0; i < labels.size(); ++i) {
            tmp[labels[i]] = i;
            m_renum[labels[i]] = i;
        }
 
        this->private_renumber(tmp);
        m_new_label = xt::amax(m_label)() + 1;
        std::iota(m_renum.begin(), m_renum.begin() + m_new_label, 0);
    }
 
private:
    void merge_detail(ptrdiff_t a, ptrdiff_t b)
    {
        // -> head[list(b)] = head[a]
        ptrdiff_t i = m_renum[b];
        ptrdiff_t target = m_renum[a];
        m_renum[b] = target;
        while (true) {
            i = m_next[i];
            if (i == -1) {
                break;
            }
            m_renum[i] = target;
        }
        // -> list(head[a]).append(list(b))
        while (m_next[a] != -1) {
            a = m_next[a];
        }
        m_next[a] = b;
    }
 
    ptrdiff_t merge(ptrdiff_t* labels, size_t nlabels)
    {
        std::sort(labels, labels + nlabels);
        nlabels = std::unique(labels, labels + nlabels) - labels;
        ptrdiff_t target = labels[0];
        for (size_t i = 1; i < nlabels; ++i) {
            this->merge_detail(target, labels[i]);
        }
        return target;
    }
 
public:
    void advance(size_t n = 1)
    {
        size_t nmerge = 0;
        size_t nconnected;
        ptrdiff_t iidx;
        size_t nsteps = std::min(m_step + n, m_idx.size());
        std::array<ptrdiff_t, Dim> index;
        std::array<ptrdiff_t, Dim> shape;
        std::array<ptrdiff_t, Dim> strides;
        for (size_t i = 0; i < Dim; ++i) {
            shape[i] = static_cast<ptrdiff_t>(m_shape[i]);
            strides[i] = static_cast<ptrdiff_t>(m_label.strides()[i]);
        }
 
        for (; m_step < nsteps; ++m_step) {
            auto idx = static_cast<ptrdiff_t>(m_idx(m_step));
            m_s.flat(idx) += 1;
 
            // if the block was labelled before there is nothing further to do
            if (m_s.flat(idx) > 1) {
                continue;
            }
 
            nconnected = 0;
            for (size_t j = 0; j < m_dx.shape(0); ++j) {
                if constexpr (Dim == 1) {
                    ptrdiff_t nn = shape[0];
                    ptrdiff_t ii = (nn + idx + m_dx(j)) % nn; // index corrected for periodicity
                    iidx = ii;
                }
                else if constexpr (Dim == 2) {
                    detail::unravel_index(idx, strides, index);
                    ptrdiff_t nn = shape[0];
                    ptrdiff_t mm = shape[1];
                    ptrdiff_t ii = (nn + index[0] + m_dx(j, 0)) % nn;
                    ptrdiff_t jj = (mm + index[1] + m_dx(j, 1)) % mm;
                    iidx = ii * mm + jj;
                }
                if (m_label.flat(iidx) != 0) {
                    m_connected[nconnected] = m_renum[m_label.flat(iidx)];
                    nconnected++;
                }
            }
            if (nconnected == 0) {
                m_label.flat(idx) = m_new_label;
                m_new_label += 1;
            }
            else if (nconnected == 1) {
                m_label.flat(idx) = m_connected[0];
            }
            else {
                // mark all labels in the list for merging
                // `m_label` is not yet updated to avoid looping over all blocks too frequently
                // the new label can be read by `m_renum[lab]` (as done above)
                m_label.flat(idx) = this->merge(&m_connected[0], nconnected);
                nmerge++;
 
                // every so often: apply the renumbering to `m_label`
                // the linked labels in `m_next` can be released
                if (nmerge > 100) {
                    this->private_renumber(m_renum);
                    this->clean_next();
                    nmerge = 0;
                }
            }
        }
 
        // apply the renumbering to `m_label`
        // the linked labels in `m_next` can be released
        if (nmerge > 0) {
            this->private_renumber(m_renum);
            this->clean_next();
        }
    }
 
    void advance_to(double t, bool floor = false)
    {
        if (t >= m_time.back()) {
            return this->advance(m_time.size() - m_step);
        }
        if (t <= m_time(m_step)) {
            return;
        }
 
        for (size_t i = 0; i < m_time.size() - m_step; ++i) {
            if (m_time(m_step + i) > t) {
                if (floor) {
                    return this->advance(i - 1);
                }
                else {
                    if (m_time(m_step + i) - t < t - m_time(m_step + i - 1)) {
                        return this->advance(i);
                    }
                    else {
                        return this->advance(i - 1);
                    }
                }
            }
        }
    }
 
    std::string repr() const
    {
        return detail::get_namespace() + "AvalancheSegmenter" + std::to_string(Dim) + " " +
               detail::shape_to_string(m_shape);
    }
 
    const auto& shape() const
    {
        return m_label.shape();
    }
 
    auto size() const
    {
        return m_label.size();
    }
 
    const auto& s() const
    {
        return m_s;
    }
 
    // todo: allow resetting a cluster map ?
 
    const auto& labels() const
    {
        return m_label;
    }
 
    size_t nstep()
    {
        return m_step;
    }
 
    auto t() const
    {
        if (m_step == 0) {
            return 0.0;
        }
        return m_time(m_step - 1);
    }
};
 
class Avalanche {
private:
    array_type::tensor<size_t, 1> m_idx; 
    array_type::tensor<double, 1> m_x; 
    array_type::tensor<double, 1> m_t; 
 
public:
    Avalanche() = default;
 
    std::string repr() const
    {
        return detail::get_namespace() + "Avalanche" + " " +
               detail::shape_to_string(std::array<size_t, 1>{m_idx.size()});
    }
 
    const array_type::tensor<size_t, 1>& idx() const
    {
        return m_idx;
    }
 
    const array_type::tensor<double, 1>& x() const
    {
        return m_x;
    }
 
    const array_type::tensor<double, 1>& t() const
    {
        return m_t;
    }
 
    template <class System, class Func>
    void makeFailureSteps(System& system, size_t n, Func func)
    {
        m_idx = xt::empty<size_t>({n});
        m_x = xt::empty<double>({n});
        m_t = xt::empty<double>({n});
 
        for (size_t i = 0; i < n; ++i) {
            func();
            m_idx(i) = system.failure_idx();
            m_x(i) = system.failure_x();
            m_t(i) = system.t();
        }
    }
 
    template <class System, class Func>
    size_t makeFailureSteps_chrono(System& system, size_t n, double elapsed, Func func)
    {
        std::chrono::steady_clock::time_point tic = std::chrono::steady_clock::now();
        auto dt = std::chrono::duration_cast<std::chrono::milliseconds>(tic - tic).count();
        dt = static_cast<decltype(dt)>(elapsed * 1e3);
 
        std::array<size_t, 1> shape = {n};
        m_idx.resize(shape);
        m_x.resize(shape);
        m_t.resize(shape);
 
        for (size_t i = 0; i < n; ++i) {
            func();
            m_idx(i) = system.failure_idx();
            m_x(i) = system.failure_x();
            m_t(i) = system.t();
 
            std::chrono::steady_clock::time_point toc = std::chrono::steady_clock::now();
            if (std::chrono::duration_cast<std::chrono::milliseconds>(toc - tic).count() > dt) {
                // todo: remove copy if xtensor.resize is item preserving
                // https://github.com/xtensor-stack/xtensor/issues/2745
                size_t nret = i + 1;
                shape[0] = nret;
                {
                    std::vector<size_t> tmp(nret);
                    std::copy(m_idx.begin(), m_idx.begin() + nret, tmp.begin());
                    m_idx.resize(shape);
                    std::copy(tmp.begin(), tmp.end(), m_idx.begin());
                }
                {
                    std::vector<double> tmp(nret);
                    std::copy(m_x.begin(), m_x.begin() + nret, tmp.begin());
                    m_x.resize(shape);
                    std::copy(tmp.begin(), tmp.end(), m_x.begin());
                }
                {
                    std::vector<double> tmp(nret);
                    std::copy(m_t.begin(), m_t.begin() + nret, tmp.begin());
                    m_t.resize(shape);
                    std::copy(tmp.begin(), tmp.end(), m_t.begin());
                }
                return nret;
            }
        }
 
        return n;
    }
 
    template <class System>
    void makeAthermalFailureSteps(System& system, size_t n)
    {
        this->makeFailureSteps(system, n, [&system]() { return system.makeAthermalFailureStep(); });
    }
 
    template <class System>
    void makeWeakestFailureSteps(System& system, size_t n, bool allow_stable = false)
    {
        this->makeFailureSteps(system, n, [&system, allow_stable]() {
            return system.makeWeakestFailureStep(allow_stable);
        });
    }
 
    template <class System>
    void makeThermalFailureSteps(System& system, size_t n)
    {
        this->makeFailureSteps(system, n, [&system]() { return system.makeThermalFailureStep(); });
    }
 
    template <class System>
    size_t makeAthermalFailureSteps_chrono(System& system, size_t n, double elapsed)
    {
        return this->makeFailureSteps_chrono(system, n, elapsed, [&system]() {
            return system.makeAthermalFailureStep();
        });
    }
 
    template <class System>
    size_t makeWeakestFailureSteps_chrono(
        System& system,
        size_t n,
        double elapsed,
        bool allow_stable = false
    )
    {
        return this->makeFailureSteps_chrono(system, n, elapsed, [&system, allow_stable]() {
            return system.makeWeakestFailureStep(allow_stable);
        });
    }
 
    template <class System>
    size_t makeThermalFailureSteps_chrono(System& system, size_t n, double elapsed)
    {
        return this->makeFailureSteps_chrono(system, n, elapsed, [&system]() {
            return system.makeThermalFailureStep();
        });
    }
};
 
} // namespace GooseEPM
 
#endif