From 3a55b7905eb8c44d9799e79ab8402bea488b4db4 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Thu, 26 Sep 2019 11:30:40 +1000
Subject: [PATCH 01/20] Initial framework for fnirs raw to optical density
 conversion

---
 mne/preprocessing/__init__.py                 |  1 +
 mne/preprocessing/optical_density.py          | 25 +++++++++++++++++++
 .../tests/test_optical_density.py             | 21 ++++++++++++++++
 3 files changed, 47 insertions(+)
 create mode 100644 mne/preprocessing/optical_density.py
 create mode 100644 mne/preprocessing/tests/test_optical_density.py

diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py
index 0c59396a91d..0672b1900c1 100644
--- a/mne/preprocessing/__init__.py
+++ b/mne/preprocessing/__init__.py
@@ -21,3 +21,4 @@
 from .stim import fix_stim_artifact
 from .maxwell import maxwell_filter
 from .xdawn import Xdawn
+from .optical_density import optical_density
diff --git a/mne/preprocessing/optical_density.py b/mne/preprocessing/optical_density.py
new file mode 100644
index 00000000000..ac5e6f36325
--- /dev/null
+++ b/mne/preprocessing/optical_density.py
@@ -0,0 +1,25 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+from ..io import BaseRaw
+from ..utils import _validate_type
+
+
+def optical_density(raw, verbose=None):
+    r"""Convert NIRS raw data to optical density.
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        The raw data.
+    %(verbose)s
+
+    Returns
+    -------
+    raw : instance of Raw
+        The modified raw instance. Operates in place.
+
+    """
+    _validate_type(raw, BaseRaw, 'raw')
+    return raw
diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
new file mode 100644
index 00000000000..7da129d257c
--- /dev/null
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -0,0 +1,21 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import os.path as op
+
+from mne.datasets.testing import data_path
+from mne.io import read_raw_nirx, BaseRaw
+from mne.preprocessing import optical_density
+from mne.utils import _validate_type
+
+fname_nirx = op.join(data_path(download=False),
+                     'NIRx', 'nirx_15_2_recording_w_short')
+
+
+def test_optical_density():
+    """Test fix stim artifact."""
+    raw = read_raw_nirx(fname_nirx, preload=True)
+    od = optical_density(raw, verbose=None)
+    _validate_type(od, BaseRaw, 'raw')
+    print(od)

From e30713987ef55d62ccc47486d08e810f46c7650e Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Thu, 26 Sep 2019 13:22:19 +1000
Subject: [PATCH 02/20] Add optical density to docs

---
 doc/python_reference.rst             | 1 +
 mne/preprocessing/optical_density.py | 3 +--
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/python_reference.rst b/doc/python_reference.rst
index e0a6566365e..73a731861c3 100644
--- a/doc/python_reference.rst
+++ b/doc/python_reference.rst
@@ -292,6 +292,7 @@ Projections:
    compute_proj_raw
    read_proj
    write_proj
+   optical_density
 
 :py:mod:`mne.channels`:
 
diff --git a/mne/preprocessing/optical_density.py b/mne/preprocessing/optical_density.py
index ac5e6f36325..83695372ed9 100644
--- a/mne/preprocessing/optical_density.py
+++ b/mne/preprocessing/optical_density.py
@@ -6,14 +6,13 @@
 from ..utils import _validate_type
 
 
