Error when reading e.g. stresses or strains from binout of model with eroded elements with numpy version >= 2.0

[BerSchneider]

🐛 Describe the bug
The code np.array(data) in method _get_variable in file binout.py leads to a ValueError if tuples in data have different lengths and numpy version >= 2.0, e.g.
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (26,) + inhomogeneous part.
This problem was already raised in issue #62, but now the dependency is numpy>2.2.

🔢 To Reproduce
Read e.g. stresses or strains from binout of model with eroded elements with numpy version >= 2.0

🖥️ Setup

• lasso-python version: development
• numpy version: 2.3.4

Possible Solutions
Would the numpy 1.x behavior be replicated if one used np.array(data, dtype=object) if the tuples have different lengths (either by explicitly checking or using a try except block)?
Another possibility would be to fill the missing values with e.g. np.nans but one would probably need to check the element IDs to fill in the correct columns.
In issue #62 there is a reply containing some code showing how the issue could be handled. If I understand it correctly, it just pads the missing values with zeros but does not take into account which elements were deleted.

[marcodesiderio]

I have implemented the solution under #62 in my 'local' lasso version and has worked great so far. Working in simulations where shells, solid & beam elements frequently get eroded and I have never seen anything wrong with the data. I think the basic assumption that any missing values are associated with an element that has been eroded is quite reasonable and robust.

[BerSchneider]

The reason I think the implementation in #62 might be confusing (at least if my understanding is correct) is the following. If no elements are eroded and I want the data of a particular element over all output time steps I get that by looking at the column associated with that element (at least if one can safely assume that the order does not change over the time steps). But if elements are eroded and the missing values are just padded with zeros at the end that is no longer the case. E.g. if I look at the last column, and the element associated with that column in the first time step is not eroded but others are, then there are zeros in that column for later time steps. With the behavior of numpy 1.x it was at least clear that further processing steps are necessary.
I guess it all comes down to what you do with the data. If you are only using, e.g., the strains of all available elements in an assessment then there should be no issue. However, if you are accessing strains of elements in individual parts of your simulation model (in my case, individual bones, e.g., the ribs), then by padding the array with zeros at the end, you will no longer be able to correctly access the strains of the elements in question or at least one has to perform further processing steps.

[marcodesiderio]

You are totally right! I actually thought I did fix this in my code, but eventually forgot to add that.
My idea is to use the idx elout output to match the ids with the actual output.

[marcodesiderio]

So in my program, just as context, I am reading the binout data and then storing it in a big dictionary.

With this in mind, the code I am pasting below will be the workaround I have AFTER extracting the binout data with lasso (so, binout.read('elout', element_type, output_type) -> remap data to correct data structure -> store data in dictionary).

This is currently working for me:

def build_mapping(ids_all: np.ndarray, nip_all: np.ndarray):
    """
    Construct a mapping array that expands each element ID into one label per
    integration point. Each label is formatted as "element_id/iptN", where N 
    is the integration-point index for that element.

    Parameters
    ----------
    ids_all : np.ndarray
        2D array of shape (n_timesteps, n_elements). Each row contains the
        element IDs for that timestep.
    nip_all : np.ndarray
        2D array of shape (n_timesteps, n_elements) specifying the number of
        integration points per element. If None, empty, or not 2D, all elements
        are assumed to have one integration point.

    Returns
    -------
    mapping : np.ndarray
        2D object array of shape (n_timesteps, nip_total), where nip_total is the
        total number of integration points for the first timestep. Each entry is
        a string of the form "element_id/iptN". Rows are padded with zeros if a
        timestep contains fewer total integration points.
    """
    
    if nip_all is None or nip_all.size == 0 or nip_all.ndim < 2:
        # => All elements have 1 integration point
        # Build a new nip_all with the same shape as ids_all but filled with 1s
        nip_all = np.ones_like(ids_all, dtype=int)

    nip_total = np.sum(nip_all[0, :])
    n_timesteps = ids_all.shape[0]

    mapping = np.zeros((n_timesteps, nip_total), dtype=object)

    for i in range(n_timesteps):
        start_idx = 0
        for idx, eid in enumerate(ids_all[i, :]):
            new_mapping = np.array([f"{eid}/ipt{j+1}" for j in range(nip_all[i, idx])])
            mapping[i, start_idx:start_idx+new_mapping.shape[0]] = new_mapping
            start_idx += new_mapping.shape[0]

    return mapping

def align_data(data: np.ndarray, mapping: np.ndarray) -> np.ndarray:
    """
    Reorder `data` so that each column index corresponds to the same label
    as in mapping[0, :], for all timesteps.

    Assumes:
    - data.shape == mapping.shape
    - Each row of `mapping` contains the same labels as row 0, plus possibly
      padded zeros (0 or '0') in some columns.
    """
    assert data.shape == mapping.shape, "Data and mapping must have the same shape"
    n_timesteps, n_cols = data.shape

    # Canonical labels from first row (assume no padding in row 0)
    headers = mapping[0, :]

    # Map label -> canonical column index, ignoring any zero-like padding just in case
    label_to_idx = {
        label: j
        for j, label in enumerate(headers)
        if not (label == 0 or label == "0" or label is None)
    }

    # Initialize aligned data with zeros (same behavior as your original)
    new_data = np.zeros((n_timesteps, n_cols), dtype=np.float64)

    for t in range(n_timesteps):
        row_labels = mapping[t, :]

        for src_col, label in enumerate(row_labels):
            # Skip padding: if label not in dict, we leave new_data[t, *] as 0
            dest_col = label_to_idx.get(label, None)
            if dest_col is None:
                continue

            new_data[t, dest_col] = data[t, src_col]

    return new_data

