{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f1ede658-cfee-4dbe-9c28-961fff6d8933",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-05-13T05:27:58.424027Z",
     "start_time": "2024-05-13T05:27:58.417893Z"
    },
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "MNE: Close-up on forward solutions\n",
    "=================================="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b946432d-6c67-4f32-9c95-0ce7cc1e7280",
   "metadata": {},
   "source": [
    "In previous notebooks we computed the forward solution operator given the raw data, transformed data, BEM model and source space.\n",
    "In this notebook we investigate deeper the forward solution operator. Let us import the forward solution operator we previously computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2fe1ecc4-4231-4cf0-a73d-01433677e9f3",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2024-05-14T08:30:13.101825Z"
    },
    "jupyter": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading forward solution from C:\\Users\\hz3752\\PycharmProjects\\mne_bids_pipeline\\data\\meg\\Sub-0037\\sub-0037-fwd.fif...\n",
      "    Reading a source space...\n",
      "    Computing patch statistics...\n",
      "    Patch information added...\n",
      "    Distance information added...\n",
      "    [done]\n",
      "    Reading a source space...\n",
      "    Computing patch statistics...\n",
      "    Patch information added...\n",
      "    Distance information added...\n",
      "    [done]\n",
      "    2 source spaces read\n",
      "    Desired named matrix (kind = 3523) not available\n",
      "    Read MEG forward solution (5124 sources, 207 channels, free orientations)\n",
      "    Source spaces transformed to the forward solution coordinate frame\n"
     ]
    }
   ],
   "source": [
    "import mne\n",
    "\n",
    "subjects_dir = r'C:\\Users\\hz3752\\PycharmProjects\\mne_bids_pipeline\\data\\meg\\Sub-0037\\sub-0037-fwd.fif'\n",
    "\n",
    "fwd = mne.read_forward_solution(subjects_dir, include=(), exclude=(), ordered=True, verbose=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70dd7b67-9fef-4c8e-a0cd-3023e10ae5cc",
   "metadata": {},
   "source": [
    "Two possible source orientations exist in MNE, documented in `source_ori` in `fwd`:\n",
    "- FIFF.FIFFV_MNE_FIXED_ORI: Fixed orientation mode, the orientation of the source currents is constrained to be perpendicular to the cortical surface. This means that the direction of the current is assumed to be known a priori and is fixed during the inverse solution process. The primary advantage of this approach is its simplicity and reduced computational load. By constraining the direction of the sources, the number of variables in the inverse problem is effectively halved, which can make the solutions more stable under conditions of noisy data or when the data are limited.\n",
    "- When to use fixed orientation? Fixed orientation is often used when a strong prior knowledge of the orientation of the cortical currents exists (such as orientation aligned with the local sulcal or gyral geometry) or when computational simplicity is required.\n",
    "\n",
    "\n",
    "- FIFF.FIFFV_MNE_FREE_ORI: Free Orientation: In contrast, free orientation allows the orientation of the source currents to vary and be estimated as part of the source localization process. This mode does not impose any constraints on the direction of the currents, thereby enabling the estimation of the most probable orientation based on the data itself. This can potentially lead to more accurate reconstructions of the underlying neural activity, particularly in complex scenarios where the orientation of the sources cannot be easily predicted.\n",
    "- When to use free orientation?: Free orientation is particularly useful in research settings where the precise modeling of neural activity is crucial, and where there is enough high-quality data to support the increased complexity of the inverse solution.\n",
    "For the data we red, let us print which source orientation was used:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4c4411e0-a3ce-4c02-b8cf-177db0a41498",
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 (FIFFV_MNE_FREE_ORI)\n"
     ]
    }
   ],
   "source": [
    "print(fwd['source_ori'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "294b6549-8d72-4b56-884e-42c078d1500e",
   "metadata": {},
   "source": [
    "The location of the sources are in `source_rr`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "084d08d2-7fd3-4f53-9504-aa8532e8c3f3",
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.0210264  -0.06770435  0.06570835]\n",
      " [-0.01708252 -0.06216683  0.07929588]\n",
      " [-0.01315619 -0.06514898  0.07275735]\n",
      " ...\n",
      " [ 0.00288066  0.02885441  0.04047859]\n",
      " [ 0.01054338  0.0624362   0.0432792 ]\n",
      " [ 0.01025996  0.08508232  0.04505657]]\n"
     ]
    }
   ],
   "source": [
    "print(fwd['source_rr'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be507e64-a985-46b8-beaa-91d4543936e0",
   "metadata": {},
   "source": [
    "The forward solution operator are in `sol` they are of shape (n_channels, n_sources * n_orientation) where n_orientation is the number of possible orientations of the source (3 orientations for the free-orientation setup and 1 for the fixed orientation setup), let us print the shape of the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a678487d-9ebb-463c-abbc-2b3f32dc30ed",
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(207, 15372)\n"
     ]
    }
   ],
   "source": [
    "print(fwd['sol']['data'].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bfda8a3-e062-4bda-95c3-9449e61f024c",
   "metadata": {},
   "source": [
    "This means that we used 15372/3 = 5124 sources and the solutions are computed for 207 MEG gradiometers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ea8ffd20-a08b-4b07-955e-52aea6f1c9c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening raw data file C:\\Users\\hz3752\\PycharmProjects\\mne_bids_pipeline\\data\\meg\\Sub-0037\\sub-01_01-eyes-closed-raw.fif...\n",
      "    Range : 0 ... 321999 =      0.000 ...   321.999 secs\n",
      "Ready.\n",
      "Reading 0 ... 321999  =      0.000 ...   321.999 secs...\n"
     ]
    }
   ],
   "source": [
    "import mne\n",
    "# Load your raw data\n",
    "raw = mne.io.read_raw_fif(r'C:\\Users\\hz3752\\PycharmProjects\\mne_bids_pipeline\\data\\meg\\Sub-0037\\sub-01_01-eyes-closed-raw.fif', preload=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "856dd1b1-02cc-4e16-b517-db3df13d2c03",
   "metadata": {
    "jupyter": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using up to 1610 segments\n",
      "Number of samples used : 322000\n",
      "[done]\n"
     ]
    }
   ],
   "source": [
    "#Add covariance matrix computation\n",
    "\n",
    "#cov = mne.compute_covariance\n",
    "\n",
    "cov = mne.compute_raw_covariance(raw)\n",
    "cov.save(r'C:\\Users\\hz3752\\PycharmProjects\\mne_bids_pipeline\\data\\meg\\Sub-0037\\sub-01-cov.fif')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "161a90e3-108d-409a-9032-45d0e6218616",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converting forward solution to surface orientation\n",
      "    Average patch normals will be employed in the rotation to the local surface coordinates....\n",
      "    Converting to surface-based source orientations...\n",
      "    [done]\n",
      "Computing inverse operator with 207 channels.\n",
      "    207 out of 207 channels remain after picking\n",
      "Selected 207 channels\n",
      "Applying loose dipole orientations to surface source spaces: 0.2\n",
      "Whitening the forward solution.\n",
      "Computing rank from covariance with rank=None\n",
      "    Using tolerance 8.2e-13 (2.2e-16 eps * 207 dim * 18  max singular value)\n",
      "    Estimated rank (mag): 207\n",
      "    MAG: rank 207 computed from 207 data channels with 0 projectors\n",
      "    Setting small MAG eigenvalues to zero (without PCA)\n",
      "Creating the source covariance matrix\n",
      "Adjusting source covariance matrix.\n",
      "Computing SVD of whitened and weighted lead field matrix.\n",
      "    largest singular value = 14.1305\n",
      "    scaling factor to adjust the trace = 7.29474e+27 (nchan = 207 nzero = 3)\n"
     ]
    }
   ],
   "source": [
    " inv = mne.minimum_norm.make_inverse_operator(raw.info, fwd, cov, depth=None, loose='auto', fixed=False) #fixed=False: Ignoring dipole direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de52ad8f-3a30-4fe6-a897-94e9eb22efc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "subj='sub-01'\n",
    "        #--------------------------STCs--------------------------------#\n",
    "\n",
    "print '%s: Creating STCs...'%subj\n",
    "        os.makedirs('STC/%s' %subj)\n",
    "        for ev in evoked:\n",
    "            stc = mne.minimum_norm.apply_inverse(ev, inv, lambda2=lambda2, method='dSPM')\n",
    "            # mophing stcs to the fsaverage using precomputed matrix method:\n",
    "            vertices_to = mne.grade_to_vertices('fsaverage', grade=4, subjects_dir=subjects_dir) #fsaverage's source space\n",
    "            morph_mat = mne.compute_morph_matrix(subject_from=subj, subject_to='fsaverage', vertices_from=stc.vertices, vertices_to=vertices_to, subjects_dir=subjects_dir)\n",
    "            stc_morph = mne.morph_data_precomputed(subject_from=subj, subject_to='fsaverage', stc_from=stc, vertices_to=vertices_to, morph_mat=morph_mat)\n",
    "            stc_morph.save('STC/%s/%s_%s_dSPM' %(subj,subj,ev.comment))\n",
    "            del stc, stc_morph\n",
    "        print '>> DONE CREATING STCS FOR SUBJ=%s'%subj\n",
    "        print '-----------------------------------------\\n'\n",
    "\n",
    "        #deleting variables\n",
    "        del epochs_rej, evoked, info, trans, src, fwd, cov, inv"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}