-def optical_density(raw, verbose=None):
+def optical_density(raw):
     r"""Convert NIRS raw data to optical density.
 
     Parameters
     ----------
     raw : instance of Raw
         The raw data.
-    %(verbose)s
 
     Returns
     -------

From bd6479165994a56f4395eb2fa1e22b722f42adff Mon Sep 17 00:00:00 2001
From: Robert Luke <748691+rob-luke@users.noreply.github.com>
Date: Fri, 27 Sep 2019 12:21:16 +1000
Subject: [PATCH 03/20] Update mne/preprocessing/tests/test_optical_density.py

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>
---
 mne/preprocessing/tests/test_optical_density.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 7da129d257c..b32bb1347a1 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -13,6 +13,7 @@
                      'NIRx', 'nirx_15_2_recording_w_short')
 
 
+@testing.requires_testing_data
 def test_optical_density():
     """Test fix stim artifact."""
     raw = read_raw_nirx(fname_nirx, preload=True)

From 7337b5e15d30381bba41de821f306d7166b42bd0 Mon Sep 17 00:00:00 2001
From: Robert Luke <748691+rob-luke@users.noreply.github.com>
Date: Fri, 27 Sep 2019 12:21:37 +1000
Subject: [PATCH 04/20] Update mne/preprocessing/tests/test_optical_density.py

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>
---
 mne/preprocessing/tests/test_optical_density.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index b32bb1347a1..47804311ea8 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -17,6 +17,6 @@
 def test_optical_density():
     """Test fix stim artifact."""
     raw = read_raw_nirx(fname_nirx, preload=True)
-    od = optical_density(raw, verbose=None)
+    od = optical_density(raw)
     _validate_type(od, BaseRaw, 'raw')
     print(od)

From e6a5a1bc487d5af9ad0b90a947dfbfe0687bfe74 Mon Sep 17 00:00:00 2001
From: Robert Luke <748691+rob-luke@users.noreply.github.com>
Date: Fri, 27 Sep 2019 12:36:00 +1000
Subject: [PATCH 05/20] Import testing

---
 mne/preprocessing/tests/test_optical_density.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 47804311ea8..3ad58c41cb3 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -8,6 +8,7 @@
 from mne.io import read_raw_nirx, BaseRaw
 from mne.preprocessing import optical_density
 from mne.utils import _validate_type
+from mne.datasets import testing
 
 fname_nirx = op.join(data_path(download=False),
                      'NIRx', 'nirx_15_2_recording_w_short')

From 59a96c830429f443784c1db68bf7643f49594121 Mon Sep 17 00:00:00 2001
From: Robert Luke <748691+rob-luke@users.noreply.github.com>
Date: Sat, 28 Sep 2019 07:38:05 +1000
Subject: [PATCH 06/20] Put optical density docs in preprocessing section

---
 doc/python_reference.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/python_reference.rst b/doc/python_reference.rst
index 73a731861c3..eceb55a148f 100644
--- a/doc/python_reference.rst
+++ b/doc/python_reference.rst
@@ -292,7 +292,6 @@ Projections:
    compute_proj_raw
    read_proj
    write_proj
-   optical_density
 
 :py:mod:`mne.channels`:
 
@@ -360,6 +359,7 @@ Projections:
    read_ica
    run_ica
    corrmap
+   optical_density
 
 EEG referencing:
 

From 97e0e1d84921488ae4cc51e6491064afd60af309 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Mon, 30 Sep 2019 15:25:44 +1000
Subject: [PATCH 07/20] Implement optical density conversion

---
 mne/preprocessing/optical_density.py            | 12 ++++++++++++
 mne/preprocessing/tests/test_optical_density.py |  1 -
 2 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/mne/preprocessing/optical_density.py b/mne/preprocessing/optical_density.py
index 83695372ed9..77e496bcabb 100644
--- a/mne/preprocessing/optical_density.py
+++ b/mne/preprocessing/optical_density.py
@@ -3,7 +3,11 @@
 # License: BSD (3-clause)
 
 from ..io import BaseRaw
+from ..io.constants import FIFF
 from ..utils import _validate_type
+from ..io.pick import _picks_to_idx
+
+from numpy import mean, divide, log, ones, concatenate
 
 
 def optical_density(raw):
@@ -21,4 +25,12 @@ def optical_density(raw):
 
     """
     _validate_type(raw, BaseRaw, 'raw')
