Speed up preprocess_strains_with_limits by removing per-call np.isclose

[magnusfjeldolsen]

ConstitutiveLaw.preprocess_strains_with_limits snaps strains that fall within a
small tolerance of an ultimate-strain limit onto that exact limit. It runs at the
top of every constitutive law's get_stress/get_tangent, so the integrators
call it every time they evaluate stress or tangential stiffness — for the fiber
integrator, once per fiber on each iteration of the equilibrium solver (typically
tens of thousands of strain values per moment curvature calculation).

On every call it does more work than needed: it runs np.isclose twice and
allocates two np.zeros_like(eps) arrays, even though it is only comparing each
strain against the limit.

idxs = np.isclose(eps, np.zeros_like(eps) + eps_max, atol=1e-6)
eps[idxs] = eps_max
idxs = np.isclose(eps, np.zeros_like(eps) + eps_min, atol=1e-6)
eps[idxs] = eps_min

preprocess_strains_with_limits reads the limits once per call via
get_ultimate_strain(), so within a call the closeness threshold is a single
scalar shared by every value in eps. np.isclose(a, b) tests
|a - b| <= atol + rtol * |b|
(numpy docs);
here b is the limit, atol=1e-6 is passed and rtol keeps its default 1e-5,
so the threshold is just 1e-6 + 1e-5 * |limit|. Precomputing that scalar once
per call and comparing with a plain abs(...) <= tol gives bit-identical masks
for every law — it removes only the per-element np.isclose machinery and the
zeros_like allocations, not any behaviour.

Fiber benefits most — every fiber goes through this method. For Marin it usually
only benefits the rebars, except when Marin falls back to discretizing a
constitutive law whose stress–strain curve isn't a polynomial.

Example: 40 randomly generated sections, fiber integrator

Each section is randomly rectangular or circular (50/50), with parameters drawn
from:

parameter	rectangular	circular
outer size [mm]	width 200–500, height 300–800	diameter 300–800
cover [mm]	30–50	30–50
bar diameter [mm]	16 / 20 / 25 / 32	16 / 20 / 25 / 32
number of bars	2–5 per layer (top + bottom)	6–12 around perimeter
concrete class	C25–C45	C25–C45

Each section runs bending capacity + moment–curvature + N–M interaction domain,
with the current np.isclose implementation timed against the
precomputed-threshold version (identical results on every section):

metric	per-section speedup
calculation time	median 1.82×, range 1.56–2.15× (p10 1.69, p90 2.04)

So roughly a ~1.8× faster fiber calculation.

[talledodiego]

Hi, thanks for this very promising improvement. I will look at the PR ASAP.

It runs at the top of every constitutive law's get_stress/get_tangent, so the integrators call it every time they evaluate stress or tangential stiffness — for the fiber integrator, once per fiber on each iteration of the equilibrium solver (typically tens of thousands of strain values per moment curvature calculation).

Regarding this statement I am not sure. The function get_stress is not called once per fiber, but once per iteration. With fiber integrator get_stress is called once passing an array with the strains of all fibers. Indeed I profiled the function calculate_bending_strength() and preprocess_strains_with_limits is called 68 times with a total time of 0.001 s. The cumulated time is 0.007 s. The total time of calculate_bending_strength was 0.106 seconds, so that function alone is responsible for less than 1% of the total time.
For calculate_nm_interaction_domain() the preprocess_strain_with_limits is called 70 times with a total time of 0.001 s and a cumtime of 0.007 seconds. The total time for the function was 0.117 s, so again less than 1% of the total time.
For calculate_moment_curvature() the weight is higher: the preprocess_strain_with_limits is called 1294 times, with a total time of 0.018 s and a cumtime of 0.16 s. The total time for the function was 0.41s, so the weight of preprocess if of about 40%. What I would like to think more is if it is correct that the number of calls is so much higher for moment-curvature analysis, but I will think more about it in the coming days.

I was under the impression that doing explicitly numpy.isclose or the equation reported should take more or less the same time. Some profiling shows a difference of about 25% faster with the equation reported. I think anyhow I would prefer for clarity and transparency using numpy.isclose even if a bit less efficient (considering what said above it should not be the real bottleneck), even if we can think about this if also @mortenengen agrees on this.

I have the suspect anyhow the more time consuming and useless call is the zeros_like that has to create an array in memory. Numpy already handles scalar comparison. So these calls are equivalent:

a = np.linspace(-0.003000001, 0.003000001, 1_000_000)
lim = 0.003

def preprocess(eps, eps_max):
    eps[np.isclose(eps, eps_max, atol=1e-6)] = eps_max
    eps[np.isclose(eps, -eps_max, atol=1e-6)] = -eps_max
    return eps

def preprocess_old(eps, eps_max):
    idxs = np.isclose(eps, np.zeros_like(eps) + eps_max, atol=1e-6)
    eps[idxs] = eps_max
    idx = np.isclose(eps, np.zeros_like(eps) - eps_max, atol=1e-6)
    eps[idx] = -eps_max
    return eps

A quick profile of these two functions shows a bit of improvement (1.3x) even if this is not inside our algorithms and therefore the conclusions could be hard to extend.

[mortenengen]

Thanks for bringing this up @magnusfjeldolsen, and thanks for investigating @talledodiego. I agree with the observation that the useless call to zeros_like should be removed. If indeed we can demonstrate that np.isclose is significantly slower than the bare equations we should also refine this since the function is called many times throughout our algorithms, and since the get_stress and get_tangent methods are so fundamental to the library.