init gDCNN v0.2

01_gDCNN_label_ica.py

+"""
+gDCNN label_ica
+Authors:
+ - f.boers@fz-juelich.de
+ - j.dammers@fz-juelich.de
+
+Motivation:
+    In contrast to [1], this routine is capable of processing
+    4D-Neuroimaging, CTF and Neuromag data.
+    Performs ICA on chunks of MEG data (2-3 minutes) and apply
+    an automatic labeling (cardiac or ocular) for all
+    ICA components in all chops of data.
+    The information is used for training a generalised version
+    of the DCNN model (gDCNN)
+
+    Remarks - about data:
+    - Data should not be filtered (recommendation).
+    - !It is assumed that bad channels are marked as bad. It is highly recommended
+      to replace bad channels with an interpolated signal.
+    - !MEG system specific noise reduction should be applied (e.g., sss) prior to
+      this analysis
+    - !ECG and EOG channels must be of good quality. Horizontal EOG can be used
+      in addition to vertical EOG.
+    - Raw continuous data should be used (no epochs)
+    - Besides resting state also task data should be used for training
+    - equal numbers of female/male would be nice
+
+    Remarks - about config:
+    - !MEG system specs must be defined in config file
+    - !ECG and EOG channel names must be provided and fit to all data to be processed.
+      In cases of changes in AUX names for different raw files you must use a different config file
+    - for different data sets you may want to use different directories
+
+
+Workflow
+    1. Compute ICA on filtered and down sampled data.
+       Note, cleaning will be applied on unfiltered data later.
+       (or raw data that has been filtered differently)
+    2. Identification of ECG and EOG components using standard methods
+       as implemented in MNE-Python (CTPS and correlation)
+    3. Results will be save in as dict in numpy format. Only ICA components
+       together relevant information about the data (DCNN object) is stored
+    4. Data is used to train the gDCNN model (in Juelich)
+    5. gDCNN-based artifact rejection ist tested against standard methods for all MEG systems
+
+[1] Ahmad Hasasneh, Nikolas Kampel, Praveen Sripad, N. Jon Shah, and Juergen Dammers
+    "Deep Learning Approach for Automatic Classification of Ocular and Cardiac
+    Artifacts in MEG Data"
+    Journal of Engineering, vol. 2018, Article ID 1350692,10 pages, 2018.
+    https://doi.org/10.1155/2018/1350692
+"""
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# import libraries
+#
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+from dcnn_base import DCNN_CONFIG
+from dcnn_main import DCNN
+from dcnn_utils import find_files
+
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# select your config file
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# fnconfig = 'config_4D_CAU.yaml'
+# fnconfig = 'config_4D_INTEXT.yaml'
+
+fnconfig = 'config_MEGIN.yaml'
+# fnconfig = 'config_CTF_Paris.yaml'
+# fnconfig = 'config_CTF_Philly.yaml'
+
+#  --- MacOSX -----
+# fnconfig = 'config_4D_mac.yaml'
+# fnconfig = 'config_MEGIN_mac.yaml'
+# fnconfig = 'config_CTF_Paris_mac.yaml'
+# fnconfig = 'config_CTF_Philly_mac.yaml'
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# load config details and file list to process
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+cfg = DCNN_CONFIG(verbose=True)
+cfg.load(fname=fnconfig)
+
+# get file list to process
+path_in = cfg.config['path']['data_meg']
+fnames = find_files(path_in, pattern='*-raw.fif')
+
+# ++++++++++++++++++++++++++++++++++++++++++++++
+#
+# loop across files
+#
+#++++++++++++++++++++++++++++++++++++++++++++++
+for fnraw in fnames:
+
+    # init DCNN obejct with config details
+    dcnn = DCNN(**cfg.config)
+    dcnn.verbose = True
+
+    # read raw data and apply noise reduction and downsample data
+    dcnn.meg.update(fname=fnraw)
+
+    # chop filtered data, apply ICA on chops, auto-label ICs and save results to disk
+    fgdcnn = dcnn.label_ica(save=True)
+
+print ('>>> finished ICA auto-labeling')

02_gDCNN_label_check.py