+    picks = _picks_to_idx(raw.info, 'fnirs_raw')
+    data_means = mean(raw.get_data(), 1)
+    for ii in picks:
+        data = -1.0 * log(divide(raw.get_data(ii).T,
+                           ones((len(raw), 1)) * data_means[ii]))
+        raw._data[ii, :] = data.T
+        raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
+
     return raw
diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 3ad58c41cb3..cc96af8d83e 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -20,4 +20,3 @@ def test_optical_density():
     raw = read_raw_nirx(fname_nirx, preload=True)
     od = optical_density(raw)
     _validate_type(od, BaseRaw, 'raw')
-    print(od)

From acad961e40f9b2406cbb0427cc7c70a28f886680 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Mon, 30 Sep 2019 15:28:54 +1000
Subject: [PATCH 08/20] Fix optical density code style issues

---
 mne/preprocessing/optical_density.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mne/preprocessing/optical_density.py b/mne/preprocessing/optical_density.py
index 77e496bcabb..2733437f172 100644
--- a/mne/preprocessing/optical_density.py
+++ b/mne/preprocessing/optical_density.py
@@ -7,7 +7,7 @@
 from ..utils import _validate_type
 from ..io.pick import _picks_to_idx
 
-from numpy import mean, divide, log, ones, concatenate
+from numpy import mean, divide, log, ones
 
 
 def optical_density(raw):
@@ -29,7 +29,7 @@ def optical_density(raw):
     data_means = mean(raw.get_data(), 1)
     for ii in picks:
         data = -1.0 * log(divide(raw.get_data(ii).T,
-                           ones((len(raw), 1)) * data_means[ii]))
+                          ones((len(raw), 1)) * data_means[ii]))
         raw._data[ii, :] = data.T
         raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
 

From d0ed6f45e1de522bba6debdeb4890fa6ddbabaed Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Tue, 1 Oct 2019 08:58:45 +1000
Subject: [PATCH 09/20] Add plotting of optical density data

---
 mne/viz/raw.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/viz/raw.py b/mne/viz/raw.py
