code: def run_glm(
      Y,
      X,
      noise_model="ar1",
      bins=100,
      n_jobs=1,
      verbose=0,
      random_state=None,
  ):
      """:term:`GLM` fit for an :term:`fMRI` data matrix.
  
      Parameters
      ----------
      Y : array of shape (n_time_points, n_voxels)
          The :term:`fMRI` data.
  
      X : array of shape (n_time_points, n_regressors)
          The design matrix.
  
      noise_model : {'ar(N)', 'ols'}, default='ar1'
          The temporal variance model.
          To specify the order of an autoregressive model place the
          order after the characters `ar`, for example to specify a third order
          model use `ar3`.
  
      bins : int, default=100
          Maximum number of discrete bins for the AR coef histogram.
          If an autoregressive model with order greater than one is specified
          then adaptive quantification is performed and the coefficients
          will be clustered via K-means with `bins` number of clusters.
  
      n_jobs : int, default=1
          The number of CPUs to use to do the computation. -1 means
          'all CPUs'.
  
      verbose : int, default=0
          The verbosity level.
  
      random_state : int or numpy.random.RandomState, default=None
          Random state seed to sklearn.cluster.KMeans for autoregressive models
          of order at least 2 ('ar(N)' with n >= 2).
  
          .. versionadded:: 0.9.1
  
      Returns
      -------
      labels : array of shape (n_voxels,),
          A map of values on voxels used to identify the corresponding model.
  
      results : dict,
          Keys correspond to the different labels values
          values are RegressionResults instances corresponding to the voxels.
  
      """
      acceptable_noise_models = ["ols", "arN"]
      if (noise_model[:2] != "ar") and (noise_model != "ols"):
          raise ValueError(
              f"Acceptable noise models are {acceptable_noise_models}. "
              f"You provided 'noise_model={noise_model}'."
          )
      if Y.shape[0] != X.shape[0]:
          raise ValueError(
              "The number of rows of Y "
              "should match the number of rows of X.\n"
              f"You provided X with shape {X.shape} "
              f"and Y with shape {Y.shape}."
          )
  
      # Create the model
      ols_result = OLSModel(X).fit(Y)
  
      if noise_model[:2] == "ar":
          err_msg = (
              "AR order must be a positive integer specified as arN, "
              "where N is an integer. E.g. ar3. "
              f"You provided {noise_model}."
          )
          try:
              ar_order = int(noise_model[2:])
          except ValueError:
              raise ValueError(err_msg)
  
          # compute the AR coefficients
          ar_coef_ = _yule_walker(ols_result.residuals.T, ar_order)
          del ols_result
          if len(ar_coef_[0]) == 1:
              ar_coef_ = ar_coef_[:, 0]
  
          # Either bin the AR1 coefs or cluster ARN coefs
          if ar_order == 1:
              for idx in range(len(ar_coef_)):
                  ar_coef_[idx] = (
                      (ar_coef_[idx] * bins).astype(int) * 1.0 / bins
                  )
              labels = np.array([str(val) for val in ar_coef_])
          else:  # AR(N>1) case
              n_clusters = np.min([bins, Y.shape[1]])
              kmeans = KMeans(
                  n_clusters=n_clusters,
                  n_init=10,
                  random_state=random_state,
              ).fit(ar_coef_)
              ar_coef_ = kmeans.cluster_centers_[kmeans.labels_]
  
              # Create a set of rounded values for the labels with _ between
              # each coefficient
              cluster_labels = kmeans.cluster_centers_.copy()
              cluster_labels = np.array(
                  [
                      "_".join(map(str, np.round(a, 2)))
                      for a in cluster_labels
                  ]
              )
              # Create labels and coef per voxel
              labels = np.array(
                  [cluster_labels[i] for i in kmeans.labels_]
              )
  
          unique_labels = np.unique(labels)
          results = {}
  
          # Fit the AR model according to current AR(N) estimates
          ar_result = Parallel(n_jobs=n_jobs, verbose=verbose)(
              delayed(_ar_model_fit)(
                  X, ar_coef_[labels == val][0], Y[:, labels == val]
              )
              for val in unique_labels
          )
  
          # Converting the key to a string is required for AR(N>1) cases
          for val, result in zip(unique_labels, ar_result):
              results[val] = result
          del unique_labels
          del ar_result
  
      else:
          labels = np.zeros(Y.shape[1])
          results = {0.0: ols_result}
  
      return labels, results
  
imports: from __future__ import annotations
  import csv
  import glob
  import os
  import sys
  import time
  from pathlib import Path
  from warnings import warn
  import numpy
  import pandas
  from joblib import Memory, Parallel, delayed
  from nibabel import Nifti1Image
  from sklearn.base import clone
  from sklearn.cluster import KMeans
  from nilearn._utils import fill_doc, stringify_path
  from nilearn._utils.niimg_conversions import check_niimg
  from nilearn._utils.param_validation import check_run_sample_masks
  from nilearn.glm._base import BaseGLM
  from nilearn.glm.contrasts import compute_fixed_effect_contrast, expression_to_contrast_vector
  from nilearn.glm.first_level.design_matrix import make_first_level_design_matrix
  from nilearn.glm.regression import ARModel, OLSModel, RegressionResults, SimpleRegressionResults
  from nilearn.image import get_data
  from nilearn.interfaces.bids import get_bids_files, parse_bids_filename
  from nilearn.interfaces.bids.query import infer_repetition_time_from_dataset, infer_slice_timing_start_time_from_dataset
  from nilearn.interfaces.bids.utils import bids_entities, check_bids_label
  from nilearn.interfaces.fmriprep.load_confounds import load_confounds
label_desc: Use image node_id 4 for calling a global function with example usage: run_glm(Y, X, noise_model, bins, n_jobs, verbose, random_state) and returns: labels,  results