+"""
+gDCNN label_check
+
+    check ECG and EOG labeled ICA components and apply correction if applicable.
+
+"""
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# import libraries
+#
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+from dcnn_base import DCNN_CONFIG
+from dcnn_main import DCNN
+from dcnn_utils import find_files
+
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# select your config file
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# fnconfig = 'config_4D_CAU.yaml'
+# fnconfig = 'config_4D_INTEXT.yaml'
+
+# fnconfig = 'config_MEGIN.yaml'
+fnconfig = 'config_CTF_Paris.yaml'
+# fnconfig = 'config_CTF_Philly.yaml'
+
+#  --- MacOSX -----
+# fnconfig = 'config_4D_mac.yaml'
+# fnconfig = 'config_MEGIN_mac.yaml'
+# fnconfig = 'config_CTF_Paris_mac.yaml'
+
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# load config details and file list to process
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+cfg = DCNN_CONFIG(verbose=True)
+cfg.load(fname=fnconfig)
+
+# get file list to process
+path_in = cfg.config['path']['data_train']
+fnames = find_files(path_in, pattern='*.npz')
+
+# ++++++++++++++++++++++++++++++++++++++++++++++
+#
+# loop across files to process
+#
+#++++++++++++++++++++++++++++++++++++++++++++++
+for fname in fnames:
+
+    # init object with config details
+    dcnn = DCNN(cfg.config)
+    dcnn.load_gdcnn(fname)
+
+    # check IC labels (and apply corrections)
+    dcnn.check_labels(save=True)

dcnn_report.py