index c8707189422..7c484179f52 100644
--- a/mne/viz/raw.py
+++ b/mne/viz/raw.py
@@ -332,7 +332,7 @@ def plot_raw(raw, events=None, duration=10.0, start=0.0, n_channels=20,
     for t in ['grad', 'mag']:
         inds += [pick_types(info, meg=t, ref_meg=False, exclude=[])]
         types += [t] * len(inds[-1])
-    for t in ['hbo', 'hbr', 'fnirs_raw']:
+    for t in ['hbo', 'hbr', 'fnirs_raw', 'fnirs_od']:
         inds += [pick_types(info, meg=False, ref_meg=False, fnirs=t,
                             exclude=[])]
         types += [t] * len(inds[-1])

From 9efc70098148bc52599613de8df08bd73647eb47 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Tue, 1 Oct 2019 09:08:50 +1000
Subject: [PATCH 10/20] Add defaults for optical density plotting

---
 mne/defaults.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/mne/defaults.py b/mne/defaults.py
index 9d362a988ac..8f3f265e480 100644
--- a/mne/defaults.py
+++ b/mne/defaults.py
@@ -11,20 +11,22 @@
                ref_meg='steelblue', misc='k', stim='k', resp='k', chpi='k',
                exci='k', ias='k', syst='k', seeg='saddlebrown', dipole='k',
                gof='k', bio='k', ecog='k', hbo='darkblue', hbr='b',
-               fnirs_raw='k'),
+               fnirs_raw='k', fnirs_od='k'),
     units=dict(mag='fT', grad='fT/cm', eeg='uV', eog='uV', ecg='uV', emg='uV',
                misc='AU', seeg='mV', dipole='nAm', gof='GOF', bio='uV',
-               ecog='uV', hbo='uM', hbr='uM', ref_meg='fT', fnirs_raw='V'),
+               ecog='uV', hbo='uM', hbr='uM', ref_meg='fT', fnirs_raw='V',
+               fnirs_od='V'),
     # scalings for the units
     scalings=dict(mag=1e15, grad=1e13, eeg=1e6, eog=1e6, emg=1e6, ecg=1e6,
                   misc=1.0, seeg=1e3, dipole=1e9, gof=1.0, bio=1e6, ecog=1e6,
-                  hbo=1e6, hbr=1e6, ref_meg=1e15, fnirs_raw=1),
+                  hbo=1e6, hbr=1e6, ref_meg=1e15, fnirs_raw=1.0, fnirs_od=1.0),
     # rough guess for a good plot
     scalings_plot_raw=dict(mag=1e-12, grad=4e-11, eeg=20e-6, eog=150e-6,
                            ecg=5e-4, emg=1e-3, ref_meg=1e-12, misc='auto',
                            stim=1, resp=1, chpi=1e-4, exci=1, ias=1, syst=1,
                            seeg=1e-4, bio=1e-6, ecog=1e-4, hbo=10e-6,
-                           hbr=10e-6, whitened=10., fnirs_raw=2e-2),
+                           hbr=10e-6, whitened=10., fnirs_raw=2e-2,
+                           fnirs_od=2e-2),
     scalings_cov_rank=dict(mag=1e12, grad=1e11, eeg=1e5,  # ~100x scalings
                            seeg=1e1, ecog=1e4, hbo=1e4, hbr=1e4),
     ylim=dict(mag=(-600., 600.), grad=(-200., 200.), eeg=(-200., 200.),
@@ -35,7 +37,8 @@
                 ecg='ECG', emg='EMG', misc='misc', seeg='sEEG', bio='BIO',
                 dipole='Dipole', ecog='ECoG', hbo='Oxyhemoglobin',
                 ref_meg='Reference Magnetometers', fnirs_raw='fNIRS',
-                hbr='Deoxyhemoglobin', gof='Goodness of fit'),
+                fnirs_od='fNIRS', hbr='Deoxyhemoglobin',
+                gof='Goodness of fit'),
     mask_params=dict(marker='o',
                      markerfacecolor='w',
                      markeredgecolor='k',

From f42cd3ea1862042d89674cbd46d538f48993cecf Mon Sep 17 00:00:00 2001
From: Robert Luke <748691+rob-luke@users.noreply.github.com>
Date: Thu, 3 Oct 2019 06:46:06 +1000
Subject: [PATCH 11/20] Update mne/preprocessing/optical_density.py

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>
---
 mne/preprocessing/optical_density.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/mne/preprocessing/optical_density.py b/mne/preprocessing/optical_density.py
index 2733437f172..cb1f81894f4 100644
--- a/mne/preprocessing/optical_density.py
+++ b/mne/preprocessing/optical_density.py
@@ -30,7 +30,8 @@ def optical_density(raw):
     for ii in picks:
         data = -1.0 * log(divide(raw.get_data(ii).T,
                           ones((len(raw), 1)) * data_means[ii]))
-        raw._data[ii, :] = data.T
+        raw._data[ii] /= means[ii]
+        np.log(raw._data[ii], out=raw._data[ii])
         raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
 
     return raw

From 687b63f680808bdf83b2959887263e6634929ee1 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Thu, 3 Oct 2019 08:07:10 +1000
Subject: [PATCH 12/20] Handle negative intensities in optical density
 conversion

---
 mne/preprocessing/__init__.py                  |  2 +-
 ...{optical_density.py => _optical_density.py} | 18 +++++++++++++-----
 .../tests/test_optical_density.py              | 13 ++++++++++++-
 3 files changed, 26 insertions(+), 7 deletions(-)
 rename mne/preprocessing/{optical_density.py => _optical_density.py} (51%)

diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py
index 0672b1900c1..583202476b4 100644
--- a/mne/preprocessing/__init__.py
+++ b/mne/preprocessing/__init__.py
@@ -21,4 +21,4 @@
 from .stim import fix_stim_artifact
 from .maxwell import maxwell_filter
 from .xdawn import Xdawn
-from .optical_density import optical_density
+from ._optical_density import optical_density
diff --git a/mne/preprocessing/optical_density.py b/mne/preprocessing/_optical_density.py
similarity index 51%
rename from mne/preprocessing/optical_density.py
rename to mne/preprocessing/_optical_density.py
index cb1f81894f4..e401d14648e 100644
--- a/mne/preprocessing/optical_density.py
+++ b/mne/preprocessing/_optical_density.py
@@ -7,7 +7,7 @@
 from ..utils import _validate_type
 from ..io.pick import _picks_to_idx
 
-from numpy import mean, divide, log, ones
+from numpy import mean, abs, log, where
 
 
 def optical_density(raw):
@@ -27,11 +27,19 @@ def optical_density(raw):
     _validate_type(raw, BaseRaw, 'raw')
     picks = _picks_to_idx(raw.info, 'fnirs_raw')
     data_means = mean(raw.get_data(), 1)
+
+    # The devices measure light intensity. Negative light intensities should
+    # not occur. If they do it is likely due to hardware or movement issues.
+    # Set all negative values to abs(x), this also has the benefit of ensuring
+    # that the means are all greater than zero for the division below.
+    if len(where(raw._data[picks] <= 0)[0]) > 0:
+        # TODO How to throw warning without exiting tests
+        print("    Negative intensities encountered. Setting to abs(x)")
+        raw._data[picks] = abs(raw._data[picks])
+
     for ii in picks:
-        data = -1.0 * log(divide(raw.get_data(ii).T,
-                          ones((len(raw), 1)) * data_means[ii]))
-        raw._data[ii] /= means[ii]
-        np.log(raw._data[ii], out=raw._data[ii])
+        raw._data[ii] /= data_means[ii]
+        log(raw._data[ii], out=raw._data[ii])
         raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
 
     return raw
diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index cc96af8d83e..440443644f8 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -10,13 +10,24 @@
 from mne.utils import _validate_type
 from mne.datasets import testing
 
+from numpy import mean
+
 fname_nirx = op.join(data_path(download=False),
                      'NIRx', 'nirx_15_2_recording_w_short')
 
 
 @testing.requires_testing_data
 def test_optical_density():
-    """Test fix stim artifact."""
+    """Test that the optical density conversion returns BaseRaw."""
+    raw = read_raw_nirx(fname_nirx, preload=True)
+    od = optical_density(raw)
+    _validate_type(od, BaseRaw, 'raw')
+
+
+@testing.requires_testing_data
+def test_optical_density_zeromean():
+    """Test that optical density can process zero mean data."""
     raw = read_raw_nirx(fname_nirx, preload=True)
+    raw._data[4] -= mean(raw._data[4])
     od = optical_density(raw)
     _validate_type(od, BaseRaw, 'raw')

From f80c0de19abc3158844eeec502010e3713b1eb71 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Thu, 3 Oct 2019 15:30:37 +1000
Subject: [PATCH 13/20] Add in -1* that was removed in previous commit

---
 mne/preprocessing/_optical_density.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mne/preprocessing/_optical_density.py b/mne/preprocessing/_optical_density.py
index e401d14648e..13e516304c9 100644
--- a/mne/preprocessing/_optical_density.py
+++ b/mne/preprocessing/_optical_density.py
@@ -40,6 +40,7 @@ def optical_density(raw):
     for ii in picks:
         raw._data[ii] /= data_means[ii]
         log(raw._data[ii], out=raw._data[ii])
+        raw._data[ii] *= -1
         raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
 
     return raw

From 0971ef72dd4fa1dc0f71116c0b2da30d4d21d5e0 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Fri, 4 Oct 2019 09:37:28 +1000
Subject: [PATCH 14/20] Add tests and fix imports for optical density

---
 mne/preprocessing/_optical_density.py           | 16 ++++++++--------
 mne/preprocessing/tests/test_optical_density.py | 15 ++++++++++-----
 2 files changed, 18 insertions(+), 13 deletions(-)

diff --git a/mne/preprocessing/_optical_density.py b/mne/preprocessing/_optical_density.py
index 13e516304c9..5d302f904a6 100644
--- a/mne/preprocessing/_optical_density.py
+++ b/mne/preprocessing/_optical_density.py
@@ -1,13 +1,14 @@
 # Authors: Robert Luke <mail@robertluke.net>
+#          Eric Larson <larson.eric.d@gmail.com>
 #
 # License: BSD (3-clause)
 
 from ..io import BaseRaw
 from ..io.constants import FIFF
-from ..utils import _validate_type
+from ..utils import _validate_type, warn
 from ..io.pick import _picks_to_idx
 
-from numpy import mean, abs, log, where
+import numpy as np
 
 
 def optical_density(raw):
@@ -26,20 +27,19 @@ def optical_density(raw):
     """
     _validate_type(raw, BaseRaw, 'raw')
     picks = _picks_to_idx(raw.info, 'fnirs_raw')
-    data_means = mean(raw.get_data(), 1)
+    data_means = np.mean(raw.get_data(), 1)
 
     # The devices measure light intensity. Negative light intensities should
     # not occur. If they do it is likely due to hardware or movement issues.
     # Set all negative values to abs(x), this also has the benefit of ensuring
     # that the means are all greater than zero for the division below.
-    if len(where(raw._data[picks] <= 0)[0]) > 0:
-        # TODO How to throw warning without exiting tests
-        print("    Negative intensities encountered. Setting to abs(x)")
-        raw._data[picks] = abs(raw._data[picks])
+    if len(np.where(raw._data[picks] <= 0)[0]) > 0:
+        warn("Negative intensities encountered. Setting to abs(x)")
+        raw._data[picks] = np.abs(raw._data[picks])
 
     for ii in picks:
         raw._data[ii] /= data_means[ii]
-        log(raw._data[ii], out=raw._data[ii])
+        np.log(raw._data[ii], out=raw._data[ii])
         raw._data[ii] *= -1
         raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
 
diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 440443644f8..02c029e375d 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -10,7 +10,8 @@
 from mne.utils import _validate_type
 from mne.datasets import testing
 
-from numpy import mean
+import numpy as np
+import pytest as pytest
 
 fname_nirx = op.join(data_path(download=False),
                      'NIRx', 'nirx_15_2_recording_w_short')
@@ -18,16 +19,20 @@
 
 @testing.requires_testing_data
 def test_optical_density():
-    """Test that the optical density conversion returns BaseRaw."""
+    """Test return type for optical density."""
     raw = read_raw_nirx(fname_nirx, preload=True)
+    assert 'fnirs_raw' in raw
+    assert 'fnirs_od' not in raw
     od = optical_density(raw)
     _validate_type(od, BaseRaw, 'raw')
+    assert 'fnirs_raw' not in raw
+    assert 'fnirs_od' in raw
 
 
 @testing.requires_testing_data
 def test_optical_density_zeromean():
     """Test that optical density can process zero mean data."""
     raw = read_raw_nirx(fname_nirx, preload=True)
-    raw._data[4] -= mean(raw._data[4])
-    od = optical_density(raw)
-    _validate_type(od, BaseRaw, 'raw')
+    raw._data[4] -= np.mean(raw._data[4])
+    with pytest.warns(RuntimeWarning, match='Negative'):
+        od = optical_density(raw)

From 91ece6475c6ede682bb08b08a553b9edcbe55818 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Fri, 4 Oct 2019 11:30:52 +1000
Subject: [PATCH 15/20] Add tests for optical density conversion

---
 .../tests/test_optical_density.py             | 26 ++++++++++++++++---
 1 file changed, 23 insertions(+), 3 deletions(-)

diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 02c029e375d..3594850e68f 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -1,4 +1,5 @@
 # Authors: Robert Luke <mail@robertluke.net>
+#          Eric Larson <larson.eric.d@gmail.com>
 #
 # License: BSD (3-clause)
 
@@ -23,8 +24,8 @@ def test_optical_density():
     raw = read_raw_nirx(fname_nirx, preload=True)
     assert 'fnirs_raw' in raw
     assert 'fnirs_od' not in raw
-    od = optical_density(raw)
-    _validate_type(od, BaseRaw, 'raw')
+    raw = optical_density(raw)
+    _validate_type(raw, BaseRaw, 'raw')
     assert 'fnirs_raw' not in raw
     assert 'fnirs_od' in raw
 
@@ -35,4 +36,23 @@ def test_optical_density_zeromean():
     raw = read_raw_nirx(fname_nirx, preload=True)
     raw._data[4] -= np.mean(raw._data[4])
     with pytest.warns(RuntimeWarning, match='Negative'):
-        od = optical_density(raw)
+        raw = optical_density(raw)
+    assert 'fnirs_od' in raw
+
+
+@testing.requires_testing_data
+def test_optical_density_manual():
+    """Test optical density on known values"""
+    test_tol = 0.01
+    raw = read_raw_nirx(fname_nirx, preload=True)
+    # log(1) = 0
+    raw._data[4] = np.ones((145))
+    # log(0.5)/-1 = 0.69
+    # log(1.5)/-1 = -0.40
+    test_data = np.tile([0.5, 1.5], 73)[:145]
+    raw._data[5] = test_data
+
+    od = optical_density(raw)
+    assert np.sum(od.get_data()[4] - np.zeros(145)) == 0
+    assert np.abs(od.get_data()[5][0] - 0.69) < test_tol
+    assert np.abs(od.get_data()[5][1] - -0.40) < test_tol

From 9575d4708f4cf0f6925100d15de74a5e7da7f0f3 Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Fri, 4 Oct 2019 13:55:10 +1000
Subject: [PATCH 16/20] Use assert_allclose for testing optical density

---
 mne/preprocessing/tests/test_optical_density.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 3594850e68f..782a714f5a1 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -11,8 +11,10 @@
 from mne.utils import _validate_type
 from mne.datasets import testing
 
-import numpy as np
 import pytest as pytest
+import numpy as np
+from numpy.testing import assert_allclose
+
 
 fname_nirx = op.join(data_path(download=False),
                      'NIRx', 'nirx_15_2_recording_w_short')
@@ -53,6 +55,5 @@ def test_optical_density_manual():
     raw._data[5] = test_data
 
     od = optical_density(raw)
-    assert np.sum(od.get_data()[4] - np.zeros(145)) == 0
-    assert np.abs(od.get_data()[5][0] - 0.69) < test_tol
-    assert np.abs(od.get_data()[5][1] - -0.40) < test_tol
+    assert_allclose(od.get_data([4]), 0.)
+    assert_allclose(od.get_data([5])[0, :2], [0.69, -0.4], atol=test_tol)

From f3fe457fb1dc4f26c8d1c3a86c218d313c44892c Mon Sep 17 00:00:00 2001
From: Eric Larson <larson.eric.d@gmail.com>
Date: Fri, 4 Oct 2019 10:00:11 -0400
Subject: [PATCH 17/20] Update test_optical_density.py

---
 mne/preprocessing/tests/test_optical_density.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 782a714f5a1..7ba8671dd47 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -44,7 +44,7 @@ def test_optical_density_zeromean():
 
 @testing.requires_testing_data
 def test_optical_density_manual():
-    """Test optical density on known values"""
+    """Test optical density on known values."""
     test_tol = 0.01
     raw = read_raw_nirx(fname_nirx, preload=True)
     # log(1) = 0

From 2509f1edfb4c573a8c6a9e96cec57718153e38dd Mon Sep 17 00:00:00 2001
From: Robert Luke <748691+rob-luke@users.noreply.github.com>
Date: Sat, 5 Oct 2019 07:47:25 +1000
Subject: [PATCH 18/20] Update mne/preprocessing/_optical_density.py

Co-Authored-By: Alexandre Gramfort <alexandre.gramfort@m4x.org>
---
 mne/preprocessing/_optical_density.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/preprocessing/_optical_density.py b/mne/preprocessing/_optical_density.py
index 5d302f904a6..e319c7bf580 100644
--- a/mne/preprocessing/_optical_density.py
+++ b/mne/preprocessing/_optical_density.py
@@ -27,7 +27,7 @@ def optical_density(raw):
     """
     _validate_type(raw, BaseRaw, 'raw')
     picks = _picks_to_idx(raw.info, 'fnirs_raw')
-    data_means = np.mean(raw.get_data(), 1)
+    data_means = np.mean(raw.get_data(), axis=1)
 
     # The devices measure light intensity. Negative light intensities should
     # not occur. If they do it is likely due to hardware or movement issues.

From 79a9e766f934a15e778b063c2265c4a5f357857c Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Sat, 5 Oct 2019 11:33:42 +1000
Subject: [PATCH 19/20] Do not operate in place for optical density

---
 mne/preprocessing/_optical_density.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/mne/preprocessing/_optical_density.py b/mne/preprocessing/_optical_density.py
index e319c7bf580..3ba8005df4b 100644
--- a/mne/preprocessing/_optical_density.py
+++ b/mne/preprocessing/_optical_density.py
@@ -22,9 +22,10 @@ def optical_density(raw):
     Returns
     -------
     raw : instance of Raw
-        The modified raw instance. Operates in place.
+        The modified raw instance.
 
     """
+    raw = raw.copy()
     _validate_type(raw, BaseRaw, 'raw')
     picks = _picks_to_idx(raw.info, 'fnirs_raw')
     data_means = np.mean(raw.get_data(), axis=1)

From b2940ff1817d0eb98c97d7c9f462f76317111c2e Mon Sep 17 00:00:00 2001
From: Robert Luke <hello@robertluke.net>
Date: Sat, 5 Oct 2019 12:01:46 +1000
Subject: [PATCH 20/20] Force data to load and modify test to include
 preload=False condition

---
 mne/preprocessing/_optical_density.py           | 2 +-
 mne/preprocessing/tests/test_optical_density.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/mne/preprocessing/_optical_density.py b/mne/preprocessing/_optical_density.py
index 3ba8005df4b..ea96835ed6b 100644
--- a/mne/preprocessing/_optical_density.py
+++ b/mne/preprocessing/_optical_density.py
@@ -25,7 +25,7 @@ def optical_density(raw):
         The modified raw instance.
 
     """
-    raw = raw.copy()
+    raw = raw.copy().load_data()
     _validate_type(raw, BaseRaw, 'raw')
     picks = _picks_to_idx(raw.info, 'fnirs_raw')
     data_means = np.mean(raw.get_data(), axis=1)
diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py
index 7ba8671dd47..1883039e1b1 100644
--- a/mne/preprocessing/tests/test_optical_density.py
+++ b/mne/preprocessing/tests/test_optical_density.py
@@ -23,7 +23,7 @@
 @testing.requires_testing_data
 def test_optical_density():
     """Test return type for optical density."""
-    raw = read_raw_nirx(fname_nirx, preload=True)
+    raw = read_raw_nirx(fname_nirx, preload=False)
     assert 'fnirs_raw' in raw
     assert 'fnirs_od' not in raw
     raw = optical_density(raw)