Essentially

1. build a mapping matrix similar to the ['elout', element_type, 'ids'] output you get from lasso, but just expanded to take into account that different elements will have different number of integration points. If you always have 1 ipt, then you do not need this (simply, raw_data[np.where(ids_all == requested_id)] will return the correct array)

2. use that mapping matrix rearrange the outputs so that one column now corresponds to one element.

I think a similar logic can be implemented directly into lasso. It's just easier for me to fix my own code however.

If there is a easier or faster way of achieving this, please share! - The 'mapping matrix' takes a bit to be generated, align_data runs quite fast.

[marcodesiderio]

** update **
I added the optional argument nodes_all=None to allow mapping also when reading eloutdet, where outputs are provided also at element nodes locations.

with

nodes_all = binout.read(['eloutdet', element_type, 'nqt'])

since eloutdet output (stress, strain) is in the form of EID_n_IPT_i, ..., EID_n_NODE_j - then simply add some extra headers for the nodes.

def build_mapping(ids_all: np.ndarray, nip_all: np.ndarray, nodes_all=None):
    """
    Construct a mapping array that expands each element ID into one label per
    integration point. Each label is formatted as "element_id/iptN", where N
    is the integration-point index for that element.

    Parameters
    ----------
    ids_all : np.ndarray
        2D array of shape (n_timesteps, n_elements). Each row contains the
        element IDs for that timestep.
    nip_all : np.ndarray
        2D array of shape (n_timesteps, n_elements) specifying the number of
        integration points per element. If None, empty, or not 2D, all elements
        are assumed to have one integration point.
    nodes_all : np.ndarray (optional)
        2D array of shape (n_timesteps, n_elements) specifying the number of
        nodes per element. If None, no nodes are included in the mapping.

    Returns
    -------
    mapping : np.ndarray
        2D object array of shape (n_timesteps, nip_total), where nip_total is the
        total number of integration points for the first timestep. Each entry is
        a string of the form "element_id/iptN". Rows are padded with zeros if a
        timestep contains fewer total integration points.
    """

    ids_all = ids_all.astype(int)
    nip_all = nip_all.astype(int)

    if nodes_all is not None:
        nodes_all = nodes_all.astype(int)
        nodes_total = np.sum(nodes_all[0, :])
    else:
        nodes_total = 0

    if nip_all is None or nip_all.size == 0 or nip_all.ndim < 2:
        # => All elements have 1 integration point
        # Build a new nip_all with the same shape as ids_all but filled with 1s
        nip_all = np.ones_like(ids_all, dtype=int)

    nip_total = np.sum(nip_all[0, :])

    n_timesteps = ids_all.shape[0]

    n_rows = n_timesteps
    n_cols = nip_total + nodes_total

    mapping = np.zeros((n_rows, n_cols), dtype=object)

    for i in range(n_timesteps):
        start_idx = 0
        for idx, eid in enumerate(ids_all[i, :]):
            new_mapping = [f"{eid}/ipt{j+1}" for j in range(nip_all[i, idx])]
            if nodes_all is not None:
                new_mapping += [f"{eid}/N{j+1}" for j in range(nodes_all[i, idx])]
            new_mapping = np.array(new_mapping)
            mapping[i, start_idx:start_idx+new_mapping.shape[0]] = new_mapping
            start_idx += new_mapping.shape[0]

    return mapping

def align_data(data: np.ndarray, mapping: np.ndarray) -> np.ndarray:
    """
    Reorder `data` so that each column index corresponds to the same label
    as in mapping[0, :], for all timesteps.

    Assumes:
    - data.shape == mapping.shape
    - Each row of `mapping` contains the same labels as row 0, plus possibly
      padded zeros (0 or '0') in some columns.
    """
    assert data.shape == mapping.shape, "Data and mapping must have the same shape"
    n_timesteps, n_cols = data.shape

    # Canonical labels from first row (assume no padding in row 0)
    headers = mapping[0, :]

    # Map label -> canonical column index, ignoring any zero-like padding just in case
    label_to_idx = {
        label: j
        for j, label in enumerate(headers)
        if not (label == 0 or label == "0" or label is None)
    }

    # Initialize aligned data with zeros (same behavior as your original)
    new_data = np.zeros((n_timesteps, n_cols), dtype=np.float64)

    for t in range(n_timesteps):
        row_labels = mapping[t, :]

        for src_col, label in enumerate(row_labels):
            # Skip padding: if label not in dict, we leave new_data[t, *] as 0
            dest_col = label_to_idx.get(label, None)
            if dest_col is None:
                continue

            new_data[t, dest_col] = data[t, src_col]

    return new_data

[re-felix]

This looks really good, @marcodesiderio. Thanks for sharing.
Just to clarify, this would still work if the one possible solution brought by @BerSchneider (replicating the numpy 1.x behaviour by using np.array(data, dtype=object)) would be implemented into lasso, right?

[marcodesiderio]

I believe it should but would need to be tested - I am not sure how the missing values would be treated in that case (and therefore how the 'ids' output would look like).