+#!/usr/bin/env python3
+# -+-coding: utf-8 -+-
+
+"""
+"""
+
+#--------------------------------------------
+# Authors: Frank Boers <f.boers@fz-juelich.de> 
+#
+#-------------------------------------------- 
+# Date: 16.12.19
+#-------------------------------------------- 
+# License: BSD (3-clause)
+#--------------------------------------------
+# Updates
+#--------------------------------------------
+
+import os
+from distutils.dir_util import mkpath
+
+#import pandas as pd
+#from PIL import Image
+
+import warnings
+import mne
+from mne.report import Report
+
+from jumeg.base.jumeg_base         import JUMEG_SLOTS
+from jumeg.base.jumeg_base         import jumeg_base as jb
+
+from jumeg.base.jumeg_base_config  import JuMEG_CONFIG as jCFG
+from jumeg.base                    import jumeg_logger
+
+
+__version__= "2020.05.05.001"
+
+logger = jumeg_logger.get_logger()
+
+class MNE_REPORT(JUMEG_SLOTS):
+    """
+    saving
+     noise_reducer fig
+     ica figs in HDF
+     epocher figs
+
+    and show in MNE Reports
+    """
+    __slots__ = { "postfix","html_extention","h5_extention","info_name","experiment","subject_id",
+                  "title","open_browser","write_hdf5","raw_psd","image_format","overwrite","verbose","debug",
+                  "_isOpen","_path","_fhdf","_MNE_REPORT"}
+    
+    def __init__(self,**kwargs):
+        super().__init__()
+        self._init()
+        
+        self._path        = "."
+       #---
+        self.image_format   = "png"
+        self.postfix        = "report"
+        self.html_extention = ".html"
+        self.h5_extention   = ".h5"
+        self.raw_psd        = False
+        self.title          = "JuMEG Preprocessing"
+        self.subject_id     = "0815"
+        
+        self._update_from_kwargs(**kwargs)
+        
+    @property
+    def path(self):
+        return self._path
+
+    @path.setter
+    def path(self,v):
+        if v:
+           self._path = jb.expandvars(v)
+           mkpath(self._path,mode=0o770)
+        else:
+            self._path=v
+    
+    @property
+    def fname(self): return self.experiment+"_"+self.subject_id+"-"+self.postfix
+    
+    @property
+    def fullname(self):
+        return os.path.join(self.path,self.fname)
+
+    @property
+    def html_name(self): return self.fullname + self.html_extention
+    @property
+    def hdf_name(self): return self.fullname + self.h5_extention
+    
+    @property
+    def image_extention(self): return "."+self.image_format
+    
+    @property
+    def isOpen(self): return self._isOpen
+
+    @property
+    def MNEreport(self): return self._MNE_REPORT
+    
+    @property
+    def fhdf(self): return self._fhdf
+    @fhdf.setter
+    def fhdf(self,v):
+        self._fhdf = jb.expandvars(v)
+        mkpath(os.path.dirname(self._fhdf),mode=0o770)
+    
+    def _update_from_kwargs(self,**kwargs):
+        super()._update_from_kwargs(**kwargs)
+        if "path" in kwargs:
+           self.path = kwargs.get("path",self._path)
+        
+    def open(self,**kwargs):
+        """
+
+        :param kwargs:
+        :return:
+        """
+        self._update_from_kwargs(**kwargs)
+        self._isOpen = False
+        # logger.info("report path (stage): {}".format(self.fullname))
+       
+        if self.overwrite:
+           try:
+              for fext in [ self.h5_extention,self.html_extention ]:
+                  fname = self.fullname+fext
+                  if os.path.isfile(fname):
+                     os.remove(self.fname)
+           except:
+              logger.exception("ERROR: can not overwrite report ile: {}".format(self.fullname))
+              return False
+       
+       # https://mne.tools/dev/auto_tutorials/misc/plot_report.html
+        try:
+            if os.path.isfile(self.hdf_name):
+               self._MNE_REPORT = mne.open_report(self.hdf_name)
+            else:
+               self._MNE_REPORT = mne.Report(info_fname=self.info_name,title=self.title,image_format=self.image_format,
+                                             raw_psd=self.raw_psd,verbose=self.verbose)
+               logger.info("open Report (h5): \n   -> {}".format(self._MNE_REPORT))
+               
+            self._isOpen = True
+        except:
+            logger.exception("ERROR: can not open or create MNE Report {}".format(self.hdf_name))
+        return self._isOpen
+
+    def save(self,overwrite=True):
+        if not self.isOpen:
+            logger.exception("ERROR in saving JuMEG MNE report: {}\n ---> Report not open\n".format(self.fullname))
+            return self.isOpen
+        
+        # mkpath( self.path,mode=0o770)
+      
+        if overwrite:
+           if os.path.isfile(self.html_name):
+              os.remove(self.html_name)
+           if os.path.isfile(self.hdf_name):
+              os.remove(self.hdf_name)
+       #--- html
+        self.MNEreport.save(self.html_name,overwrite=overwrite,open_browser=self.open_browser)
+        logger.info("DONE saving JuMEG MNE report [overwrite: {}] : HTML: {}\n".format(overwrite,self.html_name))
+       #--- h5
+        self.MNEreport.save( self.hdf_name, overwrite=overwrite,open_browser=False)
+        logger.info("DONE saving JuMEG MNE report [overwrite: {}] : HDF5: {}\n".format(overwrite,self.hdf_name))
+     
+        return self.isOpen
+    
+    def  update_report(self,path=None,data=None,section=None,prefix=None,replace=True):
+        """
+        load img from list
+        add to report
+        :param path   : report image path
+        :param data   : dict or list of pngs
+        :param section: section in report e.g.: Noise reducer, ICA
+        :param prefix : prefix for caption e.g.: ICA-0815_
+        :return:
+        """
+        if not data: return False
+ 
+        fimgs    = []
+        captions = []
+        section  = section
+    
+        if isinstance(data,(dict)):
+           for k in data.keys():
+               files    = []
+               fimgs    = []
+               captions = []
+               section  = k
+               
+               if data.get(k):
+                  files.extend( data.get(k) )
+               for f in files:
+                   fimgs.append( os.path.join(path,f.rstrip()) )
+                   captions.append( prefix+"-" + os.path.basename(f.rstrip().rsplit(self.image_extention,1)[0]) )
+               if self.debug:
+                  logger.debug("update MNE report: {}\n counts: {} ->\n {}".format(section,len(fimgs),fimgs) )
+               self.MNEreport.add_images_to_section(fimgs,captions=captions,section=section,replace=replace)
+                   
+           return True
+        
+        if isinstance(data,(list)):#
+           for f in data:
+               fimgs.append( os.path.join(path,f.rstrip()) )
+               captions.append( prefix+"-" + os.path.basename(f.rstrip().rsplit(self.image_extention,1)[0]) )
+        else:
+           fimgs.append(os.path.join(path,data.rstrip()))
+           captions.append(prefix + "-" + os.path.basename(data.rstrip().rsplit(self.image_extention,1)[0]))
+        if self.debug:
+           logger.debug("update MNE report: {}\n counts: {} ->\n {}".format(section,len(fimgs),fimgs))
+        self.MNEreport.add_images_to_section(fimgs,captions=captions,section=section,replace=replace)
+
+        if self.verbose:
+           logger.info("update MNE report: {}\n counts: {} ->\n {}".format(section,len(fimgs),fimgs) )
+
+        return True
+
+class JuMEG_REPORT(JUMEG_SLOTS):
+    """
+    saving
+     noise_reducer fig
+     ica figs in HDF
+     epocher figs
+    
+    and show in MNE Reports
+    """
+    __slots__= {"stage","_path","fname","report_cfg_extention","_REPORT_CFG","_REPORT","_jumeg_cfg"}
+    
+    def __init__(self,**kwargs):
+        super().__init__()
+        self.init()
+        self._jumeg_cfg = {"run":True,"save":True,"overwrite":False,"path":"report",
+                           "noise_reducer":{"run":True,"extention":"nr-raw.png"},
+                           "ica":{"run":True,"extention":"ar.png"}}
+        self._path = "."
+        self.report_cfg_extention = "-report.yaml"
+        self._REPORT = MNE_REPORT(**kwargs)
+        self._REPORT_CFG = jCFG()
+        
+    @property
+    def jumeg_cfg(self): return self._jumeg_cfg
+
+    @property
+    def report_cfg(self): return self._REPORT_CFG
+
+    @property
+    def Report(self):
+        return self._REPORT
+
+    @property
+    def verbose(self): return self._REPORT.verbose
+    @verbose.setter
+    def verbose(self,v):
+        self._REPORT.verbose=v
+    
+    @property
+    def path(self): return self._path
+    @path.setter
+    def path(self,v):
+        pext = self._jumeg_cfg.get("path","report")
+        if not v.endswith(pext):
+           self._path = jb.expandvars( os.path.join(v,pext) )
+        else:
+           self._path = jb.expandvars(v)
+      
+    def update_report_cfg(self):
+       #--- ['0815_TEST_20200412_1001_3', 'c']
+        f    = self.fname.split(",",1)[0].rsplit("_",1)[0] + self.report_cfg_extention
+        fcfg = os.path.join(self.path,f)
+        self._REPORT_CFG.update(config=fcfg)
+        
+    def run(self,**kwargs):
+        """
+        :param kwargs:
+        :return:
+        stage=jpl.stage,subject_id=subject_id,fname=raw_fname,config=jpl.config.get("report")
+        jReport.run(path=report_path,report_config_file=report_config_file,
+                    experiment="QUARTERS",subject_id=210857,config=config)
+            
+        open/read cfg /reprt/fname-report.yaml
+        MNEreport.open
+         update NR
+         update ICa
+         save HDf5
+         save htlm
+         
+        report config as dict
+        report_path=report_path,path=path,fname=raw_fname,subject_id=210857,config=config
+       
+        report:
+         run: True
+         save: True
+         overwrite: False
+         noise_reducer:
+          run: True
+         ica:
+          run: True
+        """
+        # logger.info(kwargs)
+        self._update_from_kwargs(**kwargs)
+        self.Report._update_from_kwargs(**kwargs)
+        
+       #--- try from jumeg config <report>  config=config.get("report")
+        if "config" in kwargs:
+           self._jumeg_cfg = kwargs.get("config")
+        
+        #logger.info("jumeg config report: {}".format( self.jumeg_cfg))
+        
+       #--- update from kwargs
+        self.stage = jb.expandvars(self.stage)
+        self.path  = kwargs.get("path",self.path) #--- report image path / image yaml file
+        
+        if self.stage.endswith("mne"):
+           rp = self.stage.replace("mne","report")
+        else:
+           rp = self.stage
+       #--- setup & open MNE report
+        if not self.Report.open(path=rp):  return False
+      
+       #--- report image config
+        self.update_report_cfg()
+        
+       #--- noise reducer
+        if self.jumeg_cfg.get("noise_reducer",False):
+           cfg = self.report_cfg.GetDataDict("noise_reducer")
+           if cfg:
+              self.Report.update_report(data=cfg.get("files"), path=self.path,section="Noise Reducer",prefix="NR")
+       #--- ica
+        if self.jumeg_cfg.get("ica",False):
+           cfg = self.report_cfg.GetDataDict("ica")
+           if cfg:
+              self.Report.update_report(data=cfg,path=self.path,section="ICA",prefix="ICA")
+       #--- save
+        if self.jumeg_cfg.get("save",False):
+           self.Report.save(overwrite=True)
+     
+if __name__ == "__main__":
+  # jumeg_logger.setup_script_logging(logger=logger)
+   pass
\ No newline at end of file

