code: def run_vanilla_fpi(
      A,
      obs,
      num_obs,
      num_states,
      prior=None,
      num_iter=10,
      dF=1.0,
      dF_tol=0.001,
  ):
      """
      Update marginal posterior beliefs over hidden states using mean-field variational inference, via
      fixed point iteration.
  
      Parameters
      ----------
      A: ``numpy.ndarray`` of dtype object
          Sensory likelihood mapping or 'observation model', mapping from hidden states to observations. Each element ``A[m]`` of
          stores an ``np.ndarray`` multidimensional array for observation modality ``m``, whose entries ``A[m][i, j, k, ...]`` store
          the probability of observation level ``i`` given hidden state levels ``j, k, ...``
      obs: numpy 1D array or numpy ndarray of dtype object
          The observation (generated by the environment). If single modality, this should be a 1D ``np.ndarray``
          (one-hot vector representation). If multi-modality, this should be ``np.ndarray`` of dtype object whose entries are 1D one-hot vectors.
      num_obs: list of ints
          List of dimensionalities of each observation modality
      num_states: list of ints
          List of dimensionalities of each observation modality
      prior: numpy ndarray of dtype object, default None
          Prior over hidden states. If absent, prior is set to be the log uniform distribution over hidden states (identical to the
          initialisation of the posterior)
      num_iter: int, default 10
          Number of variational fixed-point iterations to run until convergence.
      dF: float, default 1.0
          Initial free energy gradient (dF/dt) before updating in the course of gradient descent.
      dF_tol: float, default 0.001
          Threshold value of the time derivative of the variational free energy (dF/dt), to be checked at
          each iteration. If dF <= dF_tol, the iterations are halted pre-emptively and the final
          marginal posterior belief(s) is(are) returned
  
      Returns
      ----------
      qs: numpy 1D array, numpy ndarray of dtype object, optional
          Marginal posterior beliefs over hidden states at current timepoint
      """
  
      # get model dimensions
      n_modalities = len(num_obs)
      n_factors = len(num_states)
  
      """
      =========== Step 1 ===========
          Loop over the observation modalities and use assumption of independence 
          among observation modalitiesto multiply each modality-specific likelihood 
          onto a single joint likelihood over hidden factors [size num_states]
      """
  
      likelihood = get_joint_likelihood(A, obs, num_states)
  
      likelihood = spm_log_single(likelihood)
  
      """
      =========== Step 2 ===========
          Create a flat posterior (and prior if necessary)
      """
  
      qs = np.empty(n_factors, dtype=object)
      for factor in range(n_factors):
          qs[factor] = np.ones(num_states[factor]) / num_states[factor]
  
      """
      If prior is not provided, initialise prior to be identical to posterior 
      (namely, a flat categorical distribution). Take the logarithm of it (required for 
      FPI algorithm below).
      """
      if prior is None:
          prior = obj_array_uniform(num_states)
  
      prior = spm_log_obj_array(prior)  # log the prior
  
      """
      =========== Step 3 ===========
          Initialize initial free energy
      """
      prev_vfe = calc_free_energy(qs, prior, n_factors)
  
      """
      =========== Step 4 ===========
          If we have a single factor, we can just add prior and likelihood because there is a unique FE minimum that can reached instantaneously,
          otherwise we run fixed point iteration
      """
  
      if n_factors == 1:
          qL = spm_dot(likelihood, qs, [0])
  
          return to_obj_array(softmax(qL + prior[0]))
  
      else:
          """
          =========== Step 5 ===========
          Run the FPI scheme
          """
  
          curr_iter = 0
          while curr_iter < num_iter and dF >= dF_tol:
              # Initialise variational free energy
              vfe = 0
  
              # arg_list = [likelihood, list(range(n_factors))]
              # arg_list = arg_list + list(chain(*([qs_i,[i]] for i, qs_i in enumerate(qs)))) + [list(range(n_factors))]
              # LL_tensor = np.einsum(*arg_list)
  
              qs_all = qs[0]
              for factor in range(n_factors - 1):
                  qs_all = qs_all[..., None] * qs[factor + 1]
              LL_tensor = likelihood * qs_all
  
              for factor, qs_i in enumerate(qs):
                  # qL = np.einsum(LL_tensor, list(range(n_factors)), 1.0/qs_i, [factor], [factor])
                  qL = (
                      np.einsum(
                          LL_tensor, list(range(n_factors)), [factor]
                      )
                      / qs_i
                  )
                  qs[factor] = softmax(qL + prior[factor])
  
              # List of orders in which marginal posteriors are sequentially multiplied into the joint likelihood:
              # First order loops over factors starting at index = 0, second order goes in reverse
              # factor_orders = [range(n_factors), range((n_factors - 1), -1, -1)]
  
              # iteratively marginalize out each posterior marginal from the joint log-likelihood
              # except for the one associated with a given factor
              # for factor_order in factor_orders:
              #     for factor in factor_order:
              #         qL = spm_dot(likelihood, qs, [factor])
              #         qs[factor] = softmax(qL + prior[factor])
  
              # calculate new free energy
              vfe = calc_free_energy(qs, prior, n_factors, likelihood)
  
              # stopping condition - time derivative of free energy
              dF = np.abs(prev_vfe - vfe)
              prev_vfe = vfe
  
              curr_iter += 1
  
          return qs
  
imports: import numpy
  from pymdp.maths import spm_dot, get_joint_likelihood, softmax, calc_free_energy, spm_log_single, spm_log_obj_array
  from pymdp.utils import to_obj_array, obj_array_uniform
  from itertools import chain
label_desc: Use image node_id 1 for calling a global function with example usage: run_vanilla_fpi(A, obs, num_obs, num_states, prior, num_iter, dF, dF_tol) and returns: to_obj_array, qs