Inconsistent scalar vs. array types for prior / likelihood / posterior in Link cause to_inference_data to fail

[marie-cloet2000]

Description

I was adapting the code from this example notebook to my own problem. The setup follows the notebook closely; I only modified some hyperparameters, the prior, and the forward models.
When calling

idata = tda.to_inference_data(my_chain, level='fine', burnin=burnin)

I encounter the following error:

File ".../tinyDA/diagnostics.py", line 153, in <lambda>
    getattribute = lambda link, attribute: np.array([link.prior, link.likelihood, link.posterior])
                                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

Diagnosis

After inspecting the values of the Link object, I found that:

• link.prior is a NumPy array of shape (1,)
• link.likelihood is a Python float
• link.posterior is a NumPy array of shape (1,)

This mixture of types causes np.array([prior, likelihood, posterior]) to fail.

Looking at the Link class in tinyDA/link.py, both prior and likelihood are documented as floats, and posterior is computed as their sum. However, no normalization or type checking is applied in Link.__init__, so prior may silently become a 1‑element array depending on the prior distribution implementation.

Proposed fix

I plan to submit a PR that normalizes prior and likelihood to Python scalars (e.g. via explicit scalar extraction) in the Link initializer, so that downstream diagnostics can safely assume consistent scalar types.

Please let me know if you have any suggestions regarding the implementation approach.

By the way, I have no clue why the prior values are NumPy arrays. Could this be an issue as well?

genAI statement: as a non-native English speaker, I used Microsoft Copilot to polish my Issue Description.

---

Appendix

Where my code differs from the example notebook:

# prior distribution
my_prior = stats.uniform(loc=-1e2, scale=102)

class LinearODEModel:
    def __init__(self, datapoints, method, y0):
        # set the timesteps, where we are collecting the model output
        self.datapoints = datapoints

        # set the span of the integration
        self.t_span = [0, self.datapoints[-1]]

        # set the integration method
        self.method = method

        # constant parameters 
        self.p = -1 * np.ones(4)
        self.A = np.eye(4)
        self.q = np.ones(4)
        self.alpha = 1
        self.epsilon = 0.1

        # set the initial condition
        self.y0 = y0

    def dydx(self, t, y, p1):
        """ODE system to generate data"""
        # modify the parameter vector
        _p = self.p.copy()
        _p[0] = p1
        
        # compute rhs of ODE
        rhs = np.hstack((self.alpha*y[0] + _p @ y[1:], (self.A @ self.q * y[0] - y[1:]) / self.epsilon))

        return rhs

    def __call__(self, parameters):
        p1 = parameters[0]
        
        # solve the initial value problem
        self.y = solve_ivp(lambda t, y: self.dydx(t, y, p1), self.t_span, self.y0, method=self.method, t_eval=self.datapoints)

        # return the results, only if the integration succeeded
        if self.y.success:
            return self.y.y[0,:], True
        else:
            return np.zeros_like(self.y.y[0,:]), False