dcnn_run.py

+#!/usr/bin/env python
+from utils import find_files,get_args,dict2str
+
+#################   FROM FRANK ###################################
+import os,sys
+from dcnn_base   import DCNN_CONFIG
+from dcnn_main   import DCNN
+from dcnn_logger import setup_script_logging,init_logfile
+__version__= "2020-06-05-001"
+
+####################################################
+
+#==============================================
+# USER SETTINGS 4 Magic
+# ==============================================
+basedir  =  None
+data_meg = 'meg_rawdata/Philly'
+fnconfig = 'config_CTF_Philly.yaml'
+pattern  ='*-raw.fif'
+verbose  = True
+
+do_label_ica   = False
+do_label_check = True
+do_log2file    = False
+
+
+def run(fnconfig=None,basedir=None,data_meg=None,pattern='-raw.fif',
+        verbose=False,do_label_ica=False,do_label_check=False,do_log2file=False):
+    
+    # ToDo run for list of files or search under dir
+    
+    # ToDo check if this make sense?
+    #  if __DEFAULTS__: # if defaults define as global see on top of the script
+    #     cfg = DCNN_CONFIG(defaults=__DEFAULTS__)
+    
+    
+    cfg = DCNN_CONFIG(verbose=verbose)
+    cfg.load(fname=fnconfig)
+    
+    dcnn = DCNN(**cfg.config)  # init object with config details
+    dcnn.verbose = True
+    
+    if basedir: # FB test
+       cfg.config['path']['basedir'] = basedir
+    if data_meg:
+       cfg.config['path']['data_meg']= data_meg  # input directory
+    
+    # ==========================================================
+    # run ICA auto labelling
+    # ==========================================================
+    if do_label_ica:
+        # get file list to process
+        path_in = os.path.join( cfg.config['path']['basedir'],cfg.config['path']['data_meg'])
+        
+        fnames = find_files(path_in, pattern=pattern)
+        if not fnames:
+            logger.exception("ERROR No files found:\n  -> path: {}\n  -> pattern: {}".format(path_in,pattern))
+            import sys
+            sys.exit()
+    
+        fnraw = fnames[0]  # use first file for testing
+    
+       # -- init logfile
+        if do_log2file:
+           flog = os.path.splitext( fnraw )[0]
+           Hlog = init_logfile(logger=logger,fname=flog +".log",overwrite=True,name=flog)
+           msg= ["DCNN : {}".format(__version__),
+                 "  -> raw file: {}".format(fnraw),
+                 "  -> logfile : {}".format(Hlog.filename)]
+           logger.info("\n".join(msg))
+           if verbose:
+              cfg.info()
+           
+        # -- read raw data and apply noise reduction
+        dcnn.meg.update(fname=fnraw)
+        
+        # -- apply ICA on chops, label ICs save results to disk
+        # TODO: store chop-times, ECG,EOG in raw.annotations
+        #   - chop:  make use of annotations in get_chop_times_indices()
+        fgdcnn = dcnn.label_ica(save=True)
+        
+        if verbose:
+           dcnn.get_info()
+        
+    # ==========================================================
+    # check ICA auto labelling
+    # ==========================================================
+    if do_label_check:
+        path_in = cfg.config['path']['data_train']
+        fnames  = find_files(path_in, pattern= '*.npz')
+        fname   = fnames[0]  # use first file for testing
+    
+        dcnn.load_gdcnn(fname)
+        
+        if verbose:
+           dcnn.get_info()
+        
+           logger.info("dcnn ica chop  dump.\n{}\n{}\n\n".format( dcnn.ica.chop, dict2str(dcnn.ica.chop.dump())  ))
+           logger.info("dcnn ica n_chop.\n{}\n\n".format(dcnn.ica.chop.n_chop))
+
+
+
+
+if __name__ == "__main__":
+    # --  get parameter / flags from cmd line
+    argv = sys.argv
+    opt, parser = get_args(argv,version=__version__)
+    if len(argv) < 2:
+       parser.print_help()
+       sys.exit(-1)