Saccade

README.md

@@ -3,6 +3,9 @@
 
 A package to support neuroscientists in analyzing MEAs.
 
+## For Docker 1.17:
+`sudo docker run --network=eyecandy_default -v $(pwd):/data tbenst/glia analyze -e "http://eyecandy:3000" /data/R2_E1 convert`
+
 
 ## How to analyze data:
 1. Convert *.mcd files into *.voltages

glia/eyecandy.py

@@ -27,6 +27,7 @@
 Seconds = float
 ms = float
 UnitSpikeTrains = List[Dict[str,np.ndarray]]
+
 
 # TODO thread with async/await
 
@@ -96,11 +97,15 @@ def create_epl_gen_v2(program, epl, window_width, window_height, seed,
         raise
 
     stimuli = []
-    for json in response['stimulusList']:
-        value = json["value"]
-        value["stimulusIndex"] = json["stimulusIndex"]
+    for j in response['stimulusList']:
+        value = j["value"]
+        logger.debug(f"value in stimulusList: {str(j.keys())}")
+        value["stimulusIndex"] = j["stimulusIndex"]
         stimuli.append(value)
-    metadata = demjson.decode(response["metadata"])
+    logger.debug(f"metadata: {type(response['metadata'])}")
+    # logger.debug(f"stimuli: {stimuli}")
+    metadata = response["metadata"]
+    # metadata = demjson.decode(response["metadata"])
     assert type(metadata)==dict
     return (metadata,
         iter(sorted(stimuli, key=lambda x: x['stimulusIndex'])))
@@ -385,6 +390,7 @@ def get_stimulus_index_start_times(filtered,sampling_rate, stimulus_gen, percent
     next_stimulus = next(stimulus_gen)
     stimulus_list = []
     for i,state in np.ndenumerate(filtered):
+        # logger.debug(f"# of stimuli: {len(stimulus_list)}")
         i = i[0]
 
         #

glia/files.py

@@ -25,6 +25,18 @@ def match_filename(start,ext='txt'):
     else:
         raise(ValueError("Could not match file, try specifying full filename"))
 
+
+def find_notebook(directory):
+    notebooks = glob(os.path.join(directory, 'lab*.yml')) + \
+        glob(os.path.join(directory, 'lab*.yaml'))
+    if len(notebooks)==0:
+        raise ValueError("no lab notebooks (.yml) were found. Either add to directory," \
+            "or specify file path with -n.")
+    elif len(notebooks)>1:
+        logger.warning(f"""Found multiple possible lab notebooks.
+        Using {notebooks[0]}. If wrong, try manually specifying""")
+    return notebooks[0]
+
 def read_raw_voltage(raw_filename):
     """Read in a raw file exported from MCS datatool."""
     header, offset = get_header(raw_filename)

glia/machine_learning.py

@@ -1,13 +1,16 @@
 from collections import namedtuple
 from .pipeline import get_unit
 from .types import Unit
-from .functional import f_map, pmap, flatten, compose
+from .functional import f_map, pmap, flatten, compose, f_filter, group_by
 from .pipeline import group_dict_to_list
+from .eyecandy import checkerboard_contrast
 import numpy as np
 from scipy.sparse import csr_matrix
 from tqdm import tqdm
 from functools import partial
 import logging
+from sklearn import datasets, svm, metrics, neighbors
+import pandas as pd
 logger = logging.getLogger('glia')
 
 
@@ -173,18 +176,95 @@ def experiments_to_ndarrays(experiments, get_class=lambda x: x['metadata']['clas
     for i,array in enumerate(arrays):
         data[i] = array
 
-    # for idx in indices:
-        # data[idx] = 1
-    # for i,e in gen:
-    #     for unit_id, spikes in e['units'].items():
-    #         (row, column, unit_num) = key_map[unit_id]
-    #         for spike in spikes:
-    #             s = int(np.floor(spike*1000))
-    #             if s>1000:
-    #                 print('>1000',spikes, e)
-    #                 raise ValueError()
-    #             data[i,s,row,column,unit_num] = 1
-    #     classes[i] = get_class(e)
+    return (data, classes)
 
 
-    return (data, classes)
+def bin_100ms(data):
+    # turn the data in a (samples, feature) matrix from 100ms time bins:
+    (nconditions, nsizes, n_training, timesteps, n_x, n_y, n_units) = data.shape
+    new_steps = int(timesteps/100)
+    return np.sum(data.reshape(
+                        (nconditions, nsizes,
+                            n_training, new_steps, 100, n_x,n_y,n_units)),
+                    axis=4).reshape(
+                        (nconditions, nsizes, \
+                            n_training, new_steps*n_x*n_y*n_units))
+
+
+
+def bin_sum(data):
+    (nconditions, nsizes, n_training, timesteps, n_x, n_y, n_units) = data.shape
+    return np.sum(data,axis=3).reshape(
+        (nconditions, nsizes, n_training, n_x*n_y*n_units))
+
+letter_map = {'K': 4, 'C': 1, 'V': 9, 'N': 5, 'R': 7, 'H': 3, 'O': 6, 'Z': 10, 'D': 2, 'S': 8, 'BLANK': 0}
+letter_classes = list(map(lambda x: x[0],
+                   sorted(list(letter_map.items()),
+                          key=lambda x: x[1])))
+
+def classifier_helper(classifier, training, validation, classes=letter_classes):
+    training_data, training_target = training
+    validation_data, validation_target = validation
+
+    classifier.fit(training_data, training_target)
+    predicted = classifier.predict(validation_data)
+    expected = validation_target
+
+    report = metrics.classification_report(expected, predicted)
+    confusion = confusion_matrix(expected, predicted, classes)
+    return (report, confusion)
+
+
+def confusion_matrix(expected, predicted, classes=letter_classes):
+    m = metrics.confusion_matrix(expected, predicted)
+
+    return pd.DataFrame(data=m, index=classes,
+                              columns=classes)
+
+
+def px_to_logmar(px,px_per_deg=12.524):
+    minutes = px/px_per_deg*60
+    return np.log10(minutes)
+
+def get_stimulus_parameters(stimulus_list, stimulus_type, parameter):
+    f = compose(
+        f_filter(lambda x: x["stimulus"]['stimulusType']==stimulus_type),
+        partial(group_by,
+                key=lambda x: x["stimulus"][parameter])
+        )
+    parameters = sorted(list(f(stimulus_list).keys()))
+    logger.debug(f"Parameters: {parameters}")
+    assert len(parameters)>0
+    return parameters
+
+def get_image_parameters(stimulus_list):
+    f = compose(
+        f_filter(lambda x: x["stimulus"]['stimulusType']=='IMAGE'),
+        partial(group_by,
+                key=lambda x: x["stimulus"]["metadata"]["parameter"])
+        )
+    parameters = sorted(list(f(stimulus_list).keys()))
+    logger.debug(f"Parameters: {parameters}")
+    assert len(parameters)>0
+    return parameters
+
+
+def get_checkerboard_contrasts(stimulus_list):
+    f = compose(
+        f_filter(lambda x: x["stimulus"]['stimulusType']=='CHECKERBOARD'),
+        partial(group_by,
+                key=lambda x: checkerboard_contrast(x["stimulus"]))
+        )
+    contrasts = [float(x) for x in sorted(list(f(stimulus_list).keys()))]
+    assert len(contrasts)>0
+    return contrasts
+
+def svm_helper(training_data, training_target, validation_data, validation_target):
+    # Create a classifier: a support vector classifier
+    classifier = svm.SVC()
+    classifier.fit(training_data, training_target)
+
+    predicted = classifier.predict(validation_data)
+    expected = validation_target
+
+    return metrics.accuracy_score(expected, predicted)

glia/plot.py

@@ -18,6 +18,7 @@
 from .pipeline import get_unit
 from .functional import zip_dictionaries
 from .types import Unit
+from .machine_learning import tvt_by_percentage, svm_helper
 from glia.config import logger
 
 Seconds = float
@@ -492,3 +493,24 @@ def raster_group(fig, axis_gen, data):
     ax.set_ylabel("trials")
     ax.set_xlim((0,longest_group))
     ax.set_ylim((0,trial))
+
+
+def error_bars(data, target, ndraws=20):
+    n = data.shape[0]
+    accuracy = np.full((ndraws,), 0)
+    (ntrain, nvalid, _) = tvt_by_percentage(n,60,40,0)
+    indices = np.arange(n)
+    for i in range(ndraws):
+        np.random.shuffle(indices)
+        training_ind = indices[0:ntrain]
+        validation_ind = indices[ntrain:]
+
+        training_data = data[training_ind]
+        training_target = target[training_ind]
+        validation_data = data[validation_ind]
+        validation_target = target[validation_ind]
+
+        accuracy[i] = svm_helper(training_data, training_target,
+            validation_data, validation_target)
+    std = np.std(accuracy)
+    return (np.mean(accuracy),std,std)