code: def test_mapmri_signal_fitting(radial_order=6):
      gtab = get_gtab_taiwan_dsi()
      l1, l2, l3 = [0.0015, 0.0003, 0.0003]
      S, _ = generate_signal_crossing(gtab, l1, l2, l3)
  
      mapm = MapmriModel(
          gtab, radial_order=radial_order, laplacian_weighting=0.02
      )
      mapfit = mapm.fit(S)
      S_reconst = mapfit.predict(gtab, 1.0)
  
      # test the signal reconstruction
      S = S / S[0]
      nmse_signal = np.sqrt(np.sum((S - S_reconst) ** 2)) / (S.sum())
      assert_almost_equal(nmse_signal, 0.0, 3)
  
      # Test with multidimensional signals:
      mapm = MapmriModel(
          gtab, radial_order=radial_order, laplacian_weighting=0.02
      )
      # Each voxel is identical:
      mapfit = mapm.fit(
          S[:, None, None].T * np.ones((3, 3, 3, S.shape[0]))
      )
  
      # Predict back with an array of ones or a single value:
      for S0 in [S[0], np.ones((3, 3, 3, 203))]:
          S_reconst = mapfit.predict(gtab, S0=S0)
          # test the signal reconstruction for one voxel:
          nmse_signal = np.sqrt(
              np.sum((S - S_reconst[0, 0, 0]) ** 2)
          ) / (S.sum())
          assert_almost_equal(nmse_signal, 0.0, 3)
  
      # do the same for isotropic implementation
      with warnings.catch_warnings():
          warnings.filterwarnings(
              "ignore",
              message=descoteaux07_legacy_msg,
              category=PendingDeprecationWarning,
          )
          mapm = MapmriModel(
              gtab,
              radial_order=radial_order,
              laplacian_weighting=0.0001,
              anisotropic_scaling=False,
          )
      mapfit = mapm.fit(S)
      with warnings.catch_warnings():
          warnings.filterwarnings(
              "ignore",
              message=descoteaux07_legacy_msg,
              category=PendingDeprecationWarning,
          )
          S_reconst = mapfit.predict(gtab, 1.0)
  
      # test the signal reconstruction
      S = S / S[0]
      nmse_signal = np.sqrt(np.sum((S - S_reconst) ** 2)) / (S.sum())
      assert_almost_equal(nmse_signal, 0.0, 3)
  
      # do the same without the positivity constraint:
      with warnings.catch_warnings():
          warnings.filterwarnings(
              "ignore",
              message=descoteaux07_legacy_msg,
              category=PendingDeprecationWarning,
          )
          mapm = MapmriModel(
              gtab,
              radial_order=radial_order,
              laplacian_weighting=0.0001,
              positivity_constraint=False,
              anisotropic_scaling=False,
          )
  
      mapfit = mapm.fit(S)
      with warnings.catch_warnings():
          warnings.filterwarnings(
              "ignore",
              message=descoteaux07_legacy_msg,
              category=PendingDeprecationWarning,
          )
          S_reconst = mapfit.predict(gtab, 1.0)
  
      # test the signal reconstruction
      S = S / S[0]
      nmse_signal = np.sqrt(np.sum((S - S_reconst) ** 2)) / (S.sum())
      assert_almost_equal(nmse_signal, 0.0, 3)
  
      # Repeat with a gtab with big_delta and small_delta:
      gtab.big_delta = 5
      gtab.small_delta = 3
      with warnings.catch_warnings():
          warnings.filterwarnings(
              "ignore",
              message=descoteaux07_legacy_msg,
              category=PendingDeprecationWarning,
          )
          mapm = MapmriModel(
              gtab,
              radial_order=radial_order,
              laplacian_weighting=0.0001,
              positivity_constraint=False,
              anisotropic_scaling=False,
          )
  
      mapfit = mapm.fit(S)
      with warnings.catch_warnings():
          warnings.filterwarnings(
              "ignore",
              message=descoteaux07_legacy_msg,
              category=PendingDeprecationWarning,
          )
          S_reconst = mapfit.predict(gtab, 1.0)
  
      # test the signal reconstruction
      S = S / S[0]
      nmse_signal = np.sqrt(np.sum((S - S_reconst) ** 2)) / (S.sum())
      assert_almost_equal(nmse_signal, 0.0, 3)
  
      if mapmri.have_cvxpy:
          # Positivity constraint and anisotropic scaling:
          with warnings.catch_warnings():
              warnings.filterwarnings(
                  "ignore",
                  message=descoteaux07_legacy_msg,
                  category=PendingDeprecationWarning,
              )
              mapm = MapmriModel(
                  gtab,
                  radial_order=radial_order,
                  laplacian_weighting=0.0001,
                  positivity_constraint=True,
                  anisotropic_scaling=False,
                  pos_radius=2,
              )
  
          mapfit = mapm.fit(S)
          with warnings.catch_warnings():
              warnings.filterwarnings(
                  "ignore",
                  message=descoteaux07_legacy_msg,
                  category=PendingDeprecationWarning,
              )
              S_reconst = mapfit.predict(gtab, 1.0)
  
          # test the signal reconstruction
          S = S / S[0]
          nmse_signal = np.sqrt(np.sum((S - S_reconst) ** 2)) / (
              S.sum()
          )
          assert_almost_equal(nmse_signal, 0.0, 3)
  
          # Positivity constraint and anisotropic scaling:
          with warnings.catch_warnings():
              warnings.filterwarnings(
                  "ignore",
                  message=descoteaux07_legacy_msg,
                  category=PendingDeprecationWarning,
              )
              mapm = MapmriModel(
                  gtab,
                  radial_order=radial_order,
                  laplacian_weighting=None,
                  positivity_constraint=True,
                  anisotropic_scaling=False,
                  pos_radius=2,
              )
  
          mapfit = mapm.fit(S)
          with warnings.catch_warnings():
              warnings.filterwarnings(
                  "ignore",
                  message=descoteaux07_legacy_msg,
                  category=PendingDeprecationWarning,
              )
              S_reconst = mapfit.predict(gtab, 1.0)
  
          # test the signal reconstruction
          S = S / S[0]
          nmse_signal = np.sqrt(np.sum((S - S_reconst) ** 2)) / (
              S.sum()
          )
          assert_almost_equal(nmse_signal, 0.0, 2)
  
imports: import platform
  import time
  from math import factorial
  from scipy.special import gamma
  import scipy.integrate
  import warnings
  import numpy
  from numpy.testing import assert_almost_equal, assert_array_almost_equal, assert_equal, assert_, assert_raises
  import pytest
  from dipy.core.sphere_stats import angular_similarity
  from dipy.core.subdivide_octahedron import create_unit_sphere
  from dipy.data import get_gtab_taiwan_dsi, default_sphere
  from dipy.direction.peaks import peak_directions
  from dipy.reconst.mapmri import MapmriModel, mapmri_index_matrix
  from dipy.reconst import dti, mapmri
  from dipy.reconst.odf import gfa
  from dipy.reconst.tests.test_dsi import sticks_and_ball_dummies
  from dipy.reconst.shm import sh_to_sf, descoteaux07_legacy_msg
  from dipy.sims.voxel import multi_tensor, multi_tensor_pdf, add_noise, single_tensor, cylinders_and_ball_soderman
  from dipy.testing.decorators import set_random_number_generator
label_desc: Use image node_id 8 for calling a global function with example usage: test_mapmri_signal_fitting(radial_order) without return types