Gridsearch nkg

invokers/run_gridsearch.py

@@ -1,95 +1,130 @@
 from pathlib import Path
 import argparse
+import time
+
+from kymata.datasets.data_root import data_root_path
 from kymata.gridsearch.plain import do_gridsearch
 from kymata.io.functions import load_function
 from kymata.io.mne import load_emeg_pack
+from kymata.io.nkg import save_expression_set
+from kymata.plot.plot import expression_plot
+
+_default_output_dir = Path(data_root_path(), "output")
 
 
 def main():
 
+    _default_output_dir.mkdir(exist_ok=True, parents=False)
+
     parser = argparse.ArgumentParser(description='Gridsearch Params')
-    parser.add_argument('--emeg_sample_rate', type=int, default=1000,
+    parser.add_argument('--emeg-sample-rate', type=int, default=1000,
                         help='sampling rate of the emeg machine (not implemented yet)')
-    parser.add_argument('--snr', type=float, default=3,
-                        help='inverse solution snr')
-    parser.add_argument('--downsample_rate', type=int, default=5,
-                        help='downsample_rate')
-    parser.add_argument('--base_dir', type=str, default="/imaging/projects/cbu/kymata/data/dataset_4-english-narratives/",
-                        help='base data directory')
-    parser.add_argument('--data_path', type=str, default="intrim_preprocessing_files/3_trialwise_sensorspace/evoked_data",
-                        help='data path after base dir')
-    parser.add_argument('--function_path', type=str, default="predicted_function_contours/GMSloudness/stimulisig",
-                        help='snr')
-    parser.add_argument('--function_name', type=str, default="d_IL2",
-                        help='function name in stimulisig')
-    parser.add_argument('--emeg_file', type=str, default="participant_01-ave",
-                        help='emeg_file_name')
-    parser.add_argument('--ave_mode', type=str, default="ave",
+    parser.add_argument('--snr', type=float, default=3, help='inverse solution snr')
+    parser.add_argument('--downsample-rate', type=int, default=5, help='downsample_rate')
+    parser.add_argument('--base-dir', type=str, default='/imaging/projects/cbu/kymata/data/dataset_4-english-narratives/', help='base data directory')
+    parser.add_argument('--data-path', type=str, default='intrim_preprocessing_files/3_trialwise_sensorspace/evoked_data', help='data path after base dir')
+    parser.add_argument('--function-path', type=str, default='predicted_function_contours/GMSloudness/stimulisig', help='location of function stimulisig')
+    parser.add_argument('--save-expression-set-location', type=Path, default=Path(_default_output_dir),
+                        help="Save the results of the gridsearch into an ExpressionSet .nkg file")
+    parser.add_argument('--save-plot-location', type=Path, default=Path(_default_output_dir),
+                        help="Save an expression plots, and other plots, in this location")
+    parser.add_argument('--overwrite', action="store_true", help="Silently overwrite existing files.")
+    parser.add_argument('--function-name', type=str, default="IL", help='function name in stimulisig')
+    parser.add_argument('--emeg-file', type=str, default="participant_01-ave", help='emeg_file_name')
+    parser.add_argument('--ave-mode', type=str, default="ave",
                         help='either ave or add, either average over the list of repetitions or treat them as extra data')
-    parser.add_argument('--inverse_operator', type=str, default="intrim_preprocessing_files/4_hexel_current_reconstruction/inverse-operators",
-                        help='inverse solution path')
-    parser.add_argument('--seconds_per_split', type=float, default=0.5,
+    parser.add_argument('--inverse-operator-dir', type=str, default=None, help='inverse solution path')
+    parser.add_argument('--inverse-operator-name', type=str, default="participant_01_ico5-3L-loose02-cps-nodepth-fusion.fif",
+                        help='inverse solution name')
+    parser.add_argument('--seconds-per-split', type=float, default=0.5,
                         help='seconds in each split of the recording, also maximum range of latencies being checked')
-    parser.add_argument('--n_splits', type=int, default=800,
+    parser.add_argument('--n-splits', type=int, default=800,
                         help='number of splits to split the recording into, (set to 400/seconds_per_split for full file)')
-    parser.add_argument('--n_derangements', type=int, default=1,
-                        help='inverse solution snr')
-    parser.add_argument('--start_latency', type=float, default=-100,
+    parser.add_argument('--n-derangements', type=int, default=1,
+                        help='number of deragements for the null distribution')
+    parser.add_argument('--start-latency', type=float, default=-100,
                         help='earliest latency to check in cross correlation')
-    parser.add_argument('--emeg_t_start', type=float, default=-200,
+    parser.add_argument('--emeg-t-start', type=float, default=-200,
                         help='start of the emeg evoked files relative to the start of the function')
-    parser.add_argument('--audio_shift_correction', type=float, default=0.000_537_5,
+    parser.add_argument('--audio-shift-correction', type=float, default=0.000_537_5,
                         help='audio shift correction, for every second of function, add this number of seconds (to the start of the emeg split) per seconds of emeg seen')
     args = parser.parse_args()
     args.base_dir = Path(args.base_dir)
 
-
     emeg_dir = Path(args.base_dir, args.data_path)
     emeg_paths = [Path(emeg_dir, args.emeg_file)]
 
-    participants = ['participant_01',
+    participants = ['pilot_01',
+                    'pilot_02',
+                    'participant_01',
                     'participant_01b',
                     'participant_02',
                     'participant_03',
                     'participant_04',
                     'participant_05',
-                    'pilot_01',
-                    'pilot_02']
+                    'participant_07',
+                    'participant_08',
+                    'participant_09',
+                    'participant_10',
+                    'participant_11',
+                    'participant_12',
+                    'participant_13',
+                    'participant_14',
+                    'participant_15',
+                    'participant_16',
+                    'participant_17'
+                    ]
 
     reps = [f'_rep{i}' for i in range(8)] + ['-ave']
 
     # emeg_paths = [Path(emeg_dir, p + r) for p in participants[:2] for r in reps[-1:]]
 
-    inverse_operator = Path(args.base_dir, args.inverse_operator, f"{participants[0]}_ico5-3L-loose02-cps-nodepth.fif")
+    start = time.time()
+
+    if args.inverse_operator_dir is None:
+        inverse_operator = None
+    else:
+        inverse_operator = Path(args.base_dir, args.inverse_operator_dir, args.inverse_operator_name)
 
     # Load data
-    emeg, ch_names = load_emeg_pack(emeg_paths,
-                                    need_names=False,
-                                    ave_mode=args.ave_mode,
-                                    inverse_operator=None, #inverse_operator, # set to None/inverse_operator if you want to run on sensor space/source space
-                                    p_tshift=None,
-                                    snr=args.snr)
+    emeg_values, ch_names = load_emeg_pack(emeg_paths,
+                                           need_names=True,
+                                           ave_mode=args.ave_mode,
+                                           inverse_operator=inverse_operator,
+                                           p_tshift=None,
+                                           snr=args.snr)
 
     func = load_function(Path(args.base_dir, args.function_path),
                          func_name=args.function_name,
                          bruce_neurons=(5, 10))
     func = func.downsampled(args.downsample_rate)
 
+    channel_space = "source" if inverse_operator is not None else "sensor"
+
     es = do_gridsearch(
-        emeg_values=emeg,
-        sensor_names=ch_names,
+        emeg_values=emeg_values,
+        channel_names=ch_names,
+        channel_space=channel_space,
         function=func,
         seconds_per_split=args.seconds_per_split,
         n_derangements=args.n_derangements,
         n_splits=args.n_splits,
         start_latency=args.start_latency,
+        plot_location=args.save_plot_location,
         emeg_t_start=args.emeg_t_start,
         emeg_sample_rate=args.emeg_sample_rate,
         audio_shift_correction=args.audio_shift_correction,
         ave_mode=args.ave_mode,
+        overwrite=args.overwrite,
     )
 
-    # expression_plot(es)
+    if args.save_expression_set_location is not None:
+        save_expression_set(es, to_path_or_file = Path(args.save_expression_set_location, args.function_name + '_gridsearch.nkg'), overwrite=args.overwrite)
+
+    expression_plot(es, paired_axes=channel_space == "source", save_to=Path(args.save_plot_location, args.function_name + '_gridsearch.png'), overwrite=args.overwrite)
+
+    print(f'Time taken for code to run: {time.time() - start:.4f}')
+
 
 if __name__ == '__main__':
     main()

kymata/config/dataset4.yaml

@@ -31,6 +31,7 @@ supress_excessive_plots_and_prompts: True
 # Inverse operator
 eeg: True
 meg: True
+inverse_operator: "intrim_preprocessing_files/4_hexel_current_reconstruction/inverse-operators"
 
 # Method to estimate noise covariance matrix
 cov_method: 'grand_ave'            # grand_ave | empty_room | run_start

kymata/datasets/data_root.py

@@ -0,0 +1,73 @@
+from os import getcwd, getenv
+from pathlib import Path
+from typing import Optional
+
+from kymata.io.file import path_type
+
+
+_DATA_PATH_ENVIRONMENT_VAR_NAME = "KYMATA_DATA_ROOT"
+DATA_DIR_NAME = "kymata-toolbox-data"
+
+# Places downloaded datasets could go, in order of preference
+_preferred_default_data_locations = [
+    Path(Path(__file__).parent.parent.parent),  # kymata/../data_dir (next to kymata dir)
+    Path(getcwd()),                             # <cwd>/data_dir
+    Path(Path.home(), "Documents"),             # ~/Documents/data_dir
+    Path(Path.home()),                          # ~/data_dir
+]
+
+
+def data_root_path(data_root: Optional[path_type] = None) -> Path:
+
+    # Check if the data root has been specified
+
+    # Might be in an environmental variable
+    if data_root is None:
+        data_root: path_type | None = getenv(_DATA_PATH_ENVIRONMENT_VAR_NAME, default=None)
+
+    # Might have been supplied as an argument
+    if data_root is not None:
+        if isinstance(data_root, str):
+            data_root = Path(data_root)
+        # Data root specified
+        if not data_root.exists():
+            raise FileNotFoundError(f"data_root {str(data_root)} specified but does not exist")
+        if not data_root.is_dir():
+            raise NotADirectoryError(f"Please specify a directory ({str(data_root)} is not a directory)")
+
+        return data_root
+
+    else:
+        # Data root not specified
+
+        # Check if the data root already exists
+        for loc in _preferred_default_data_locations:
+            if (here := Path(loc, DATA_DIR_NAME)).exists():
+                data_root = here
+                break
+
+        # If not, attempt to create it
+        if data_root is None:
+            here: Path | None = None
+            for loc in _preferred_default_data_locations:
+                here = Path(loc, DATA_DIR_NAME)
+                try:
+                    here.mkdir()
+                    break
+                # If it fails for sensible reasons, no sweat, we'll fall through to the next option
+                except (FileNotFoundError, OSError):
+                    # Parent didn't exist, not writeable, etc
+                    pass
+            # Did we make it?
+            if here is not None and here.exists():
+                data_root = here
+            else:
+                raise FileNotFoundError("Failed to create data root directory")
+
+        # Data root location has been derived, rather than prespecified, so feed that back to the user to avoid a
+        # different location somehow being derived next time
+        print(f"Data root set at {str(data_root)}.")
+        print(f"Consider setting this as environmental variable {_DATA_PATH_ENVIRONMENT_VAR_NAME} to ensure it's reused"
+              f" next time.")
+        print(f"Hint: $> {_DATA_PATH_ENVIRONMENT_VAR_NAME}=\"{str(data_root)}\"")
+        return data_root

kymata/datasets/sample.py

@@ -1,23 +1,15 @@
 from abc import ABC, abstractmethod
-from os import getenv, getcwd, remove, rmdir
+from os import remove, rmdir
 from pathlib import Path
 from typing import Optional
 from urllib import request
 
+from kymata.datasets.data_root import data_root_path
 from kymata.entities.expression import HexelExpressionSet, SensorExpressionSet
 from kymata.io.file import path_type
 from kymata.io.nkg import load_expression_set
 
-_DATA_PATH_ENVIRONMENT_VAR_NAME = "KYMATA_DATA_ROOT"
-_DATA_DIR_NAME = "kymata-toolbox-data/tutorial_nkg_data"
-
-# Places downloaded datasets could go, in order of preference
-_preferred_default_data_locations = [
-    Path(Path(__file__).parent.parent.parent),  # kymata/../data_dir (next to kymata dir)
-    Path(getcwd()),                             # <cwd>/data_dir
-    Path(Path.home(), "Documents"),             # ~/Documents/data_dir
-    Path(Path.home()),                          # ~/data_dir
-]
+_SAMPLE_DATA_DIR_NAME = "tutorial_nkg_data"
 
 
 class SampleDataset(ABC):
@@ -36,9 +28,13 @@ def __init__(self,
                  download: bool):
         self.name: str = name
         self.filenames: list[str] = filenames
-        self.data_root: Path = data_root_path(data_root)
+        self.data_root: Path = Path(data_root_path(data_root), _SAMPLE_DATA_DIR_NAME)
         self.remote_root: str = remote_root
 
+        # Create the default location, if it's being used
+        if data_root is None:
+            self.data_root.mkdir(exist_ok=True)
+
         if download:
             self.download()
 
@@ -141,62 +137,6 @@ def to_expressionset(self) -> SensorExpressionSet:
         return es
 
 
-def data_root_path(data_root: Optional[path_type] = None) -> Path:
-
-    # Check if the data root has been specified
-
-    # Might be in an environmental variable
-    if data_root is None:
-        data_root: path_type | None = getenv(_DATA_PATH_ENVIRONMENT_VAR_NAME, default=None)
-
-    # Might have been supplied as an argument
-    if data_root is not None:
-        if isinstance(data_root, str):
-            data_root = Path(data_root)
-        # Data root specified
-        if not data_root.exists():
-            raise FileNotFoundError(f"data_root {str(data_root)} specified but does not exist")
-        if not data_root.is_dir():
-            raise NotADirectoryError(f"Please specify a directory ({str(data_root)} is not a directory)")
-
-        return data_root
-
-    else:
-        # Data root not specified
-
-        # Check if the data root already exists
-        for loc in _preferred_default_data_locations:
-            if (here := Path(loc, _DATA_DIR_NAME)).exists():
-                data_root = here
-                break
-
-        # If not, attempt to create it
-        if data_root is None:
-            here: Path | None = None
-            for loc in _preferred_default_data_locations:
-                here = Path(loc, _DATA_DIR_NAME)
-                try:
-                    here.mkdir()
-                    break
-                # If it fails for sensible reasons, no sweat, we'll fall through to the next option
-                except (FileNotFoundError, OSError):
-                    # Parent didn't exist, not writeable, etc
-                    pass
-            # Did we make it?
-            if here is not None and here.exists():
-                data_root = here
-            else:
-                raise FileNotFoundError("Failed to create data root directory")
-
-        # Data root location has been derived, rather than prespecified, so feed that back to the user to avoid a
-        # different location somehow being derived next time
-        print(f"Data root set at {str(data_root)}.")
-        print(f"Consider setting this as environmental variable {_DATA_PATH_ENVIRONMENT_VAR_NAME} to ensure it's reused"
-              f" next time.")
-        print(f"Hint: $> {_DATA_PATH_ENVIRONMENT_VAR_NAME}=\"{str(data_root)}\"")
-        return data_root
-
-
 def delete_dataset(local_dataset: SampleDataset):
     # Make sure it's not silent
     print(f"Deleting dataset {local_dataset.name}")

kymata/entities/expression.py

@@ -92,8 +92,8 @@ def __init__(self,
                 data = data[i]
                 data = self._init_prep_data(data)
                 # Check validity of input data dimensions
-                assert len(channels) == data.shape[0], f"{channel_coord_name} mismatch for {f}"
-                assert len(latencies) == data.shape[1], f"Latencies mismatch for {f}"
+                assert len(channels) == data.shape[0], f"{channel_coord_name} mismatch for {f}: {len(channels)} {channel_coord_name} versus data shape {data.shape}"
+                assert len(latencies) == data.shape[1], f"Latencies mismatch for {f}: {len(latencies)} latencies versus data shape {data.shape}"
                 dataset_dict[layer] = DataArray(
                     data=data,
                     dims=self._dims,
@@ -403,6 +403,9 @@ def best_functions(self) -> DataFrame:
         return super()._best_functions_for_layer(LAYER_SCALP)
 
 
+log_base = 10
+
+
 def p_to_logp(arraylike: ArrayLike) -> ArrayLike:
     """The one-stop-shop for converting from p-values to log p-values."""
     return log10(arraylike)