diff --git a/.flake8 b/.flake8 new file mode 100644 index 000000000..a40a644a0 --- /dev/null +++ b/.flake8 @@ -0,0 +1,2 @@ +[flake8] +extend-ignore = E203 diff --git a/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_neuropixels_project_cloud_api.py b/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_neuropixels_project_cloud_api.py index 11a7398cd..96b588fa4 100644 --- a/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_neuropixels_project_cloud_api.py +++ b/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_neuropixels_project_cloud_api.py @@ -6,27 +6,28 @@ from allensdk.brain_observatory.behavior.behavior_session import ( BehaviorSession, ) +from allensdk.brain_observatory.ecephys._probe import ProbeWithLFPMeta from allensdk.brain_observatory.ecephys.behavior_ecephys_session import ( BehaviorEcephysSession, ) -from allensdk.brain_observatory.ecephys._probe import ProbeWithLFPMeta from allensdk.core.dataframe_utils import ( enforce_df_int_typing, - return_one_dataframe_row_only + return_one_dataframe_row_only, ) INTEGER_COLUMNS = [ - "prior_exposures_to_image_set", "ecephys_session_id", "unit_count", - "probe_count", "channel_count" + "prior_exposures_to_image_set", + "ecephys_session_id", + "unit_count", + "probe_count", + "channel_count", ] class VisualBehaviorNeuropixelsProjectCloudApi(ProjectCloudApiBase): - MANIFEST_COMPATIBILITY = ["0.1.0", "10.0.0"] def _load_manifest_tables(self): - self._get_ecephys_session_table() self._get_behavior_session_table() self._get_unit_table() @@ -56,7 +57,7 @@ def get_behavior_session( row = return_one_dataframe_row_only( input_table=self._behavior_session_table, index_value=behavior_session_id, - table_name="behavior_session_table" + table_name="behavior_session_table", ) row = row.squeeze() ecephys_session_id = row.ecephys_session_id @@ -68,7 +69,7 @@ def get_behavior_session( row = return_one_dataframe_row_only( input_table=self._ecephys_session_table, index_value=ecephys_session_id, - table_name="ecephys_session_table" + table_name="ecephys_session_table", ) file_id = str(int(row[self.cache.file_id_column])) @@ -79,7 +80,6 @@ def get_behavior_session( def get_ecephys_session( self, ecephys_session_id: int ) -> BehaviorEcephysSession: - """get a BehaviorEcephysSession by specifying ecephys_session_id Parameters @@ -95,7 +95,7 @@ def get_ecephys_session( session_meta = return_one_dataframe_row_only( input_table=self._ecephys_session_table, index_value=ecephys_session_id, - table_name="ecephys_session_table" + table_name="ecephys_session_table", ) probes_meta = self._probe_table[ (self._probe_table["ecephys_session_id"] == ecephys_session_id) @@ -122,10 +122,9 @@ def f(): probe_meta = { p.name: ProbeWithLFPMeta( lfp_csd_filepath=make_lazy_load_filepath_function( - file_id=str(int(getattr( - p, self.cache.file_id_column))) - ), - lfp_sampling_rate=p.lfp_sampling_rate + file_id=str(int(getattr(p, self.cache.file_id_column))) + ), + lfp_sampling_rate=p.lfp_sampling_rate, ) for p in probes_meta.itertuples(index=False) } diff --git a/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_project_cloud_api.py b/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_project_cloud_api.py index 1b566e225..5ace2a643 100644 --- a/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_project_cloud_api.py +++ b/allensdk/brain_observatory/behavior/behavior_project_cache/project_apis/data_io/behavior_project_cloud_api.py @@ -13,17 +13,23 @@ from allensdk.brain_observatory.behavior.behavior_session import ( BehaviorSession, ) -from allensdk.core.utilities import literal_col_eval from allensdk.core.dataframe_utils import ( enforce_df_int_typing, - return_one_dataframe_row_only + return_one_dataframe_row_only, ) +from allensdk.core.utilities import literal_col_eval COL_EVAL_LIST = ["ophys_experiment_id", "ophys_container_id", "driver_line"] -INTEGER_COLUMNS = ["session_number", "prior_exposures_to_image_set", - "ophys_session_id", "imaging_plane_group_count", - "imaging_plane_group", "targeted_areas", - "num_depths_per_area", "num_targeted_structures"] +INTEGER_COLUMNS = [ + "session_number", + "prior_exposures_to_image_set", + "ophys_session_id", + "imaging_plane_group_count", + "imaging_plane_group", + "targeted_areas", + "num_depths_per_area", + "num_targeted_structures", +] def sanitize_data_columns( @@ -107,17 +113,19 @@ def get_behavior_session( row = return_one_dataframe_row_only( input_table=self._behavior_session_table, index_value=behavior_session_id, - table_name="behavior_session_table" + table_name="behavior_session_table", ) row = row.squeeze() - has_file_id = (not pd.isna(row[self.cache.file_id_column]) - and row[self.cache.file_id_column] > 0) + has_file_id = ( + not pd.isna(row[self.cache.file_id_column]) + and row[self.cache.file_id_column] > 0 + ) if not has_file_id: oeid = row.ophys_experiment_id[0] row = return_one_dataframe_row_only( input_table=self._ophys_experiment_table, index_value=oeid, - table_name="ophys_experiment_table" + table_name="ophys_experiment_table", ) file_id = str(int(row[self.cache.file_id_column])) data_path = self._get_data_path(file_id=file_id) @@ -141,8 +149,7 @@ def get_behavior_ophys_experiment( row = return_one_dataframe_row_only( input_table=self._ophys_experiment_table, index_value=ophys_experiment_id, - table_name="ophys_experiment_table" - + table_name="ophys_experiment_table", ) file_id = str(int(row[self.cache.file_id_column])) data_path = self._get_data_path(file_id=file_id) diff --git a/allensdk/brain_observatory/behavior/behavior_project_cache/tables/util/prior_exposure_processing.py b/allensdk/brain_observatory/behavior/behavior_project_cache/tables/util/prior_exposure_processing.py index e813228a1..a355304f4 100644 --- a/allensdk/brain_observatory/behavior/behavior_project_cache/tables/util/prior_exposure_processing.py +++ b/allensdk/brain_observatory/behavior/behavior_project_cache/tables/util/prior_exposure_processing.py @@ -171,6 +171,7 @@ def __get_prior_exposure_count( elif agg_method == "cumsum": df["to"] = to df_index_name = df.index.name + def cumsum(x): return x.cumsum().shift(fill_value=0).astype("int64") @@ -184,8 +185,7 @@ def cumsum(x): return counts.reindex(index) -def add_experience_level_ophys( - input_df: pd.DataFrame) -> pd.DataFrame: +def add_experience_level_ophys(input_df: pd.DataFrame) -> pd.DataFrame: """ adds a column to ophys tables that contains a string indicating whether a session had exposure level of Familiar, @@ -210,36 +210,34 @@ def add_experience_level_ophys( # do not modify in place table = input_df.copy(deep=True) - session_number = 'session_number' \ - if 'session_number' in table.columns else 'session' + session_number = ( + "session_number" if "session_number" in table.columns else "session" + ) # add experience_level column with strings indicating relevant conditions - table['experience_level'] = 'None' + table["experience_level"] = "None" - session_training = table.session_type.str.startswith('TRAINING') + session_training = table.session_type.str.startswith("TRAINING") train_indices = table[session_training].index.values - table.loc[train_indices, 'experience_level'] = 'Training' + table.loc[train_indices, "experience_level"] = "Training" session_0123 = table[session_number].isin([0, 1, 2, 3]) familiar_indices = table[session_0123].index.values - table.loc[familiar_indices, 'experience_level'] = 'Familiar' + table.loc[familiar_indices, "experience_level"] = "Familiar" session_456 = table[session_number].isin([4, 5, 6]) - zero_prior_exp = (table.prior_exposures_to_image_set == 0) + zero_prior_exp = table.prior_exposures_to_image_set == 0 - novel_indices = table[session_456 - & zero_prior_exp].index.values + novel_indices = table[session_456 & zero_prior_exp].index.values - table.loc[novel_indices, 'experience_level'] = 'Novel 1' + table.loc[novel_indices, "experience_level"] = "Novel 1" session_456 = table[session_number].isin([4, 5, 6]) - nonzero_prior_exp = (table.prior_exposures_to_image_set != 0) - novel_gt_1_indices = table[ - session_456 - & nonzero_prior_exp].index.values + nonzero_prior_exp = table.prior_exposures_to_image_set != 0 + novel_gt_1_indices = table[session_456 & nonzero_prior_exp].index.values - table.loc[novel_gt_1_indices, 'experience_level'] = 'Novel >1' + table.loc[novel_gt_1_indices, "experience_level"] = "Novel >1" return table diff --git a/allensdk/brain_observatory/behavior/dprime.py b/allensdk/brain_observatory/behavior/dprime.py index d731a8b98..0c6e41ec3 100644 --- a/allensdk/brain_observatory/behavior/dprime.py +++ b/allensdk/brain_observatory/behavior/dprime.py @@ -1,38 +1,57 @@ -import pandas as pd import numpy as np -from scipy.stats import norm - +import pandas as pd from allensdk import one - +from scipy.stats import norm SLIDING_WINDOW = 100 + def get_go_responses(hit=None, miss=None, aborted=None): assert len(hit) == len(miss) == len(aborted) not_aborted = np.logical_not(np.array(aborted, dtype=bool)) hit = np.array(hit, dtype=bool)[not_aborted] miss = np.array(miss, dtype=bool)[not_aborted] - # Go responses are nan when catch (aborted are masked out); 0 for miss, 1 for hit - # This allows pd.Series.rolling to ignore non-go trial data - go_responses = np.empty_like(hit, dtype='float') - go_responses.fill(float('nan')) + # Go responses are nan when catch (aborted are masked out); 0 for miss, 1 + # for hit. This allows pd.Series.rolling to ignore non-go trial data + go_responses = np.empty_like(hit, dtype="float") + go_responses.fill(float("nan")) go_responses[hit] = 1 go_responses[miss] = 0 return go_responses -def get_hit_rate(hit=None, miss=None, aborted=None, sliding_window=SLIDING_WINDOW): +def get_hit_rate( + hit=None, miss=None, aborted=None, sliding_window=SLIDING_WINDOW +): go_responses = get_go_responses(hit=hit, miss=miss, aborted=aborted) - hit_rate = pd.Series(go_responses).rolling(window=sliding_window, min_periods=0).mean().values + hit_rate = ( + pd.Series(go_responses) + .rolling(window=sliding_window, min_periods=0) + .mean() + .values + ) return hit_rate -def get_trial_count_corrected_hit_rate(hit=None, miss=None, aborted=None, sliding_window=SLIDING_WINDOW): +def get_trial_count_corrected_hit_rate( + hit=None, miss=None, aborted=None, sliding_window=SLIDING_WINDOW +): go_responses = get_go_responses(hit=hit, miss=miss, aborted=aborted) - go_responses_count = pd.Series(go_responses).rolling(window=sliding_window, min_periods=0).count() - hit_rate = pd.Series(go_responses).rolling(window=sliding_window, min_periods=0).mean().values - trial_count_corrected_hit_rate = np.vectorize(trial_number_limit)(hit_rate, go_responses_count) + go_responses_count = ( + pd.Series(go_responses) + .rolling(window=sliding_window, min_periods=0) + .count() + ) + hit_rate = ( + pd.Series(go_responses) + .rolling(window=sliding_window, min_periods=0) + .mean() + .values + ) + trial_count_corrected_hit_rate = np.vectorize(trial_number_limit)( + hit_rate, go_responses_count + ) return trial_count_corrected_hit_rate @@ -42,35 +61,73 @@ def get_catch_responses(correct_reject=None, false_alarm=None, aborted=None): correct_reject = np.array(correct_reject, dtype=bool)[not_aborted] false_alarm = np.array(false_alarm, dtype=bool)[not_aborted] - # Catch responses are nan when go (aborted are masked out); 0 for correct-rejection, 1 for false-alarm - # This allows pd.Series.rolling to ignore non-catch trial data - catch_responses = np.empty_like(correct_reject, dtype='float') - catch_responses.fill(float('nan')) + # Catch responses are nan when go (aborted are masked out); 0 for + # correct-rejection, 1 for false-alarm This allows pd.Series.rolling to + # ignore non-catch trial data + catch_responses = np.empty_like(correct_reject, dtype="float") + catch_responses.fill(float("nan")) catch_responses[false_alarm] = 1 catch_responses[correct_reject] = 0 return catch_responses -def get_false_alarm_rate(correct_reject=None, false_alarm=None, aborted=None, sliding_window=SLIDING_WINDOW): - catch_responses = get_catch_responses(correct_reject=correct_reject, false_alarm=false_alarm, aborted=aborted) - false_alarm_rate = pd.Series(catch_responses).rolling(window=sliding_window, min_periods=0).mean().values +def get_false_alarm_rate( + correct_reject=None, + false_alarm=None, + aborted=None, + sliding_window=SLIDING_WINDOW, +): + catch_responses = get_catch_responses( + correct_reject=correct_reject, false_alarm=false_alarm, aborted=aborted + ) + false_alarm_rate = ( + pd.Series(catch_responses) + .rolling(window=sliding_window, min_periods=0) + .mean() + .values + ) return false_alarm_rate -def get_trial_count_corrected_false_alarm_rate(correct_reject=None, false_alarm=None, aborted=None, sliding_window=SLIDING_WINDOW): - catch_responses = get_catch_responses(correct_reject=correct_reject, false_alarm=false_alarm, aborted=aborted) - catch_responses_count = pd.Series(catch_responses).rolling(window=sliding_window, min_periods=0).count() - false_alarm_rate = pd.Series(catch_responses).rolling(window=sliding_window, min_periods=0).mean().values - trial_count_corrected_false_alarm_rate = np.vectorize(trial_number_limit)(false_alarm_rate, catch_responses_count) +def get_trial_count_corrected_false_alarm_rate( + correct_reject=None, + false_alarm=None, + aborted=None, + sliding_window=SLIDING_WINDOW, +): + catch_responses = get_catch_responses( + correct_reject=correct_reject, false_alarm=false_alarm, aborted=aborted + ) + catch_responses_count = ( + pd.Series(catch_responses) + .rolling(window=sliding_window, min_periods=0) + .count() + ) + false_alarm_rate = ( + pd.Series(catch_responses) + .rolling(window=sliding_window, min_periods=0) + .mean() + .values + ) + trial_count_corrected_false_alarm_rate = np.vectorize(trial_number_limit)( + false_alarm_rate, catch_responses_count + ) return trial_count_corrected_false_alarm_rate -def get_rolling_dprime(rolling_hit_rate, rolling_fa_rate, sliding_window=SLIDING_WINDOW): - return np.array([get_dprime(hr, far, sliding_window=SLIDING_WINDOW) for hr, far in zip(rolling_hit_rate, rolling_fa_rate)]) +def get_rolling_dprime( + rolling_hit_rate, rolling_fa_rate, sliding_window=SLIDING_WINDOW +): + return np.array( + [ + get_dprime(hr, far, sliding_window=SLIDING_WINDOW) + for hr, far in zip(rolling_hit_rate, rolling_fa_rate) + ] + ) def get_dprime(hit_rate, fa_rate, sliding_window=SLIDING_WINDOW): - """ calculates the d-prime for a given hit rate and false alarm rate + """calculates the d-prime for a given hit rate and false alarm rate https://en.wikipedia.org/wiki/Sensitivity_index Parameters ---------- @@ -85,9 +142,9 @@ def get_dprime(hit_rate, fa_rate, sliding_window=SLIDING_WINDOW): d_prime """ - limits = (1/SLIDING_WINDOW, 1 - 1/SLIDING_WINDOW) - assert limits[0] > 0.0, 'limits[0] must be greater than 0.0' - assert limits[1] < 1.0, 'limits[1] must be less than 1.0' + limits = (1 / SLIDING_WINDOW, 1 - 1 / SLIDING_WINDOW) + assert limits[0] > 0.0, "limits[0] must be greater than 0.0" + assert limits[1] < 1.0, "limits[1] must be less than 1.0" Z = norm.ppf # Limit values in order to avoid d' infinity @@ -101,6 +158,6 @@ def trial_number_limit(p, N): if N == 0: return np.nan if not pd.isnull(p): - p = np.max((p, 1. / (2 * N))) - p = np.min((p, 1 - 1. / (2 * N))) + p = np.max((p, 1.0 / (2 * N))) + p = np.min((p, 1 - 1.0 / (2 * N))) return p diff --git a/allensdk/brain_observatory/chisquare_categorical.py b/allensdk/brain_observatory/chisquare_categorical.py index 9b0725f5c..6ad254171 100644 --- a/allensdk/brain_observatory/chisquare_categorical.py +++ b/allensdk/brain_observatory/chisquare_categorical.py @@ -7,172 +7,183 @@ """ # TODO: Fix header -import numpy as np import warnings +import numpy as np + -def chisq_from_stim_table(stim_table, - columns, - mean_sweep_events, - num_shuffles=1000, - verbose=False): +def chisq_from_stim_table( + stim_table, columns, mean_sweep_events, num_shuffles=1000, verbose=False +): # stim_table is a pandas DataFrame with len = num_sweeps - # columns is a list of column names that define the categories (e.g. ['Ori','Contrast']) - # mean_sweep_events is a numpy array with shape (num_sweeps,num_cells) - - sweep_categories = stim_table_to_categories(stim_table,columns,verbose=verbose) - p_vals = compute_chi_shuffle(mean_sweep_events,sweep_categories,num_shuffles=num_shuffles) - + # columns is a list of column names that define the categories (e.g. + # ['Ori','Contrast']) mean_sweep_events is a numpy array with shape + # (num_sweeps,num_cells) + + sweep_categories = stim_table_to_categories( + stim_table, columns, verbose=verbose + ) + p_vals = compute_chi_shuffle( + mean_sweep_events, sweep_categories, num_shuffles=num_shuffles + ) + return p_vals -def compute_chi_shuffle(mean_sweep_events, - sweep_categories, - num_shuffles=1000): +def compute_chi_shuffle( + mean_sweep_events, sweep_categories, num_shuffles=1000 +): # mean_sweep_events is a numpy array with shape (num_sweeps,num_cells) # sweep_conditions is a numpy array with shape (num_sweeps) # sweep_conditions gives the category label for each sweep - - (num_sweeps,num_cells) = np.shape(mean_sweep_events) - + + (num_sweeps, num_cells) = np.shape(mean_sweep_events) + if len(sweep_categories) != num_sweeps: - warnings.warn('sweep_categories and num_sweeps do not match') + warnings.warn("sweep_categories and num_sweeps do not match") return np.nan - sweep_categories_dummy = make_category_dummy(sweep_categories) - - expected = compute_expected(mean_sweep_events,sweep_categories_dummy) - observed = compute_observed(mean_sweep_events,sweep_categories_dummy) - chi_actual = compute_chi(observed,expected) - - chi_shuffle = np.zeros((num_cells,num_shuffles)) + + expected = compute_expected(mean_sweep_events, sweep_categories_dummy) + observed = compute_observed(mean_sweep_events, sweep_categories_dummy) + chi_actual = compute_chi(observed, expected) + + chi_shuffle = np.zeros((num_cells, num_shuffles)) for ns in range(num_shuffles): - shuffle_sweeps = np.random.choice(num_sweeps,size=(num_sweeps,)) + shuffle_sweeps = np.random.choice(num_sweeps, size=(num_sweeps,)) shuffle_sweep_events = mean_sweep_events[shuffle_sweeps] - - shuffle_expected = compute_expected(shuffle_sweep_events,sweep_categories_dummy) - shuffle_observed = compute_observed(shuffle_sweep_events,sweep_categories_dummy) - - chi_shuffle[:,ns] = compute_chi(shuffle_observed,shuffle_expected) - - p_vals = np.mean(chi_actual.reshape(num_cells,1)0,chi,0.0) - return np.sum(chi,axis=1) +def compute_chi(observed, expected): + chi = (observed - expected) ** 2 / expected + chi = np.where(expected > 0, chi, 0.0) + return np.sum(chi, axis=1) + -# %% \ No newline at end of file +# %% diff --git a/allensdk/brain_observatory/ecephys/ecephys_project_api/ecephys_project_warehouse_api.py b/allensdk/brain_observatory/ecephys/ecephys_project_api/ecephys_project_warehouse_api.py index c219a5de2..08ffa7e7a 100644 --- a/allensdk/brain_observatory/ecephys/ecephys_project_api/ecephys_project_warehouse_api.py +++ b/allensdk/brain_observatory/ecephys/ecephys_project_api/ecephys_project_warehouse_api.py @@ -1,17 +1,16 @@ -import re -import json import ast +import json +import re -import pandas as pd import numpy as np +import pandas as pd -from .rma_engine import RmaEngine, AsyncRmaEngine from .ecephys_project_api import EcephysProjectApi -from .utilities import rma_macros, build_and_execute +from .rma_engine import AsyncRmaEngine, RmaEngine +from .utilities import build_and_execute, rma_macros class EcephysProjectWarehouseApi(EcephysProjectApi): - movie_re = re.compile(r".*natural_movie_(?P\d+).npy") scene_re = re.compile(r".*/(?P\d+).tiff") @@ -28,27 +27,38 @@ def get_session_data(self, session_id, **kwargs): "[attachable_type$eq'EcephysSession']" r"[attachable_id$eq{{session_id}}]" ), - engine=self.rma_engine.get_rma_tabular, session_id=session_id + engine=self.rma_engine.get_rma_tabular, + session_id=session_id, ) if well_known_files.shape[0] != 1: - raise ValueError(f"expected exactly 1 nwb file for session {session_id}, found: {well_known_files}") + raise ValueError( + f"expected exactly 1 nwb file for session {session_id}, found: {well_known_files}" # noqa: E501 + ) download_link = well_known_files.iloc[0]["download_link"] return self.rma_engine.stream(download_link) def get_natural_movie_template(self, number): - well_known_files = self.stimulus_templates[self.stimulus_templates["movie_number"] == number] + well_known_files = self.stimulus_templates[ + self.stimulus_templates["movie_number"] == number + ] if well_known_files.shape[0] != 1: - raise ValueError(f"expected exactly one natural movie template with number {number}, found {well_known_files}") + raise ValueError( + f"expected exactly one natural movie template with number {number}, found {well_known_files}" # noqa: E501 + ) download_link = well_known_files.iloc[0]["download_link"] return self.rma_engine.stream(download_link) def get_natural_scene_template(self, number): - well_known_files = self.stimulus_templates[self.stimulus_templates["scene_number"] == number] + well_known_files = self.stimulus_templates[ + self.stimulus_templates["scene_number"] == number + ] if well_known_files.shape[0] != 1: - raise ValueError(f"expected exactly one natural scene template with number {number}, found {well_known_files}") + raise ValueError( + f"expected exactly one natural scene template with number {number}, found {well_known_files}" # noqa: E501 + ) download_link = well_known_files.iloc[0]["download_link"] return self.rma_engine.stream(download_link) @@ -68,7 +78,7 @@ def _list_stimulus_templates(self, ecephys_product_id=714914585): r"[attachable_id$eq{{ecephys_product_id}}]" ), engine=self.rma_engine.get_rma_tabular, - ecephys_product_id=ecephys_product_id + ecephys_product_id=ecephys_product_id, ) scene_number = [] @@ -97,24 +107,29 @@ def get_probe_lfp_data(self, probe_id): "[attachable_type$eq'EcephysProbe']" r"[attachable_id$eq{{probe_id}}]" ), - engine=self.rma_engine.get_rma_tabular, probe_id=probe_id + engine=self.rma_engine.get_rma_tabular, + probe_id=probe_id, ) if well_known_files.shape[0] != 1: - raise ValueError(f"expected exactly 1 LFP NWB file for probe {probe_id}, found: {well_known_files}") + raise ValueError( + f"expected exactly 1 LFP NWB file for probe {probe_id}, found: {well_known_files}" # noqa: E501 + ) download_link = well_known_files.loc[0, "download_link"] return self.rma_engine.stream(download_link) - def get_sessions(self, session_ids=None, has_eye_tracking=None, stimulus_names=None): + def get_sessions( + self, session_ids=None, has_eye_tracking=None, stimulus_names=None + ): response = build_and_execute( ( "{% import 'rma_macros' as rm %}" "{% import 'macros' as m %}" "criteria=model::EcephysSession" r"{{rm.optional_contains('id',session_ids)}}" - r"{%if has_eye_tracking is not none%}[fail_eye_tracking$eq{{m.str(not has_eye_tracking).lower()}}]{%endif%}" - r"{{rm.optional_contains('stimulus_name',stimulus_names,True)}}" + r"{%if has_eye_tracking is not none%}[fail_eye_tracking$eq{{m.str(not has_eye_tracking).lower()}}]{%endif%}" # noqa: E501 + r"{{rm.optional_contains('stimulus_name',stimulus_names,True)}}" # noqa: E501 ",rma::include,specimen(donor(age))" ",well_known_files(well_known_file_type)" ), @@ -122,7 +137,7 @@ def get_sessions(self, session_ids=None, has_eye_tracking=None, stimulus_names=N engine=self.rma_engine.get_rma_tabular, session_ids=session_ids, has_eye_tracking=has_eye_tracking, - stimulus_names=stimulus_names + stimulus_names=stimulus_names, ) response.set_index("id", inplace=True) @@ -152,8 +167,13 @@ def get_sessions(self, session_ids=None, has_eye_tracking=None, stimulus_names=N response["genotype"] = genotype response["has_nwb"] = has_nwb - response.drop(columns=["specimen", "fail_eye_tracking", "well_known_files"], inplace=True) - response.rename(columns={"stimulus_name": "session_type"}, inplace=True) + response.drop( + columns=["specimen", "fail_eye_tracking", "well_known_files"], + inplace=True, + ) + response.rename( + columns={"stimulus_name": "session_type"}, inplace=True + ) return response @@ -169,7 +189,7 @@ def get_probes(self, probe_ids=None, session_ids=None): base=rma_macros(), engine=self.rma_engine.get_rma_tabular, session_ids=session_ids, - probe_ids=probe_ids + probe_ids=probe_ids, ) response.set_index("id", inplace=True) # Clarify name for external users @@ -186,50 +206,62 @@ def get_channels(self, channel_ids=None, probe_ids=None): r"{{rm.optional_contains('ecephys_probe_id',probe_ids)}}" ",rma::include,structure" ",rma::options[tabular$eq'" - "ecephys_channels.id" - ",ecephys_probe_id" - ",local_index" - ",probe_horizontal_position" - ",probe_vertical_position" - ",anterior_posterior_ccf_coordinate" - ",dorsal_ventral_ccf_coordinate" - ",left_right_ccf_coordinate" - ",structures.id as ecephys_structure_id" - ",structures.acronym as ecephys_structure_acronym" + "ecephys_channels.id" + ",ecephys_probe_id" + ",local_index" + ",probe_horizontal_position" + ",probe_vertical_position" + ",anterior_posterior_ccf_coordinate" + ",dorsal_ventral_ccf_coordinate" + ",left_right_ccf_coordinate" + ",structures.id as ecephys_structure_id" + ",structures.acronym as ecephys_structure_acronym" "']" ), base=rma_macros(), engine=self.rma_engine.get_rma_tabular, probe_ids=probe_ids, - channel_ids=channel_ids + channel_ids=channel_ids, ) response.set_index("id", inplace=True) return response - def get_units(self, unit_ids=None, channel_ids=None, probe_ids=None, session_ids=None, *a, **k): + def get_units( + self, + unit_ids=None, + channel_ids=None, + probe_ids=None, + session_ids=None, + *a, + **k, + ): response = build_and_execute( ( "{% import 'macros' as m %}" "criteria=model::EcephysUnit" - r"{% if unit_ids is not none %},rma::criteria[id$in{{m.comma_sep(unit_ids)}}]{% endif %}" - r"{% if channel_ids is not none %},rma::criteria[ecephys_channel_id$in{{m.comma_sep(channel_ids)}}]{% endif %}" - r"{% if probe_ids is not none %},rma::criteria,ecephys_channel(ecephys_probe[id$in{{m.comma_sep(probe_ids)}}]){% endif %}" - r"{% if session_ids is not none %},rma::criteria,ecephys_channel(ecephys_probe(ecephys_session[id$in{{m.comma_sep(session_ids)}}])){% endif %}" + r"{% if unit_ids is not none %},rma::criteria[id$in{{m.comma_sep(unit_ids)}}]{% endif %}" # noqa: E501 + r"{% if channel_ids is not none %},rma::criteria[ecephys_channel_id$in{{m.comma_sep(channel_ids)}}]{% endif %}" # noqa: E501 + r"{% if probe_ids is not none %},rma::criteria,ecephys_channel(ecephys_probe[id$in{{m.comma_sep(probe_ids)}}]){% endif %}" # noqa: E501 + r"{% if session_ids is not none %},rma::criteria,ecephys_channel(ecephys_probe(ecephys_session[id$in{{m.comma_sep(session_ids)}}])){% endif %}" # noqa: E501 ), - base=rma_macros(), engine=self.rma_engine.get_rma_tabular, + base=rma_macros(), + engine=self.rma_engine.get_rma_tabular, session_ids=session_ids, probe_ids=probe_ids, channel_ids=channel_ids, - unit_ids=unit_ids + unit_ids=unit_ids, ) response.set_index("id", inplace=True) return response - def get_unit_analysis_metrics(self, unit_ids=None, ecephys_session_ids=None, session_types=None): - """ Download analysis metrics - precalculated descriptions of unitwise responses to visual stimulation. + def get_unit_analysis_metrics( + self, unit_ids=None, ecephys_session_ids=None, session_types=None + ): + """Download analysis metrics - precalculated descriptions of unitwise + responses to visual stimulation. Parameters ---------- @@ -255,15 +287,15 @@ def get_unit_analysis_metrics(self, unit_ids=None, ecephys_session_ids=None, ses ( "{% import 'macros' as m %}" "criteria=model::EcephysUnitMetricBundle" - r"{% if unit_ids is not none %},rma::criteria[ecephys_unit_id$in{{m.comma_sep(unit_ids)}}]{% endif %}" - r"{% if session_ids is not none %},rma::criteria,ecephys_unit(ecephys_channel(ecephys_probe(ecephys_session[id$in{{m.comma_sep(session_ids)}}]))){% endif %}" - r"{% if session_types is not none %},rma::criteria,ecephys_unit(ecephys_channel(ecephys_probe(ecephys_session[stimulus_name$in{{m.comma_sep(session_types, True)}}]))){% endif %}" + r"{% if unit_ids is not none %},rma::criteria[ecephys_unit_id$in{{m.comma_sep(unit_ids)}}]{% endif %}" # noqa: E501 + r"{% if session_ids is not none %},rma::criteria,ecephys_unit(ecephys_channel(ecephys_probe(ecephys_session[id$in{{m.comma_sep(session_ids)}}]))){% endif %}" # noqa: E501 + r"{% if session_types is not none %},rma::criteria,ecephys_unit(ecephys_channel(ecephys_probe(ecephys_session[stimulus_name$in{{m.comma_sep(session_types, True)}}]))){% endif %}" # noqa: E501 ), base=rma_macros(), engine=self.rma_engine.get_rma_list, session_ids=ecephys_session_ids, unit_ids=unit_ids, - session_types=session_types + session_types=session_types, ) output = [] @@ -278,7 +310,9 @@ def get_unit_analysis_metrics(self, unit_ids=None, ecephys_session_ids=None, ses for colname in output.columns: try: - output[colname] = output.apply(lambda row: ast.literal_eval(str(row[colname])), axis=1) + output[colname] = output.apply( + lambda row: ast.literal_eval(str(row[colname])), axis=1 + ) except ValueError: pass @@ -287,9 +321,10 @@ def get_unit_analysis_metrics(self, unit_ids=None, ecephys_session_ids=None, ses # but switched with one another. This snippet unswitches them. columns = set(output.columns.values.tolist()) if "p_value_rf" in columns and "on_screen_rf" in columns: - pv_is_bool = np.issubdtype(output["p_value_rf"].values[0], bool) - on_screen_is_float = np.issubdtype(output["on_screen_rf"].values[0].dtype, np.floating) + on_screen_is_float = np.issubdtype( + output["on_screen_rf"].values[0].dtype, np.floating + ) # this is not a good test, but it avoids the case where we fix # these in the data for a future release, but diff --git a/allensdk/brain_observatory/ecephys/stimulus_table/naming_utilities.py b/allensdk/brain_observatory/ecephys/stimulus_table/naming_utilities.py index bb49953b6..adf75b8f0 100644 --- a/allensdk/brain_observatory/ecephys/stimulus_table/naming_utilities.py +++ b/allensdk/brain_observatory/ecephys/stimulus_table/naming_utilities.py @@ -3,14 +3,14 @@ import numpy as np - -GABOR_DIAMETER_RE = \ - re.compile(r"gabor_(\d*\.{0,1}\d*)_{0,1}deg(?:_\d+ms){0,1}") +GABOR_DIAMETER_RE = re.compile( + r"gabor_(\d*\.{0,1}\d*)_{0,1}deg(?:_\d+ms){0,1}" +) GENERIC_MOVIE_RE = re.compile( - r"natural_movie_" + - r"(?P\d+|one|two|three|four|five|six|seven|eight|nine)" + - r"(_shuffled){0,1}(_more_repeats){0,1}" + r"natural_movie_" + + r"(?P\d+|one|two|three|four|five|six|seven|eight|nine)" + + r"(_shuffled){0,1}(_more_repeats){0,1}" ) DIGIT_NAMES = { "1": "one", @@ -28,8 +28,7 @@ def drop_empty_columns(table): - """ Remove from the stimulus table columns whose values are all nan - """ + """Remove from the stimulus table columns whose values are all nan""" to_drop = [] @@ -42,7 +41,7 @@ def drop_empty_columns(table): def collapse_columns(table): - """ merge, where possible, columns that describe the same parameter. This + """merge, where possible, columns that describe the same parameter. This is pretty conservative - it only matches columns by capitalization and it only overrides nans. """ @@ -53,7 +52,6 @@ def collapse_columns(table): for col in table.columns: for transformed in (col.upper(), col.capitalize()): if transformed in colnames and col != transformed: - col_notna = ~(table[col].isna()) trans_notna = ~(table[transformed].isna()) if (col_notna & trans_notna).sum() != 0: @@ -77,24 +75,23 @@ def add_number_to_shuffled_movie( template="natural_movie_{}_shuffled", tmp_colname="__movie_number__", ): - """ - """ + """ """ if not table[stim_colname].str.contains(SHUFFLED_MOVIE_RE).any(): return table table = table.copy() - table[tmp_colname] = \ - table[stim_colname].str.extract(natural_movie_re, - expand=True)["number"] + table[tmp_colname] = table[stim_colname].str.extract( + natural_movie_re, expand=True + )["number"] unique_numbers = [ item for item in table[tmp_colname].dropna(inplace=False).unique() ] if len(unique_numbers) != 1: raise ValueError( - "unable to uniquely determine a movie number for this session. " + - f"Candidates: {unique_numbers}" + "unable to uniquely determine a movie number for this session. " + + f"Candidates: {unique_numbers}" ) movie_number = unique_numbers[0] @@ -119,7 +116,7 @@ def standardize_movie_numbers( digit_names=DIGIT_NAMES, stim_colname="stimulus_name", ): - """ Natural movie stimuli in visual coding are numbered using words, like + """Natural movie stimuli in visual coding are numbered using words, like "natural_movie_two" rather than "natural_movie_2". This function ensures that all of the natural movie stimuli in an experiment are named by that convention. @@ -153,14 +150,15 @@ def replace(match_obj): warnings.filterwarnings("ignore", category=UserWarning) movie_rows = table[stim_colname].str.contains(movie_re, na=False) - table.loc[movie_rows, stim_colname] = \ - table.loc[movie_rows, stim_colname].str.replace(numeral_re, replace, regex=True) + table.loc[movie_rows, stim_colname] = table.loc[ + movie_rows, stim_colname + ].str.replace(numeral_re, replace, regex=True) return table def map_stimulus_names(table, name_map=None, stim_colname="stimulus_name"): - """ Applies a mappting to the stimulus names in a stimulus table + """Applies a mappting to the stimulus names in a stimulus table Parameters ---------- @@ -188,7 +186,6 @@ def map_stimulus_names(table, name_map=None, stim_colname="stimulus_name"): def map_column_names(table, name_map=None, ignore_case=True): - if ignore_case and name_map is not None: name_map = {key.lower(): value for key, value in name_map.items()} diff --git a/allensdk/brain_observatory/receptive_field_analysis/chisquarerf.py b/allensdk/brain_observatory/receptive_field_analysis/chisquarerf.py index 8fa1094c7..dc5fd0b28 100644 --- a/allensdk/brain_observatory/receptive_field_analysis/chisquarerf.py +++ b/allensdk/brain_observatory/receptive_field_analysis/chisquarerf.py @@ -38,7 +38,6 @@ import scipy.ndimage.filters as filt import scipy.stats as stats - ON_LUMINANCE = 255 OFF_LUMINANCE = 0 @@ -49,15 +48,15 @@ def chi_square_binary(events, LSN_template): # *****INPUT***** # events: 2D numpy bool with shape (num_trials,num_cells) for presence # or absence of response on a given trial - # LSN_template: 3D numpy int8 with shape (num_trials,num_y_pixels,num_x_pixels) - # for luminance at each pixel location + # LSN_template: 3D numpy int8 with shape (num_trials,num_y_pixels, + # num_x_pixels) for luminance at each pixel location # # *****OUTPUT***** - # chi_square_grid_NLL: 3D numpy float with shape (num_cells,num_y_pixels,num_x_pixels) - # that gives the p value for the hypothesis that a receptive field is contained - # within a 7x7 pixel mask centered on a given pixel location as measured - # by a chi-square test for the responses among the pixels (both on and off) - # that fall within the mask. + # chi_square_grid_NLL: 3D numpy float with shape (num_cells,num_y_pixels, + # num_x_pixels) that gives the p value for the hypothesis that a + # receptive field is contained within a 7x7 pixel mask centered on a + # given pixel location as measured by a chi-square test for the responses + # among the pixels (both on and off) that fall within the mask. num_trials = np.shape(events)[0] num_cells = np.shape(events)[1] @@ -77,40 +76,44 @@ def chi_square_binary(events, LSN_template): # trials_per_pixel has shape (num_y,num_x,2) trials_per_pixel = np.sum(trial_matrix, axis=3) - # get the sum of the number of events across all trials that each pixel is active - # events_per_pixel has shape (num_cells,num_y,num_x,2) + # get the sum of the number of events across all trials that each pixel is + # active events_per_pixel has shape (num_cells,num_y,num_x,2) events_per_pixel = get_events_per_pixel(events, trial_matrix) # smooth stimulus-triggered average spatially with a gaussian for n in range(num_cells): for on_off in range(2): - events_per_pixel[n, :, :, on_off] = smooth_STA(events_per_pixel[n, :, :, on_off]) + events_per_pixel[n, :, :, on_off] = smooth_STA( + events_per_pixel[n, :, :, on_off] + ) # calculate the p_value for each exclusion region chi_square_grid = np.zeros((num_cells, num_y, num_x)) for y in range(num_y): for x in range(num_x): - exclusion_mask = np.ones((num_y, num_x, 2)) * disc_masks[y, x, :, :].reshape(num_y, num_x, 1) - p_vals, __ = chi_square_within_mask(exclusion_mask, events_per_pixel, trials_per_pixel) + exclusion_mask = np.ones((num_y, num_x, 2)) * disc_masks[ + y, x, :, : + ].reshape(num_y, num_x, 1) + p_vals, __ = chi_square_within_mask( + exclusion_mask, events_per_pixel, trials_per_pixel + ) chi_square_grid[:, y, x] = p_vals return chi_square_grid -def get_peak_significance(chi_square_grid_NLL, - LSN_template, - alpha=0.05): +def get_peak_significance(chi_square_grid_NLL, LSN_template, alpha=0.05): # ****INPUT***** # chi_square_grid_NLL: result of chi_square_binary(events,LSN_template) - # LSN_template: 3D numpy int8 with shape (num_trials,num_y_pixels,num_x_pixels) - # for luminance at each pixel location + # LSN_template: 3D numpy int8 with shape (num_trials,num_y_pixels, + # num_x_pixels) for luminance at each pixel location # # *****OUTPUT***** # significant_cells: 1D numpy bool with shape (num_cells,) that indicates # whether or not each cell has a location with a significant chance of a # true receptive field - # best_exclusion_region_list: list of 2D numpy bool with len (num_cells). For each cell, - # the array is a mask for the best pixel, or all false + # best_exclusion_region_list: list of 2D numpy bool with len (num_cells). + # For each cell, the array is a mask for the best pixel, or all false num_cells = np.shape(chi_square_grid_NLL)[0] num_y = np.shape(chi_square_grid_NLL)[1] @@ -125,36 +128,55 @@ def get_peak_significance(chi_square_grid_NLL, # find the smallest p-value and determine if it's significant significant_cells = np.zeros((num_cells)).astype(bool) best_p = np.zeros((num_cells)) - p_value_correction_factor_per_pixel = (1.0 * num_y * num_x / pixels_per_mask_per_pixel) + p_value_correction_factor_per_pixel = ( + 1.0 * num_y * num_x / pixels_per_mask_per_pixel + ) best_exclusion_region_list = [] corrected_p_value_array_list = [] for n in range(num_cells): - # Sidak correction: - p_value_corrected_per_pixel = 1-np.power((1-chi_square_grid[n, :,:]), p_value_correction_factor_per_pixel) + p_value_corrected_per_pixel = 1 - np.power( + (1 - chi_square_grid[n, :, :]), p_value_correction_factor_per_pixel + ) corrected_p_value_array_list.append(p_value_corrected_per_pixel) - y, x = np.unravel_index(p_value_corrected_per_pixel.argmin(), (num_y, num_x)) + y, x = np.unravel_index( + p_value_corrected_per_pixel.argmin(), (num_y, num_x) + ) # if more than one p-value that maxes out, use the median location if np.sum(p_value_corrected_per_pixel == 0.0) > 1: - - y, x = np.unravel_index(np.argwhere(p_value_corrected_per_pixel.flatten() == 0.0)[:, 0], (num_y, num_x)) + y, x = np.unravel_index( + np.argwhere(p_value_corrected_per_pixel.flatten() == 0.0)[ + :, 0 + ], + (num_y, num_x), + ) center_y, center_x = locate_median(y, x) - best_p[n] = p_value_corrected_per_pixel[y,x] + best_p[n] = p_value_corrected_per_pixel[y, x] if best_p[n] < alpha: significant_cells[n] = True - best_exclusion_region_list.append(disc_masks[y, x, :,:].astype(bool)) + best_exclusion_region_list.append( + disc_masks[y, x, :, :].astype(bool) + ) else: - best_exclusion_region_list.append(np.zeros((disc_masks.shape[0], disc_masks.shape[1]), dtype=bool)) + best_exclusion_region_list.append( + np.zeros( + (disc_masks.shape[0], disc_masks.shape[1]), dtype=bool + ) + ) - return significant_cells, best_p, corrected_p_value_array_list, best_exclusion_region_list + return ( + significant_cells, + best_p, + corrected_p_value_array_list, + best_exclusion_region_list, + ) def locate_median(y, x): - med_x = np.median(x) med_y = np.median(y) center_x = x[0] @@ -166,7 +188,7 @@ def locate_median(y, x): dc_x = center_x - med_x dc_y = center_y - med_y - if np.sqrt(dx ** 2 + dy ** 2) < np.sqrt(dc_x ** 2 + dc_y ** 2): + if np.sqrt(dx**2 + dy**2) < np.sqrt(dc_x**2 + dc_y**2): center_x = x[i] center_y = y[i] @@ -178,29 +200,29 @@ def pvalue_to_NLL(p_values, max_NLL=10.0): def NLL_to_pvalue(NLLs, log_base=10.0): - return (log_base ** (-NLLs)) + return log_base ** (-NLLs) def get_events_per_pixel(responses_np, trial_matrix): - '''Obtain a matrix linking cellular responses to pixel activity. + """Obtain a matrix linking cellular responses to pixel activity. Parameters ---------- responses_np : np.ndarray - Dimensions are (nTrials, nCells). Boolean values indicate presence/absence - of a response on a given trial. + Dimensions are (nTrials, nCells). Boolean values indicate + presence/absence of a response on a given trial. trial_matrix : np.ndarray - Dimensions are (nYPixels, nXPixels, {on, off}, nTrials). Boolean values + Dimensions are (nYPixels, nXPixels, {on, off}, nTrials). Boolean values indicate that a pixel was on/off on a particular trial. Returns ------- events_per_pixel : np.ndarray - Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Values for each - cell, pixel, and on/off state are the sum of events for that cell across - all trials where the pixel was in the on/off state. + Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Values for each + cell, pixel, and on/off state are the sum of events for that cell + across all trials where the pixel was in the on/off state. - ''' + """ num_cells = np.shape(responses_np)[1] num_y = np.shape(trial_matrix)[0] @@ -211,20 +233,22 @@ def get_events_per_pixel(responses_np, trial_matrix): for x in range(num_x): for on_off in range(2): frames = np.argwhere(trial_matrix[y, x, on_off, :])[:, 0] - events_per_pixel[:, y, x, on_off] = np.sum(responses_np[frames, :], axis=0) + events_per_pixel[:, y, x, on_off] = np.sum( + responses_np[frames, :], axis=0 + ) return events_per_pixel def smooth_STA(STA, gauss_std=0.75, total_degrees=64): - '''Smooth an image by convolution with a gaussian kernel + """Smooth an image by convolution with a gaussian kernel Parameters ---------- STA : np.ndarray Input image gauss_std : numeric, optional - Standard deviation of the gaussian kernel. Will be applied to the + Standard deviation of the gaussian kernel. Will be applied to the upsampled image, so units are visual degrees. Default is 0.75 total_degrees : int, optional Size in visual degrees of the input image along its zeroth (row) axis. @@ -234,19 +258,23 @@ def smooth_STA(STA, gauss_std=0.75, total_degrees=64): ------- STA_smoothed : np.ndarray Smoothed image - ''' + """ deg_per_pnt = total_degrees // STA.shape[0] STA_interpolated = interpolate_RF(STA, deg_per_pnt) - STA_interpolated_smoothed = filt.gaussian_filter(STA_interpolated, gauss_std) - STA_smoothed = deinterpolate_RF(STA_interpolated_smoothed, STA.shape[1], STA.shape[0], deg_per_pnt) + STA_interpolated_smoothed = filt.gaussian_filter( + STA_interpolated, gauss_std + ) + STA_smoothed = deinterpolate_RF( + STA_interpolated_smoothed, STA.shape[1], STA.shape[0], deg_per_pnt + ) return STA_smoothed def interpolate_RF(rf_map, deg_per_pnt): - '''Upsample an image - + """Upsample an image + Parameters ---------- rf_map : np.ndarray @@ -258,16 +286,32 @@ def interpolate_RF(rf_map, deg_per_pnt): ------- interpolated : np.ndarray Upsampled image - ''' + """ x_pnts = np.shape(rf_map)[1] y_pnts = np.shape(rf_map)[0] - x_coor = np.arange(-(x_pnts - 1) * deg_per_pnt / 2, (x_pnts + 1) * deg_per_pnt / 2, deg_per_pnt) - y_coor = np.arange(-(y_pnts - 1) * deg_per_pnt / 2, (y_pnts + 1) * deg_per_pnt / 2, deg_per_pnt) - - x_interpolated = np.arange(-(x_pnts - 1) * deg_per_pnt / 2, deg_per_pnt / 2 + (x_pnts / 2 - 1) * deg_per_pnt + 1, 1) - y_interpolated = np.arange(-(y_pnts - 1) * deg_per_pnt / 2, deg_per_pnt / 2 + (y_pnts / 2 - 1) * deg_per_pnt + 1, 1) + x_coor = np.arange( + -(x_pnts - 1) * deg_per_pnt / 2, + (x_pnts + 1) * deg_per_pnt / 2, + deg_per_pnt, + ) + y_coor = np.arange( + -(y_pnts - 1) * deg_per_pnt / 2, + (y_pnts + 1) * deg_per_pnt / 2, + deg_per_pnt, + ) + + x_interpolated = np.arange( + -(x_pnts - 1) * deg_per_pnt / 2, + deg_per_pnt / 2 + (x_pnts / 2 - 1) * deg_per_pnt + 1, + 1, + ) + y_interpolated = np.arange( + -(y_pnts - 1) * deg_per_pnt / 2, + deg_per_pnt / 2 + (y_pnts / 2 - 1) * deg_per_pnt + 1, + 1, + ) interpolated = si.interp2d(x_coor, y_coor, rf_map) interpolated = interpolated(x_interpolated, y_interpolated) @@ -276,7 +320,7 @@ def interpolate_RF(rf_map, deg_per_pnt): def deinterpolate_RF(rf_map, x_pnts, y_pnts, deg_per_pnt): - '''Downsample an image + """Downsample an image Parameters ---------- @@ -288,18 +332,26 @@ def deinterpolate_RF(rf_map, x_pnts, y_pnts, deg_per_pnt): Count of sample points along the zeroth (row) axis deg_per_pnt : numeric scale factor - + Returns ------- sampled_yx : np.ndarray Downsampled image - ''' + """ # x_pnts = 28 # y_pnts = 16 - x_interpolated = np.arange(-(x_pnts - 1) * deg_per_pnt / 2, deg_per_pnt / 2 + (x_pnts / 2 - 1) * deg_per_pnt + 1, 1) - y_interpolated = np.arange(-(y_pnts - 1) * deg_per_pnt / 2, deg_per_pnt / 2 + (y_pnts / 2 - 1) * deg_per_pnt + 1, 1) + x_interpolated = np.arange( + -(x_pnts - 1) * deg_per_pnt / 2, + deg_per_pnt / 2 + (x_pnts / 2 - 1) * deg_per_pnt + 1, + 1, + ) + y_interpolated = np.arange( + -(y_pnts - 1) * deg_per_pnt / 2, + deg_per_pnt / 2 + (y_pnts / 2 - 1) * deg_per_pnt + 1, + 1, + ) x_deinterpolate = np.arange(0, len(x_interpolated), deg_per_pnt) y_deinterpolate = np.arange(0, len(y_interpolated), deg_per_pnt) @@ -311,31 +363,31 @@ def deinterpolate_RF(rf_map, x_pnts, y_pnts, deg_per_pnt): def chi_square_within_mask(exclusion_mask, events_per_pixel, trials_per_pixel): - '''Determine if cells respond preferentially to on/off pixels in a mask using - a chi2 test. - + """Determine if cells respond preferentially to on/off pixels in a mask + using a chi2 test. + Parameters ---------- exclusion_mask : np.ndarray - Dimensions are (nYPixels, nXPixels, {on, off}). Integer indicator for INCLUSION (!) - of a pixel within the testing region. + Dimensions are (nYPixels, nXPixels, {on, off}). Integer indicator for + INCLUSION (!) of a pixel within the testing region. events_per_pixel : np.ndarray - Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Integer values + Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Integer values are response counts by cell to on/off luminance at each pixel. trials_per_pixel : np.ndarray - Dimensions are (nYPixels, nXPixels, {on, off}). Integer values are + Dimensions are (nYPixels, nXPixels, {on, off}). Integer values are counts of trials where a pixel is on/off. Returns ------- - p_vals : np.ndarray - One-dimensional, of length nCells. Float values are p-values - for the hypothesis that a given cell has a receptive field within the + p_vals : np.ndarray + One-dimensional, of length nCells. Float values are p-values + for the hypothesis that a given cell has a receptive field within the exclusion mask. chi : np.ndarray - Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Values (float) + Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Values (float) are squared residual event counts divided by expected event counts. - ''' + """ num_y = np.shape(exclusion_mask)[0] num_x = np.shape(exclusion_mask)[1] @@ -344,12 +396,16 @@ def chi_square_within_mask(exclusion_mask, events_per_pixel, trials_per_pixel): degrees_of_freedom = int(np.sum(exclusion_mask)) - 1 # observed_by_pixel has shape (num_cells,num_y,num_x,2) - expected_by_pixel = get_expected_events_by_pixel(exclusion_mask, events_per_pixel, trials_per_pixel) - observed_by_pixel = (events_per_pixel * exclusion_mask.reshape(1, num_y, num_x, 2)).astype(float) + expected_by_pixel = get_expected_events_by_pixel( + exclusion_mask, events_per_pixel, trials_per_pixel + ) + observed_by_pixel = ( + events_per_pixel * exclusion_mask.reshape(1, num_y, num_x, 2) + ).astype(float) # calculate test statistic given observed and expected residual_by_pixel = observed_by_pixel - expected_by_pixel - chi = (residual_by_pixel ** 2) / expected_by_pixel + chi = (residual_by_pixel**2) / expected_by_pixel chi_sum = np.nansum(chi, axis=(1, 2, 3)) # get p-value given test statistic and degrees of freedom @@ -358,28 +414,30 @@ def chi_square_within_mask(exclusion_mask, events_per_pixel, trials_per_pixel): return p_vals, chi -def get_expected_events_by_pixel(exclusion_mask, events_per_pixel, trials_per_pixel): - '''Calculate expected number of events per pixel +def get_expected_events_by_pixel( + exclusion_mask, events_per_pixel, trials_per_pixel +): + """Calculate expected number of events per pixel Parameters ---------- exclusion_mask : np.ndarray - Dimensions are (nYPixels, nXPixels, {on, off}). Integer indicator for INCLUSION (!) - of a pixel within the testing region. + Dimensions are (nYPixels, nXPixels, {on, off}). Integer indicator for + INCLUSION (!) of a pixel within the testing region. events_per_pixel : np.ndarray - Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Integer values + Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Integer values are response counts by cell to on/off luminance at each pixel. trials_per_pixel : np.ndarray - Dimensions are (nYPixels, nXPixels, {on, off}). Integer values are + Dimensions are (nYPixels, nXPixels, {on, off}). Integer values are counts of trials where a pixel is on/off. Returns ------- np.ndarray : - Dimensions (nCells, nYPixels, nXPixels, {on, off}). Float values are + Dimensions (nCells, nYPixels, nXPixels, {on, off}). Float values are pixelwise counts of events expected if events are evenly distributed in mask across trials. - ''' + """ num_y = np.shape(exclusion_mask)[0] num_x = np.shape(exclusion_mask)[1] @@ -393,79 +451,88 @@ def get_expected_events_by_pixel(exclusion_mask, events_per_pixel, trials_per_pi total_trials = np.sum(masked_trials).astype(float) total_events_by_cell = np.sum(masked_events, axis=(1, 2, 3)).astype(float) - + expected_by_cell_per_trial = total_events_by_cell / total_trials - return masked_trials * expected_by_cell_per_trial.reshape(num_cells, 1, 1, 1) + return masked_trials * expected_by_cell_per_trial.reshape( + num_cells, 1, 1, 1 + ) -def build_trial_matrix(LSN_template, - num_trials, - on_off_luminance=(ON_LUMINANCE, OFF_LUMINANCE)): - '''Construct indicator arrays for on/off pixels across trials. +def build_trial_matrix( + LSN_template, num_trials, on_off_luminance=(ON_LUMINANCE, OFF_LUMINANCE) +): + """Construct indicator arrays for on/off pixels across trials. Parameters ---------- LSN_template : np.ndarray - Dimensions are (nTrials, nYPixels, nXPixels). Luminance values per pixel - and trial. The size of the first dimension may be larger than the num_trials - argument (in which case only the first num_trials slices will be used) - but may not be smaller. + Dimensions are (nTrials, nYPixels, nXPixels). Luminance values per + pixel and trial. The size of the first dimension may be larger than the + num_trials argument (in which case only the first num_trials slices + will be used) but may not be smaller. num_trials : int The number of trials (left-justified) to build indicators for. on_off_luminance : array-like, optional - The zeroth element is the luminance value of a pixel when on, the first when off. - Defaults are [255, 0]. + The zeroth element is the luminance value of a pixel when on, the first + when off. Defaults are [255, 0]. Returns ------- trial_mat : np.ndarray - Dimensions are (nYPixels, nXPixels, {on, off}, nTrials). Boolean values + Dimensions are (nYPixels, nXPixels, {on, off}, nTrials). Boolean values indicate that a pixel was on/off on a particular trial. - ''' + """ _, num_y, num_x = np.shape(LSN_template) - trial_mat = np.zeros( (num_y, num_x, 2, num_trials), dtype=bool ) + trial_mat = np.zeros((num_y, num_x, 2, num_trials), dtype=bool) for y in range(num_y): for x in range(num_x): for oo, on_off in enumerate(on_off_luminance): - - frame = np.argwhere( LSN_template[:num_trials, y, x] == on_off )[:, 0] + frame = np.argwhere(LSN_template[:num_trials, y, x] == on_off)[ + :, 0 + ] trial_mat[y, x, oo, frame] = True return trial_mat -def get_disc_masks(LSN_template, radius=3, on_luminance=ON_LUMINANCE, off_luminance=OFF_LUMINANCE): - '''Obtain an indicator mask surrounding each pixel. The mask is a square, excluding pixels which - are coactive on any trial with the main pixel. +def get_disc_masks( + LSN_template, + radius=3, + on_luminance=ON_LUMINANCE, + off_luminance=OFF_LUMINANCE, +): + """Obtain an indicator mask surrounding each pixel. The mask is a square, + excluding pixels which are coactive on any trial with the main pixel. Parameters ---------- LSN_template : np.ndarray - Dimensions are (nTrials, nYPixels, nXPixels). Luminance values per pixel - and trial. + Dimensions are (nTrials, nYPixels, nXPixels). Luminance values per + pixel and trial. radius : int The base mask will be a box whose sides are 2 * radius + 1 in length. on_luminance : int, optional The value of the luminance for on trials. Default is 255 off_luminance : int, optional The value of the luminance for off trials. Default is 0 - + Returns ------- masks : np.ndarray - Dimensions are (nYPixels, nXPixels, nYPixels, nXPixels). The first 2 - dimensions describe the pixel from which the mask was computed. The last - 2 serve as the dimensions of the mask images themselves. Masks are binary - arrays of type float, with 1 indicating inside, 0 outside. - ''' + Dimensions are (nYPixels, nXPixels, nYPixels, nXPixels). The first 2 + dimensions describe the pixel from which the mask was computed. The + last 2 serve as the dimensions of the mask images themselves. Masks are + binary arrays of type float, with 1 indicating inside, 0 outside. + """ num_y = np.shape(LSN_template)[1] num_x = np.shape(LSN_template)[2] - # convert template to true on trials a pixel is not gray and false when gray + # convert template to true on trials a pixel is not gray and false when + # gray LSN_binary = np.where(LSN_template == off_luminance, 1, LSN_template) LSN_binary = np.where(LSN_binary == on_luminance, 1, LSN_binary) LSN_binary = np.where(LSN_binary == 1, 1.0, 0.0) @@ -476,10 +543,14 @@ def get_disc_masks(LSN_template, radius=3, on_luminance=ON_LUMINANCE, off_lumina masks = np.zeros((num_y, num_x, num_y, num_x)) for y in range(num_y): for x in range(num_x): - trials_not_gray = np.argwhere( LSN_binary[:, y, x] > 0 )[:, 0] - raw_mask = np.divide( LSN_binary[trials_not_gray, :, :].sum(axis=0), on_trials ) + trials_not_gray = np.argwhere(LSN_binary[:, y, x] > 0)[:, 0] + raw_mask = np.divide( + LSN_binary[trials_not_gray, :, :].sum(axis=0), on_trials + ) - center_y, center_x = np.unravel_index( raw_mask.argmax(), (num_y, num_x) ) + center_y, center_x = np.unravel_index( + raw_mask.argmax(), (num_y, num_x) + ) # include center pixel in mask raw_mask[center_y, center_x] = 0.0 @@ -500,11 +571,12 @@ def get_disc_masks(LSN_template, radius=3, on_luminance=ON_LUMINANCE, off_lumina # don't include far away pixels that just happen # to not have any trials in common with center pixel clean_mask = np.ones(np.shape(raw_mask)) - clean_mask[y_min:y_max, x_min:x_max] = raw_mask[y_min:y_max, x_min:x_max] + clean_mask[y_min:y_max, x_min:x_max] = raw_mask[ + y_min:y_max, x_min:x_max + ] masks[y, x, :, :] = clean_mask masks = np.where(masks > 0, 0.0, 1.0) return masks - diff --git a/allensdk/brain_observatory/receptive_field_analysis/eventdetection.py b/allensdk/brain_observatory/receptive_field_analysis/eventdetection.py index b7b25f50a..a63f9145f 100644 --- a/allensdk/brain_observatory/receptive_field_analysis/eventdetection.py +++ b/allensdk/brain_observatory/receptive_field_analysis/eventdetection.py @@ -33,46 +33,48 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # -from .utilities import smooth import numpy as np import scipy.stats as sps -def detect_events(data, cell_index, stimulus, debug_plots=False): - +from .utilities import smooth +def detect_events(data, cell_index, stimulus, debug_plots=False): stimulus_table = data.get_stimulus_table(stimulus) dff_trace = data.get_dff_traces()[1][cell_index, :] - - k_min = 0 k_max = 10 delta = 3 dff_trace = smooth(dff_trace, 5) - - var_dict = {} debug_dict = {} - for ii, fi in enumerate(stimulus_table['start'].values): - - if ii > 0 and stimulus_table.iloc[ii].start == stimulus_table.iloc[ii-1].end: + for ii, fi in enumerate(stimulus_table["start"].values): + if ( + ii > 0 + and stimulus_table.iloc[ii].start + == stimulus_table.iloc[ii - 1].end + ): offset = 1 else: offset = 0 if fi + k_min >= 0 and fi + k_max <= len(dff_trace): - trace = dff_trace[fi + k_min+1+offset:fi + k_max+1+offset] + trace = dff_trace[ + fi + k_min + 1 + offset : fi + k_max + 1 + offset + ] xx = (trace - trace[0])[delta] - (trace - trace[0])[0] - yy = max((trace - trace[0])[delta + 2] - (trace - trace[0])[0 + 2], - (trace - trace[0])[delta + 3] - (trace - trace[0])[0 + 3], - (trace - trace[0])[delta + 4] - (trace - trace[0])[0 + 4]) + yy = max( + (trace - trace[0])[delta + 2] - (trace - trace[0])[0 + 2], + (trace - trace[0])[delta + 3] - (trace - trace[0])[0 + 3], + (trace - trace[0])[delta + 4] - (trace - trace[0])[0 + 4], + ) var_dict[ii] = (trace[0], trace[-1], xx, yy) - debug_dict[fi + k_min+1+offset] = (ii, trace) + debug_dict[fi + k_min + 1 + offset] = (ii, trace) xx_list, yy_list = [], [] for _, _, xx, yy in var_dict.values(): @@ -82,17 +84,31 @@ def detect_events(data, cell_index, stimulus, debug_plots=False): mu_x = np.median(xx_list) mu_y = np.median(yy_list) - xx_centered = np.array(xx_list)-mu_x - yy_centered = np.array(yy_list)-mu_y + xx_centered = np.array(xx_list) - mu_x + yy_centered = np.array(yy_list) - mu_y std_factor = 1 - std_x = 1./std_factor*np.percentile(np.abs(xx_centered), [100*(1-2*(1-sps.norm.cdf(std_factor)))]) - std_y = 1./std_factor*np.percentile(np.abs(yy_centered), [100*(1-2*(1-sps.norm.cdf(std_factor)))]) + std_x = ( + 1.0 + / std_factor + * np.percentile( + np.abs(xx_centered), + [100 * (1 - 2 * (1 - sps.norm.cdf(std_factor)))], + ) + ) + std_y = ( + 1.0 + / std_factor + * np.percentile( + np.abs(yy_centered), + [100 * (1 - 2 * (1 - sps.norm.cdf(std_factor)))], + ) + ) curr_inds = [] allowed_sigma = 4 for ii, (xi, yi) in enumerate(zip(xx_centered, yy_centered)): - if np.sqrt(((xi)/std_x)**2+((yi)/std_y)**2) < allowed_sigma: + if np.sqrt(((xi) / std_x) ** 2 + ((yi) / std_y) ** 2) < allowed_sigma: curr_inds.append(True) else: curr_inds.append(False) @@ -104,15 +120,15 @@ def detect_events(data, cell_index, stimulus, debug_plots=False): Cov_Factor = np.linalg.cholesky(Cov) Cov_Factor_Inv = np.linalg.inv(Cov_Factor) - #=================================================================================================================== + # ========================================================================= - noise_threshold = max(allowed_sigma * std_x + mu_x, allowed_sigma * std_y + mu_y) + noise_threshold = max( + allowed_sigma * std_x + mu_x, allowed_sigma * std_y + mu_y + ) mu_array = np.array([mu_x, mu_y]) yes_set, no_set = set(), set() for ii, (t0, tf, xx, yy) in var_dict.items(): - - - xi_z, yi_z = Cov_Factor_Inv.dot((np.array([xx,yy]) - mu_array)) + xi_z, yi_z = Cov_Factor_Inv.dot((np.array([xx, yy]) - mu_array)) # Conditions in order: # 1) Outside noise blob @@ -120,38 +136,44 @@ def detect_events(data, cell_index, stimulus, debug_plots=False): # 3) Change evoked by this trial, not previous # 4) At end of trace, ended up outside of noise floor - if np.sqrt(xi_z**2 + yi_z**2) > 4 and yy > .05 and xx < yy and tf > noise_threshold/2: + if ( + np.sqrt(xi_z**2 + yi_z**2) > 4 + and yy > 0.05 + and xx < yy + and tf > noise_threshold / 2 + ): yes_set.add(ii) else: no_set.add(ii) - - assert len(var_dict) == len(stimulus_table) b = np.zeros(len(stimulus_table), dtype=bool) for yi in yes_set: b[yi] = True - if debug_plots == True: + if debug_plots: import matplotlib.pyplot as plt - fig, ax = plt.subplots(1,2) + + fig, ax = plt.subplots(1, 2) # ax[0].plot(dff_trace) for key, val in debug_dict.items(): ti, trace = val if ti in no_set: - ax[0].plot(np.arange(key, key+len(trace)), trace, 'b') + ax[0].plot(np.arange(key, key + len(trace)), trace, "b") elif ti in yes_set: - ax[0].plot(np.arange(key, key + len(trace)), trace, 'r', linewidth=2) + ax[0].plot( + np.arange(key, key + len(trace)), trace, "r", linewidth=2 + ) else: raise Exception for ii in yes_set: - ax[1].plot([var_dict[ii][2]], [var_dict[ii][3]], 'r.') + ax[1].plot([var_dict[ii][2]], [var_dict[ii][3]], "r.") for ii in no_set: - ax[1].plot([var_dict[ii][2]], [var_dict[ii][3]], 'b.') + ax[1].plot([var_dict[ii][2]], [var_dict[ii][3]], "b.") - print('number_of_events: %d' % b.sum()) + print("number_of_events: %d" % b.sum()) plt.show() return b diff --git a/allensdk/brain_observatory/receptive_field_analysis/postprocessing.py b/allensdk/brain_observatory/receptive_field_analysis/postprocessing.py index aacab6d25..de12f29d9 100644 --- a/allensdk/brain_observatory/receptive_field_analysis/postprocessing.py +++ b/allensdk/brain_observatory/receptive_field_analysis/postprocessing.py @@ -33,105 +33,157 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # -from .fit_parameters import get_gaussian_fit_single_channel, compute_distance, compute_overlap -from .chisquarerf import chi_square_binary, get_peak_significance, pvalue_to_NLL -from .utilities import upsample_image_to_degrees import collections + import numpy as np -import sys -def get_gaussian_fit(rf): +from .chisquarerf import ( + chi_square_binary, + get_peak_significance, + pvalue_to_NLL, +) +from .fit_parameters import ( + compute_distance, + compute_overlap, + get_gaussian_fit_single_channel, +) +from .utilities import upsample_image_to_degrees - fit_parameters_dict_combined = {'on':collections.defaultdict(list), 'off':collections.defaultdict(list)} - counter = {'on':0, 'off':0} - for on_off_key in ['on', 'off']: + +def get_gaussian_fit(rf): + fit_parameters_dict_combined = { + "on": collections.defaultdict(list), + "off": collections.defaultdict(list), + } + counter = {"on": 0, "off": 0} + for on_off_key in ["on", "off"]: fit_parameters_dict = fit_parameters_dict_combined[on_off_key] - for ci in range(rf[on_off_key]['fdr_mask']['attrs']['number_of_components']): - curr_component_mask = upsample_image_to_degrees(np.logical_not(rf[on_off_key]['fdr_mask']['data'][ci,:,:])) > .5 - rf_response = upsample_image_to_degrees(rf[on_off_key]['rts_convolution']['data'].copy()) + for ci in range( + rf[on_off_key]["fdr_mask"]["attrs"]["number_of_components"] + ): + curr_component_mask = ( + upsample_image_to_degrees( + np.logical_not( + rf[on_off_key]["fdr_mask"]["data"][ci, :, :] + ) + ) + > 0.5 + ) + rf_response = upsample_image_to_degrees( + rf[on_off_key]["rts_convolution"]["data"].copy() + ) rf_response[curr_component_mask] = 0 if rf_response.sum() > 0: - get_gaussian_fit_single_channel(rf_response, fit_parameters_dict) + get_gaussian_fit_single_channel( + rf_response, fit_parameters_dict + ) counter[on_off_key] += 1 - for ii_off in range(counter['on']): - fit_parameters_dict_combined['on']['distance'].append([None]*counter['off']) - fit_parameters_dict_combined['on']['overlap'].append([None] * counter['off']) - - for ii_off in range(counter['off']): - fit_parameters_dict_combined['off']['distance'].append([None] * counter['on']) - fit_parameters_dict_combined['off']['overlap'].append([None]*counter['on']) - - - for ii_on in range(counter['on']): - for ii_off in range(counter['off']): - center_on = fit_parameters_dict_combined['on']['center_x'][ii_on], fit_parameters_dict_combined['on']['center_y'][ii_on] - center_off = fit_parameters_dict_combined['off']['center_x'][ii_off], fit_parameters_dict_combined['off']['center_y'][ii_off] + for ii_off in range(counter["on"]): + fit_parameters_dict_combined["on"]["distance"].append( + [None] * counter["off"] + ) + fit_parameters_dict_combined["on"]["overlap"].append( + [None] * counter["off"] + ) + + for ii_off in range(counter["off"]): + fit_parameters_dict_combined["off"]["distance"].append( + [None] * counter["on"] + ) + fit_parameters_dict_combined["off"]["overlap"].append( + [None] * counter["on"] + ) + + for ii_on in range(counter["on"]): + for ii_off in range(counter["off"]): + center_on = ( + fit_parameters_dict_combined["on"]["center_x"][ii_on], + fit_parameters_dict_combined["on"]["center_y"][ii_on], + ) + center_off = ( + fit_parameters_dict_combined["off"]["center_x"][ii_off], + fit_parameters_dict_combined["off"]["center_y"][ii_off], + ) curr_distance = compute_distance(center_on, center_off) - fit_parameters_dict_combined['on']['distance'][ii_on][ii_off] = curr_distance - fit_parameters_dict_combined['off']['distance'][ii_off][ii_on] = curr_distance - - data_on = fit_parameters_dict_combined['on']['data'][ii_on] - data_off = fit_parameters_dict_combined['off']['data'][ii_off] + fit_parameters_dict_combined["on"]["distance"][ii_on][ + ii_off + ] = curr_distance + fit_parameters_dict_combined["off"]["distance"][ii_off][ + ii_on + ] = curr_distance + + data_on = fit_parameters_dict_combined["on"]["data"][ii_on] + data_off = fit_parameters_dict_combined["off"]["data"][ii_off] curr_overlap = compute_overlap(data_on, data_off) - fit_parameters_dict_combined['on']['overlap'][ii_on][ii_off] = curr_overlap - fit_parameters_dict_combined['off']['overlap'][ii_off][ii_on] = curr_overlap + fit_parameters_dict_combined["on"]["overlap"][ii_on][ + ii_off + ] = curr_overlap + fit_parameters_dict_combined["off"]["overlap"][ii_off][ + ii_on + ] = curr_overlap return fit_parameters_dict_combined, counter -def run_postprocessing(data, rf): - stimulus = rf['attrs']['stimulus'] +def run_postprocessing(data, rf): + stimulus = rf["attrs"]["stimulus"] # Gaussian fit postprocessing: fit_parameters_dict_combined, counter = get_gaussian_fit(rf) - for on_off_key in ['on', 'off']: - + for on_off_key in ["on", "off"]: if counter[on_off_key] > 0: - - rf[on_off_key]['gaussian_fit'] = {} - rf[on_off_key]['gaussian_fit']['attrs'] = {} + rf[on_off_key]["gaussian_fit"] = {} + rf[on_off_key]["gaussian_fit"]["attrs"] = {} fit_parameters_dict = fit_parameters_dict_combined[on_off_key] for key, val in fit_parameters_dict.items(): - - if key == 'data': - rf[on_off_key]['gaussian_fit']['data'] = np.array(val) + if key == "data": + rf[on_off_key]["gaussian_fit"]["data"] = np.array(val) else: - rf[on_off_key]['gaussian_fit']['attrs'][key] = np.array(val) + rf[on_off_key]["gaussian_fit"]["attrs"][key] = np.array( + val + ) # Chi squared test statistic postprocessing: - cell_index = rf['attrs']['cell_index'] + # cell_index = rf["attrs"]["cell_index"] locally_sparse_noise_template = data.get_stimulus_template(stimulus) - event_array = np.zeros((rf['event_vector']['data'].shape[0], 1), dtype=bool) - event_array[:,0] = rf['event_vector']['data'] + event_array = np.zeros( + (rf["event_vector"]["data"].shape[0], 1), dtype=bool + ) + event_array[:, 0] = rf["event_vector"]["data"] - chi_squared_grid = chi_square_binary(event_array, locally_sparse_noise_template) - alpha = rf['on']['fdr_mask']['attrs']['alpha'] - assert rf['off']['fdr_mask']['attrs']['alpha'] == alpha + chi_squared_grid = chi_square_binary( + event_array, locally_sparse_noise_template + ) + alpha = rf["on"]["fdr_mask"]["attrs"]["alpha"] + assert rf["off"]["fdr_mask"]["attrs"]["alpha"] == alpha chi_square_grid_NLL = pvalue_to_NLL(chi_squared_grid) - peak_significance = get_peak_significance(chi_square_grid_NLL, locally_sparse_noise_template, alpha=alpha) + peak_significance = get_peak_significance( + chi_square_grid_NLL, locally_sparse_noise_template, alpha=alpha + ) significant = peak_significance[0][0] min_p = peak_significance[1][0] pvalues_chi_square = peak_significance[2][0] best_exclusion_region_mask = peak_significance[3][0] chi_squared_grid_dict = { - 'best_exclusion_region_mask':{'data':best_exclusion_region_mask}, - 'attrs':{'significant':significant, 'alpha': alpha, 'min_p':min_p}, - 'pvalues':{'data':pvalues_chi_square} - } + "best_exclusion_region_mask": {"data": best_exclusion_region_mask}, + "attrs": {"significant": significant, "alpha": alpha, "min_p": min_p}, + "pvalues": {"data": pvalues_chi_square}, + } - rf['chi_squared_analysis'] = chi_squared_grid_dict + rf["chi_squared_analysis"] = chi_squared_grid_dict return rf + if __name__ == "__main__": # csid = 517472416 # triple! - csid = 517526760 # two ON + csid = 517526760 # two ON # csid = 539917553 # csid = 540988186 diff --git a/allensdk/brain_observatory/receptive_field_analysis/utilities.py b/allensdk/brain_observatory/receptive_field_analysis/utilities.py index 187af810d..1a2a6ef13 100644 --- a/allensdk/brain_observatory/receptive_field_analysis/utilities.py +++ b/allensdk/brain_observatory/receptive_field_analysis/utilities.py @@ -33,137 +33,174 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # -from scipy.ndimage.filters import gaussian_filter +import warnings + import numpy as np import scipy.interpolate as spinterp -from .tools import dict_generator from allensdk.api.warehouse_cache.cache import memoize -import os -import warnings +from scipy.ndimage.filters import gaussian_filter from skimage.measure import block_reduce -def upsample_image_to_degrees(img): +from .tools import dict_generator + - upsample = 74.4/img.shape[0] +def upsample_image_to_degrees(img): + upsample = 74.4 / img.shape[0] x = np.arange(img.shape[0]) y = np.arange(img.shape[1]) - g = spinterp.interp2d(y, x, img, kind='linear') - ZZ_on = g(np.arange(0, img.shape[1], 1. / upsample), np.arange(0, img.shape[0], 1. / upsample)) + g = spinterp.interp2d(y, x, img, kind="linear") + ZZ_on = g( + np.arange(0, img.shape[1], 1.0 / upsample), + np.arange(0, img.shape[0], 1.0 / upsample), + ) return ZZ_on + def convolve(img, sigma=4): - ''' + """ 2D Gaussian convolution - ''' + """ if img.sum() == 0: return img img_pad = np.zeros((3 * img.shape[0], 3 * img.shape[1])) - img_pad[img.shape[0]:2 * img.shape[0], img.shape[1]:2 * img.shape[1]] = img + img_pad[ + img.shape[0] : 2 * img.shape[0], img.shape[1] : 2 * img.shape[1] + ] = img x = np.arange(3 * img.shape[0]) y = np.arange(3 * img.shape[1]) - g = spinterp.interp2d(y, x, img_pad, kind='linear') + g = spinterp.interp2d(y, x, img_pad, kind="linear") if img.shape[0] == 16: upsample = 4 - offset = -(1 - .625) + offset = -(1 - 0.625) elif img.shape[0] == 8: upsample = 8 - offset = -(1 - .5625) + offset = -(1 - 0.5625) else: raise NotImplementedError - ZZ_on = g(offset + np.arange(0, img.shape[1] * 3, 1. / upsample), offset + np.arange(0, img.shape[0] * 3, 1. / upsample)) - ZZ_on_f = gaussian_filter(ZZ_on, float(sigma), mode='constant') + ZZ_on = g( + offset + np.arange(0, img.shape[1] * 3, 1.0 / upsample), + offset + np.arange(0, img.shape[0] * 3, 1.0 / upsample), + ) + ZZ_on_f = gaussian_filter(ZZ_on, float(sigma), mode="constant") z_on_new = block_reduce(ZZ_on_f, (upsample, upsample)) z_on_new = z_on_new / z_on_new.sum() * img.sum() - z_on_new = z_on_new[img.shape[0]:2 * img.shape[0], img.shape[1]:2 * img.shape[1]] + z_on_new = z_on_new[ + img.shape[0] : 2 * img.shape[0], img.shape[1] : 2 * img.shape[1] + ] return z_on_new - @memoize def get_A(data, stimulus): - stimulus_table = data.get_stimulus_table(stimulus) - stimulus_template = data.get_stimulus_template(stimulus)[stimulus_table['frame'].values, :,:] + stimulus_template = data.get_stimulus_template(stimulus)[ + stimulus_table["frame"].values, :, : + ] - number_of_pixels = stimulus_template.shape[1]*stimulus_template.shape[2] + number_of_pixels = stimulus_template.shape[1] * stimulus_template.shape[2] - A = np.zeros((2*number_of_pixels, stimulus_template.shape[0])) + A = np.zeros((2 * number_of_pixels, stimulus_template.shape[0])) for fi in range(stimulus_template.shape[0]): - A[:number_of_pixels, fi] = (stimulus_template[fi,:,:].flatten() > 127).astype(float) - A[number_of_pixels:, fi] = (stimulus_template[fi, :, :].flatten() < 127).astype(float) + A[:number_of_pixels, fi] = ( + stimulus_template[fi, :, :].flatten() > 127 + ).astype(float) + A[number_of_pixels:, fi] = ( + stimulus_template[fi, :, :].flatten() < 127 + ).astype(float) return A + @memoize def get_A_blur(data, stimulus): - stimulus_table = data.get_stimulus_table(stimulus) - stimulus_template = data.get_stimulus_template(stimulus)[stimulus_table['frame'].values, :, :] + stimulus_template = data.get_stimulus_template(stimulus)[ + stimulus_table["frame"].values, :, : + ] A = get_A(data, stimulus).copy() - number_of_pixels = A.shape[0] // 2 for fi in range(A.shape[1]): - A[:number_of_pixels,fi] = convolve(A[:number_of_pixels, fi].reshape(stimulus_template.shape[1], stimulus_template.shape[2])).flatten() - A[number_of_pixels:,fi] = convolve(A[number_of_pixels:, fi].reshape(stimulus_template.shape[1], stimulus_template.shape[2])).flatten() - + A[:number_of_pixels, fi] = convolve( + A[:number_of_pixels, fi].reshape( + stimulus_template.shape[1], stimulus_template.shape[2] + ) + ).flatten() + A[number_of_pixels:, fi] = convolve( + A[number_of_pixels:, fi].reshape( + stimulus_template.shape[1], stimulus_template.shape[2] + ) + ).flatten() return A -def get_shuffle_matrix(data, event_vector, A, number_of_shuffles=5000, response_detection_error_std_dev=.1): +def get_shuffle_matrix( + data, + event_vector, + A, + number_of_shuffles=5000, + response_detection_error_std_dev=0.1, +): number_of_events = event_vector.sum() number_of_pixels = A.shape[0] // 2 - shuffle_data = np.zeros((2*number_of_pixels, number_of_shuffles)) + shuffle_data = np.zeros((2 * number_of_pixels, number_of_shuffles)) evr = range(len(event_vector)) for ii in range(number_of_shuffles): - - size = number_of_events + int(np.round(response_detection_error_std_dev*number_of_events*np.random.randn())) + size = number_of_events + int( + np.round( + response_detection_error_std_dev + * number_of_events + * np.random.randn() + ) + ) shuffled_event_inds = np.random.choice(evr, size=size, replace=False) b_tmp = np.zeros(len(event_vector), dtype=bool) b_tmp[shuffled_event_inds] = True - shuffle_data[:, ii] = A[:,b_tmp].sum(axis=1)/float(size) + shuffle_data[:, ii] = A[:, b_tmp].sum(axis=1) / float(size) return shuffle_data -def get_sparse_noise_epoch_mask_list(st, number_of_acquisition_frames, threshold=7): - delta = (st.start.values[1:] - st.end.values[:-1]) +def get_sparse_noise_epoch_mask_list( + st, number_of_acquisition_frames, threshold=7 +): + delta = st.start.values[1:] - st.end.values[:-1] cut_inds = np.where(delta > threshold)[0] + 1 epoch_mask_list = [] if len(cut_inds) > 2: - warnings.warn('more than 2 epochs cut') - print(' %d %d' % (len(delta), cut_inds)) - - for ii in range(len(cut_inds)+1): + warnings.warn("more than 2 epochs cut") + print(" %d %d" % (len(delta), cut_inds)) + for ii in range(len(cut_inds) + 1): if ii == 0: first_ind = st.iloc[0].start else: - first_ind = st.iloc[cut_inds[ii-1]].start + first_ind = st.iloc[cut_inds[ii - 1]].start if ii == len(cut_inds): last_ind_inclusive = st.iloc[-1].end else: - last_ind_inclusive = st.iloc[cut_inds[ii]-1].end + last_ind_inclusive = st.iloc[cut_inds[ii] - 1].end - epoch_mask_list.append((first_ind,last_ind_inclusive)) + epoch_mask_list.append((first_ind, last_ind_inclusive)) return epoch_mask_list -def smooth(x,window_len=11,window='hanning', mode='valid'): + +def smooth(x, window_len=11, window="hanning", mode="valid"): """smooth the data using a window with requested size. This method is based on the convolution of a scaled window with the signal. @@ -173,9 +210,11 @@ def smooth(x,window_len=11,window='hanning', mode='valid'): input: x: the input signal - window_len: the dimension of the smoothing window; should be an odd integer - window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman' - flat window will produce a moving average smoothing. + window_len: the dimension of the smoothing window; should be an odd + integer + window: the type of window from 'flat', 'hanning', 'hamming', + 'bartlett', 'blackman' flat window will produce a moving average + smoothing. output: the smoothed signal @@ -188,11 +227,13 @@ def smooth(x,window_len=11,window='hanning', mode='valid'): see also: - numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve - scipy.signal.lfilter + numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, + numpy.convolve scipy.signal.lfilter - TODO: the window parameter could be the window itself if an array instead of a string - NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y. + TODO: the window parameter could be the window itself if an array instead + of a string + NOTE: length(output) != length(input), to correct this: + return y[(window_len/2-1):-(window_len/2)] instead of just y. """ if x.ndim != 1: @@ -201,37 +242,38 @@ def smooth(x,window_len=11,window='hanning', mode='valid'): if x.size < window_len: raise ValueError("Input vector needs to be bigger than window size.") - - if window_len<3: + if window_len < 3: return x + if window not in ["flat", "hanning", "hamming", "bartlett", "blackman"]: + raise ValueError( + "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', " + "'blackman'" + ) - if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: - raise ValueError("Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'") - - - s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]] - #print(len(s)) - if window == 'flat': #moving average - w=np.ones(window_len,'d') + s = np.r_[x[window_len - 1 : 0 : -1], x, x[-1:-window_len:-1]] + # print(len(s)) + if window == "flat": # moving average + w = np.ones(window_len, "d") else: - w=eval('np.'+window+'(window_len)') + w = eval("np." + window + "(window_len)") - y=np.convolve(w/w.sum(),s,mode=mode) + y = np.convolve(w / w.sum(), s, mode=mode) return y def get_components(receptive_field_data): - s1, s2 = receptive_field_data.shape - candidate_pixel_list = np.where(receptive_field_data.flatten()==True)[0] - pixel_coord_dict = dict((px, (int(px/s2), (px - s2 * int(px/s2)), px% (s1 * s2) == px)) for px in candidate_pixel_list) + candidate_pixel_list = np.where(receptive_field_data.flatten())[0] + pixel_coord_dict = dict( + (px, (int(px / s2), (px - s2 * int(px / s2)), px % (s1 * s2) == px)) + for px in candidate_pixel_list + ) component_list = [] for curr_pixel in candidate_pixel_list: - curr_x, curr_y, curr_frame = pixel_coord_dict[curr_pixel] component_list.append([curr_pixel]) @@ -239,14 +281,13 @@ def get_components(receptive_field_data): for ii, curr_component in enumerate(component_list): dist_to_component_dict[ii] = np.inf for other_pixel in curr_component: - other_x, other_y, other_frame = pixel_coord_dict[other_pixel] if other_frame == curr_frame: x_dist = np.abs(curr_x - other_x) y_dist = np.abs(curr_y - other_y) curr_dist = max(x_dist, y_dist) - if curr_dist < dist_to_component_dict[ii]: + if curr_dist < dist_to_component_dict[ii]: dist_to_component_dict[ii] = curr_dist # Merge all components with a distance leq 1 to current pixel @@ -261,33 +302,43 @@ def get_components(receptive_field_data): new_component_list.append(tmp) component_list = new_component_list - if len(component_list) == 0: - return np.zeros((1, receptive_field_data.shape[0], receptive_field_data.shape[1])), len(component_list) + return np.zeros( + (1, receptive_field_data.shape[0], receptive_field_data.shape[1]) + ), len(component_list) elif len(component_list) == 1: - return_array = np.zeros((1,receptive_field_data.shape[0], receptive_field_data.shape[1])) + return_array = np.zeros( + (1, receptive_field_data.shape[0], receptive_field_data.shape[1]) + ) else: - return_array = np.zeros((len(component_list), receptive_field_data.shape[0], receptive_field_data.shape[1])) + return_array = np.zeros( + ( + len(component_list), + receptive_field_data.shape[0], + receptive_field_data.shape[1], + ) + ) for ii, component in enumerate(component_list): - curr_component_mask = np.zeros_like(receptive_field_data, dtype=bool).flatten() + curr_component_mask = np.zeros_like( + receptive_field_data, dtype=bool + ).flatten() curr_component_mask[component] = True - return_array[ii,:,:] = curr_component_mask.reshape(receptive_field_data.shape) - - + return_array[ii, :, :] = curr_component_mask.reshape( + receptive_field_data.shape + ) return return_array, len(component_list) -def get_attribute_dict(rf): +def get_attribute_dict(rf): attribute_dict = {} for x in dict_generator(rf): - if x[-3] == 'attrs': + if x[-3] == "attrs": if len(x[:-3]) == 0: key = x[-2] else: - key = '/'.join(['/'.join(x[:-3]), x[-2]]) + key = "/".join(["/".join(x[:-3]), x[-2]]) attribute_dict[key] = x[-1] return attribute_dict - diff --git a/allensdk/brain_observatory/stimulus_info.py b/allensdk/brain_observatory/stimulus_info.py index e37f2a0d9..c257028b6 100755 --- a/allensdk/brain_observatory/stimulus_info.py +++ b/allensdk/brain_observatory/stimulus_info.py @@ -33,84 +33,120 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # -import six +import itertools + import numpy as np import scipy.ndimage.interpolation as spndi -from PIL import Image +import six from allensdk.api.warehouse_cache.cache import memoize -import itertools +from PIL import Image # some handles for stimulus types -DRIFTING_GRATINGS = 'drifting_gratings' -DRIFTING_GRATINGS_SHORT = 'dg' -DRIFTING_GRATINGS_COLOR = '#a6cee3' +DRIFTING_GRATINGS = "drifting_gratings" +DRIFTING_GRATINGS_SHORT = "dg" +DRIFTING_GRATINGS_COLOR = "#a6cee3" -STATIC_GRATINGS = 'static_gratings' -STATIC_GRATINGS_SHORT = 'sg' -STATIC_GRATINGS_COLOR = '#1f78b4' +STATIC_GRATINGS = "static_gratings" +STATIC_GRATINGS_SHORT = "sg" +STATIC_GRATINGS_COLOR = "#1f78b4" -NATURAL_MOVIE_ONE = 'natural_movie_one' -NATURAL_MOVIE_ONE_SHORT = 'nm1' -NATURAL_MOVIE_ONE_COLOR = '#b2df8a' +NATURAL_MOVIE_ONE = "natural_movie_one" +NATURAL_MOVIE_ONE_SHORT = "nm1" +NATURAL_MOVIE_ONE_COLOR = "#b2df8a" -NATURAL_MOVIE_TWO = 'natural_movie_two' -NATURAL_MOVIE_TWO_SHORT = 'nm2' -NATURAL_MOVIE_TWO_COLOR = '#33a02c' +NATURAL_MOVIE_TWO = "natural_movie_two" +NATURAL_MOVIE_TWO_SHORT = "nm2" +NATURAL_MOVIE_TWO_COLOR = "#33a02c" -NATURAL_MOVIE_THREE = 'natural_movie_three' -NATURAL_MOVIE_THREE_SHORT = 'nm3' -NATURAL_MOVIE_THREE_COLOR = '#fb9a99' +NATURAL_MOVIE_THREE = "natural_movie_three" +NATURAL_MOVIE_THREE_SHORT = "nm3" +NATURAL_MOVIE_THREE_COLOR = "#fb9a99" -NATURAL_SCENES = 'natural_scenes' -NATURAL_SCENES_SHORT = 'ns' -NATURAL_SCENES_COLOR = '#e31a1c' +NATURAL_SCENES = "natural_scenes" +NATURAL_SCENES_SHORT = "ns" +NATURAL_SCENES_COLOR = "#e31a1c" -# note that this stimulus is equivalent to LOCALLY_SPARSE_NOISE_4DEG in session C2 files -LOCALLY_SPARSE_NOISE = 'locally_sparse_noise' -LOCALLY_SPARSE_NOISE_SHORT = 'lsn' -LOCALLY_SPARSE_NOISE_COLOR = '#fdbf6f' +# note that this stimulus is equivalent to LOCALLY_SPARSE_NOISE_4DEG in session +# C2 files +LOCALLY_SPARSE_NOISE = "locally_sparse_noise" +LOCALLY_SPARSE_NOISE_SHORT = "lsn" +LOCALLY_SPARSE_NOISE_COLOR = "#fdbf6f" -LOCALLY_SPARSE_NOISE_4DEG = 'locally_sparse_noise_4deg' -LOCALLY_SPARSE_NOISE_4DEG_SHORT = 'lsn4' -LOCALLY_SPARSE_NOISE_4DEG_COLOR = '#fdbf6f' +LOCALLY_SPARSE_NOISE_4DEG = "locally_sparse_noise_4deg" +LOCALLY_SPARSE_NOISE_4DEG_SHORT = "lsn4" +LOCALLY_SPARSE_NOISE_4DEG_COLOR = "#fdbf6f" -LOCALLY_SPARSE_NOISE_8DEG = 'locally_sparse_noise_8deg' -LOCALLY_SPARSE_NOISE_8DEG_SHORT = 'lsn8' -LOCALLY_SPARSE_NOISE_8DEG_COLOR = '#ff7f00' +LOCALLY_SPARSE_NOISE_8DEG = "locally_sparse_noise_8deg" +LOCALLY_SPARSE_NOISE_8DEG_SHORT = "lsn8" +LOCALLY_SPARSE_NOISE_8DEG_COLOR = "#ff7f00" -SPONTANEOUS_ACTIVITY = 'spontaneous' -SPONTANEOUS_ACTIVITY_SHORT = 'sp' -SPONTANEOUS_ACTIVITY_COLOR = '#cab2d6' +SPONTANEOUS_ACTIVITY = "spontaneous" +SPONTANEOUS_ACTIVITY_SHORT = "sp" +SPONTANEOUS_ACTIVITY_COLOR = "#cab2d6" # handles for stimulus names -THREE_SESSION_A = 'three_session_A' -THREE_SESSION_B = 'three_session_B' -THREE_SESSION_C = 'three_session_C' -THREE_SESSION_C2 = 'three_session_C2' - -SESSION_LIST = [THREE_SESSION_A, THREE_SESSION_B, THREE_SESSION_C, THREE_SESSION_C2] +THREE_SESSION_A = "three_session_A" +THREE_SESSION_B = "three_session_B" +THREE_SESSION_C = "three_session_C" +THREE_SESSION_C2 = "three_session_C2" + +SESSION_LIST = [ + THREE_SESSION_A, + THREE_SESSION_B, + THREE_SESSION_C, + THREE_SESSION_C2, +] SESSION_STIMULUS_MAP = { - THREE_SESSION_A: [DRIFTING_GRATINGS, NATURAL_MOVIE_ONE, NATURAL_MOVIE_THREE, SPONTANEOUS_ACTIVITY], - THREE_SESSION_B: [STATIC_GRATINGS, NATURAL_SCENES, NATURAL_MOVIE_ONE, SPONTANEOUS_ACTIVITY], - THREE_SESSION_C: [LOCALLY_SPARSE_NOISE, NATURAL_MOVIE_ONE, NATURAL_MOVIE_TWO, SPONTANEOUS_ACTIVITY], - THREE_SESSION_C2: [LOCALLY_SPARSE_NOISE_4DEG, LOCALLY_SPARSE_NOISE_8DEG, NATURAL_MOVIE_ONE, NATURAL_MOVIE_TWO, SPONTANEOUS_ACTIVITY] + THREE_SESSION_A: [ + DRIFTING_GRATINGS, + NATURAL_MOVIE_ONE, + NATURAL_MOVIE_THREE, + SPONTANEOUS_ACTIVITY, + ], + THREE_SESSION_B: [ + STATIC_GRATINGS, + NATURAL_SCENES, + NATURAL_MOVIE_ONE, + SPONTANEOUS_ACTIVITY, + ], + THREE_SESSION_C: [ + LOCALLY_SPARSE_NOISE, + NATURAL_MOVIE_ONE, + NATURAL_MOVIE_TWO, + SPONTANEOUS_ACTIVITY, + ], + THREE_SESSION_C2: [ + LOCALLY_SPARSE_NOISE_4DEG, + LOCALLY_SPARSE_NOISE_8DEG, + NATURAL_MOVIE_ONE, + NATURAL_MOVIE_TWO, + SPONTANEOUS_ACTIVITY, + ], } -LOCALLY_SPARSE_NOISE_STIMULUS_TYPES = [LOCALLY_SPARSE_NOISE, LOCALLY_SPARSE_NOISE_4DEG, LOCALLY_SPARSE_NOISE_8DEG] -NATURAL_MOVIE_STIMULUS_TYPES = [NATURAL_MOVIE_ONE, NATURAL_MOVIE_TWO, NATURAL_MOVIE_THREE] +LOCALLY_SPARSE_NOISE_STIMULUS_TYPES = [ + LOCALLY_SPARSE_NOISE, + LOCALLY_SPARSE_NOISE_4DEG, + LOCALLY_SPARSE_NOISE_8DEG, +] +NATURAL_MOVIE_STIMULUS_TYPES = [ + NATURAL_MOVIE_ONE, + NATURAL_MOVIE_TWO, + NATURAL_MOVIE_THREE, +] LOCALLY_SPARSE_NOISE_DIMENSIONS = { - LOCALLY_SPARSE_NOISE: [ 16, 28 ], - LOCALLY_SPARSE_NOISE_4DEG: [ 16, 28 ], - LOCALLY_SPARSE_NOISE_8DEG: [ 8, 14 ], - } + LOCALLY_SPARSE_NOISE: [16, 28], + LOCALLY_SPARSE_NOISE_4DEG: [16, 28], + LOCALLY_SPARSE_NOISE_8DEG: [8, 14], +} LOCALLY_SPARSE_NOISE_PIXELS = { LOCALLY_SPARSE_NOISE: 45, LOCALLY_SPARSE_NOISE_4DEG: 45, LOCALLY_SPARSE_NOISE_8DEG: 90, - } +} NATURAL_SCENES_PIXELS = (918, 1174) NATURAL_MOVIE_PIXELS = (1080, 1920) @@ -126,13 +162,14 @@ LOCALLY_SPARSE_NOISE_PIXEL_SIZE = { LOCALLY_SPARSE_NOISE: 4.65, LOCALLY_SPARSE_NOISE_4DEG: 4.65, - LOCALLY_SPARSE_NOISE_8DEG: 9.3 + LOCALLY_SPARSE_NOISE_8DEG: 9.3, } RADIANS_TO_DEGREES = 57.2958 + def sessions_with_stimulus(stimulus): - """ Return the names of the sessions that contain a given stimulus. """ + """Return the names of the sessions that contain a given stimulus.""" sessions = set() for session, session_stimuli in six.iteritems(SESSION_STIMULUS_MAP): @@ -143,12 +180,17 @@ def sessions_with_stimulus(stimulus): def stimuli_in_session(session, allow_unknown=True): - """ Return a list what stimuli are available in a given session. + """Return a list what stimuli are available in a given session. Parameters ---------- session: string - Must be one of: [stimulus_info.THREE_SESSION_A, stimulus_info.THREE_SESSION_B, stimulus_info.THREE_SESSION_C, stimulus_info.THREE_SESSION_C2] + Must be one of: [ + stimulus_info.THREE_SESSION_A, + stimulus_info.THREE_SESSION_B, + stimulus_info.THREE_SESSION_C, + stimulus_info.THREE_SESSION_C2 + ] """ try: return SESSION_STIMULUS_MAP[session] @@ -156,37 +198,39 @@ def stimuli_in_session(session, allow_unknown=True): if allow_unknown: return [] else: - raise + raise e def all_stimuli(): - """ Return a list of all stimuli in the data set """ - return set([v for k, vl in six.iteritems(SESSION_STIMULUS_MAP) for v in vl]) + """Return a list of all stimuli in the data set""" + return set( + [v for k, vl in six.iteritems(SESSION_STIMULUS_MAP) for v in vl] + ) -class BinaryIntervalSearchTree(object): +class BinaryIntervalSearchTree(object): @staticmethod def from_df(input_df): - search_list = input_df.to_dict('records') - - + search_list = input_df.to_dict("records") new_list = [] for x in search_list: - if x['start'] == x['end']: - new_list.append((x['start'], x['end'], x)) + if x["start"] == x["end"]: + new_list.append((x["start"], x["end"], x)) else: - # -.01 prevents endpoint-overlapping intervals; assigns ties to intervals that start at requested index - new_list.append((x['start'], x['end'] - .01, x)) + # -.01 prevents endpoint-overlapping intervals; assigns ties to + # intervals that start at requested index + new_list.append((x["start"], x["end"] - 0.01, x)) return BinaryIntervalSearchTree(new_list) - def __init__(self, search_list): - """Create a binary tree to search for a point within a list of intervals. Assumes that the intervals are - non-overlapping. If two intervals share an endpoint, the left-side wins the tie. + """Create a binary tree to search for a point within a list of + intervals. Assumes that the intervals are non-overlapping. If two + intervals share an endpoint, the left-side wins the tie. - :param search_list: list of interval tuples; in the tuple, first element is interval start, then interval - end (inclusive), then the return value for the lookup + :param search_list: list of interval tuples; in the tuple, first + element is interval start, then interval end (inclusive), then the + return value for the lookup Example: bist = BinaryIntervalSearchTree([(0,.5,'A'), (1,2,'B')]) @@ -194,13 +238,13 @@ def __init__(self, search_list): """ # Double-check that the list is sorted - search_list = sorted(search_list, key=lambda x:x[0]) + search_list = sorted(search_list, key=lambda x: x[0]) - # Check that the intervals are non-overlapping (except potentially at the end point) + # Check that the intervals are non-overlapping (except potentially at + # the end point) for x, y in zip(search_list[:-1], search_list[1:]): assert x[1] <= y[0] - self.data = {} self.add(search_list) @@ -211,15 +255,17 @@ def add(self, input_list, tmp=None): if len(input_list) == 1: self.data[tuple(tmp)] = input_list[0] else: - self.add(input_list[:int(len(input_list)/2)], tmp=tmp+[0]) - self.add(input_list[int(len(input_list)/2):], tmp=tmp+[1]) - self.data[tuple(tmp)] = input_list[int(len(input_list)/2)-1] + self.add(input_list[: int(len(input_list) / 2)], tmp=tmp + [0]) + self.add(input_list[int(len(input_list) / 2) :], tmp=tmp + [1]) + self.data[tuple(tmp)] = input_list[int(len(input_list) / 2) - 1] def search(self, fi, tmp=None): if tmp is None: tmp = [] - if (self.data[tuple(tmp)][0] <= fi) and (fi <= self.data[tuple(tmp)][1]): + if (self.data[tuple(tmp)][0] <= fi) and ( + fi <= self.data[tuple(tmp)][1] + ): return_val = self.data[tuple(tmp)] elif fi < self.data[tuple(tmp)][1]: return_val = self.search(fi, tmp=tmp + [0]) @@ -229,67 +275,73 @@ def search(self, fi, tmp=None): assert (return_val[0] <= fi) and (fi <= return_val[1]) return return_val -class StimulusSearch(object): +class StimulusSearch(object): def __init__(self, nwb_dataset): - self.nwb_data = nwb_dataset self.epoch_df = nwb_dataset.get_stimulus_epoch_table() - self.master_df = nwb_dataset.get_stimulus_table('master') + self.master_df = nwb_dataset.get_stimulus_table("master") self.epoch_bst = BinaryIntervalSearchTree.from_df(self.epoch_df) self.master_bst = BinaryIntervalSearchTree.from_df(self.master_df) @memoize def search(self, fi): - try: - # Look in fine-grain tree: search_result = self.master_bst.search(fi) return search_result except KeyError: - # Current frame not found in a fine-grain interval; # see if it is unregistered to a coarse-grain epoch: try: - # THis will thow KeyError if not in coarse-grain epoch self.epoch_bst.search(fi) - # Frame is in a coarse-grain epoch, but not a fine grain interval; - # look backwards to find most recent find nearest matching interval - if fi < self.epoch_df.iloc[0]['start']: + # Frame is in a coarse-grain epoch, but not a fine grain + # interval; look backwards to find most recent find nearest + # matching interval + if fi < self.epoch_df.iloc[0]["start"]: return None else: - return self.search(fi-1) + return self.search(fi - 1) except KeyError: - # Frame is unregistered at the coarse level; return None return None + def rotate(X, Y, theta): x = np.array([X, Y]) - M = np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta), np.cos(theta)]]) - if len(x.shape) in [1,2]: + M = np.array( + [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] + ) + if len(x.shape) in [1, 2]: assert x.shape[0] == 2 return M.dot(x) elif len(x.shape) == 3: M2 = M[:, :, np.newaxis, np.newaxis] x2 = x[np.newaxis, :, :] - return (M2*x2).sum(axis=1) + return (M2 * x2).sum(axis=1) else: raise NotImplementedError -def get_spatial_grating(height=None, aspect_ratio=None, ori=None, pix_per_cycle=None, phase=None, p2p_amp=2, baseline=0): +def get_spatial_grating( + height=None, + aspect_ratio=None, + ori=None, + pix_per_cycle=None, + phase=None, + p2p_amp=2, + baseline=0, +): aspect_ratio = float(aspect_ratio) _height_prime = 100 - sf = 1./(float(pix_per_cycle)/(height/float(_height_prime))) + sf = 1.0 / (float(pix_per_cycle) / (height / float(_height_prime))) # Final height set by zoom below: - y, x = (_height_prime,_height_prime*aspect_ratio) + y, x = (_height_prime, _height_prime * aspect_ratio) theta = ori * np.pi / 180.0 # convert to radians @@ -302,19 +354,25 @@ def get_spatial_grating(height=None, aspect_ratio=None, ori=None, pix_per_cycle= img = np.cos(2.0 * np.pi * Xp * sf + ph) - return (p2p_amp/2.)*spndi.zoom(img, height/float(_height_prime)) + baseline + return (p2p_amp / 2.0) * spndi.zoom( + img, height / float(_height_prime) + ) + baseline # def grating_to_screen(self, phase, spatial_frequency, orientation, **kwargs): -def get_spatio_temporal_grating(t, temporal_frequency=None, **kwargs): - kwargs['phase'] = kwargs.pop('phase', 0) + (float(t)*temporal_frequency)%1 +def get_spatio_temporal_grating(t, temporal_frequency=None, **kwargs): + kwargs["phase"] = ( + kwargs.pop("phase", 0) + (float(t) * temporal_frequency) % 1 + ) return get_spatial_grating(**kwargs) -def map_template_coordinate_to_monitor_coordinate(template_coord, monitor_shape, template_shape): +def map_template_coordinate_to_monitor_coordinate( + template_coord, monitor_shape, template_shape +): rx, cx = template_coord n_pixels_r, n_pixels_c = monitor_shape tr, tc = template_shape @@ -324,113 +382,184 @@ def map_template_coordinate_to_monitor_coordinate(template_coord, monitor_shape, return rx_new, cx_new -def map_monitor_coordinate_to_template_coordinate(monitor_coord, monitor_shape, template_shape): +def map_monitor_coordinate_to_template_coordinate( + monitor_coord, monitor_shape, template_shape +): rx, cx = monitor_coord n_pixels_r, n_pixels_c = monitor_shape tr, tc = template_shape - rx_new = rx - float((n_pixels_r - tr) / 2) - cx_new = cx - float((n_pixels_c - tc) / 2) + rx_new = rx - float((n_pixels_r - tr) / 2) + cx_new = cx - float((n_pixels_c - tc) / 2) return rx_new, cx_new -def lsn_coordinate_to_monitor_coordinate(lsn_coordinate, monitor_shape, stimulus_type): +def lsn_coordinate_to_monitor_coordinate( + lsn_coordinate, monitor_shape, stimulus_type +): template_shape = LOCALLY_SPARSE_NOISE_DIMENSIONS[stimulus_type] pixels_per_patch = LOCALLY_SPARSE_NOISE_PIXELS[stimulus_type] rx, cx = lsn_coordinate tr, tc = template_shape - return map_template_coordinate_to_monitor_coordinate((rx*pixels_per_patch, cx*pixels_per_patch), - monitor_shape, - (tr*pixels_per_patch, tc*pixels_per_patch)) + return map_template_coordinate_to_monitor_coordinate( + (rx * pixels_per_patch, cx * pixels_per_patch), + monitor_shape, + (tr * pixels_per_patch, tc * pixels_per_patch), + ) -def monitor_coordinate_to_lsn_coordinate(monitor_coordinate, monitor_shape, stimulus_type): +def monitor_coordinate_to_lsn_coordinate( + monitor_coordinate, monitor_shape, stimulus_type +): pixels_per_patch = LOCALLY_SPARSE_NOISE_PIXELS[stimulus_type] tr, tc = LOCALLY_SPARSE_NOISE_DIMENSIONS[stimulus_type] - rx, cx = map_monitor_coordinate_to_template_coordinate(monitor_coordinate, monitor_shape, (tr*pixels_per_patch, tc*pixels_per_patch)) - - return (rx/pixels_per_patch, cx/pixels_per_patch) - -def natural_scene_coordinate_to_monitor_coordinate(natural_scene_coordinate, monitor_shape): - - return map_template_coordinate_to_monitor_coordinate(natural_scene_coordinate, monitor_shape, NATURAL_SCENES_PIXELS) - -def natural_movie_coordinate_to_monitor_coordinate(natural_movie_coordinate, monitor_shape): - - local_y = 1.*NATURAL_MOVIE_PIXELS[0]*natural_movie_coordinate[0]/NATURAL_MOVIE_DIMENSIONS[0] - local_x = 1. * NATURAL_MOVIE_PIXELS[1] * natural_movie_coordinate[1] / NATURAL_MOVIE_DIMENSIONS[1] - - return map_template_coordinate_to_monitor_coordinate((local_y, local_x), monitor_shape, NATURAL_MOVIE_PIXELS) - -def map_stimulus_coordinate_to_monitor_coordinate(template_coordinate, monitor_shape, stimulus_type): - + rx, cx = map_monitor_coordinate_to_template_coordinate( + monitor_coordinate, + monitor_shape, + (tr * pixels_per_patch, tc * pixels_per_patch), + ) + + return (rx / pixels_per_patch, cx / pixels_per_patch) + + +def natural_scene_coordinate_to_monitor_coordinate( + natural_scene_coordinate, monitor_shape +): + return map_template_coordinate_to_monitor_coordinate( + natural_scene_coordinate, monitor_shape, NATURAL_SCENES_PIXELS + ) + + +def natural_movie_coordinate_to_monitor_coordinate( + natural_movie_coordinate, monitor_shape +): + local_y = ( + 1.0 + * NATURAL_MOVIE_PIXELS[0] + * natural_movie_coordinate[0] + / NATURAL_MOVIE_DIMENSIONS[0] + ) + local_x = ( + 1.0 + * NATURAL_MOVIE_PIXELS[1] + * natural_movie_coordinate[1] + / NATURAL_MOVIE_DIMENSIONS[1] + ) + + return map_template_coordinate_to_monitor_coordinate( + (local_y, local_x), monitor_shape, NATURAL_MOVIE_PIXELS + ) + + +def map_stimulus_coordinate_to_monitor_coordinate( + template_coordinate, monitor_shape, stimulus_type +): if stimulus_type in LOCALLY_SPARSE_NOISE_STIMULUS_TYPES: - return lsn_coordinate_to_monitor_coordinate(template_coordinate, monitor_shape, stimulus_type) + return lsn_coordinate_to_monitor_coordinate( + template_coordinate, monitor_shape, stimulus_type + ) elif stimulus_type in NATURAL_MOVIE_STIMULUS_TYPES: - return natural_movie_coordinate_to_monitor_coordinate(template_coordinate, monitor_shape) + return natural_movie_coordinate_to_monitor_coordinate( + template_coordinate, monitor_shape + ) elif stimulus_type == NATURAL_SCENES: - return natural_scene_coordinate_to_monitor_coordinate(template_coordinate, monitor_shape) - elif stimulus_type in [DRIFTING_GRATINGS, STATIC_GRATINGS, SPONTANEOUS_ACTIVITY]: + return natural_scene_coordinate_to_monitor_coordinate( + template_coordinate, monitor_shape + ) + elif stimulus_type in [ + DRIFTING_GRATINGS, + STATIC_GRATINGS, + SPONTANEOUS_ACTIVITY, + ]: return template_coordinate else: - raise NotImplementedError # pragma: no cover + raise NotImplementedError # pragma: no cover -def monitor_coordinate_to_natural_movie_coordinate(monitor_coordinate, monitor_shape): - local_y, local_x = map_monitor_coordinate_to_template_coordinate(monitor_coordinate, monitor_shape, NATURAL_MOVIE_PIXELS) +def monitor_coordinate_to_natural_movie_coordinate( + monitor_coordinate, monitor_shape +): + local_y, local_x = map_monitor_coordinate_to_template_coordinate( + monitor_coordinate, monitor_shape, NATURAL_MOVIE_PIXELS + ) - return float(NATURAL_MOVIE_DIMENSIONS[0])*local_y/NATURAL_MOVIE_PIXELS[0], float(NATURAL_MOVIE_DIMENSIONS[1])*local_x/NATURAL_MOVIE_PIXELS[1] + return ( + float(NATURAL_MOVIE_DIMENSIONS[0]) * local_y / NATURAL_MOVIE_PIXELS[0], + float(NATURAL_MOVIE_DIMENSIONS[1]) * local_x / NATURAL_MOVIE_PIXELS[1], + ) -def map_monitor_coordinate_to_stimulus_coordinate(monitor_coordinate, monitor_shape, stimulus_type): +def map_monitor_coordinate_to_stimulus_coordinate( + monitor_coordinate, monitor_shape, stimulus_type +): if stimulus_type in LOCALLY_SPARSE_NOISE_STIMULUS_TYPES: - return monitor_coordinate_to_lsn_coordinate(monitor_coordinate, monitor_shape, stimulus_type) + return monitor_coordinate_to_lsn_coordinate( + monitor_coordinate, monitor_shape, stimulus_type + ) elif stimulus_type == NATURAL_SCENES: - return map_monitor_coordinate_to_template_coordinate(monitor_coordinate, monitor_shape, NATURAL_SCENES_PIXELS) + return map_monitor_coordinate_to_template_coordinate( + monitor_coordinate, monitor_shape, NATURAL_SCENES_PIXELS + ) elif stimulus_type in NATURAL_MOVIE_STIMULUS_TYPES: - return monitor_coordinate_to_natural_movie_coordinate(monitor_coordinate, monitor_shape) - elif stimulus_type in [DRIFTING_GRATINGS, STATIC_GRATINGS, SPONTANEOUS_ACTIVITY]: + return monitor_coordinate_to_natural_movie_coordinate( + monitor_coordinate, monitor_shape + ) + elif stimulus_type in [ + DRIFTING_GRATINGS, + STATIC_GRATINGS, + SPONTANEOUS_ACTIVITY, + ]: return monitor_coordinate else: - raise NotImplementedError # pragma: no cover + raise NotImplementedError # pragma: no cover + -def map_stimulus(source_stimulus_coordinate, source_stimulus_type, target_stimulus_type, monitor_shape): - mc = map_stimulus_coordinate_to_monitor_coordinate(source_stimulus_coordinate, monitor_shape, source_stimulus_type) - return map_monitor_coordinate_to_stimulus_coordinate(mc, monitor_shape, target_stimulus_type) +def map_stimulus( + source_stimulus_coordinate, + source_stimulus_type, + target_stimulus_type, + monitor_shape, +): + mc = map_stimulus_coordinate_to_monitor_coordinate( + source_stimulus_coordinate, monitor_shape, source_stimulus_type + ) + return map_monitor_coordinate_to_stimulus_coordinate( + mc, monitor_shape, target_stimulus_type + ) -def translate_image_and_fill(img, translation=(0,0)): + +def translate_image_and_fill(img, translation=(0, 0)): # first coordinate is horizontal, second is vertical roll = (int(translation[0]), -int(translation[1])) - im2 = np.roll(img, roll, (1,0)) + im2 = np.roll(img, roll, (1, 0)) if roll[1] >= 0: - im2[:roll[1],:] = STIMULUS_GRAY + im2[: roll[1], :] = STIMULUS_GRAY else: - im2[roll[1]:,:] = STIMULUS_GRAY + im2[roll[1] :, :] = STIMULUS_GRAY if roll[0] >= 0: - im2[:,:roll[0]] = STIMULUS_GRAY + im2[:, : roll[0]] = STIMULUS_GRAY else: - im2[:,roll[0]:] = STIMULUS_GRAY + im2[:, roll[0] :] = STIMULUS_GRAY return im2 -class Monitor(object): +class Monitor(object): def __init__(self, n_pixels_r, n_pixels_c, panel_size, spatial_unit): - self.spatial_unit = spatial_unit - if spatial_unit == 'cm': - self.spatial_conversion_factor = 1. + if spatial_unit == "cm": + self.spatial_conversion_factor = 1.0 else: - raise NotImplementedError # pragma: no cover + raise NotImplementedError # pragma: no cover self._panel_size = panel_size self.n_pixels_r = n_pixels_r @@ -445,195 +574,290 @@ def mask(self): @property def panel_size(self): - return self._panel_size*self.spatial_conversion_factor + return self._panel_size * self.spatial_conversion_factor @property def aspect_ratio(self): - return float(self.n_pixels_c)/self.n_pixels_r + return float(self.n_pixels_c) / self.n_pixels_r @property def height(self): - return self.spatial_conversion_factor*np.sqrt(self.panel_size**2/(1+self.aspect_ratio**2)) + return self.spatial_conversion_factor * np.sqrt( + self.panel_size**2 / (1 + self.aspect_ratio**2) + ) @property def width(self): - return self.height*self.aspect_ratio + return self.height * self.aspect_ratio def set_spatial_unit(self, new_unit): if new_unit == self.spatial_unit: pass - elif new_unit == 'inch' and self.spatial_unit == 'cm': - self.spatial_conversion_factor *= .393701 - elif new_unit == 'cm' and self.spatial_unit == 'inch': - self.spatial_conversion_factor *= 1./.393701 + elif new_unit == "inch" and self.spatial_unit == "cm": + self.spatial_conversion_factor *= 0.393701 + elif new_unit == "cm" and self.spatial_unit == "inch": + self.spatial_conversion_factor *= 1.0 / 0.393701 else: - raise NotImplementedError # pragma: no cover + raise NotImplementedError # pragma: no cover self.spatial_unit = new_unit @property def pixel_size(self): - return float(self.width)/self.n_pixels_c - - def pixels_to_visual_degrees(self, n, distance_from_monitor, small_angle_approximation=True): - - - if small_angle_approximation == True: - return n*self.pixel_size/distance_from_monitor*RADIANS_TO_DEGREES # radians to degrees + return float(self.width) / self.n_pixels_c + + def pixels_to_visual_degrees( + self, n, distance_from_monitor, small_angle_approximation=True + ): + if small_angle_approximation: + return ( + n + * self.pixel_size + / distance_from_monitor + * RADIANS_TO_DEGREES + ) # radians to degrees else: - return 2*np.arctan(n*1./2*self.pixel_size / distance_from_monitor) * RADIANS_TO_DEGREES # radians to degrees - - def visual_degrees_to_pixels(self, vd, distance_from_monitor, small_angle_approximation=True): - - if small_angle_approximation == True: - return vd*(distance_from_monitor/self.pixel_size/RADIANS_TO_DEGREES) + return ( + 2 + * np.arctan( + n * 1.0 / 2 * self.pixel_size / distance_from_monitor + ) + * RADIANS_TO_DEGREES + ) # radians to degrees + + def visual_degrees_to_pixels( + self, vd, distance_from_monitor, small_angle_approximation=True + ): + if small_angle_approximation: + return vd * ( + distance_from_monitor / self.pixel_size / RADIANS_TO_DEGREES + ) else: raise NotImplementedError - - def lsn_image_to_screen(self, img, stimulus_type, origin='lower', background_color=STIMULUS_GRAY, translation=(0,0)): - + def lsn_image_to_screen( + self, + img, + stimulus_type, + origin="lower", + background_color=STIMULUS_GRAY, + translation=(0, 0), + ): # assert img.dtype == np.uint8 - - full_image = np.full((self.n_pixels_r, self.n_pixels_c), background_color, dtype=np.uint8) + full_image = np.full( + (self.n_pixels_r, self.n_pixels_c), + background_color, + dtype=np.uint8, + ) pixels_per_patch = float(LOCALLY_SPARSE_NOISE_PIXELS[stimulus_type]) - target_size = tuple( int(pixels_per_patch * dimsize) for dimsize in img.shape[::-1] ) - img_full_res = np.array(Image.fromarray(img).resize(target_size, 0)) # 0 -> nearest neighbor interpolator - - mr, mc = lsn_coordinate_to_monitor_coordinate((0, 0), (self.n_pixels_r, self.n_pixels_c), stimulus_type) - Mr, Mc = lsn_coordinate_to_monitor_coordinate(img.shape, (self.n_pixels_r, self.n_pixels_c), stimulus_type) - full_image[int(mr):int(Mr), int(mc):int(Mc)] = img_full_res - - full_image = translate_image_and_fill(full_image, translation=translation) - - if origin == 'lower': + target_size = tuple( + int(pixels_per_patch * dimsize) for dimsize in img.shape[::-1] + ) + img_full_res = np.array( + Image.fromarray(img).resize(target_size, 0) + ) # 0 -> nearest neighbor interpolator + + mr, mc = lsn_coordinate_to_monitor_coordinate( + (0, 0), (self.n_pixels_r, self.n_pixels_c), stimulus_type + ) + Mr, Mc = lsn_coordinate_to_monitor_coordinate( + img.shape, (self.n_pixels_r, self.n_pixels_c), stimulus_type + ) + full_image[int(mr) : int(Mr), int(mc) : int(Mc)] = img_full_res + + full_image = translate_image_and_fill( + full_image, translation=translation + ) + + if origin == "lower": return full_image - elif origin == 'upper': + elif origin == "upper": return np.flipud(full_image) else: raise Exception return full_image - def natural_scene_image_to_screen(self, img, origin='lower', translation=(0,0)): - - full_image = np.full((self.n_pixels_r, self.n_pixels_c), 127, dtype=np.uint8) - mr, mc = natural_scene_coordinate_to_monitor_coordinate((0, 0), (self.n_pixels_r, self.n_pixels_c)) - Mr, Mc = natural_scene_coordinate_to_monitor_coordinate((img.shape[0], img.shape[1]), (self.n_pixels_r, self.n_pixels_c)) - full_image[int(mr):int(Mr), int(mc):int(Mc)] = img - - full_image = translate_image_and_fill(full_image, translation=translation) - - - if origin == 'lower': + def natural_scene_image_to_screen( + self, img, origin="lower", translation=(0, 0) + ): + full_image = np.full( + (self.n_pixels_r, self.n_pixels_c), 127, dtype=np.uint8 + ) + mr, mc = natural_scene_coordinate_to_monitor_coordinate( + (0, 0), (self.n_pixels_r, self.n_pixels_c) + ) + Mr, Mc = natural_scene_coordinate_to_monitor_coordinate( + (img.shape[0], img.shape[1]), (self.n_pixels_r, self.n_pixels_c) + ) + full_image[int(mr) : int(Mr), int(mc) : int(Mc)] = img + + full_image = translate_image_and_fill( + full_image, translation=translation + ) + + if origin == "lower": return np.flipud(full_image) - elif origin == 'upper': + elif origin == "upper": return full_image else: raise Exception - def natural_movie_image_to_screen(self, img, origin='lower', translation=(0,0)): - - img = np.array(Image.fromarray(img).resize(NATURAL_MOVIE_PIXELS[::-1], 2)).astype(np.uint8) # 2 -> bilinear interpolator + def natural_movie_image_to_screen( + self, img, origin="lower", translation=(0, 0) + ): + img = np.array( + Image.fromarray(img).resize(NATURAL_MOVIE_PIXELS[::-1], 2) + ).astype( + np.uint8 + ) # 2 -> bilinear interpolator assert img.dtype == np.uint8 - full_image = np.full((self.n_pixels_r, self.n_pixels_c), 127, dtype=np.uint8) - mr, mc = map_template_coordinate_to_monitor_coordinate((0, 0), (self.n_pixels_r, self.n_pixels_c), NATURAL_MOVIE_PIXELS) - Mr, Mc = map_template_coordinate_to_monitor_coordinate((img.shape[0], img.shape[1]), (self.n_pixels_r, self.n_pixels_c), NATURAL_MOVIE_PIXELS) - - full_image[int(mr):int(Mr), int(mc):int(Mc)] = img - - full_image = translate_image_and_fill(full_image, translation=translation) - - if origin == 'lower': + full_image = np.full( + (self.n_pixels_r, self.n_pixels_c), 127, dtype=np.uint8 + ) + mr, mc = map_template_coordinate_to_monitor_coordinate( + (0, 0), (self.n_pixels_r, self.n_pixels_c), NATURAL_MOVIE_PIXELS + ) + Mr, Mc = map_template_coordinate_to_monitor_coordinate( + (img.shape[0], img.shape[1]), + (self.n_pixels_r, self.n_pixels_c), + NATURAL_MOVIE_PIXELS, + ) + + full_image[int(mr) : int(Mr), int(mc) : int(Mc)] = img + + full_image = translate_image_and_fill( + full_image, translation=translation + ) + + if origin == "lower": return np.flipud(full_image) - elif origin == 'upper': + elif origin == "upper": return full_image else: raise Exception - def spatial_frequency_to_pix_per_cycle(self, spatial_frequency, distance_from_monitor): - + def spatial_frequency_to_pix_per_cycle( + self, spatial_frequency, distance_from_monitor + ): # How many cycles do I want to see post warp: - number_of_cycles = spatial_frequency*2*np.degrees(np.arctan(self.width/2./distance_from_monitor)) + number_of_cycles = ( + spatial_frequency + * 2 + * np.degrees(np.arctan(self.width / 2.0 / distance_from_monitor)) + ) # How many pixels to I have pre-warp to place my cycles on: _, m_col = np.where(self.mask != 0) - number_of_pixels = (m_col.max() - m_col.min()) - - return float(number_of_pixels)/number_of_cycles - - - def grating_to_screen(self, phase, spatial_frequency, orientation, distance_from_monitor, p2p_amp=256, baseline=127, translation=(0,0)): - - pix_per_cycle = self.spatial_frequency_to_pix_per_cycle(spatial_frequency, distance_from_monitor) - - full_image = get_spatial_grating(height=self.n_pixels_r, - aspect_ratio=self.aspect_ratio, - ori=orientation, - pix_per_cycle=pix_per_cycle, - phase=phase, - p2p_amp=p2p_amp, - baseline=baseline) - - full_image = translate_image_and_fill(full_image, translation=translation) + number_of_pixels = m_col.max() - m_col.min() + + return float(number_of_pixels) / number_of_cycles + + def grating_to_screen( + self, + phase, + spatial_frequency, + orientation, + distance_from_monitor, + p2p_amp=256, + baseline=127, + translation=(0, 0), + ): + pix_per_cycle = self.spatial_frequency_to_pix_per_cycle( + spatial_frequency, distance_from_monitor + ) + + full_image = get_spatial_grating( + height=self.n_pixels_r, + aspect_ratio=self.aspect_ratio, + ori=orientation, + pix_per_cycle=pix_per_cycle, + phase=phase, + p2p_amp=p2p_amp, + baseline=baseline, + ) + + full_image = translate_image_and_fill( + full_image, translation=translation + ) return full_image def get_mask(self): - - mask = make_display_mask(display_shape=(self.n_pixels_c, self.n_pixels_r)).T + mask = make_display_mask( + display_shape=(self.n_pixels_c, self.n_pixels_r) + ).T assert mask.shape[0] == self.n_pixels_r assert mask.shape[1] == self.n_pixels_c return mask - def show_image(self, img, ax=None, show=True, mask=False, warp=False, origin='lower'): + def show_image( + self, img, ax=None, show=True, mask=False, warp=False, origin="lower" + ): import matplotlib.pyplot as plt - assert img.shape == (self.n_pixels_r, self.n_pixels_c) or img.shape == (self.n_pixels_r, self.n_pixels_c, 4) + + assert img.shape == ( + self.n_pixels_r, + self.n_pixels_c, + ) or img.shape == (self.n_pixels_r, self.n_pixels_c, 4) if ax is None: fig, ax = plt.subplots(1, 1) - if warp == True: + if warp: img = self.warp_image(img) - if warp == True: - assert mask == False + if warp: + assert mask is False - ax.imshow(img, origin=origin, cmap=plt.cm.gray, interpolation='none') + ax.imshow(img, origin=origin, cmap=plt.cm.gray, interpolation="none") - if mask == True: - mask = make_display_mask(display_shape=(self.n_pixels_c, self.n_pixels_r)).T + if mask: + mask = make_display_mask( + display_shape=(self.n_pixels_c, self.n_pixels_r) + ).T alpha_mask = np.zeros((mask.shape[0], mask.shape[1], 4)) alpha_mask[:, :, 2] = 1 - mask - alpha_mask[:, :, 3] = .4 - ax.imshow(alpha_mask, origin=origin, interpolation='none') + alpha_mask[:, :, 3] = 0.4 + ax.imshow(alpha_mask, origin=origin, interpolation="none") ax.axes.get_xaxis().set_visible(False) ax.axes.get_yaxis().set_visible(False) - if origin == 'upper': + if origin == "upper": ax.set_ylim((img.shape[0], 0)) - elif origin == 'lower': + elif origin == "lower": ax.set_ylim((0, img.shape[0])) else: raise Exception ax.set_xlim((0, img.shape[1])) - if show == True: + if show: plt.show() - def map_stimulus(self, source_stimulus_coordinate, source_stimulus_type, target_stimulus_type): + def map_stimulus( + self, + source_stimulus_coordinate, + source_stimulus_type, + target_stimulus_type, + ): monitor_shape = (self.n_pixels_r, self.n_pixels_c) - return map_stimulus(source_stimulus_coordinate, source_stimulus_type, target_stimulus_type, monitor_shape) - -class ExperimentGeometry(object): + return map_stimulus( + source_stimulus_coordinate, + source_stimulus_type, + target_stimulus_type, + monitor_shape, + ) - def __init__(self, distance, mon_height_cm, mon_width_cm, mon_res, eyepoint): +class ExperimentGeometry(object): + def __init__( + self, distance, mon_height_cm, mon_width_cm, mon_res, eyepoint + ): self.distance = distance self.mon_height_cm = mon_height_cm self.mon_width_cm = mon_width_cm @@ -650,86 +874,121 @@ def warp_coordinates(self): return self._warp_coordinates def generate_warp_coordinates(self): - - display_shape=self.mon_res + display_shape = self.mon_res x = np.array(range(display_shape[0])) - display_shape[0] / 2 y = np.array(range(display_shape[1])) - display_shape[1] / 2 display_coords = np.array(list(itertools.product(y, x))) - warp_coorinates = warp_stimulus_coords(display_coords, - distance=self.distance, - mon_height_cm=self.mon_height_cm, - mon_width_cm=self.mon_width_cm, - mon_res=self.mon_res, - eyepoint=self.eyepoint) + warp_coorinates = warp_stimulus_coords( + display_coords, + distance=self.distance, + mon_height_cm=self.mon_height_cm, + mon_width_cm=self.mon_width_cm, + mon_res=self.mon_res, + eyepoint=self.eyepoint, + ) warp_coorinates[:, 0] += display_shape[1] / 2 warp_coorinates[:, 1] += display_shape[0] / 2 return warp_coorinates + class BrainObservatoryMonitor(Monitor): - ''' - http://help.brain-map.org/display/observatory/Documentation?preview=/10616846/10813485/VisualCoding_VisualStimuli.pdf + """ + http://help.brain-map.org/display/observatory/Documentation?preview=/10616846/10813485/VisualCoding_VisualStimuli.pdf # noqa: E501 https://www.cnet.com/products/asus-pa248q/specs/ - ''' + """ def __init__(self, experiment_geometry=None): - height, width = MONITOR_DIMENSIONS - super(BrainObservatoryMonitor, self).__init__(height, width, 61.214, 'cm') + super(BrainObservatoryMonitor, self).__init__( + height, width, 61.214, "cm" + ) if experiment_geometry is None: - self.experiment_geometry = ExperimentGeometry(distance=float(MONITOR_DISTANCE), mon_height_cm=self.height, mon_width_cm=self.width, mon_res=(self.n_pixels_c, self.n_pixels_r), eyepoint=(0.5, 0.5)) + self.experiment_geometry = ExperimentGeometry( + distance=float(MONITOR_DISTANCE), + mon_height_cm=self.height, + mon_width_cm=self.width, + mon_res=(self.n_pixels_c, self.n_pixels_r), + eyepoint=(0.5, 0.5), + ) else: self.experiment_geometry = experiment_geometry def lsn_image_to_screen(self, img, **kwargs): - - if img.shape == tuple(LOCALLY_SPARSE_NOISE_DIMENSIONS[LOCALLY_SPARSE_NOISE]): - return super(BrainObservatoryMonitor, self).lsn_image_to_screen(img, LOCALLY_SPARSE_NOISE, **kwargs) - elif img.shape == tuple(LOCALLY_SPARSE_NOISE_DIMENSIONS[LOCALLY_SPARSE_NOISE_4DEG]): - return super(BrainObservatoryMonitor, self).lsn_image_to_screen(img, LOCALLY_SPARSE_NOISE_4DEG, **kwargs) - elif img.shape == tuple(LOCALLY_SPARSE_NOISE_DIMENSIONS[LOCALLY_SPARSE_NOISE_8DEG]): - return super(BrainObservatoryMonitor, self).lsn_image_to_screen(img, LOCALLY_SPARSE_NOISE_8DEG, **kwargs) - else: # pragma: no cover - raise RuntimeError # pragma: no cover + if img.shape == tuple( + LOCALLY_SPARSE_NOISE_DIMENSIONS[LOCALLY_SPARSE_NOISE] + ): + return super(BrainObservatoryMonitor, self).lsn_image_to_screen( + img, LOCALLY_SPARSE_NOISE, **kwargs + ) + elif img.shape == tuple( + LOCALLY_SPARSE_NOISE_DIMENSIONS[LOCALLY_SPARSE_NOISE_4DEG] + ): + return super(BrainObservatoryMonitor, self).lsn_image_to_screen( + img, LOCALLY_SPARSE_NOISE_4DEG, **kwargs + ) + elif img.shape == tuple( + LOCALLY_SPARSE_NOISE_DIMENSIONS[LOCALLY_SPARSE_NOISE_8DEG] + ): + return super(BrainObservatoryMonitor, self).lsn_image_to_screen( + img, LOCALLY_SPARSE_NOISE_8DEG, **kwargs + ) + else: # pragma: no cover + raise RuntimeError # pragma: no cover def warp_image(self, img, **kwargs): - assert img.shape == (self.n_pixels_r, self.n_pixels_c) - assert self.spatial_unit == 'cm' - - return spndi.map_coordinates(img, self.experiment_geometry.warp_coordinates.T).reshape((self.n_pixels_r, self.n_pixels_c)) - - def grating_to_screen(self, phase, spatial_frequency, orientation, **kwargs): - - return super(BrainObservatoryMonitor, self).grating_to_screen(phase, spatial_frequency, orientation, - self.experiment_geometry.distance, - p2p_amp = 256, baseline = 127, **kwargs) + assert self.spatial_unit == "cm" + + return spndi.map_coordinates( + img, self.experiment_geometry.warp_coordinates.T + ).reshape((self.n_pixels_r, self.n_pixels_c)) + + def grating_to_screen( + self, phase, spatial_frequency, orientation, **kwargs + ): + return super(BrainObservatoryMonitor, self).grating_to_screen( + phase, + spatial_frequency, + orientation, + self.experiment_geometry.distance, + p2p_amp=256, + baseline=127, + **kwargs, + ) def pixels_to_visual_degrees(self, n, **kwargs): - - return super(BrainObservatoryMonitor, self).pixels_to_visual_degrees(n, self.experiment_geometry.distance, **kwargs) + return super(BrainObservatoryMonitor, self).pixels_to_visual_degrees( + n, self.experiment_geometry.distance, **kwargs + ) def visual_degrees_to_pixels(self, vd, **kwargs): - - return super(BrainObservatoryMonitor, self).visual_degrees_to_pixels(vd, self.experiment_geometry.distance, **kwargs) - -def warp_stimulus_coords(vertices, - distance=15.0, - mon_height_cm=32.5, - mon_width_cm=51.0, - mon_res=(1920, 1200), - eyepoint=(0.5, 0.5)): - ''' - For a list of screen vertices, provides a corresponding list of texture coordinates. + return super(BrainObservatoryMonitor, self).visual_degrees_to_pixels( + vd, self.experiment_geometry.distance, **kwargs + ) + + +def warp_stimulus_coords( + vertices, + distance=15.0, + mon_height_cm=32.5, + mon_width_cm=51.0, + mon_res=(1920, 1200), + eyepoint=(0.5, 0.5), +): + """ + For a list of screen vertices, provides a corresponding list of texture + coordinates. Parameters ---------- vertices: numpy.ndarray - [[x0,y0], [x1,y1], ...] A set of vertices to convert to texture positions. + [[x0,y0], [x1,y1], ...] A set of vertices to convert to texture + positions. distance: float distance from the monitor in cm. mon_height_cm: float @@ -743,17 +1002,18 @@ def warp_stimulus_coords(vertices, Returns ------- np.ndarray - x,y coordinates shaped like the input that describe what pixel coordinates - are displayed an the input coordinates after warping the stimulus. + x,y coordinates shaped like the input that describe what pixel + coordinates are displayed an the input coordinates after warping the + stimulus. - ''' + """ mon_width_cm = float(mon_width_cm) mon_height_cm = float(mon_height_cm) distance = float(distance) mon_res_x, mon_res_y = float(mon_res[0]), float(mon_res[1]) - vertices = vertices.astype('float') + vertices = vertices.astype("float") # from pixels (-1920/2 -> 1920/2) to stimulus space (-0.5->0.5) vertices[:, 0] = vertices[:, 0] / mon_res_x @@ -806,22 +1066,30 @@ def warp_stimulus_coords(vertices, def make_display_mask(display_shape=(1920, 1200)): - ''' Build a display-shaped mask that indicates which pixels are on screen after warping the stimulus. ''' + """Build a display-shaped mask that indicates which pixels are on screen + after warping the stimulus. + """ x = np.array(range(display_shape[0])) - display_shape[0] / 2 y = np.array(range(display_shape[1])) - display_shape[1] / 2 display_coords = np.array(list(itertools.product(x, y))) warped_coords = warp_stimulus_coords(display_coords).astype(int) - off_warped_coords = np.array([warped_coords[:, 0] + display_shape[0] / 2, - warped_coords[:, 1] + display_shape[1] / 2]) + off_warped_coords = np.array( + [ + warped_coords[:, 0] + display_shape[0] / 2, + warped_coords[:, 1] + display_shape[1] / 2, + ] + ) used_coords = set() for i in range(off_warped_coords.shape[1]): used_coords.add((off_warped_coords[0, i], off_warped_coords[1, i])) - used_coords = (np.array([x for (x, y) in used_coords]).astype(int), - np.array([y for (x, y) in used_coords]).astype(int)) + used_coords = ( + np.array([x for (x, y) in used_coords]).astype(int), + np.array([y for (x, y) in used_coords]).astype(int), + ) mask = np.zeros(display_shape) @@ -830,9 +1098,12 @@ def make_display_mask(display_shape=(1920, 1200)): return mask -def mask_stimulus_template(template_display_coords, template_shape, display_mask=None, threshold=1.0): - ''' Build a mask for a stimulus template of a given shape and display coordinates that indicates - which part of the template is on screen after warping. +def mask_stimulus_template( + template_display_coords, template_shape, display_mask=None, threshold=1.0 +): + """Build a mask for a stimulus template of a given shape and display + coordinates that indicates which part of the template is on screen after + warping. Parameters ---------- @@ -843,16 +1114,17 @@ def mask_stimulus_template(template_display_coords, template_shape, display_mask (width,height) of the display template display_mask: np.ndarray - boolean 2D mask indicating which display coordinates are on screen after warping. + boolean 2D mask indicating which display coordinates are on screen + after warping. threshold: float - Fraction of pixels associated with a template display coordinate that should remain - on screen to count as belonging to the mask. + Fraction of pixels associated with a template display coordinate that + should remain on screen to count as belonging to the mask. Returns ------- tuple: (template mask, pixel fraction) - ''' + """ if display_mask is None: display_mask = make_display_mask() @@ -860,8 +1132,10 @@ def mask_stimulus_template(template_display_coords, template_shape, display_mask mask = np.zeros(template_shape, dtype=bool) for y in range(template_shape[1]): for x in range(template_shape[0]): - tdcm = np.where((template_display_coords[0, :, :] == x) & ( - template_display_coords[1, :, :] == y)) + tdcm = np.where( + (template_display_coords[0, :, :] == x) + & (template_display_coords[1, :, :] == y) + ) v = display_mask[tdcm] f = np.sum(v) / len(v) frac[x, y] = f diff --git a/allensdk/brain_observatory/vbn_2022/metadata_writer/metadata_writer.py b/allensdk/brain_observatory/vbn_2022/metadata_writer/metadata_writer.py index bde2ea72b..7112d261f 100644 --- a/allensdk/brain_observatory/vbn_2022/metadata_writer/metadata_writer.py +++ b/allensdk/brain_observatory/vbn_2022/metadata_writer/metadata_writer.py @@ -4,13 +4,16 @@ import allensdk import argschema import pandas as pd -from allensdk.brain_observatory.data_release_utils.metadata_utils.id_generator import ( # noqa +from allensdk.brain_observatory.behavior.behavior_project_cache.project_metadata_writer.schemas import ( # noqa: E501 + DataReleaseToolsInputSchema, +) +from allensdk.brain_observatory.data_release_utils.metadata_utils.id_generator import ( # noqa: E501 FileIDGenerator, ) -from allensdk.brain_observatory.data_release_utils.metadata_utils.utils import ( # noqa +from allensdk.brain_observatory.data_release_utils.metadata_utils.utils import ( # noqa: E501 add_file_paths_to_metadata_table, ) -from allensdk.brain_observatory.vbn_2022.metadata_writer.dataframe_manipulations import ( # noqa +from allensdk.brain_observatory.vbn_2022.metadata_writer.dataframe_manipulations import ( # noqa: E501 strip_substructure_acronym_df, ) from allensdk.brain_observatory.vbn_2022.metadata_writer.lims_queries import ( @@ -20,7 +23,6 @@ session_tables_from_ecephys_session_id_list, units_table_from_ecephys_session_id_list, ) -from allensdk.brain_observatory.behavior.behavior_project_cache.project_metadata_writer.schemas import DataReleaseToolsInputSchema # noqa: E501 from allensdk.brain_observatory.vbn_2022.metadata_writer.schemas import ( VBN2022MetadataWriterInputSchema, ) @@ -181,7 +183,7 @@ def run(self): ecephys_session_id_list=session_id_list, failed_ecephys_session_id_list=failed_session_list, probe_ids_to_skip=probe_ids_to_skip, - n_workers=self.args["n_workers"] + n_workers=self.args["n_workers"], ) ecephys_nwb_dir = pathlib.Path(self.args["ecephys_nwb_dir"]) diff --git a/allensdk/core/dataframe_utils.py b/allensdk/core/dataframe_utils.py index 2b674ea97..8bf63814b 100644 --- a/allensdk/core/dataframe_utils.py +++ b/allensdk/core/dataframe_utils.py @@ -1,4 +1,5 @@ from typing import List + import pandas as pd # Null value fill for integer Pandas.Series objects. NWB currently doesn't @@ -10,10 +11,11 @@ def patch_df_from_other( - target_df: pd.DataFrame, - source_df: pd.DataFrame, - columns_to_patch: List[str], - index_column: str) -> pd.DataFrame: + target_df: pd.DataFrame, + source_df: pd.DataFrame, + columns_to_patch: List[str], + index_column: str, +) -> pd.DataFrame: """ Overwrite column values in target_df from column values in source_df in rows where the two dataframes @@ -70,9 +72,11 @@ def patch_df_from_other( target_df[column] = None if index_column in columns_to_patch: - msg += (f"{index_column} is in the list of " - f"columns to patch {columns_to_patch}; " - "unsure how to handle that case\n") + msg += ( + f"{index_column} is in the list of " + f"columns to patch {columns_to_patch}; " + "unsure how to handle that case\n" + ) if len(msg) > 0: msg = f"failures in patch_df_from_other:\n{msg}" @@ -83,10 +87,7 @@ def patch_df_from_other( patch_df = source_df[columns_to_patch + [index_column]] patch_df = patch_df.set_index(index_column) - target_df.update( - patch_df, - join='left', - overwrite=True) + target_df.update(patch_df, join="left", overwrite=True) target_df = target_df.reset_index() if original_index is not None: @@ -95,8 +96,8 @@ def patch_df_from_other( def enforce_df_column_order( - input_df: pd.DataFrame, - column_order: List[str]) -> pd.DataFrame: + input_df: pd.DataFrame, column_order: List[str] +) -> pd.DataFrame: """Return the data frame but with columns ordered. Parameters @@ -121,13 +122,14 @@ def enforce_df_column_order( # Get the full list of columns in the data frame with our ordered columns # first. pruned_order.extend( - list(set(input_df.columns).difference(set(pruned_order)))) + list(set(input_df.columns).difference(set(pruned_order))) + ) return input_df[pruned_order] -def enforce_df_int_typing(input_df: pd.DataFrame, - int_columns: List[str], - use_pandas_type=False) -> pd.DataFrame: +def enforce_df_int_typing( + input_df: pd.DataFrame, int_columns: List[str], use_pandas_type=False +) -> pd.DataFrame: """Enforce integer typing for columns that may have lost int typing when combined into the final DataFrame. @@ -153,16 +155,15 @@ def enforce_df_int_typing(input_df: pd.DataFrame, for col in int_columns: if col in input_df.columns: if use_pandas_type: - input_df[col] = input_df[col].astype('Int64') + input_df[col] = input_df[col].astype("Int64") else: - input_df[col] = \ - input_df[col].fillna(INT_NULL).astype(int) + input_df[col] = input_df[col].fillna(INT_NULL).astype(int) return input_df -def return_one_dataframe_row_only(input_table: pd.DataFrame, - index_value: int, - table_name: str) -> pd.Series: +def return_one_dataframe_row_only( + input_table: pd.DataFrame, index_value: int, table_name: str +) -> pd.Series: """Lookup and return one and only one row from the DataFrame returning an informative error if no or multiple rows are returned for a given index. @@ -204,5 +205,5 @@ def return_one_dataframe_row_only(input_table: pd.DataFrame, f"{input_table.index.name}. For " f"{index_value} " f" there are {len(row)} entries." - ) + ) return row diff --git a/allensdk/core/json_utilities.py b/allensdk/core/json_utilities.py index e85649516..6c26196f0 100644 --- a/allensdk/core/json_utilities.py +++ b/allensdk/core/json_utilities.py @@ -33,11 +33,11 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # +import logging +import re + import numpy as np import simplejson as json -import math -import re -import logging ju_logger = logging.getLogger(__name__) @@ -52,45 +52,51 @@ def read(file_name): - """ Shortcut reading JSON from a file. """ - with open(file_name, 'rb') as f: - json_string = f.read().decode('utf-8') - if len(json_string)==0: # If empty file - json_string='{}' # Create a string that will give an empty JSON object instead of an error + """Shortcut reading JSON from a file.""" + with open(file_name, "rb") as f: + json_string = f.read().decode("utf-8") + if len(json_string) == 0: # If empty file + # Create a string that will give an empty JSON object instead of an + # error + json_string = "{}" json_obj = json.loads(json_string) return json_obj def write(file_name, obj): - """ Shortcut for writing JSON to a file. This also takes care of serializing numpy and data types. """ - with open(file_name, 'wb') as f: + """Shortcut for writing JSON to a file. This also takes care of + serializing numpy and data types.""" + with open(file_name, "wb") as f: try: - f.write(write_string(obj)) # Python 2.7 + f.write(write_string(obj)) # Python 2.7 except TypeError: - f.write(bytes(write_string(obj), 'utf-8')) # Python 3 + f.write(bytes(write_string(obj), "utf-8")) # Python 3 def write_string(obj): - """ Shortcut for writing JSON to a string. This also takes care of serializing numpy and data types. """ - return json.dumps(obj, - indent=2, - ignore_nan=True, - default=json_handler, - iterable_as_array=True) - - -def read_url(url, method='POST'): - if method == 'GET': + """Shortcut for writing JSON to a string. This also takes care of + serializing numpy and data types.""" + return json.dumps( + obj, + indent=2, + ignore_nan=True, + default=json_handler, + iterable_as_array=True, + ) + + +def read_url(url, method="POST"): + if method == "GET": return read_url_get(url) - elif method == 'POST': + elif method == "POST": return read_url_post(url) else: - raise Exception('Unknown request method: (%s)' % method) + raise Exception("Unknown request method: (%s)" % method) def read_url_get(url): - '''Transform a JSON contained in a file into an equivalent + """Transform a JSON contained in a file into an equivalent nested python dict. Parameters @@ -105,15 +111,15 @@ def read_url_get(url): Note: if the input is a bare array or literal, for example, the output will be of the corresponding type. - ''' + """ response = urllib_request.urlopen(url) - json_string = response.read().decode('utf-8') + json_string = response.read().decode("utf-8") return json.loads(json_string) def read_url_post(url): - '''Transform a JSON contained in a file into an equivalent + """Transform a JSON contained in a file into an equivalent nested python dict. Parameters @@ -128,18 +134,19 @@ def read_url_post(url): Note: if the input is a bare array or literal, for example, the output will be of the corresponding type. - ''' + """ urlp = urlparse.urlparse(url) main_url = urlparse.urlunsplit( - (urlp.scheme, urlp.netloc, urlp.path, '', '')) + (urlp.scheme, urlp.netloc, urlp.path, "", "") + ) data = json.dumps(dict(urlparse.parse_qsl(urlp.query))) handler = urllib_request.HTTPHandler() opener = urllib_request.build_opener(handler) request = urllib_request.Request(main_url, data) - request.add_header("Content-Type", 'application/json') - request.get_method = lambda: 'POST' + request.add_header("Content-Type", "application/json") + request.get_method = lambda: "POST" try: response = opener.open(request) @@ -155,8 +162,9 @@ def read_url_post(url): def json_handler(obj): - """ Used by write_json convert a few non-standard types to things that the json package can handle. """ - if hasattr(obj, 'to_dict'): + """Used by write_json convert a few non-standard types to things that the + json package can handle.""" + if hasattr(obj, "to_dict"): return obj.to_dict() elif isinstance(obj, np.ndarray): return obj.tolist() @@ -164,26 +172,21 @@ def json_handler(obj): return float(obj) elif isinstance(obj, np.integer): return int(obj) - elif (isinstance(obj, bool) or - isinstance(obj, np.bool_)): + elif isinstance(obj, bool) or isinstance(obj, np.bool_): return bool(obj) - elif hasattr(obj, 'isoformat'): + elif hasattr(obj, "isoformat"): return obj.isoformat() else: raise TypeError( - 'Object of type %s with value of %s is not JSON serializable' % - (type(obj), repr(obj))) + "Object of type %s with value of %s is not JSON serializable" + % (type(obj), repr(obj)) + ) class JsonComments(object): - _oneline_comment = re.compile(r"\/\/.*$", - re.MULTILINE) - _multiline_comment_start = re.compile(r"\/\*", - re.MULTILINE | - re.DOTALL) - _multiline_comment_end = re.compile(r"\*\/", - re.MULTILINE | - re.DOTALL) + _oneline_comment = re.compile(r"\/\/.*$", re.MULTILINE) + _multiline_comment_start = re.compile(r"\/\*", re.MULTILINE | re.DOTALL) + _multiline_comment_end = re.compile(r"\*\/", re.MULTILINE | re.DOTALL) _blank_line = re.compile(r"\n?^\s*$", re.MULTILINE) _carriage_return = re.compile(r"\r$", re.MULTILINE) @@ -202,12 +205,13 @@ def read_file(cls, file_name): return json_object except ValueError: ju_logger.error( - "Could not load json object from file: %s" % (file_name)) + "Could not load json object from file: %s" % (file_name) + ) raise @classmethod def remove_comments(cls, json_string): - '''Strip single and multiline javascript-style comments. + """Strip single and multiline javascript-style comments. Parameters ---------- @@ -220,17 +224,17 @@ def remove_comments(cls, json_string): Copy of the input with comments removed. Note: A JSON decoder MAY accept and ignore comments. - ''' - json_string = JsonComments._oneline_comment.sub('', json_string) - json_string = JsonComments._carriage_return.sub('', json_string) + """ + json_string = JsonComments._oneline_comment.sub("", json_string) + json_string = JsonComments._carriage_return.sub("", json_string) json_string = JsonComments.remove_multiline_comments(json_string) - json_string = JsonComments._blank_line.sub('', json_string) + json_string = JsonComments._blank_line.sub("", json_string) return json_string @classmethod def remove_multiline_comments(cls, json_string): - '''Rebuild input without substrings matching /*...*/. + """Rebuild input without substrings matching /*...*/. Parameters ---------- @@ -241,22 +245,24 @@ def remove_multiline_comments(cls, json_string): ------- string Copy of the input without the comments. - ''' + """ new_json = [] start_iter = JsonComments._multiline_comment_start.finditer( - json_string) + json_string + ) json_slice_start = 0 for comment_start in start_iter: json_slice_end = comment_start.start() new_json.append(json_string[json_slice_start:json_slice_end]) search_start = comment_start.end() - comment_end = JsonComments._multiline_comment_end.search(json_string[ - search_start:]) + comment_end = JsonComments._multiline_comment_end.search( + json_string[search_start:] + ) if comment_end is None: break else: json_slice_start = search_start + comment_end.end() new_json.append(json_string[json_slice_start:]) - return ''.join(new_json) + return "".join(new_json) diff --git a/allensdk/core/mouse_connectivity_cache.py b/allensdk/core/mouse_connectivity_cache.py index 7a431ee36..fd4c6901b 100644 --- a/allensdk/core/mouse_connectivity_cache.py +++ b/allensdk/core/mouse_connectivity_cache.py @@ -33,23 +33,23 @@ # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # -from allensdk.api.warehouse_cache.cache import Cache, get_default_manifest_file -from allensdk.api.queries.mouse_connectivity_api import MouseConnectivityApi - -from .reference_space_cache import ReferenceSpaceCache +import functools +import operator as op +import os +import re +import warnings import nrrd -import re -import os -import SimpleITK as sitk -import pandas as pd import numpy as np +import pandas as pd +import SimpleITK as sitk +from allensdk.api.queries.mouse_connectivity_api import MouseConnectivityApi +from allensdk.api.warehouse_cache.cache import Cache, get_default_manifest_file from allensdk.config.manifest import Manifest -import warnings -import operator as op -import functools from six.moves import reduce +from .reference_space_cache import ReferenceSpaceCache + class MouseConnectivityCache(ReferenceSpaceCache): """ @@ -478,6 +478,7 @@ def col_rn(x): return pd.DataFrame(x).rename( columns={"section_data_set_id": "experiment_id"} ) + return self.api.get_structure_unionizes( [experiment_id], path=file_name, @@ -572,11 +573,17 @@ def rank_structures( :, output_keys ] - this_experiment_unionizes = unionizes[unionizes['experiment_id'] == eid] - this_experiment_unionizes = this_experiment_unionizes.sort_values(by=rank_on, ascending=False) - this_experiment_unionizes = this_experiment_unionizes.loc[:, output_keys] + this_experiment_unionizes = unionizes[ + unionizes["experiment_id"] == eid + ] + this_experiment_unionizes = this_experiment_unionizes.sort_values( + by=rank_on, ascending=False + ) + this_experiment_unionizes = this_experiment_unionizes.loc[ + :, output_keys + ] - records = this_experiment_unionizes.to_dict('records') + records = this_experiment_unionizes.to_dict("records") if len(records) > n: records = records[:n] results.append(records) @@ -867,9 +874,7 @@ def get_affine_parameters( if direction not in ("trv", "tvr"): raise ValueError( "invalid direction: {}. direction must be one of tvr," - "trv".format( - direction - ) + "trv".format(direction) ) file_name = self.get_cache_path(file_name, self.ALIGNMENT3D_KEY) @@ -878,7 +883,7 @@ def get_affine_parameters( strategy="lazy", path=file_name, section_data_set_id=section_data_set_id, - **Cache.cache_json() + **Cache.cache_json(), ) alignment_re = re.compile(r"{}_(?P\d+)".format(direction)) diff --git a/allensdk/internal/mouse_connectivity/interval_unionize/data_utilities.py b/allensdk/internal/mouse_connectivity/interval_unionize/data_utilities.py index 9fcb57e10..aed561f3f 100755 --- a/allensdk/internal/mouse_connectivity/interval_unionize/data_utilities.py +++ b/allensdk/internal/mouse_connectivity/interval_unionize/data_utilities.py @@ -1,104 +1,119 @@ import logging -#import SimpleITK as sitk +# import SimpleITK as sitk import nrrd import numpy as np -#def read(path): -# return np.swapaxes(sitk.GetArrayFromImage(sitk.ReadImage(str(path))), 0, 2) +# def read(path): +# return np.swapaxes( +# sitk.GetArrayFromImage(sitk.ReadImage(str(path))), 0, 2) + def read(path): return np.ascontiguousarray(nrrd.read(path)[0]) - def load_annotation(annotation_path, data_mask_path=None): - '''Read data files segmenting the reference space into regions of valid + """Read data files segmenting the reference space into regions of valid and invalid data, then further among brain structures - ''' - - logging.info('getting annotation') + """ + + logging.info("getting annotation") annotation = read(annotation_path) - - logging.info('casting to signed') - # It shouldn't matter now, but there may be future structures with ids + + logging.info("casting to signed") + # It shouldn't matter now, but there may be future structures with ids # sufficiently large that we need that extra bit - logging.debug('max annotated value: {0}'.format(np.amax(annotation))) + logging.debug("max annotated value: {0}".format(np.amax(annotation))) annotation = annotation.astype(np.int32) - logging.debug('max annotated value: {0}'.format(np.amax(annotation))) - - logging.info('negating left hemisphere') - logging.debug('min annotated value: {0}'.format(np.amin(annotation))) - lr_mid = int( np.round(annotation.shape[2] / 2) ) + logging.debug("max annotated value: {0}".format(np.amax(annotation))) + + logging.info("negating left hemisphere") + logging.debug("min annotated value: {0}".format(np.amin(annotation))) + lr_mid = int(np.round(annotation.shape[2] / 2)) annotation[:, :, :lr_mid] = annotation[:, :, :lr_mid] * -1 - logging.debug('min annotated value: {0}'.format(np.amin(annotation))) + logging.debug("min annotated value: {0}".format(np.amin(annotation))) if data_mask_path is not None: - logging.info('getting_data_mask') - data_mask =read(data_mask_path) - - logging.info('applying data mask') + logging.info("getting_data_mask") + data_mask = read(data_mask_path) + + logging.info("applying data mask") annotation[np.logical_not(data_mask)] = 0 - + return annotation -def get_sum_pixels(sum_pixels_path): - logging.info('getting sum_pixels') - return {'sum_pixels': read(sum_pixels_path)} - +def get_sum_pixels(sum_pixels_path): + logging.info("getting sum_pixels") + return {"sum_pixels": read(sum_pixels_path)} -def get_sum_pixel_intensities(sum_pixel_intensities_path, injection_sum_pixel_intensities_path): - logging.info('getting sum pixel intensities') - return {'sum_pixel_intensities': read(sum_pixel_intensities_path), - 'injection_sum_pixel_intensities': read(injection_sum_pixel_intensities_path)} +def get_sum_pixel_intensities( + sum_pixel_intensities_path, injection_sum_pixel_intensities_path +): + logging.info("getting sum pixel intensities") + return { + "sum_pixel_intensities": read(sum_pixel_intensities_path), + "injection_sum_pixel_intensities": read( + injection_sum_pixel_intensities_path + ), + } -def get_cav_density(cav_density_path): - logging.info('getting cav density') - return {'cav_density': read(cav_density_path)} +def get_cav_density(cav_density_path): + logging.info("getting cav density") + return {"cav_density": read(cav_density_path)} -def get_injection_data(injection_fraction_path, injection_density_path, - injection_energy_path): - '''Read nrrd files containing injection signal data - ''' - logging.info('getting injection_fraction') +def get_injection_data( + injection_fraction_path, injection_density_path, injection_energy_path +): + """Read nrrd files containing injection signal data""" + + logging.info("getting injection_fraction") injection_fraction = read(injection_fraction_path) - - logging.info('getting injection_sum_projecting_pixels') - injection_density = read(injection_density_path) - - logging.info('getting injection_energy') + + logging.info("getting injection_sum_projecting_pixels") + injection_density = read(injection_density_path) + + logging.info("getting injection_energy") injection_energy = read(injection_energy_path) - - return {'injection_fraction': injection_fraction, - 'injection_density': injection_density, - 'injection_energy': injection_energy} - - -def get_projection_data(projection_density_path, projection_energy_path, - aav_exclusion_fraction_path=None): - '''Read nrrd files containing global signal data - ''' - - logging.info('getting projection density') - projection_density = read(projection_density_path) - - logging.info('getting projection energy') + + return { + "injection_fraction": injection_fraction, + "injection_density": injection_density, + "injection_energy": injection_energy, + } + + +def get_projection_data( + projection_density_path, + projection_energy_path, + aav_exclusion_fraction_path=None, +): + """Read nrrd files containing global signal data""" + + logging.info("getting projection density") + projection_density = read(projection_density_path) + + logging.info("getting projection energy") projection_energy = read(projection_energy_path) - + try: - logging.info('getting aav exclusion fraction') + logging.info("getting aav exclusion fraction") aav_exclusion_fraction = read(aav_exclusion_fraction_path) aav_exclusion_fraction[aav_exclusion_fraction > 0] = 1 - aav_exclusion_fraction = aav_exclusion_fraction.astype(bool, order='C') - + aav_exclusion_fraction = aav_exclusion_fraction.astype(bool, order="C") + except (IOError, OSError, RuntimeError): - logging.info('skipping aav exclusion fraction') - aav_exclusion_fraction = np.zeros(projection_density.shape, dtype=bool, order='C') - - return {'projection_density': projection_density, - 'projection_energy': projection_energy, - 'aav_exclusion_fraction': aav_exclusion_fraction} + logging.info("skipping aav exclusion fraction") + aav_exclusion_fraction = np.zeros( + projection_density.shape, dtype=bool, order="C" + ) + + return { + "projection_density": projection_density, + "projection_energy": projection_energy, + "aav_exclusion_fraction": aav_exclusion_fraction, + } diff --git a/allensdk/test/brain_observatory/behavior/behavior_project_cache/test_behavior_neuropixels_project_cloud_api.py b/allensdk/test/brain_observatory/behavior/behavior_project_cache/test_behavior_neuropixels_project_cloud_api.py index 9ec9a72d2..1dbb44794 100644 --- a/allensdk/test/brain_observatory/behavior/behavior_project_cache/test_behavior_neuropixels_project_cloud_api.py +++ b/allensdk/test/brain_observatory/behavior/behavior_project_cache/test_behavior_neuropixels_project_cloud_api.py @@ -1,29 +1,34 @@ -import pytest -import pandas as pd from pathlib import Path from unittest.mock import MagicMock, create_autospec, patch +import pandas as pd +import pytest from allensdk.api.cloud_cache.manifest import Manifest -from allensdk.brain_observatory.behavior.behavior_project_cache.project_apis.\ - data_io import behavior_neuropixels_project_cloud_api as cloudapi # noqa: E501 - -from allensdk.brain_observatory.behavior.behavior_project_cache.project_apis.\ - data_io import project_cloud_api_base as cloudapibase # noqa: E501 - -from allensdk.brain_observatory.behavior.behavior_project_cache.\ - utils import version_check, BehaviorCloudCacheVersionException -from allensdk.brain_observatory.ecephys.behavior_ecephys_session import \ - BehaviorEcephysSession - - -class MockCache(): - def __init__(self, - behavior_session_table, - ecephys_session_table, - probe_table, - channel_table, - unit_table, - cachedir): +from allensdk.brain_observatory.behavior.behavior_project_cache.project_apis.data_io import ( # noqa: E501 + behavior_neuropixels_project_cloud_api as cloudapi, +) +from allensdk.brain_observatory.behavior.behavior_project_cache.project_apis.data_io import ( # noqa: E501 + project_cloud_api_base as cloudapibase, +) +from allensdk.brain_observatory.behavior.behavior_project_cache.utils import ( + BehaviorCloudCacheVersionException, + version_check, +) +from allensdk.brain_observatory.ecephys.behavior_ecephys_session import ( + BehaviorEcephysSession, +) + + +class MockCache: + def __init__( + self, + behavior_session_table, + ecephys_session_table, + probe_table, + channel_table, + unit_table, + cachedir, + ): self.file_id_column = "file_id" self.session_table_path = cachedir / "ecephys_sessions.csv" self.behavior_session_table_path = cachedir / "behavior_sessions.csv" @@ -35,21 +40,25 @@ def __init__(self, channel_table.to_csv(self.channel_table_path, index=False) probe_table.to_csv(self.probe_table_path, index=False) unit_table.to_csv(self.unit_table_path, index=False) - behavior_session_table.to_csv(self.behavior_session_table_path, - index=False) + behavior_session_table.to_csv( + self.behavior_session_table_path, index=False + ) self._manifest = MagicMock() - self._manifest.metadata_file_names = ["behavior_sessions", - "ecephys_session", - "probes", - "channels", - "units"] + self._manifest.metadata_file_names = [ + "behavior_sessions", + "ecephys_session", + "probes", + "channels", + "units", + ] self._metadata_name_path_map = { - "behavior_sessions": self.behavior_session_table_path, - "ecephys_sessions": self.session_table_path, - "channels": self.channel_table_path, - "probes": self.probe_table_path, - "units": self.unit_table_path} + "behavior_sessions": self.behavior_session_table_path, + "ecephys_sessions": self.session_table_path, + "channels": self.channel_table_path, + "probes": self.probe_table_path, + "units": self.unit_table_path, + } def download_metadata(self, fname): return self._metadata_name_path_map[fname] @@ -59,16 +68,10 @@ def download_data(self, file_id): def metadata_path(self, fname): local_path = self._metadata_name_path_map[fname] - return { - 'local_path': local_path, - 'exists': Path(local_path).exists() - } + return {"local_path": local_path, "exists": Path(local_path).exists()} def data_path(self, file_id): - return { - 'local_path': file_id, - 'exists': True - } + return {"local_path": file_id, "exists": True} def load_last_manifest(self): return None @@ -77,55 +80,67 @@ def load_last_manifest(self): @pytest.fixture def mock_cache(tmpdir): c = { - "behavior_sessions": pd.DataFrame({ + "behavior_sessions": pd.DataFrame( + { "behavior_session_id": [1, 2, 3, 4], "ecephys_session_id": pd.Series( - [10, 11, 12, 13], - dtype='Int64'), - "mouse_id": [4, 4, 2, 1]}), - "ecephys_sessions": pd.DataFrame({ + [10, 11, 12, 13], dtype="Int64" + ), + "mouse_id": [4, 4, 2, 1], + } + ), + "ecephys_sessions": pd.DataFrame( + { "ecephys_session_id": pd.Series( - [10, 11, 12, 13], - dtype='Int64'), + [10, 11, 12, 13], dtype="Int64" + ), "behavior_session_id": [1, 2, 3, 4], - "file_id": [10, 11, 12, 13]}), - "probes": pd.DataFrame({ + "file_id": [10, 11, 12, 13], + } + ), + "probes": pd.DataFrame( + { "ecephys_probe_id": [4, 5, 6, 7], "name": ["probeA", "probeB", "probeA", "probeB"], "ecephys_session_id": [10, 10, 11, 11], "has_lfp_data": [True] * 4, "lfp_sampling_rate": [1, 2, 3, 4], - "file_id": [0, 1, 2, 3] - }), - "channels": pd.DataFrame({ + "file_id": [0, 1, 2, 3], + } + ), + "channels": pd.DataFrame( + { "ecephys_channel_id": [14, 15, 16], "ecephys_probe_id": [4, 4, 4], - "ecephys_session_id": [10, 10, 10] - }), - "units": pd.DataFrame({ + "ecephys_session_id": [10, 10, 10], + } + ), + "units": pd.DataFrame( + { "unit_id": [204, 205, 206], "ecephys_channel_id": [14, 15, 16], "ecephys_probe_id": [4, 4, 4], - "ecephys_session_id": [10, 10, 10] - }), + "ecephys_session_id": [10, 10, 10], + } + ), } # round-trip the tables through csv to pick up # pandas mods to lists fname = tmpdir / "my.csv" - c['behavior_sessions'].to_csv(fname, index=False) + c["behavior_sessions"].to_csv(fname, index=False) bst = pd.read_csv(fname) - c['ecephys_sessions'].to_csv(fname, index=False) + c["ecephys_sessions"].to_csv(fname, index=False) est = pd.read_csv(fname) - c['probes'].to_csv(fname, index=False) + c["probes"].to_csv(fname, index=False) pt = pd.read_csv(fname) - c['channels'].to_csv(fname, index=False) + c["channels"].to_csv(fname, index=False) ct = pd.read_csv(fname) - c['units'].to_csv(fname, index=False) + c["units"].to_csv(fname, index=False) ut = pd.read_csv(fname) yield (MockCache(bst, est, pt, ct, ut, tmpdir), c) @@ -133,22 +148,17 @@ def mock_cache(tmpdir): @pytest.mark.parametrize("local", [True, False]) def test_VisualBehaviorNeuropixelsProjectCloudApi( - mock_cache, - monkeypatch, - local + mock_cache, monkeypatch, local ): - mocked_cache, expected = mock_cache api = cloudapi.VisualBehaviorNeuropixelsProjectCloudApi( - mocked_cache, - skip_version_check=True, - local=False) + mocked_cache, skip_version_check=True, local=False + ) if local: api = cloudapi.VisualBehaviorNeuropixelsProjectCloudApi( - mocked_cache, - skip_version_check=True, - local=True) + mocked_cache, skip_version_check=True, local=True + ) # behavior session table as expected bst = api.get_behavior_session_table() @@ -187,63 +197,68 @@ def test_VisualBehaviorNeuropixelsProjectCloudApi( def mock_nwb(nwb_path, probe_meta): return nwb_path - monkeypatch.setattr(cloudapi.BehaviorEcephysSession, - "from_nwb_path", mock_nwb) + monkeypatch.setattr( + cloudapi.BehaviorEcephysSession, "from_nwb_path", mock_nwb + ) assert api.get_ecephys_session(12) == "12" -@pytest.mark.parametrize('has_lfp_data', [True, False]) +@pytest.mark.parametrize("has_lfp_data", [True, False]) def test_probe_meta(mock_cache, monkeypatch, has_lfp_data): """tests that probe_meta is set correctly""" mocked_cache, _ = mock_cache probe_meta_table = pd.read_csv(mocked_cache.probe_table_path) - probe_meta_table['has_lfp_data'] = has_lfp_data + probe_meta_table["has_lfp_data"] = has_lfp_data probe_meta_table.to_csv(mocked_cache.probe_table_path, index=False) api = cloudapi.VisualBehaviorNeuropixelsProjectCloudApi( - mocked_cache, - skip_version_check=True, - local=False) + mocked_cache, skip_version_check=True, local=False + ) def mock_from_nwb_path(nwb_path, probe_meta): return probe_meta ecephys_session_id = 10 - with patch.object(BehaviorEcephysSession, 'from_nwb_path', - wraps=mock_from_nwb_path): + with patch.object( + BehaviorEcephysSession, "from_nwb_path", wraps=mock_from_nwb_path + ): probe_meta = api.get_ecephys_session(ecephys_session_id) probe_meta_table = api.get_probe_table() probes = probe_meta_table.loc[ - probe_meta_table['ecephys_session_id'] == ecephys_session_id] + probe_meta_table["ecephys_session_id"] == ecephys_session_id + ] if has_lfp_data: - for probe_name in probes['name'].unique(): + for probe_name in probes["name"].unique(): obtained_probe_nwb_path = probe_meta[probe_name].lfp_csd_filepath() - expected_probe_nwb_path = str(probe_meta_table.loc[ - (probe_meta_table['name'] == probe_name) - ].iloc[0]['file_id']) + expected_probe_nwb_path = str( + probe_meta_table.loc[ + (probe_meta_table["name"] == probe_name) + ].iloc[0]["file_id"] + ) assert obtained_probe_nwb_path == expected_probe_nwb_path else: assert probe_meta is None @pytest.mark.parametrize( - "manifest_version, data_pipeline_version, cmin, cmax, exception", - [ - ("0.0.1", "2.9.0", "0.0.0", "1.0.0", False), - ("1.0.1", "2.9.0", "0.0.0", "1.0.0", True) - ]) -def test_version_check(manifest_version, data_pipeline_version, - cmin, cmax, exception): + "manifest_version, data_pipeline_version, cmin, cmax, exception", + [ + ("0.0.1", "2.9.0", "0.0.0", "1.0.0", False), + ("1.0.1", "2.9.0", "0.0.0", "1.0.0", True), + ], +) +def test_version_check( + manifest_version, data_pipeline_version, cmin, cmax, exception +): if exception: - with pytest.raises(BehaviorCloudCacheVersionException, - match=f".*{data_pipeline_version}"): - version_check( - manifest_version, - data_pipeline_version, - cmin, cmax) + with pytest.raises( + BehaviorCloudCacheVersionException, + match=f".*{data_pipeline_version}", + ): + version_check(manifest_version, data_pipeline_version, cmin, cmax) else: version_check(manifest_version, data_pipeline_version, cmin, cmax) @@ -251,21 +266,22 @@ def test_version_check(manifest_version, data_pipeline_version, def test_from_local_cache(monkeypatch): mock_manifest = create_autospec(Manifest) mock_manifest.metadata_file_names = { - 'ecephys_sessions', - 'behavior_sessions', - 'probes', - 'units', - 'channels' + "ecephys_sessions", + "behavior_sessions", + "probes", + "units", + "channels", } mock_manifest._data_pipeline = [ { "name": "AllenSDK", "version": "2.11.0", - "comment": "This is a test entry. NOT REAL." + "comment": "This is a test entry. NOT REAL.", } ] - mock_manifest.version = cloudapi \ - .VisualBehaviorNeuropixelsProjectCloudApi.MANIFEST_COMPATIBILITY[0] + mock_manifest.version = cloudapi.VisualBehaviorNeuropixelsProjectCloudApi.MANIFEST_COMPATIBILITY[ # noqa: E501 + 0 + ] mock_local_cache = create_autospec(cloudapibase.LocalCache) type(mock_local_cache.return_value)._manifest = mock_manifest @@ -273,7 +289,6 @@ def test_from_local_cache(monkeypatch): type(mock_static_local_cache.return_value)._manifest = mock_manifest with monkeypatch.context() as m: - m.setattr(cloudapibase, "LocalCache", mock_local_cache) m.setattr(cloudapibase, "StaticLocalCache", mock_static_local_cache) diff --git a/allensdk/test/brain_observatory/ecephys/stimulus_analysis/test_stimulus_analysis.py b/allensdk/test/brain_observatory/ecephys/stimulus_analysis/test_stimulus_analysis.py index 0cab53812..2cf28cd79 100644 --- a/allensdk/test/brain_observatory/ecephys/stimulus_analysis/test_stimulus_analysis.py +++ b/allensdk/test/brain_observatory/ecephys/stimulus_analysis/test_stimulus_analysis.py @@ -1,23 +1,32 @@ -import pytest -import pandas as pd import numpy as np +import pandas as pd +import pytest import xarray as xr -import warnings - -from allensdk.brain_observatory.ecephys.ecephys_session_api import EcephysSessionApi from allensdk.brain_observatory.ecephys.ecephys_session import EcephysSession -from allensdk.brain_observatory.ecephys.stimulus_analysis.stimulus_analysis import StimulusAnalysis, \ - running_modulation, lifetime_sparseness, fano_factor, overall_firing_rate, get_fr, osi, dsi - - -pd.set_option('display.max_columns', None) +from allensdk.brain_observatory.ecephys.ecephys_session_api import ( + EcephysSessionApi, +) +from allensdk.brain_observatory.ecephys.stimulus_analysis.stimulus_analysis import ( # noqa: E501 + StimulusAnalysis, + dsi, + fano_factor, + get_fr, + lifetime_sparseness, + osi, + overall_firing_rate, + running_modulation, +) + +pd.set_option("display.max_columns", None) class MockSessionApi(EcephysSessionApi): - """Mock Data to create an EcephysSession object and pass it into stimulus analysis + """Mock Data to create an EcephysSession object and pass it into stimulus + analysis # TODO: move to conftest so other tests can use data """ + def get_spike_times(self): return { 0: np.array([1, 2, 3, 4]), @@ -25,55 +34,83 @@ def get_spike_times(self): 2: np.array([1.01, 1.03, 1.02]), 3: np.array([]), 4: np.array([0.01, 1.7, 2.13, 3.19, 4.25]), - 5: np.array([1.5, 3.0, 4.5]) + 5: np.array([1.5, 3.0, 4.5]), } def get_channels(self): - return pd.DataFrame({ - 'local_index': [0, 1, 2], - 'probe_horizontal_position': [5, 10, 15], - 'probe_id': [0, 0, 1], - 'probe_vertical_position': [10, 22, 33], - 'valid_data': [False, True, True] - }, index=pd.Index(name='channel_id', data=[0, 1, 2])) + return pd.DataFrame( + { + "local_index": [0, 1, 2], + "probe_horizontal_position": [5, 10, 15], + "probe_id": [0, 0, 1], + "probe_vertical_position": [10, 22, 33], + "valid_data": [False, True, True], + }, + index=pd.Index(name="channel_id", data=[0, 1, 2]), + ) def get_units(self): - udf = pd.DataFrame({ - 'firing_rate': np.linspace(1, 3, 6), - 'isi_violations': [40, 0.5, 0.1, 0.2, 0.0, 0.1], - 'local_index': [0, 0, 1, 1, 2, 2], - 'peak_channel_id': [0, 2, 1, 1, 2, 0], - 'quality': ['good', 'good', 'good', 'bad', 'good', 'good'], - }, index=pd.Index(name='unit_id', data=np.arange(6)[::-1])) + udf = pd.DataFrame( + { + "firing_rate": np.linspace(1, 3, 6), + "isi_violations": [40, 0.5, 0.1, 0.2, 0.0, 0.1], + "local_index": [0, 0, 1, 1, 2, 2], + "peak_channel_id": [0, 2, 1, 1, 2, 0], + "quality": ["good", "good", "good", "bad", "good", "good"], + }, + index=pd.Index(name="unit_id", data=np.arange(6)[::-1]), + ) return udf def get_probes(self): - return pd.DataFrame({ - 'description': ['probeA', 'probeB'], - 'location': ['VISp', 'VISam'], - 'sampling_rate': [30000.0, 30000.0] - }, index=pd.Index(name='id', data=[0, 1])) + return pd.DataFrame( + { + "description": ["probeA", "probeB"], + "location": ["VISp", "VISam"], + "sampling_rate": [30000.0, 30000.0], + }, + index=pd.Index(name="id", data=[0, 1]), + ) def get_stimulus_presentations(self): - return pd.DataFrame({ - 'start_time': np.linspace(0.0, 4.5, 10, endpoint=True), - 'stop_time': np.linspace(0.5, 5.0, 10, endpoint=True), - 'stimulus_name': ['spontaneous'] + ['s0'] * 6 + ['spontaneous'] + ['s1'] * 2, - 'stimulus_block': [0] + [1] * 6 + [0] + [2] * 2, - 'duration': 0.5, - 'stimulus_index': [0] + [1] * 6 + [0] + [2] * 2, - 'conditions': [0, 0, 0, 0, 1, 1, 1, 0, 2, 3] # generic stimulus condition - }, index=pd.Index(name='id', data=np.arange(10))) + return pd.DataFrame( + { + "start_time": np.linspace(0.0, 4.5, 10, endpoint=True), + "stop_time": np.linspace(0.5, 5.0, 10, endpoint=True), + "stimulus_name": ["spontaneous"] + + ["s0"] * 6 + + ["spontaneous"] + + ["s1"] * 2, + "stimulus_block": [0] + [1] * 6 + [0] + [2] * 2, + "duration": 0.5, + "stimulus_index": [0] + [1] * 6 + [0] + [2] * 2, + "conditions": [ + 0, + 0, + 0, + 0, + 1, + 1, + 1, + 0, + 2, + 3, + ], # generic stimulus condition + }, + index=pd.Index(name="id", data=np.arange(10)), + ) def get_invalid_times(self): return pd.DataFrame() def get_running_speed(self): - return pd.DataFrame({ - "start_time": np.linspace(0.0, 9.9, 100), - "end_time": np.linspace(0.1, 10.0, 100), - "velocity": np.linspace(-0.1, 11.0, 100) - }) + return pd.DataFrame( + { + "start_time": np.linspace(0.0, 9.9, 100), + "end_time": np.linspace(0.1, 10.0, 100), + "velocity": np.linspace(-0.1, 11.0, 100), + } + ) @pytest.fixture @@ -84,277 +121,572 @@ def ecephys_api(): def test_unit_ids(ecephys_api): session = EcephysSession(api=ecephys_api) stim_analysis = StimulusAnalysis(ecephys_session=session) - assert(set(stim_analysis.unit_ids) == set(range(6))) - assert(stim_analysis.unit_count == 6) + assert set(stim_analysis.unit_ids) == set(range(6)) + assert stim_analysis.unit_count == 6 def test_unit_ids_filter_by_id(ecephys_api): session = EcephysSession(api=ecephys_api) stim_analysis = StimulusAnalysis(ecephys_session=session, filter=[2, 3, 1]) - assert(set(stim_analysis.unit_ids) == {1, 2, 3}) - assert(stim_analysis.unit_count == 3) + assert set(stim_analysis.unit_ids) == {1, 2, 3} + assert stim_analysis.unit_count == 3 - stim_analysis = StimulusAnalysis(ecephys_session=session, filter={'unit_id': [3, 0]}) - assert(set(stim_analysis.unit_ids) == {0, 3}) - assert(stim_analysis.unit_count == 2) + stim_analysis = StimulusAnalysis( + ecephys_session=session, filter={"unit_id": [3, 0]} + ) + assert set(stim_analysis.unit_ids) == {0, 3} + assert stim_analysis.unit_count == 2 with pytest.raises(KeyError): # If unit ids don't exists should raise an error - stim_analysis = StimulusAnalysis(ecephys_session=session, filter=[100, 200]) - units = stim_analysis.unit_ids + stim_analysis = StimulusAnalysis( + ecephys_session=session, filter=[100, 200] + ) + stim_analysis.unit_ids def test_unit_ids_filtered(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, filter={'location': 'VISp'}) - assert(set(stim_analysis.unit_ids) == {0, 2, 3, 5}) - assert(stim_analysis.unit_count == 4) + stim_analysis = StimulusAnalysis( + ecephys_session=session, filter={"location": "VISp"} + ) + assert set(stim_analysis.unit_ids) == {0, 2, 3, 5} + assert stim_analysis.unit_count == 4 - stim_analysis = StimulusAnalysis(ecephys_session=session, filter={'location': 'VISp', 'quality': 'good'}) - assert(set(stim_analysis.unit_ids) == {0, 3, 5}) - assert(stim_analysis.unit_count == 3) + stim_analysis = StimulusAnalysis( + ecephys_session=session, filter={"location": "VISp", "quality": "good"} + ) + assert set(stim_analysis.unit_ids) == {0, 3, 5} + assert stim_analysis.unit_count == 3 with pytest.raises(Exception): # No units found should raise exception - stim_analysis = StimulusAnalysis(ecephys_session=session, filter={'location': 'pSIV'}) + stim_analysis = StimulusAnalysis( + ecephys_session=session, filter={"location": "pSIV"} + ) stim_analysis.unit_ids stim_analysis.unit_count def test_stim_table(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') - assert(isinstance(stim_analysis.stim_table, pd.DataFrame)) - assert(len(stim_analysis.stim_table) == 6) - assert(stim_analysis.total_presentations == 6) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) + assert isinstance(stim_analysis.stim_table, pd.DataFrame) + assert len(stim_analysis.stim_table) == 6 + assert stim_analysis.total_presentations == 6 # Make sure certain columns exist - assert('start_time' in stim_analysis.stim_table) - assert('stop_time' in stim_analysis.stim_table) - assert('stimulus_condition_id' in stim_analysis.stim_table) - assert('stimulus_name' in stim_analysis.stim_table) - assert('duration' in stim_analysis.stim_table) + assert "start_time" in stim_analysis.stim_table + assert "stop_time" in stim_analysis.stim_table + assert "stimulus_condition_id" in stim_analysis.stim_table + assert "stimulus_name" in stim_analysis.stim_table + assert "duration" in stim_analysis.stim_table with pytest.raises(Exception): - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='0s') + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="0s" + ) stim_analysis.stim_table def test_stim_table_spontaneous(ecephys_api): - # By default table should be empty because non of the stimulus are above the duration threshold + # By default table should be empty because non of the stimulus are above + # the duration threshold session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, spontaneous_threshold=0.49) - assert(isinstance(stim_analysis.stim_table_spontaneous, pd.DataFrame)) - assert(len(stim_analysis.stim_table_spontaneous) == 2) + stim_analysis = StimulusAnalysis( + ecephys_session=session, spontaneous_threshold=0.49 + ) + assert isinstance(stim_analysis.stim_table_spontaneous, pd.DataFrame) + assert len(stim_analysis.stim_table_spontaneous) == 2 # Check that threshold is working - stim_analysis = StimulusAnalysis(ecephys_session=session, spontaneous_threshold=0.51) - assert(len(stim_analysis.stim_table_spontaneous) == 0) + stim_analysis = StimulusAnalysis( + ecephys_session=session, spontaneous_threshold=0.51 + ) + assert len(stim_analysis.stim_table_spontaneous) == 0 def test_conditionwise_psth(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0', trial_duration=0.5, - psth_resolution=0.1) - assert(isinstance(stim_analysis.conditionwise_psth, xr.DataArray)) + stim_analysis = StimulusAnalysis( + ecephys_session=session, + stimulus_key="s0", + trial_duration=0.5, + psth_resolution=0.1, + ) + assert isinstance(stim_analysis.conditionwise_psth, xr.DataArray) # assert(stim_analysis.conditionwise_psth.shape == (2, 4, 6)) - assert(stim_analysis.conditionwise_psth.coords['time_relative_to_stimulus_onset'].size == 4) # 0.5/0.1 - 1 - assert(stim_analysis.conditionwise_psth.coords['unit_id'].size == 6) - assert(stim_analysis.conditionwise_psth.coords['stimulus_condition_id'].size == 2) - assert(np.allclose(stim_analysis.conditionwise_psth[{'unit_id': 0, 'stimulus_condition_id': 1}].values, - np.array([1.0/3.0, 0.0, 0.0, 0.0]))) + assert ( + stim_analysis.conditionwise_psth.coords[ + "time_relative_to_stimulus_onset" + ].size + == 4 + ) # 0.5/0.1 - 1 + assert stim_analysis.conditionwise_psth.coords["unit_id"].size == 6 + assert ( + stim_analysis.conditionwise_psth.coords["stimulus_condition_id"].size + == 2 + ) + assert np.allclose( + stim_analysis.conditionwise_psth[ + {"unit_id": 0, "stimulus_condition_id": 1} + ].values, + np.array([1.0 / 3.0, 0.0, 0.0, 0.0]), + ) # Make sure psth doesn't fail even when all the condition_ids are unique. - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s1', trial_duration=0.5, - psth_resolution=0.1) - assert(stim_analysis.conditionwise_psth.coords['time_relative_to_stimulus_onset'].size == 4) - assert(stim_analysis.conditionwise_psth.coords['unit_id'].size == 6) - assert(stim_analysis.conditionwise_psth.coords['stimulus_condition_id'].size == 2) + stim_analysis = StimulusAnalysis( + ecephys_session=session, + stimulus_key="s1", + trial_duration=0.5, + psth_resolution=0.1, + ) + assert ( + stim_analysis.conditionwise_psth.coords[ + "time_relative_to_stimulus_onset" + ].size + == 4 + ) + assert stim_analysis.conditionwise_psth.coords["unit_id"].size == 6 + assert ( + stim_analysis.conditionwise_psth.coords["stimulus_condition_id"].size + == 2 + ) def test_conditionwise_statistics(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') - assert(len(stim_analysis.conditionwise_statistics) == 2*6) # units x condition_ids - assert(set(stim_analysis.conditionwise_statistics.index.names) == {'unit_id', 'stimulus_condition_id'}) - assert(set(stim_analysis.conditionwise_statistics.columns) == - {'spike_std', 'spike_sem', 'spike_count', 'stimulus_presentation_count', 'spike_mean'}) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) + assert ( + len(stim_analysis.conditionwise_statistics) == 2 * 6 + ) # units x condition_ids + assert set(stim_analysis.conditionwise_statistics.index.names) == { + "unit_id", + "stimulus_condition_id", + } + assert set(stim_analysis.conditionwise_statistics.columns) == { + "spike_std", + "spike_sem", + "spike_count", + "stimulus_presentation_count", + "spike_mean", + } expected = pd.Series( - [2.0, 3.0, 0.66666667, 0.57735027, 0.33333333], - ["spike_count", "stimulus_presentation_count", "spike_mean", "spike_std", "spike_sem"] + [2.0, 3.0, 0.66666667, 0.57735027, 0.33333333], + [ + "spike_count", + "stimulus_presentation_count", + "spike_mean", + "spike_std", + "spike_sem", + ], ) obtained = stim_analysis.conditionwise_statistics.loc[(0, 1)] - pd.testing.assert_series_equal(expected, obtained[expected.index], check_dtype=True, check_names=False, rtol=1e-5, atol=1e-8) + pd.testing.assert_series_equal( + expected, + obtained[expected.index], + check_dtype=True, + check_names=False, + rtol=1e-5, + atol=1e-8, + ) def test_presentationwise_spike_times(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') - assert(len(stim_analysis.presentationwise_spike_times) == 12) - assert(list(stim_analysis.presentationwise_spike_times.index.names) == ['spike_time']) - assert(set(stim_analysis.presentationwise_spike_times.columns) == {'stimulus_presentation_id', 'unit_id', 'time_since_stimulus_presentation_onset'}) - assert(stim_analysis.presentationwise_spike_times.loc[1.01]['unit_id'] == 2) - assert(stim_analysis.presentationwise_spike_times.loc[1.01]['stimulus_presentation_id'] == 2) - assert(len(stim_analysis.presentationwise_spike_times.loc[3.0]) == 2) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) + assert len(stim_analysis.presentationwise_spike_times) == 12 + assert list(stim_analysis.presentationwise_spike_times.index.names) == [ + "spike_time" + ] + assert set(stim_analysis.presentationwise_spike_times.columns) == { + "stimulus_presentation_id", + "unit_id", + "time_since_stimulus_presentation_onset", + } + assert stim_analysis.presentationwise_spike_times.loc[1.01]["unit_id"] == 2 + assert ( + stim_analysis.presentationwise_spike_times.loc[1.01][ + "stimulus_presentation_id" + ] + == 2 + ) + assert len(stim_analysis.presentationwise_spike_times.loc[3.0]) == 2 def test_presentationwise_statistics(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0', trial_duration=0.5) - assert(len(stim_analysis.presentationwise_statistics) == 6*6) # units x presentation_ids - assert(set(stim_analysis.presentationwise_statistics.index.names) == {'stimulus_presentation_id', 'unit_id'}) - assert(set(stim_analysis.presentationwise_statistics.columns) == {'spike_counts', 'stimulus_condition_id', - 'running_speed'}) - assert(stim_analysis.presentationwise_statistics.loc[1, 0]['spike_counts'] == 1.0) - assert(stim_analysis.presentationwise_statistics.loc[1, 0]['stimulus_condition_id'] == 1.0) - assert(np.isclose(stim_analysis.presentationwise_statistics.loc[1, 0]['running_speed'], 0.684848)) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0", trial_duration=0.5 + ) + assert ( + len(stim_analysis.presentationwise_statistics) == 6 * 6 + ) # units x presentation_ids + assert set(stim_analysis.presentationwise_statistics.index.names) == { + "stimulus_presentation_id", + "unit_id", + } + assert set(stim_analysis.presentationwise_statistics.columns) == { + "spike_counts", + "stimulus_condition_id", + "running_speed", + } + assert ( + stim_analysis.presentationwise_statistics.loc[1, 0]["spike_counts"] + == 1.0 + ) + assert ( + stim_analysis.presentationwise_statistics.loc[1, 0][ + "stimulus_condition_id" + ] + == 1.0 + ) + assert np.isclose( + stim_analysis.presentationwise_statistics.loc[1, 0]["running_speed"], + 0.684848, + ) def test_stimulus_conditions(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0', trial_duration=0.5) - assert(len(stim_analysis.stimulus_conditions) == 2) - assert(np.all(stim_analysis.stimulus_conditions['stimulus_name'].unique() == ['s0'])) - assert(set(stim_analysis.stimulus_conditions['conditions'].unique()) == {0, 1}) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0", trial_duration=0.5 + ) + assert len(stim_analysis.stimulus_conditions) == 2 + assert np.all( + stim_analysis.stimulus_conditions["stimulus_name"].unique() == ["s0"] + ) + assert set(stim_analysis.stimulus_conditions["conditions"].unique()) == { + 0, + 1, + } def test_running_speed(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') - assert(set(stim_analysis.running_speed.index.values) == set(range(1, 7))) - assert(np.isclose(stim_analysis.running_speed.loc[1]['running_speed'], 0.684848)) - assert(np.isclose(stim_analysis.running_speed.loc[3]['running_speed'], 1.806061)) - assert(np.isclose(stim_analysis.running_speed.loc[6]['running_speed'], 3.487879)) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) + assert set(stim_analysis.running_speed.index.values) == set(range(1, 7)) + assert np.isclose( + stim_analysis.running_speed.loc[1]["running_speed"], 0.684848 + ) + assert np.isclose( + stim_analysis.running_speed.loc[3]["running_speed"], 1.806061 + ) + assert np.isclose( + stim_analysis.running_speed.loc[6]["running_speed"], 3.487879 + ) def test_spikes(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') - assert(isinstance(stim_analysis.spikes, dict)) - assert(stim_analysis.spikes.keys() == set(range(6))) - assert(np.allclose(stim_analysis.spikes[0], [1, 2, 3, 4])) - assert(np.allclose(stim_analysis.spikes[4], [0.01, 1.7 , 2.13, 3.19, 4.25])) - assert(stim_analysis.spikes[3].size == 0) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) + assert isinstance(stim_analysis.spikes, dict) + assert stim_analysis.spikes.keys() == set(range(6)) + assert np.allclose(stim_analysis.spikes[0], [1, 2, 3, 4]) + assert np.allclose(stim_analysis.spikes[4], [0.01, 1.7, 2.13, 3.19, 4.25]) + assert stim_analysis.spikes[3].size == 0 # Check that spikes dict is filtering units session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0', filter=[0, 2]) - assert(stim_analysis.spikes.keys() == {0, 2}) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0", filter=[0, 2] + ) + assert stim_analysis.spikes.keys() == {0, 2} def test_get_preferred_condition(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') - assert(stim_analysis._get_preferred_condition(3) == 1) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) + assert stim_analysis._get_preferred_condition(3) == 1 with pytest.raises(KeyError): stim_analysis._get_preferred_condition(10) + def test_check_multiple_preferred_conditions(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0') + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0" + ) - assert(stim_analysis._check_multiple_pref_conditions(0, 'conditions', [0, 1]) is False) - assert(stim_analysis._check_multiple_pref_conditions(3, 'conditions', [0, 1]) is True) + assert ( + stim_analysis._check_multiple_pref_conditions(0, "conditions", [0, 1]) + is False + ) + assert ( + stim_analysis._check_multiple_pref_conditions(3, "conditions", [0, 1]) + is True + ) def test_get_time_to_peak(ecephys_api): session = EcephysSession(api=ecephys_api) - stim_analysis = StimulusAnalysis(ecephys_session=session, stimulus_key='s0', trial_duration=0.5) - assert(stim_analysis._get_time_to_peak(1, stim_analysis._get_preferred_condition(1)) == 0.0005) - - -@pytest.mark.parametrize('spike_counts,running_speeds,speed_threshold,expected', - [ - (np.zeros(10), np.zeros(1),1.0, (np.nan, np.nan)), # Input error, return nan - (np.zeros(5), np.full(5, 2.0), 1.0, (np.nan, np.nan)), # returns Nan, always running - (np.zeros(5), np.full(5, 2.0), 2.1, (np.nan, np.nan)), # returns Nan, always stationary - (np.zeros(5), np.array([0.0, 0.0, 2.0, 2.5, 3.0]), 1.0, (np.nan, np.nan)), # No firing, return Nans - (np.ones(5), np.array([0.0, 0.0, 2.0, 2.0, 2.0]), 1.0, (np.nan, 0.0)), # always the same fr, pval is Nan but run_mod is 0.0) - (np.array([3.0, 3.0, 1.5, 1.5, 0.9]), np.array([0.0, 0.0, 2.0, 2.5, 3.0]), 1.0, (0.013559949584378913, -0.5666666666666667)), - (np.array([3.0, 3.0, 1.5, 1.5, 1.5]), np.array([0.0, 0.0, 2.0, 2.5, 3.0]), 1.0, (0.0, -0.5)), - (np.array([3.0, 3.0, 1.5, 5.5, 2.5]), np.array([0.0, 0.0, 2.0, 2.5, 3.0]), 1.0, (0.9024099927051468, 0.052631578947368376)), - (np.array([0.0, 0.0, 4.0, 4.0]), np.array([0.0, 0.0, 2.5, 3.0]), 1.0, (0.0, 1.0)) - ]) -def test_running_modulation(spike_counts, running_speeds, speed_threshold, expected): - rm = running_modulation(spike_counts, running_speeds, speed_threshold) - assert(np.allclose(rm, expected, equal_nan=True)) + stim_analysis = StimulusAnalysis( + ecephys_session=session, stimulus_key="s0", trial_duration=0.5 + ) + assert ( + stim_analysis._get_time_to_peak( + 1, stim_analysis._get_preferred_condition(1) + ) + == 0.0005 + ) -@pytest.mark.parametrize('responses,expected', - [ - (np.array([1.0]), np.nan), # can't calculate for single point - (np.full(20, 3.2), 0.0), # always the same, sparness should be at/near 0 - (np.array([10.0, 0.0, 0.0, 0.0, 0.0, 0.0]), 1.0), # spareness should be close to 1 - (np.array([2.24, 3.6, 0.8, 2.4, 3.52, 5.68, 8.96, 0.0]), 0.43500091849856115) - ]) +@pytest.mark.parametrize( + "spike_counts,running_speeds,speed_threshold,expected", + [ + ( + np.zeros(10), + np.zeros(1), + 1.0, + (np.nan, np.nan), + ), # Input error, return nan + ( + np.zeros(5), + np.full(5, 2.0), + 1.0, + (np.nan, np.nan), + ), # returns Nan, always running + ( + np.zeros(5), + np.full(5, 2.0), + 2.1, + (np.nan, np.nan), + ), # returns Nan, always stationary + ( + np.zeros(5), + np.array([0.0, 0.0, 2.0, 2.5, 3.0]), + 1.0, + (np.nan, np.nan), + ), # No firing, return Nans + ( + np.ones(5), + np.array([0.0, 0.0, 2.0, 2.0, 2.0]), + 1.0, + (np.nan, 0.0), + ), # always the same fr, pval is Nan but run_mod is 0.0) + ( + np.array([3.0, 3.0, 1.5, 1.5, 0.9]), + np.array([0.0, 0.0, 2.0, 2.5, 3.0]), + 1.0, + (0.013559949584378913, -0.5666666666666667), + ), + ( + np.array([3.0, 3.0, 1.5, 1.5, 1.5]), + np.array([0.0, 0.0, 2.0, 2.5, 3.0]), + 1.0, + (0.0, -0.5), + ), + ( + np.array([3.0, 3.0, 1.5, 5.5, 2.5]), + np.array([0.0, 0.0, 2.0, 2.5, 3.0]), + 1.0, + (0.9024099927051468, 0.052631578947368376), + ), + ( + np.array([0.0, 0.0, 4.0, 4.0]), + np.array([0.0, 0.0, 2.5, 3.0]), + 1.0, + (0.0, 1.0), + ), + ], +) +def test_running_modulation( + spike_counts, running_speeds, speed_threshold, expected +): + rm = running_modulation(spike_counts, running_speeds, speed_threshold) + assert np.allclose(rm, expected, equal_nan=True) + + +@pytest.mark.parametrize( + "responses,expected", + [ + (np.array([1.0]), np.nan), # can't calculate for single point + ( + np.full(20, 3.2), + 0.0, + ), # always the same, sparness should be at/near 0 + ( + np.array([10.0, 0.0, 0.0, 0.0, 0.0, 0.0]), + 1.0, + ), # spareness should be close to 1 + ( + np.array([2.24, 3.6, 0.8, 2.4, 3.52, 5.68, 8.96, 0.0]), + 0.43500091849856115, + ), + ], +) def test_lifetime_sparseness(responses, expected): lts = lifetime_sparseness(responses) - assert(np.isclose(lts, expected, equal_nan=True)) - - -@pytest.mark.parametrize('spike_counts,expected', - [ - (np.zeros(10), np.nan), # mean 0.0 leads to Nan - (np.array([-1.5, 1.5]), np.nan), # mean 0 - (np.ones(5), 0.0), # no variance - (np.array([1.2, 20.0, 0.0, 36.2, 0.6]), 17.921379310344832), # High variance - (np.array([5.1, 5.3, 5.2, 5.1, 5.2]), 0.0010810810810810846), # low variance - ]) + assert np.isclose(lts, expected, equal_nan=True) + + +@pytest.mark.parametrize( + "spike_counts,expected", + [ + (np.zeros(10), np.nan), # mean 0.0 leads to Nan + (np.array([-1.5, 1.5]), np.nan), # mean 0 + (np.ones(5), 0.0), # no variance + ( + np.array([1.2, 20.0, 0.0, 36.2, 0.6]), + 17.921379310344832, + ), # High variance + ( + np.array([5.1, 5.3, 5.2, 5.1, 5.2]), + 0.0010810810810810846, + ), # low variance + ], +) def test_fano_factor(spike_counts, expected): ff = fano_factor(spike_counts) - assert(np.isclose(ff, expected, equal_nan=True)) - - -@pytest.mark.parametrize('start_times,stop_times,spike_times,expected', - [ - (np.array([0.0]), np.array([0.0]), np.linspace(0, 10.0, 10), np.nan), # nan, total time 0.0 - (np.arange(1.0, 3.0), np.arange(0.0, 2.0), np.linspace(0, 10.0, 10), np.nan), # nan, total_time negative - (np.arange(1.0, 4.0), np.arange(0.0, 2.0), np.linspace(0, 10.0, 10), np.nan), # nan, time lengths don't match - (np.array([0.0]), np.array([1.0]), np.linspace(0, 10.0, 100), 10.0), - (np.array([0.0, 9.0]), np.array([1.0, 10.0]), np.linspace(0, 10.0, 101), 10.0) # 10.0 Hz split up into blocks - ]) + assert np.isclose(ff, expected, equal_nan=True) + + +@pytest.mark.parametrize( + "start_times,stop_times,spike_times,expected", + [ + ( + np.array([0.0]), + np.array([0.0]), + np.linspace(0, 10.0, 10), + np.nan, + ), # nan, total time 0.0 + ( + np.arange(1.0, 3.0), + np.arange(0.0, 2.0), + np.linspace(0, 10.0, 10), + np.nan, + ), # nan, total_time negative + ( + np.arange(1.0, 4.0), + np.arange(0.0, 2.0), + np.linspace(0, 10.0, 10), + np.nan, + ), # nan, time lengths don't match + (np.array([0.0]), np.array([1.0]), np.linspace(0, 10.0, 100), 10.0), + ( + np.array([0.0, 9.0]), + np.array([1.0, 10.0]), + np.linspace(0, 10.0, 101), + 10.0, + ), # 10.0 Hz split up into blocks + ], +) def test_overall_firing_rate(start_times, stop_times, spike_times, expected): ofr = overall_firing_rate(start_times, stop_times, spike_times) - assert(np.isclose(ofr, expected, equal_nan=True)) - - -@pytest.mark.parametrize('spikes,sampling_freq,sweep_length,expected', - [ - (np.array([0.82764702, 0.83624702, 1.09211374]), 10, 1.5, [0.0, 0.0, 0.0, 0.0, 0.000133830, 0.004431861, 0.05412495, 0.2464033072, 0.45293459677, 0.4839428913, - 0.452934596, 0.2464033072, 0.054124958, 0.00443186162, 0.0001338306]), - (np.array([]), 10, 1.5, np.zeros(15)) - ]) + assert np.isclose(ofr, expected, equal_nan=True) + + +@pytest.mark.parametrize( + "spikes,sampling_freq,sweep_length,expected", + [ + ( + np.array([0.82764702, 0.83624702, 1.09211374]), + 10, + 1.5, + [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.000133830, + 0.004431861, + 0.05412495, + 0.2464033072, + 0.45293459677, + 0.4839428913, + 0.452934596, + 0.2464033072, + 0.054124958, + 0.00443186162, + 0.0001338306, + ], + ), + (np.array([]), 10, 1.5, np.zeros(15)), + ], +) def test_get_fr(spikes, sampling_freq, sweep_length, expected): - frs = get_fr(spikes, num_timestep_second=sampling_freq, sweep_length=sweep_length) - assert(len(frs) == int(sampling_freq*sweep_length)) - assert(np.allclose(frs, expected)) - - -@pytest.mark.parametrize('orivals,tuning,expected', - [ - (np.array([1.0]), np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]), np.nan), - (np.array([]), np.array([]), np.nan), - (np.array([0.0+0.0j, 0.52359878 + 0.j, 1.04719755 + 0.j, 1.57079633+0.j, 2.0943951+0.j, 2.61799388+0.j]), np.array([5.5, 4.44, 3.54166667, 4.10869565, 4.42, 4.55319149]), 0.07873455094232604) - ]) + frs = get_fr( + spikes, num_timestep_second=sampling_freq, sweep_length=sweep_length + ) + assert len(frs) == int(sampling_freq * sweep_length) + assert np.allclose(frs, expected) + + +@pytest.mark.parametrize( + "orivals,tuning,expected", + [ + ( + np.array([1.0]), + np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]), + np.nan, + ), + (np.array([]), np.array([]), np.nan), + ( + np.array( + [ + 0.0 + 0.0j, + 0.52359878 + 0.0j, + 1.04719755 + 0.0j, + 1.57079633 + 0.0j, + 2.0943951 + 0.0j, + 2.61799388 + 0.0j, + ] + ), + np.array([5.5, 4.44, 3.54166667, 4.10869565, 4.42, 4.55319149]), + 0.07873455094232604, + ), + ], +) def test_osi(orivals, tuning, expected): osi_val = osi(orivals, tuning) - assert(np.allclose(osi_val, expected, equal_nan=True)) - - -@pytest.mark.parametrize('orivals,tuning,expected', - [ - (np.array([1.0]), np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]), np.nan), - (np.array([]), np.array([]), np.nan), - (np.array([0.0+0.0j, 0.52359878 + 0.j, 1.04719755 + 0.j, 1.57079633+0.j, 2.0943951+0.j, 2.61799388+0.j]), np.array([5.5, 4.44, 3.54166667, 4.10869565, 4.42, 4.55319149]), 0.6126966469601506) - ]) + assert np.allclose(osi_val, expected, equal_nan=True) + + +@pytest.mark.parametrize( + "orivals,tuning,expected", + [ + ( + np.array([1.0]), + np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]), + np.nan, + ), + (np.array([]), np.array([]), np.nan), + ( + np.array( + [ + 0.0 + 0.0j, + 0.52359878 + 0.0j, + 1.04719755 + 0.0j, + 1.57079633 + 0.0j, + 2.0943951 + 0.0j, + 2.61799388 + 0.0j, + ] + ), + np.array([5.5, 4.44, 3.54166667, 4.10869565, 4.42, 4.55319149]), + 0.6126966469601506, + ), + ], +) def test_dsi(orivals, tuning, expected): dsi_val = dsi(orivals, tuning) - assert(np.allclose(dsi_val, expected, equal_nan=True)) + assert np.allclose(dsi_val, expected, equal_nan=True) -if __name__ == '__main__': +if __name__ == "__main__": # test_unit_ids() # test_unit_ids_filter_by_id() # test_unit_ids_filtered() @@ -369,12 +701,30 @@ def test_dsi(orivals, tuning, expected): # test_spikes() # test_get_preferred_condition() # test_get_time_to_peak() - # test_running_modulation(spike_counts=np.zeros(10), running_speeds=np.zeros(1), speed_threshold=1.0, + # test_running_modulation(spike_counts=np.zeros(10), + # running_speeds=np.zeros(1), + # speed_threshold=1.0, # expected=(np.nan, np.nan)) # test_lifetime_sparseness(np.array([1.0]), 1.0) # test_fano_factor([-1.5, 1.5], np.nan) # mean 0.0 leads to Nan - # test_overall_firing_rate(np.array([0.0]), np.array([0.0]), np.linspace(0, 10.0, 10)) - # test_osi(np.array([1.0]), np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]), np.nan) - test_osi(np.array([0.0+0.0j, 0.52359878 + 0.j, 1.04719755 + 0.j, 1.57079633+0.j, 2.0943951+0.j, 2.61799388+0.j]), - np.array([5.5, 4.44, 3.54166667, 4.10869565, 4.42, 4.55319149]), 0.07873455094232604) - pass \ No newline at end of file + # test_overall_firing_rate(np.array([0.0]), + # np.array([0.0]), + # np.linspace(0, 10.0, 10)) + # test_osi(np.array([1.0]), + # np.array([0.0, 30.0, 60.0, 90.0, 120.0, 150.0]), + # np.nan) + test_osi( + np.array( + [ + 0.0 + 0.0j, + 0.52359878 + 0.0j, + 1.04719755 + 0.0j, + 1.57079633 + 0.0j, + 2.0943951 + 0.0j, + 2.61799388 + 0.0j, + ] + ), + np.array([5.5, 4.44, 3.54166667, 4.10869565, 4.42, 4.55319149]), + 0.07873455094232604, + ) + pass diff --git a/allensdk/test/brain_observatory/ecephys/test_write_nwb.py b/allensdk/test/brain_observatory/ecephys/test_write_nwb.py index bc6fe4c91..d7a9f2080 100644 --- a/allensdk/test/brain_observatory/ecephys/test_write_nwb.py +++ b/allensdk/test/brain_observatory/ecephys/test_write_nwb.py @@ -1,142 +1,143 @@ +import logging import os +import platform from datetime import datetime, timezone from pathlib import Path -import logging -import platform -from unittest.mock import patch, create_autospec +from unittest.mock import create_autospec, patch -import pytest -import pynwb -import pandas as pd +import allensdk.brain_observatory.ecephys.nwb_util +import allensdk.brain_observatory.ecephys.utils +import allensdk.brain_observatory.ecephys.write_nwb.__main__ as write_nwb import numpy as np +import pandas as pd +import pynwb +import pytest import xarray as xr -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.sex import \ - Sex - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.full_genotype import \ - FullGenotype - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.age import \ - Age - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.mouse_id import \ - MouseId - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.driver_line import \ - DriverLine - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.reporter_line import \ - ReporterLine - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .behavior_metadata.stimulus_frame_rate import \ - StimulusFrameRate - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .behavior_metadata.session_type import \ - SessionType - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .behavior_metadata.equipment import \ - Equipment - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .behavior_metadata.date_of_acquisition import \ - DateOfAcquisition - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .behavior_metadata.behavior_session_uuid import \ - BehaviorSessionUUID - -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .behavior_metadata.project_code import \ - ProjectCode - -from allensdk.brain_observatory.behavior.data_objects import BehaviorSessionId -from allensdk.brain_observatory.behavior.data_objects.metadata\ - .subject_metadata.subject_metadata import \ - SubjectMetadata - -from pynwb import NWBFile, NWBHDF5IO - from allensdk.brain_observatory import dict_to_indexed_array -from allensdk.brain_observatory.behavior.data_objects.stimuli.presentations \ - import \ - Presentations -from allensdk.brain_observatory.ecephys._behavior_ecephys_metadata import \ - BehaviorEcephysMetadata -from allensdk.brain_observatory.ecephys._unit import Unit, \ - _get_filtered_and_sorted_spikes +from allensdk.brain_observatory.behavior.data_objects import BehaviorSessionId +from allensdk.brain_observatory.behavior.data_objects.metadata.behavior_metadata.behavior_session_uuid import ( # noqa: E501 + BehaviorSessionUUID, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.behavior_metadata.date_of_acquisition import ( # noqa: E501 + DateOfAcquisition, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.behavior_metadata.equipment import ( # noqa: E501 + Equipment, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.behavior_metadata.project_code import ( # noqa: E501 + ProjectCode, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.behavior_metadata.session_type import ( # noqa: E501 + SessionType, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.behavior_metadata.stimulus_frame_rate import ( # noqa: E501 + StimulusFrameRate, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.age import ( # noqa: E501 + Age, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.driver_line import ( # noqa: E501 + DriverLine, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.full_genotype import ( # noqa: E501 + FullGenotype, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.mouse_id import ( # noqa: E501 + MouseId, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.reporter_line import ( # noqa: E501 + ReporterLine, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.sex import ( # noqa: E501 + Sex, +) +from allensdk.brain_observatory.behavior.data_objects.metadata.subject_metadata.subject_metadata import ( # noqa: E501 + SubjectMetadata, +) +from allensdk.brain_observatory.behavior.data_objects.stimuli.presentations import ( # noqa: E501 + Presentations, +) +from allensdk.brain_observatory.ecephys._behavior_ecephys_metadata import ( + BehaviorEcephysMetadata, +) +from allensdk.brain_observatory.ecephys._unit import ( + Unit, + _get_filtered_and_sorted_spikes, +) from allensdk.brain_observatory.ecephys._units import Units -import allensdk.brain_observatory.ecephys.nwb_util -import allensdk.brain_observatory.ecephys.utils -from allensdk.brain_observatory.ecephys.current_source_density.__main__ \ - import write_csd_to_h5 -import allensdk.brain_observatory.ecephys.write_nwb.__main__ as write_nwb -from allensdk.brain_observatory.ecephys.ecephys_session_api \ - import EcephysNwbSessionApi +from allensdk.brain_observatory.ecephys.current_source_density.__main__ import ( # noqa: E501 + write_csd_to_h5, +) +from allensdk.brain_observatory.ecephys.ecephys_session_api import ( + EcephysNwbSessionApi, +) from allensdk.brain_observatory.ecephys.optotagging import OptotaggingTable -from allensdk.test.brain_observatory.behavior.test_eye_tracking_processing \ - import create_preload_eye_tracking_df from allensdk.brain_observatory.nwb import setup_table_for_invalid_times +from allensdk.test.brain_observatory.behavior.test_eye_tracking_processing import ( # noqa: E501 + create_preload_eye_tracking_df, +) +from pynwb import NWBHDF5IO, NWBFile @pytest.fixture def units_table(): - return pynwb.misc.Units.from_dataframe(pd.DataFrame({ - "peak_channel_id": [5, 10, 15], - "local_index": [0, 1, 2], - "quality": ["good", "good", "noise"], - "firing_rate": [0.5, 1.2, -3.14], - "snr": [1.0, 2.4, 5], - "isi_violations": [34, 39, 22] - }, index=pd.Index([11, 22, 33], name="id")), name="units") + return pynwb.misc.Units.from_dataframe( + pd.DataFrame( + { + "peak_channel_id": [5, 10, 15], + "local_index": [0, 1, 2], + "quality": ["good", "good", "noise"], + "firing_rate": [0.5, 1.2, -3.14], + "snr": [1.0, 2.4, 5], + "isi_violations": [34, 39, 22], + }, + index=pd.Index([11, 22, 33], name="id"), + ), + name="units", + ) @pytest.fixture def spike_times(): - return { - 11: [1, 2, 3, 4, 5, 6], - 22: [], - 33: [13, 4, 12] - } + return {11: [1, 2, 3, 4, 5, 6], 22: [], 33: [13, 4, 12]} @pytest.fixture def running_speed(): - return pd.DataFrame({ - "start_time": [1., 2., 3., 4., 5.], - "end_time": [2., 3., 4., 5., 6.], - "velocity": [-1., -2., -1., 0., 1.], - "net_rotation": [-np.pi, -2 * np.pi, -np.pi, 0, np.pi] - }) + return pd.DataFrame( + { + "start_time": [1.0, 2.0, 3.0, 4.0, 5.0], + "end_time": [2.0, 3.0, 4.0, 5.0, 6.0], + "velocity": [-1.0, -2.0, -1.0, 0.0, 1.0], + "net_rotation": [-np.pi, -2 * np.pi, -np.pi, 0, np.pi], + } + ) @pytest.fixture def raw_running_data(): - return pd.DataFrame({ - "frame_time": np.random.rand(4), - "dx": np.random.rand(4), - "vsig": np.random.rand(4), - "vin": np.random.rand(4), - }) + return pd.DataFrame( + { + "frame_time": np.random.rand(4), + "dx": np.random.rand(4), + "vsig": np.random.rand(4), + "vin": np.random.rand(4), + } + ) @pytest.fixture def stimulus_presentations_color(): - return pd.DataFrame({ - "alpha": [0.5, 0.4, 0.3, 0.2, 0.1], - "start_time": [1., 2., 4., 5., 6.], - "stop_time": [2., 4., 5., 6., 8.], - "stimulus_name": ['gabors', 'gabors', 'random', 'movie', 'gabors'], - "color": ["1.0", "", r"[1.0,-1.0, 10., -42.12, -.1]", "-1.0", ""] - }, index=pd.Index(name="stimulus_presentations_id", data=[0, 1, 2, 3, 4])) + return pd.DataFrame( + { + "alpha": [0.5, 0.4, 0.3, 0.2, 0.1], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + "stimulus_name": ["gabors", "gabors", "random", "movie", "gabors"], + "color": ["1.0", "", r"[1.0,-1.0, 10., -42.12, -.1]", "-1.0", ""], + }, + index=pd.Index(name="stimulus_presentations_id", data=[0, 1, 2, 3, 4]), + ) def test_roundtrip_basic_metadata(roundtripper): @@ -144,7 +145,7 @@ def test_roundtrip_basic_metadata(roundtripper): nwbfile = pynwb.NWBFile( session_description="EcephysSession", identifier="{}".format(12345), - session_start_time=dt + session_start_time=dt, ) api = roundtripper(nwbfile, EcephysNwbSessionApi) @@ -152,27 +153,34 @@ def test_roundtrip_basic_metadata(roundtripper): assert dt == api.get_session_start_time() -@pytest.mark.parametrize("metadata, expected_metadata", [ - ({ - "specimen_name": "mouse_1", - "age_in_days": 100.0, - "full_genotype": "wt", - "strain": "c57", - "sex": "F", - "stimulus_name": "brain_observatory_2.0", - "donor_id": 12345, - "species": "Mus musculus"}, - { - "specimen_name": "mouse_1", - "age_in_days": 100.0, - "age": "P100D", - "full_genotype": "wt", - "strain": "c57", - "sex": "F", - "stimulus_name": "brain_observatory_2.0", - "subject_id": "12345", - "species": "Mus musculus"}) -]) +@pytest.mark.parametrize( + "metadata, expected_metadata", + [ + ( + { + "specimen_name": "mouse_1", + "age_in_days": 100.0, + "full_genotype": "wt", + "strain": "c57", + "sex": "F", + "stimulus_name": "brain_observatory_2.0", + "donor_id": 12345, + "species": "Mus musculus", + }, + { + "specimen_name": "mouse_1", + "age_in_days": 100.0, + "age": "P100D", + "full_genotype": "wt", + "strain": "c57", + "sex": "F", + "stimulus_name": "brain_observatory_2.0", + "subject_id": "12345", + "species": "Mus musculus", + }, + ) + ], +) def test_add_metadata(nwbfile, roundtripper, metadata, expected_metadata): nwbfile = write_nwb.add_metadata_to_nwbfile(nwbfile, metadata) @@ -186,54 +194,87 @@ def test_add_metadata(nwbfile, roundtripper, metadata, expected_metadata): if obtained[key] != value: misses[key] = {"expected": value, "obtained": obtained[key]} - assert len(misses) == 0, \ - f"the following metadata items were mismatched: {misses}" - - -@pytest.mark.parametrize("presentations", [ - (pd.DataFrame({ - 'alpha': [0.5, 0.4, 0.3, 0.2, 0.1], - 'start_time': [1., 2., 4., 5., 6.], - 'stimulus_name': ['gabors', 'gabors', 'random', 'movie', 'gabors'], - 'stop_time': [2., 4., 5., 6., 8.] - }, index=pd.Index(name='stimulus_presentations_id', - data=[0, 1, 2, 3, 4]))), - - (pd.DataFrame({ - 'gabor_specific_column': [1.0, 2.0, np.nan, np.nan, 3.0], - 'mixed_column': ["a", "", "b", "", "c"], - 'movie_specific_column': [np.nan, np.nan, np.nan, 1.0, np.nan], - 'start_time': [1., 2., 4., 5., 6.], - 'stimulus_name': ['gabors', 'gabors', 'random', 'movie', 'gabors'], - 'stop_time': [2., 4., 5., 6., 8.], - 'color': [np.nan] + ['[1.0, 1.0, 1.0]'] * 4 - }, index=pd.Index(name='stimulus_presentations_id', - data=[0, 1, 2, 3, 4]))), -]) + assert ( + len(misses) == 0 + ), f"the following metadata items were mismatched: {misses}" + + +@pytest.mark.parametrize( + "presentations", + [ + ( + pd.DataFrame( + { + "alpha": [0.5, 0.4, 0.3, 0.2, 0.1], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + }, + index=pd.Index( + name="stimulus_presentations_id", data=[0, 1, 2, 3, 4] + ), + ) + ), + ( + pd.DataFrame( + { + "gabor_specific_column": [1.0, 2.0, np.nan, np.nan, 3.0], + "mixed_column": ["a", "", "b", "", "c"], + "movie_specific_column": [ + np.nan, + np.nan, + np.nan, + 1.0, + np.nan, + ], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + "color": [np.nan] + ["[1.0, 1.0, 1.0]"] * 4, + }, + index=pd.Index( + name="stimulus_presentations_id", data=[0, 1, 2, 3, 4] + ), + ) + ), + ], +) def test_add_stimulus_presentations(nwbfile, presentations, roundtripper): write_nwb.add_stimulus_timestamps(nwbfile, [0, 1]) presentations = Presentations(presentations=presentations) - presentations.to_nwb(nwbfile=nwbfile, stimulus_name_column='stimulus_name') + presentations.to_nwb(nwbfile=nwbfile, stimulus_name_column="stimulus_name") api = roundtripper(nwbfile, EcephysNwbSessionApi) obtained_stimulus_table = api.get_stimulus_presentations() - if 'color' in presentations.value: - presentations.value['color_triplet'] = [''] + ['[1.0, 1.0, 1.0]'] * 4 - presentations.value['color'] = '' + if "color" in presentations.value: + presentations.value["color_triplet"] = [""] + ["[1.0, 1.0, 1.0]"] * 4 + presentations.value["color"] = "" pd.testing.assert_frame_equal( presentations.value[sorted(presentations.value.columns)], obtained_stimulus_table[sorted(obtained_stimulus_table.columns)], - check_dtype=False) + check_dtype=False, + ) def test_add_stimulus_presentations_color( - nwbfile, - stimulus_presentations_color, - roundtripper): + nwbfile, stimulus_presentations_color, roundtripper +): write_nwb.add_stimulus_timestamps(nwbfile, [0, 1]) presentations = Presentations(presentations=stimulus_presentations_color) - presentations.to_nwb(nwbfile=nwbfile, stimulus_name_column='stimulus_name') + presentations.to_nwb(nwbfile=nwbfile, stimulus_name_column="stimulus_name") api = roundtripper(nwbfile, EcephysNwbSessionApi) obtained_stimulus_table = api.get_stimulus_presentations() @@ -246,44 +287,65 @@ def test_add_stimulus_presentations_color( if expected != obtained: mismatched = True - assert not mismatched, \ - f"expected: {expected_color}, obtained: {obtained_color}" - - -@pytest.mark.parametrize("opto_table, expected", [ - (pd.DataFrame({ - "start_time": [0., 1., 2., 3.], - "stop_time": [0.5, 1.5, 2.5, 3.5], - "level": [10., 9., 8., 7.], - "condition": ["a", "a", "b", "c"]}), - None), - - # Test for older version of optotable that used nwb reserved "name" col - (pd.DataFrame({"start_time": [0., 1., 2., 3.], - "stop_time": [0.5, 1.5, 2.5, 3.5], - "level": [10., 9., 8., 7.], - "condition": ["a", "a", "b", "c"], - "name": ["w", "x", "y", "z"]}), - pd.DataFrame({"start_time": [0., 1., 2., 3.], - "stop_time": [0.5, 1.5, 2.5, 3.5], - "level": [10., 9., 8., 7.], - "condition": ["a", "a", "b", "c"], - "stimulus_name": ["w", "x", "y", "z"], - "duration": [0.5, 0.5, 0.5, 0.5]}, - index=pd.Index(name="id", data=[0, 1, 2, 3]))), - - (pd.DataFrame({"start_time": [0., 1., 2., 3.], - "stop_time": [0.5, 1.5, 2.5, 3.5], - "level": [10., 9., 8., 7.], - "condition": ["a", "a", "b", "c"], - "stimulus_name": ["w", "x", "y", "z"]}), - None) -]) + assert ( + not mismatched + ), f"expected: {expected_color}, obtained: {obtained_color}" + + +@pytest.mark.parametrize( + "opto_table, expected", + [ + ( + pd.DataFrame( + { + "start_time": [0.0, 1.0, 2.0, 3.0], + "stop_time": [0.5, 1.5, 2.5, 3.5], + "level": [10.0, 9.0, 8.0, 7.0], + "condition": ["a", "a", "b", "c"], + } + ), + None, + ), + # Test for older version of optotable that used nwb reserved "name" col + ( + pd.DataFrame( + { + "start_time": [0.0, 1.0, 2.0, 3.0], + "stop_time": [0.5, 1.5, 2.5, 3.5], + "level": [10.0, 9.0, 8.0, 7.0], + "condition": ["a", "a", "b", "c"], + "name": ["w", "x", "y", "z"], + } + ), + pd.DataFrame( + { + "start_time": [0.0, 1.0, 2.0, 3.0], + "stop_time": [0.5, 1.5, 2.5, 3.5], + "level": [10.0, 9.0, 8.0, 7.0], + "condition": ["a", "a", "b", "c"], + "stimulus_name": ["w", "x", "y", "z"], + "duration": [0.5, 0.5, 0.5, 0.5], + }, + index=pd.Index(name="id", data=[0, 1, 2, 3]), + ), + ), + ( + pd.DataFrame( + { + "start_time": [0.0, 1.0, 2.0, 3.0], + "stop_time": [0.5, 1.5, 2.5, 3.5], + "level": [10.0, 9.0, 8.0, 7.0], + "condition": ["a", "a", "b", "c"], + "stimulus_name": ["w", "x", "y", "z"], + } + ), + None, + ), + ], +) def test_add_optotagging_table_to_nwbfile( - nwbfile, - roundtripper, - opto_table, - expected): + nwbfile, roundtripper, opto_table, expected +): opto_table["duration"] = opto_table["stop_time"] - opto_table["start_time"] opto_table = OptotaggingTable(table=opto_table) @@ -294,44 +356,57 @@ def test_add_optotagging_table_to_nwbfile( if expected is None: expected = opto_table.value - expected.index.name = 'id' + expected.index.name = "id" pd.testing.assert_frame_equal(obtained, expected, check_like=True) @pytest.mark.parametrize("roundtrip", [True, False]) -@pytest.mark.parametrize("pid,name,srate,lfp_srate,has_lfp,expected", [ +@pytest.mark.parametrize( + "pid,name,srate,lfp_srate,has_lfp,expected", [ - 12, - "a probe", - 30000.0, - 2500.0, - True, - pd.DataFrame({ - "description": ["a probe"], - "sampling_rate": [30000.0], - "lfp_sampling_rate": [2500.0], - "has_lfp_data": [True], - "location": ["See electrode locations"] - }, index=pd.Index([12], name="id")) - ] -]) + [ + 12, + "a probe", + 30000.0, + 2500.0, + True, + pd.DataFrame( + { + "description": ["a probe"], + "sampling_rate": [30000.0], + "lfp_sampling_rate": [2500.0], + "has_lfp_data": [True], + "location": ["See electrode locations"], + }, + index=pd.Index([12], name="id"), + ), + ] + ], +) def test_add_probe_to_nwbfile( + nwbfile, + roundtripper, + roundtrip, + pid, + name, + srate, + lfp_srate, + has_lfp, + expected, +): + ( + nwbfile, + _, + _, + ) = allensdk.brain_observatory.ecephys.nwb_util.add_probe_to_nwbfile( nwbfile, - roundtripper, - roundtrip, pid, - name, - srate, - lfp_srate, - has_lfp, - expected): - nwbfile, _, _ = allensdk.brain_observatory.ecephys.nwb_util \ - .add_probe_to_nwbfile(nwbfile, pid, - name=name, - sampling_rate=srate, - lfp_sampling_rate=lfp_srate, - has_lfp_data=has_lfp) + name=name, + sampling_rate=srate, + lfp_sampling_rate=lfp_srate, + has_lfp_data=has_lfp, + ) if roundtrip: obt = roundtripper(nwbfile, EcephysNwbSessionApi) else: @@ -340,21 +415,26 @@ def test_add_probe_to_nwbfile( pd.testing.assert_frame_equal(expected, obt.get_probes(), check_like=True) -@pytest.mark.parametrize("columns_to_add", [ - None, - - [("test_column_a", "description_a"), - ("test_column_b", "description_b")] -]) +@pytest.mark.parametrize( + "columns_to_add", + [ + None, + [ + ("test_column_a", "description_a"), + ("test_column_b", "description_b"), + ], + ], +) def test_add_ecephys_electrode_columns(nwbfile, columns_to_add): allensdk.brain_observatory.ecephys.nwb_util._add_ecephys_electrode_columns( - nwbfile, columns_to_add) + nwbfile, columns_to_add + ) if columns_to_add is None: - expected_columns = \ - [x[0] for x in - allensdk.brain_observatory.ecephys.nwb_util. - ELECTRODE_TABLE_DEFAULT_COLUMNS] + expected_columns = [ + x[0] + for x in allensdk.brain_observatory.ecephys.nwb_util.ELECTRODE_TABLE_DEFAULT_COLUMNS # noqa: E501 + ] else: expected_columns = [x[0] for x in columns_to_add] @@ -362,80 +442,98 @@ def test_add_ecephys_electrode_columns(nwbfile, columns_to_add): assert c in nwbfile.electrodes.colnames -@pytest.mark.parametrize(("channels, channel_number_whitelist, " - "expected_electrode_table"), [ - ([{"id": 1, - "probe_id": 1234, - "valid_data": True, - "probe_channel_number": 23, - "probe_vertical_position": 10, - "probe_horizontal_position": 10, - "anterior_posterior_ccf_coordinate": 15.0, - "dorsal_ventral_ccf_coordinate": 20.0, - "left_right_ccf_coordinate": 25.0, - "structure_acronym": "CA1", - "impedance": np.nan, - "filtering": "AP band: 500 Hz high-pass; LFP " - "band: 1000 Hz low-pass"}, - {"id": 2, - "probe_id": 1234, - "valid_data": True, - "probe_channel_number": 15, - "probe_vertical_position": 20, - "probe_horizontal_position": 20, - "anterior_posterior_ccf_coordinate": 25.0, - "dorsal_ventral_ccf_coordinate": 30.0, - "left_right_ccf_coordinate": 35.0, - "structure_acronym": "CA3", - "impedance": 42.0, - "filtering": "custom"}], - - [15, 23], - - pd.DataFrame({ - "id": [2, 1], - "probe_id": [1234, 1234], - "valid_data": [True, True], - "probe_channel_number": [15, 23], - "probe_vertical_position": [20, 10], - "probe_horizontal_position": [20, 10], - "x": [25.0, 15.0], - "y": [30.0, 20.0], - "z": [35.0, 25.0], - "location": ["CA3", "CA1"], - "impedance": [42.0, np.nan], - "filtering": ["custom", - "AP band: 500 Hz high-pass; " - "LFP band: 1000 Hz low-pass"] - }).set_index("id")) - - ]) -def test_add_ecephys_electrodes(nwbfile, channels, channel_number_whitelist, - expected_electrode_table): +@pytest.mark.parametrize( + ("channels, channel_number_whitelist, " "expected_electrode_table"), + [ + ( + [ + { + "id": 1, + "probe_id": 1234, + "valid_data": True, + "probe_channel_number": 23, + "probe_vertical_position": 10, + "probe_horizontal_position": 10, + "anterior_posterior_ccf_coordinate": 15.0, + "dorsal_ventral_ccf_coordinate": 20.0, + "left_right_ccf_coordinate": 25.0, + "structure_acronym": "CA1", + "impedance": np.nan, + "filtering": "AP band: 500 Hz high-pass; LFP " + "band: 1000 Hz low-pass", + }, + { + "id": 2, + "probe_id": 1234, + "valid_data": True, + "probe_channel_number": 15, + "probe_vertical_position": 20, + "probe_horizontal_position": 20, + "anterior_posterior_ccf_coordinate": 25.0, + "dorsal_ventral_ccf_coordinate": 30.0, + "left_right_ccf_coordinate": 35.0, + "structure_acronym": "CA3", + "impedance": 42.0, + "filtering": "custom", + }, + ], + [15, 23], + pd.DataFrame( + { + "id": [2, 1], + "probe_id": [1234, 1234], + "valid_data": [True, True], + "probe_channel_number": [15, 23], + "probe_vertical_position": [20, 10], + "probe_horizontal_position": [20, 10], + "x": [25.0, 15.0], + "y": [30.0, 20.0], + "z": [35.0, 25.0], + "location": ["CA3", "CA1"], + "impedance": [42.0, np.nan], + "filtering": [ + "custom", + "AP band: 500 Hz high-pass; " + "LFP band: 1000 Hz low-pass", + ], + } + ).set_index("id"), + ) + ], +) +def test_add_ecephys_electrodes( + nwbfile, channels, channel_number_whitelist, expected_electrode_table +): mock_device = pynwb.device.Device(name="mock_device") - mock_electrode_group = pynwb.ecephys.ElectrodeGroup(name="mock_group", - description="", - location="", - device=mock_device) + mock_electrode_group = pynwb.ecephys.ElectrodeGroup( + name="mock_group", description="", location="", device=mock_device + ) allensdk.brain_observatory.ecephys.nwb_util.add_ecephys_electrodes( - nwbfile, - channels, - mock_electrode_group, - channel_number_whitelist) + nwbfile, channels, mock_electrode_group, channel_number_whitelist + ) - obt_electrode_table = \ - nwbfile.electrodes.to_dataframe().drop(columns=["group", "group_name"]) + obt_electrode_table = nwbfile.electrodes.to_dataframe().drop( + columns=["group", "group_name"] + ) - expected_electrode_table.rename(columns={'impedance': 'imp'}, inplace=True) - pd.testing.assert_frame_equal(obt_electrode_table, - expected_electrode_table, - check_like=True) + expected_electrode_table.rename(columns={"impedance": "imp"}, inplace=True) + pd.testing.assert_frame_equal( + obt_electrode_table, expected_electrode_table, check_like=True + ) -@pytest.mark.parametrize("dc,order,exp_idx,exp_data", [ - [{"a": [1, 2, 3], "b": [4, 5, 6]}, ["a", "b"], [3, 6], [1, 2, 3, 4, 5, 6]] -]) +@pytest.mark.parametrize( + "dc,order,exp_idx,exp_data", + [ + [ + {"a": [1, 2, 3], "b": [4, 5, 6]}, + ["a", "b"], + [3, 6], + [1, 2, 3, 4, 5, 6], + ] + ], +) def test_dict_to_indexed_array(dc, order, exp_idx, exp_data): obt_idx, obt_data = dict_to_indexed_array(dc, order) assert np.allclose(exp_idx, obt_idx) @@ -443,28 +541,21 @@ def test_dict_to_indexed_array(dc, order, exp_idx, exp_data): def test_add_ragged_data_to_dynamic_table(units_table, spike_times): - allensdk.brain_observatory.ecephys.nwb_util \ - .add_ragged_data_to_dynamic_table( - table=units_table, - data=spike_times, - column_name="spike_times" - ) + allensdk.brain_observatory.ecephys.nwb_util.add_ragged_data_to_dynamic_table( # noqa: E501 + table=units_table, data=spike_times, column_name="spike_times" + ) assert np.allclose([1, 2, 3, 4, 5, 6], units_table["spike_times"][0]) assert np.allclose([], units_table["spike_times"][1]) assert np.allclose([13, 4, 12], units_table["spike_times"][2]) -@pytest.mark.parametrize("roundtrip,include_rotation", [ - [True, True], - [True, False] -]) +@pytest.mark.parametrize( + "roundtrip,include_rotation", [[True, True], [True, False]] +) def test_add_running_speed_to_nwbfile( - nwbfile, - running_speed, - roundtripper, - roundtrip, - include_rotation): + nwbfile, running_speed, roundtripper, roundtrip, include_rotation +): nwbfile = write_nwb.add_running_speed_to_nwbfile(nwbfile, running_speed) if roundtrip: api_obt = roundtripper(nwbfile, EcephysNwbSessionApi) @@ -481,13 +572,11 @@ def test_add_running_speed_to_nwbfile( @pytest.mark.parametrize("roundtrip", [[True]]) def test_add_raw_running_data_to_nwbfile( - nwbfile, - raw_running_data, - roundtripper, - roundtrip): + nwbfile, raw_running_data, roundtripper, roundtrip +): nwbfile = write_nwb.add_raw_running_data_to_nwbfile( - nwbfile, - raw_running_data) + nwbfile, raw_running_data + ) if roundtrip: api_obt = roundtripper(nwbfile, EcephysNwbSessionApi) else: @@ -496,75 +585,125 @@ def test_add_raw_running_data_to_nwbfile( obtained = api_obt.get_raw_running_data() expected = raw_running_data.rename( - columns={"dx": "net_rotation", - "vsig": "signal_voltage", - "vin": "supply_voltage"}) + columns={ + "dx": "net_rotation", + "vsig": "signal_voltage", + "vin": "supply_voltage", + } + ) pd.testing.assert_frame_equal(expected, obtained, check_like=True) @pytest.mark.parametrize( - "presentations, column_renames_map, columns_to_drop, expected", [ - (pd.DataFrame({'alpha': [0.5, 0.4, 0.3, 0.2, 0.1], - 'stimulus_name': ['gabors', - 'gabors', - 'random', - 'movie', - 'gabors'], - 'start_time': [1., 2., 4., 5., 6.], - 'stop_time': [2., 4., 5., 6., 8.]}), - {"alpha": "beta"}, - None, - pd.DataFrame({'beta': [0.5, 0.4, 0.3, 0.2, 0.1], - 'stimulus_name': ['gabors', - 'gabors', - 'random', - 'movie', - 'gabors'], - 'start_time': [1., 2., 4., 5., 6.], - 'stop_time': [2., 4., 5., 6., 8.]})), - - (pd.DataFrame({'alpha': [0.5, 0.4, 0.3, 0.2, 0.1], - 'stimulus_name': ['gabors', - 'gabors', - 'random', - 'movie', - 'gabors'], - 'start_time': [1., 2., 4., 5., 6.], - 'stop_time': [2., 4., 5., 6., 8.]}), - {"alpha": "beta"}, - ["Nonexistant_column_to_drop"], - pd.DataFrame({'beta': [0.5, 0.4, 0.3, 0.2, 0.1], - 'stimulus_name': ['gabors', - 'gabors', - 'random', - 'movie', - 'gabors'], - 'start_time': [1., 2., 4., 5., 6.], - 'stop_time': [2., 4., 5., 6., 8.]})), - - (pd.DataFrame({'alpha': [0.5, 0.4, 0.3, 0.2, 0.1], - 'stimulus_name': ['gabors', - 'gabors', - 'random', - 'movie', - 'gabors'], - 'Start': [1., 2., 4., 5., 6.], - 'End': [2., 4., 5., 6., 8.]}), - None, - ["alpha"], - pd.DataFrame({'stimulus_name': ['gabors', - 'gabors', - 'random', - 'movie', - 'gabors'], - 'start_time': [1., 2., 4., 5., 6.], - 'stop_time': [2., 4., 5., 6., 8.]})), - ]) -def test_read_stimulus_table(tmpdir_factory, presentations, - column_renames_map, columns_to_drop, expected): + "presentations, column_renames_map, columns_to_drop, expected", + [ + ( + pd.DataFrame( + { + "alpha": [0.5, 0.4, 0.3, 0.2, 0.1], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + } + ), + {"alpha": "beta"}, + None, + pd.DataFrame( + { + "beta": [0.5, 0.4, 0.3, 0.2, 0.1], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + } + ), + ), + ( + pd.DataFrame( + { + "alpha": [0.5, 0.4, 0.3, 0.2, 0.1], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + } + ), + {"alpha": "beta"}, + ["Nonexistant_column_to_drop"], + pd.DataFrame( + { + "beta": [0.5, 0.4, 0.3, 0.2, 0.1], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + } + ), + ), + ( + pd.DataFrame( + { + "alpha": [0.5, 0.4, 0.3, 0.2, 0.1], + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "Start": [1.0, 2.0, 4.0, 5.0, 6.0], + "End": [2.0, 4.0, 5.0, 6.0, 8.0], + } + ), + None, + ["alpha"], + pd.DataFrame( + { + "stimulus_name": [ + "gabors", + "gabors", + "random", + "movie", + "gabors", + ], + "start_time": [1.0, 2.0, 4.0, 5.0, 6.0], + "stop_time": [2.0, 4.0, 5.0, 6.0, 8.0], + } + ), + ), + ], +) +def test_read_stimulus_table( + tmpdir_factory, + presentations, + column_renames_map, + columns_to_drop, + expected, +): expected = expected.set_index( - pd.Index(range(expected.shape[0]), - name='stimulus_presentations_id')) + pd.Index(range(expected.shape[0]), name="stimulus_presentations_id") + ) dirname = str(tmpdir_factory.mktemp("ecephys_nwb_test")) stim_table_path = os.path.join(dirname, "stim_table.csv") @@ -576,14 +715,15 @@ def add_is_image_novel(stimulus_presentations, behavior_session_id): # not testing this for vcn return None - with patch.object(Presentations, '_add_is_image_novel', - wraps=add_is_image_novel): + with patch.object( + Presentations, "_add_is_image_novel", wraps=add_is_image_novel + ): obt = Presentations.from_path( path=stim_table_path, behavior_session_id=1, exclude_columns=columns_to_drop, columns_to_rename=column_renames_map, - sort_columns=True + sort_columns=True, ) obtained = obt.value[sorted(obt.value.columns)] @@ -604,15 +744,13 @@ def test_read_spike_times_to_dictionary(tmpdir_factory): local_to_global_unit_map = {ii: -ii for ii in spike_units} - obtained = allensdk.brain_observatory.ecephys._units \ - ._read_spike_times_to_dictionary( - spike_times_path, - spike_units_path, - local_to_global_unit_map) + obtained = allensdk.brain_observatory.ecephys._units._read_spike_times_to_dictionary( # noqa: E501 + spike_times_path, spike_units_path, local_to_global_unit_map + ) for ii in range(15): assert np.allclose( - obtained[-ii], - sorted([spike_times[ii], spike_times[15 + ii]])) + obtained[-ii], sorted([spike_times[ii], spike_times[15 + ii]]) + ) def test_read_waveforms_to_dictionary(tmpdir_factory): @@ -628,10 +766,9 @@ def test_read_waveforms_to_dictionary(tmpdir_factory): mean_waveforms = np.random.rand(nunits, nsamples, nchannels) np.save(waveforms_path, mean_waveforms, allow_pickle=False) - obtained = allensdk.brain_observatory.ecephys._units \ - ._read_waveforms_to_dictionary( - waveforms_path, - local_to_global_unit_map) + obtained = allensdk.brain_observatory.ecephys._units._read_waveforms_to_dictionary( # noqa: E501 + waveforms_path, local_to_global_unit_map + ) for ii in range(nunits): assert np.allclose(mean_waveforms[ii, :, :], obtained[-ii]) @@ -643,11 +780,10 @@ def lfp_data(): return { "data": np.arange( - total_timestamps * len(subsample_channels), - dtype=np.int16).reshape((total_timestamps, - len(subsample_channels))), + total_timestamps * len(subsample_channels), dtype=np.int16 + ).reshape((total_timestamps, len(subsample_channels))), "timestamps": np.linspace(0, 1, total_timestamps), - "subsample_channels": subsample_channels + "subsample_channels": subsample_channels, } @@ -673,7 +809,7 @@ def probe_data(): "structure_acronym": "CA1", "impedence": np.nan, "filtering": "AP band: 500 Hz high-pass; " - "LFP band: 1000 Hz low-pass" + "LFP band: 1000 Hz low-pass", }, { "id": 1, @@ -688,7 +824,7 @@ def probe_data(): "structure_acronym": "CA2", "impedence": np.nan, "filtering": "AP band: 500 Hz high-pass; " - "LFP band: 1000 Hz low-pass" + "LFP band: 1000 Hz low-pass", }, { "id": 2, @@ -703,17 +839,17 @@ def probe_data(): "structure_acronym": "CA3", "impedence": np.nan, "filtering": "AP band: 500 Hz high-pass; " - "LFP band: 1000 Hz low-pass" - } + "LFP band: 1000 Hz low-pass", + }, ], "lfp": { "input_data_path": "", "input_timestamps_path": "", "input_channels_path": "", - "output_path": "" + "output_path": "", }, "csd_path": "", - "amplitude_scale_factor": 1.0 + "amplitude_scale_factor": 1.0, } return probe_data @@ -727,7 +863,7 @@ def csd_data(): "csd_locations": np.array([[1, 2], [3, 3]]), "stimulus_name": "foo", "stimulus_index": None, - "num_trials": 1000 + "num_trials": 1000, } return csd_data @@ -745,7 +881,7 @@ def test_write_probe_lfp_file(tmpdir_factory, lfp_data, probe_data, csd_data): "input_data_path": input_data_path, "input_timestamps_path": input_timestamps_path, "input_channels_path": input_channels_path, - "output_path": output_path + "output_path": output_path, } test_session_metadata = { @@ -757,7 +893,7 @@ def test_write_probe_lfp_file(tmpdir_factory, lfp_data, probe_data, csd_data): "stimulus_name": "test_stim", "species": "Mus musculus", "donor_id": 42, - "project_code": '1234', + "project_code": "1234", } def dummy_meta_from_json(dict_repr): @@ -766,20 +902,20 @@ def dummy_meta_from_json(dict_repr): behavior_session_id=BehaviorSessionId(1), behavior_session_uuid=BehaviorSessionUUID(None), date_of_acquisition=DateOfAcquisition( - date_of_acquisition=datetime.now()), - equipment=Equipment('foo'), - session_type=SessionType('foo'), + date_of_acquisition=datetime.now() + ), + equipment=Equipment("foo"), + session_type=SessionType("foo"), stimulus_frame_rate=StimulusFrameRate(1.0), - project_code=ProjectCode('1234'), + project_code=ProjectCode("1234"), subject_metadata=SubjectMetadata( - sex=Sex(dict_repr['sex']), + sex=Sex(dict_repr["sex"]), age=Age(1), - reporter_line=ReporterLine('foo'), - driver_line=DriverLine(['foo']), - full_genotype=FullGenotype(dict_repr['full_genotype']), - mouse_id=MouseId(dict_repr['donor_id']) - - ) + reporter_line=ReporterLine("foo"), + driver_line=DriverLine(["foo"]), + full_genotype=FullGenotype(dict_repr["full_genotype"]), + mouse_id=MouseId(dict_repr["donor_id"]), + ), ) probe_data.update({"lfp": test_lfp_paths}) @@ -787,38 +923,43 @@ def dummy_meta_from_json(dict_repr): write_csd_to_h5(path=input_csd_path, **csd_data) - np.save(input_timestamps_path, - lfp_data["timestamps"], - allow_pickle=False) - np.save(input_channels_path, - lfp_data["subsample_channels"], - allow_pickle=False) + np.save(input_timestamps_path, lfp_data["timestamps"], allow_pickle=False) + np.save( + input_channels_path, lfp_data["subsample_channels"], allow_pickle=False + ) with open(input_data_path, "wb") as input_data_file: input_data_file.write(lfp_data["data"].tobytes()) - with patch.object(Units, 'from_json', wraps=lambda probe: None): - with patch.object(BehaviorEcephysMetadata, 'from_json', - wraps=dummy_meta_from_json): + with patch.object(Units, "from_json", wraps=lambda probe: None): + with patch.object( + BehaviorEcephysMetadata, "from_json", wraps=dummy_meta_from_json + ): write_nwb.write_probe_lfp_file( 4242, test_session_metadata, datetime.now(), - logging.INFO, probe_data) + logging.INFO, + probe_data, + ) - exp_electrodes = \ + exp_electrodes = ( pd.DataFrame(probe_data["channels"]).set_index("id").loc[[2, 1], :] - exp_electrodes = exp_electrodes.rename(columns={'impedance': 'imp'}) - exp_electrodes.rename(columns={"anterior_posterior_ccf_coordinate": "x", - "dorsal_ventral_ccf_coordinate": "y", - "left_right_ccf_coordinate": "z", - "structure_acronym": "location"}, - inplace=True) + ) + exp_electrodes = exp_electrodes.rename(columns={"impedance": "imp"}) + exp_electrodes.rename( + columns={ + "anterior_posterior_ccf_coordinate": "x", + "dorsal_ventral_ccf_coordinate": "y", + "left_right_ccf_coordinate": "z", + "structure_acronym": "location", + }, + inplace=True, + ) with pynwb.NWBHDF5IO(output_path, "r") as obt_io: obt_f = obt_io.read() - obt_acq = \ - obt_f.get_acquisition("probe_12345_lfp") + obt_acq = obt_f.get_acquisition("probe_12345_lfp") obt_ser = obt_acq.electrical_series["probe_12345_lfp_data"] assert np.allclose(lfp_data["data"], obt_ser.data[:]) assert np.allclose(lfp_data["timestamps"], obt_ser.timestamps[:]) @@ -826,12 +967,21 @@ def dummy_meta_from_json(dict_repr): obt_electrodes_df = obt_f.electrodes.to_dataframe() obt_electrodes = obt_electrodes_df.loc[ - :, ["probe_channel_number", - "probe_horizontal_position", - "probe_id", "probe_vertical_position", - "valid_data", "x", "y", "z", "location", "imp", - "filtering"] - ] + :, + [ + "probe_channel_number", + "probe_horizontal_position", + "probe_id", + "probe_vertical_position", + "valid_data", + "x", + "y", + "z", + "location", + "imp", + "filtering", + ], + ] assert obt_f.session_id == "4242" assert obt_f.subject.subject_id == "42" @@ -844,24 +994,30 @@ def dummy_meta_from_json(dict_repr): obt_electrodes, exp_electrodes, check_like=True, - check_dtype=False) + check_dtype=False, + ) else: pd.testing.assert_frame_equal( - obt_electrodes, - exp_electrodes, - check_like=True) + obt_electrodes, exp_electrodes, check_like=True + ) - processing_module = \ - obt_f.get_processing_module("current_source_density") + processing_module = obt_f.get_processing_module( + "current_source_density" + ) csd_series = processing_module["ecephys_csd"] assert np.allclose(csd_data["csd"], csd_series.time_series.data[:].T) - assert np.allclose(csd_data["relative_window"], - csd_series.time_series.timestamps[:]) - obt_channel_locations = \ - np.stack((csd_series.virtual_electrode_x_positions, - csd_series.virtual_electrode_y_positions), axis=1) + assert np.allclose( + csd_data["relative_window"], csd_series.time_series.timestamps[:] + ) + obt_channel_locations = np.stack( + ( + csd_series.virtual_electrode_x_positions, + csd_series.virtual_electrode_y_positions, + ), + axis=1, + ) # csd interpolated channel locations assert np.allclose([[1, 2], [3, 3]], obt_channel_locations) @@ -869,11 +1025,8 @@ def dummy_meta_from_json(dict_repr): @pytest.mark.parametrize("roundtrip", [True, False]) def test_write_probe_lfp_file_roundtrip( - tmpdir_factory, - roundtrip, - lfp_data, - probe_data, - csd_data): + tmpdir_factory, roundtrip, lfp_data, probe_data, csd_data +): expected_csd = xr.DataArray( name="CSD", data=csd_data["csd"], @@ -883,20 +1036,20 @@ def test_write_probe_lfp_file_roundtrip( "time": csd_data["relative_window"], "vertical_position": ( ("virtual_channel_index",), - csd_data["csd_locations"][:, 1] + csd_data["csd_locations"][:, 1], ), "horizontal_position": ( ("virtual_channel_index",), - csd_data["csd_locations"][:, 0] + csd_data["csd_locations"][:, 0], ), - } + }, ) expected_lfp = xr.DataArray( name="LFP", data=lfp_data["data"], dims=["time", "channel"], - coords=[lfp_data["timestamps"], [2, 1]] + coords=[lfp_data["timestamps"], [2, 1]], ) tmpdir = Path(tmpdir_factory.mktemp("probe_lfp_nwb")) @@ -910,7 +1063,7 @@ def test_write_probe_lfp_file_roundtrip( "input_data_path": input_data_path, "input_timestamps_path": input_timestamps_path, "input_channels_path": input_channels_path, - "output_path": output_path + "output_path": output_path, } probe_data.update({"lfp": test_lfp_paths}) @@ -918,29 +1071,28 @@ def test_write_probe_lfp_file_roundtrip( write_csd_to_h5(path=input_csd_path, **csd_data) - np.save(input_timestamps_path, - lfp_data["timestamps"], - allow_pickle=False) - np.save(input_channels_path, - lfp_data["subsample_channels"], - allow_pickle=False) + np.save(input_timestamps_path, lfp_data["timestamps"], allow_pickle=False) + np.save( + input_channels_path, lfp_data["subsample_channels"], allow_pickle=False + ) with open(input_data_path, "wb") as input_data_file: input_data_file.write(lfp_data["data"].tobytes()) - with patch.object(Units, 'from_json', wraps=lambda probe: None): - with patch.object(BehaviorEcephysMetadata, 'from_json', - wraps=lambda dict_repr: create_autospec( - BehaviorEcephysMetadata, instance=True)): + with patch.object(Units, "from_json", wraps=lambda probe: None): + with patch.object( + BehaviorEcephysMetadata, + "from_json", + wraps=lambda dict_repr: create_autospec( + BehaviorEcephysMetadata, instance=True + ), + ): write_nwb.write_probe_lfp_file( - 4242, - None, - datetime.now(), - logging.INFO, - probe_data) + 4242, None, datetime.now(), logging.INFO, probe_data + ) obt = EcephysNwbSessionApi( - path=None, - probe_lfp_paths={12345: NWBHDF5IO(output_path, "r").read}) + path=None, probe_lfp_paths={12345: NWBHDF5IO(output_path, "r").read} + ) obtained_lfp = obt.get_lfp(12345) obtained_csd = obt.get_current_source_density(12345) @@ -979,13 +1131,14 @@ def invalid_epochs(): def test_add_invalid_times(invalid_epochs, tmpdir_factory): - nwbfile_name = \ - str(tmpdir_factory.mktemp("test").join("test_invalid_times.nwb")) + nwbfile_name = str( + tmpdir_factory.mktemp("test").join("test_invalid_times.nwb") + ) nwbfile = NWBFile( session_description="EcephysSession", identifier="{}".format(739448407), - session_start_time=datetime.now() + session_start_time=datetime.now(), ) nwbfile = write_nwb.add_invalid_times(nwbfile, invalid_epochs) @@ -997,10 +1150,9 @@ def test_add_invalid_times(invalid_epochs, tmpdir_factory): df = nwbfile.invalid_times.to_dataframe() df_in = nwbfile_in.invalid_times.to_dataframe() - pd.testing.assert_frame_equal(df, - df_in, - check_like=True, - check_dtype=False) + pd.testing.assert_frame_equal( + df, df_in, check_like=True, check_dtype=False + ) def test_roundtrip_add_invalid_times(nwbfile, invalid_epochs, roundtripper): @@ -1041,20 +1193,12 @@ def spike_amplitudes(): @pytest.fixture def templates(): - return np.array([ + return np.array( [ - [0, 1, 2], - [0, 1, 2], - [0, 1, 2], - [10, 21, 32] - ], - [ - [0, 1, 2], - [0, 1, 2], - [0, 1, 2], - [15, 9, 4] + [[0, 1, 2], [0, 1, 2], [0, 1, 2], [10, 21, 32]], + [[0, 1, 2], [0, 1, 2], [0, 1, 2], [15, 9, 4]], ] - ]) + ) @pytest.fixture @@ -1068,36 +1212,34 @@ def expected_amplitudes(): def test_scale_amplitudes( - spike_amplitudes, - templates, - spike_templates, - expected_amplitudes): + spike_amplitudes, templates, spike_templates, expected_amplitudes +): scale_factor = 0.195 expected = expected_amplitudes * scale_factor obtained = allensdk.brain_observatory.ecephys.utils.scale_amplitudes( - spike_amplitudes, - templates, - spike_templates, - scale_factor) + spike_amplitudes, templates, spike_templates, scale_factor + ) assert np.allclose(expected, obtained) def test_read_spike_amplitudes_to_dictionary( - tmpdir_factory, - spike_amplitudes, - templates, - spike_templates, - expected_amplitudes): + tmpdir_factory, + spike_amplitudes, + templates, + spike_templates, + expected_amplitudes, +): tmpdir = str(tmpdir_factory.mktemp("spike_amps")) spike_amplitudes_path = os.path.join(tmpdir, "spike_amplitudes.npy") spike_units_path = os.path.join(tmpdir, "spike_units.npy") templates_path = os.path.join(tmpdir, "templates.npy") spike_templates_path = os.path.join(tmpdir, "spike_templates.npy") - inverse_whitening_matrix_path = \ - os.path.join(tmpdir, "inverse_whitening_matrix_path.npy") + inverse_whitening_matrix_path = os.path.join( + tmpdir, "inverse_whitening_matrix_path.npy" + ) whitening_matrix = np.diag(np.arange(3) + 1) inverse_whitening_matrix = np.linalg.inv(whitening_matrix) @@ -1105,149 +1247,213 @@ def test_read_spike_amplitudes_to_dictionary( spike_units = np.array([0, 0, 0, 1, 1]) for idx in range(templates.shape[0]): - templates[idx, :, :] = np.dot( - templates[idx, :, :], whitening_matrix - ) + templates[idx, :, :] = np.dot(templates[idx, :, :], whitening_matrix) np.save(spike_amplitudes_path, spike_amplitudes, allow_pickle=False) np.save(spike_units_path, spike_units, allow_pickle=False) np.save(templates_path, templates, allow_pickle=False) np.save(spike_templates_path, spike_templates, allow_pickle=False) - np.save(inverse_whitening_matrix_path, - inverse_whitening_matrix, - allow_pickle=False) - - obtained = allensdk.brain_observatory.ecephys._units \ - ._read_spike_amplitudes_to_dictionary( - spike_amplitudes_path, - spike_units_path, - templates_path, - spike_templates_path, - inverse_whitening_matrix_path, - scale_factor=1.0 - ) + np.save( + inverse_whitening_matrix_path, + inverse_whitening_matrix, + allow_pickle=False, + ) + + obtained = allensdk.brain_observatory.ecephys._units._read_spike_amplitudes_to_dictionary( # noqa: E501 + spike_amplitudes_path, + spike_units_path, + templates_path, + spike_templates_path, + inverse_whitening_matrix_path, + scale_factor=1.0, + ) assert np.allclose(expected_amplitudes[:3], obtained[0]) assert np.allclose(expected_amplitudes[3:], obtained[1]) -@pytest.mark.parametrize("""spike_times_mapping, - spike_amplitudes_mapping, expected""", [ - - ({12345: np.array([0, 1, 2, -1, 5, 4])}, # spike_times_mapping - - {12345: np.array([0, 1, 2, 3, 4, 5])}, # spike_amplitudes_mapping - - ({12345: np.array([0, 1, 2, 4, 5])}, # expected - {12345: np.array([0, 1, 2, 5, 4])})), - - ({12345: np.array([0, 1, 2, -1, 5, 4]), # spike_times_mapping - 54321: np.array([5, 4, 3, -1, 6])}, - - {12345: np.array([0, 1, 2, 3, 4, 5]), # spike_amplitudes_mapping - 54321: np.array([0, 1, 2, 3, 4])}, - - ({12345: np.array([0, 1, 2, 4, 5]), # expected - 54321: np.array([3, 4, 5, 6])}, - {12345: np.array([0, 1, 2, 5, 4]), - 54321: np.array([2, 1, 0, 4])})), -]) +@pytest.mark.parametrize( + """spike_times_mapping, + spike_amplitudes_mapping, expected""", + [ + ( + {12345: np.array([0, 1, 2, -1, 5, 4])}, # spike_times_mapping + {12345: np.array([0, 1, 2, 3, 4, 5])}, # spike_amplitudes_mapping + ( + {12345: np.array([0, 1, 2, 4, 5])}, # expected + {12345: np.array([0, 1, 2, 5, 4])}, + ), + ), + ( + { + 12345: np.array([0, 1, 2, -1, 5, 4]), # spike_times_mapping + 54321: np.array([5, 4, 3, -1, 6]), + }, + { + 12345: np.array( + [0, 1, 2, 3, 4, 5] + ), # spike_amplitudes_mapping + 54321: np.array([0, 1, 2, 3, 4]), + }, + ( + { + 12345: np.array([0, 1, 2, 4, 5]), # expected + 54321: np.array([3, 4, 5, 6]), + }, + { + 12345: np.array([0, 1, 2, 5, 4]), + 54321: np.array([2, 1, 0, 4]), + }, + ), + ), + ], +) def test_filter_and_sort_spikes( - spike_times_mapping, - spike_amplitudes_mapping, - expected): + spike_times_mapping, spike_amplitudes_mapping, expected +): for unit in spike_times_mapping: expected_spike_times, expected_spike_amplitudes = expected - obtained_spike_times, obtained_spike_amplitudes = \ - _get_filtered_and_sorted_spikes( - spike_times_mapping[unit], spike_amplitudes_mapping[unit]) + ( + obtained_spike_times, + obtained_spike_amplitudes, + ) = _get_filtered_and_sorted_spikes( + spike_times_mapping[unit], spike_amplitudes_mapping[unit] + ) - np.testing.assert_equal(obtained_spike_times, - expected_spike_times[unit]) - np.testing.assert_equal(obtained_spike_amplitudes, - expected_spike_amplitudes[unit]) + np.testing.assert_equal( + obtained_spike_times, expected_spike_times[unit] + ) + np.testing.assert_equal( + obtained_spike_amplitudes, expected_spike_amplitudes[unit] + ) @pytest.mark.parametrize("roundtrip", [True, False]) -@pytest.mark.parametrize("probes, parsed_probe_data", [ - ([{"id": 1234, - "name": "probeA", - "sampling_rate": 29999.9655245905, - "lfp_sampling_rate": np.nan, - "temporal_subsampling_factor": 2.0, - "lfp": None, - "spike_times_path": "/dummy_path", - "spike_clusters_files": "/dummy_path", - "mean_waveforms_path": "/dummy_path", - "channels": [{"id": 1, - "probe_id": 1234, - "valid_data": True, - "probe_channel_number": 0, - "probe_vertical_position": 10, - "probe_horizontal_position": 10, - "anterior_posterior_ccf_coordinate": 15.0, - "dorsal_ventral_ccf_coordinate": 20.0, - "left_right_ccf_coordinate": 25.0, - "structure_acronym": "CA1", - "impedence": np.nan, - "filtering": "Unknown"}, - {"id": 2, - "probe_id": 1234, - "valid_data": True, - "probe_channel_number": 1, - "probe_vertical_position": 20, - "probe_horizontal_position": 20, - "anterior_posterior_ccf_coordinate": 25.0, - "dorsal_ventral_ccf_coordinate": 30.0, - "left_right_ccf_coordinate": 35.0, - "structure_acronym": "CA3", - "impedence": np.nan, - "filtering": "Unknown"}], - - "units": [{"id": 777, - "local_index": 7, - "quality": "good", - "a": 0.5, - "b": 5}, - {"id": 778, - "local_index": 9, - "quality": "noise", - "a": 1.0, - "b": 10}]}], - - (pd.DataFrame({"id": [777, 778], "local_index": [7, 9], # units_table - "a": [0.5, 1.0], "b": [5, 10]}).set_index(keys="id", - drop=True), - {777: np.array([0., 1., 2., -1., 5., 4.]), # spike_times - 778: np.array([5., 4., 3., -1., 6.])}, - {777: np.array([0., 1., 2., 3., 4., 5.]), # spike_amplitudes - 778: np.array([0., 1., 2., 3., 4.])}, - {777: np.array([1., 2., 3., 4., 5., 6.]), # mean_waveforms - 778: np.array([1., 2., 3., 4., 5.])})) -]) -def test_add_probewise_data_to_nwbfile(monkeypatch, nwbfile, roundtripper, - roundtrip, probes, parsed_probe_data): +@pytest.mark.parametrize( + "probes, parsed_probe_data", + [ + ( + [ + { + "id": 1234, + "name": "probeA", + "sampling_rate": 29999.9655245905, + "lfp_sampling_rate": np.nan, + "temporal_subsampling_factor": 2.0, + "lfp": None, + "spike_times_path": "/dummy_path", + "spike_clusters_files": "/dummy_path", + "mean_waveforms_path": "/dummy_path", + "channels": [ + { + "id": 1, + "probe_id": 1234, + "valid_data": True, + "probe_channel_number": 0, + "probe_vertical_position": 10, + "probe_horizontal_position": 10, + "anterior_posterior_ccf_coordinate": 15.0, + "dorsal_ventral_ccf_coordinate": 20.0, + "left_right_ccf_coordinate": 25.0, + "structure_acronym": "CA1", + "impedence": np.nan, + "filtering": "Unknown", + }, + { + "id": 2, + "probe_id": 1234, + "valid_data": True, + "probe_channel_number": 1, + "probe_vertical_position": 20, + "probe_horizontal_position": 20, + "anterior_posterior_ccf_coordinate": 25.0, + "dorsal_ventral_ccf_coordinate": 30.0, + "left_right_ccf_coordinate": 35.0, + "structure_acronym": "CA3", + "impedence": np.nan, + "filtering": "Unknown", + }, + ], + "units": [ + { + "id": 777, + "local_index": 7, + "quality": "good", + "a": 0.5, + "b": 5, + }, + { + "id": 778, + "local_index": 9, + "quality": "noise", + "a": 1.0, + "b": 10, + }, + ], + } + ], + ( + pd.DataFrame( + { + "id": [777, 778], + "local_index": [7, 9], # units_table + "a": [0.5, 1.0], + "b": [5, 10], + } + ).set_index(keys="id", drop=True), + { + 777: np.array( + [0.0, 1.0, 2.0, -1.0, 5.0, 4.0] + ), # spike_times + 778: np.array([5.0, 4.0, 3.0, -1.0, 6.0]), + }, + { + 777: np.array( + [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] + ), # spike_amplitudes + 778: np.array([0.0, 1.0, 2.0, 3.0, 4.0]), + }, + { + 777: np.array( + [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] + ), # mean_waveforms + 778: np.array([1.0, 2.0, 3.0, 4.0, 5.0]), + }, + ), + ) + ], +) +def test_add_probewise_data_to_nwbfile( + monkeypatch, nwbfile, roundtripper, roundtrip, probes, parsed_probe_data +): expected_units_table = pd.read_pickle( - Path(__file__).absolute().parent / 'resources' / - 'expected_units_table.pkl') - - units = Units([Unit( - amplitude_cutoff=1.0, - cluster_id=1, - firing_rate=1.0, - id=unit['id'], - isi_violations=1.0, - local_index=unit['local_index'], - peak_channel_id=1, - presence_ratio=1.0, - quality=unit['quality'], - spike_times=parsed_probe_data[1][unit['id']], - spike_amplitudes=parsed_probe_data[2][unit['id']], - mean_waveforms=parsed_probe_data[3][unit['id']] - ) for unit in probes[0]['units']]) - - with patch.object(Units, 'from_json', return_value=units): + Path(__file__).absolute().parent + / "resources" + / "expected_units_table.pkl" + ) + + units = Units( + [ + Unit( + amplitude_cutoff=1.0, + cluster_id=1, + firing_rate=1.0, + id=unit["id"], + isi_violations=1.0, + local_index=unit["local_index"], + peak_channel_id=1, + presence_ratio=1.0, + quality=unit["quality"], + spike_times=parsed_probe_data[1][unit["id"]], + spike_amplitudes=parsed_probe_data[2][unit["id"]], + mean_waveforms=parsed_probe_data[3][unit["id"]], + ) + for unit in probes[0]["units"] + ] + ) + + with patch.object(Units, "from_json", return_value=units): nwbfile = write_nwb.add_probewise_data_to_nwbfile(nwbfile, probes) if roundtrip: @@ -1255,150 +1461,211 @@ def test_add_probewise_data_to_nwbfile(monkeypatch, nwbfile, roundtripper, else: obt = EcephysNwbSessionApi.from_nwbfile(nwbfile) - pd.testing.assert_frame_equal(obt.nwbfile.units.to_dataframe(), - expected_units_table) + pd.testing.assert_frame_equal( + obt.nwbfile.units.to_dataframe(), expected_units_table + ) @pytest.mark.parametrize("roundtrip", [True, False]) -@pytest.mark.parametrize("eye_tracking_rig_geom, expected", [ - ({"monitor_position_mm": [1., 2., 3.], - "monitor_rotation_deg": [4., 5., 6.], - "camera_position_mm": [7., 8., 9.], - "camera_rotation_deg": [10., 11., 12.], - "led_position": [13., 14., 15.], - "equipment": "test_rig"}, - - # Expected - {"geometry": pd.DataFrame({"monitor_position_mm": [1., 2., 3.], - "monitor_rotation_deg": [4., 5., 6.], - "camera_position_mm": [7., 8., 9.], - "camera_rotation_deg": [10., 11., 12.], - "led_position_mm": [13., 14., 15.]}, - index=["x", "y", "z"]), - "equipment": "test_rig"}), -]) -def test_add_eye_tracking_rig_geometry_data_to_nwbfile(nwbfile, - roundtripper, - roundtrip, - eye_tracking_rig_geom, - expected): - nwbfile = \ - write_nwb.add_eye_tracking_rig_geometry_data_to_nwbfile( - nwbfile, - eye_tracking_rig_geom) +@pytest.mark.parametrize( + "eye_tracking_rig_geom, expected", + [ + ( + { + "monitor_position_mm": [1.0, 2.0, 3.0], + "monitor_rotation_deg": [4.0, 5.0, 6.0], + "camera_position_mm": [7.0, 8.0, 9.0], + "camera_rotation_deg": [10.0, 11.0, 12.0], + "led_position": [13.0, 14.0, 15.0], + "equipment": "test_rig", + }, + # Expected + { + "geometry": pd.DataFrame( + { + "monitor_position_mm": [1.0, 2.0, 3.0], + "monitor_rotation_deg": [4.0, 5.0, 6.0], + "camera_position_mm": [7.0, 8.0, 9.0], + "camera_rotation_deg": [10.0, 11.0, 12.0], + "led_position_mm": [13.0, 14.0, 15.0], + }, + index=["x", "y", "z"], + ), + "equipment": "test_rig", + }, + ), + ], +) +def test_add_eye_tracking_rig_geometry_data_to_nwbfile( + nwbfile, roundtripper, roundtrip, eye_tracking_rig_geom, expected +): + nwbfile = write_nwb.add_eye_tracking_rig_geometry_data_to_nwbfile( + nwbfile, eye_tracking_rig_geom + ) if roundtrip: obt = roundtripper(nwbfile, EcephysNwbSessionApi) else: obt = EcephysNwbSessionApi.from_nwbfile(nwbfile) obtained_metadata = obt.get_rig_metadata() - pd.testing.assert_frame_equal(obtained_metadata["geometry"], - expected["geometry"], - check_like=True) + pd.testing.assert_frame_equal( + obtained_metadata["geometry"], expected["geometry"], check_like=True + ) assert obtained_metadata["equipment"] == expected["equipment"] @pytest.mark.parametrize("roundtrip", [True, False]) -@pytest.mark.parametrize(("eye_tracking_frame_times, eye_dlc_tracking_data, " - "eye_gaze_data, expected_pupil_data, " - "expected_gaze_data"), [ +@pytest.mark.parametrize( ( - # eye_tracking_frame_times - pd.Series([3., 4., 5., 6., 7.]), - # eye_dlc_tracking_data - {"pupil_params": create_preload_eye_tracking_df(np.full((5, 5), 1.)), - "cr_params": create_preload_eye_tracking_df(np.full((5, 5), 2.)), - "eye_params": create_preload_eye_tracking_df(np.full((5, 5), 3.))}, - # eye_gaze_data - {"raw_pupil_areas": pd.Series([2., 4., 6., 8., 10.]), - "raw_eye_areas": pd.Series([3., 5., 7., 9., 11.]), - "raw_screen_coordinates": pd.DataFrame( - {"y": [2., 4., 6., 8., 10.], - "x": [3., 5., 7., 9., 11.]}), - "raw_screen_coordinates_spherical": pd.DataFrame( - {"y": [2., 4., 6., 8., 10.], - "x": [3., 5., 7., 9., 11.]}), - "new_pupil_areas": pd.Series([2., 4., np.nan, 8., 10.]), - "new_eye_areas": pd.Series([3., 5., np.nan, 9., 11.]), - "new_screen_coordinates": pd.DataFrame( - {"y": [2., 4., np.nan, 8., 10.], - "x": [3., 5., np.nan, 9., 11.]}), - "new_screen_coordinates_spherical": pd.DataFrame({ - "y": [2., 4., np.nan, 8., 10.], - "x": [3., 5., np.nan, 9., 11.]}), - "synced_frame_timestamps": pd.Series([3., 4., 5., 6., 7.])}, - # expected_pupil_data - pd.DataFrame({"corneal_reflection_center_x": [2.] * 5, - "corneal_reflection_center_y": [2.] * 5, - "corneal_reflection_height": [4.] * 5, - "corneal_reflection_width": [4.] * 5, - "corneal_reflection_phi": [2.] * 5, - "pupil_center_x": [1.] * 5, - "pupil_center_y": [1.] * 5, - "pupil_height": [2.] * 5, - "pupil_width": [2.] * 5, - "pupil_phi": [1.] * 5, - "eye_center_x": [3.] * 5, - "eye_center_y": [3.] * 5, - "eye_height": [6.] * 5, - "eye_width": [6.] * 5, - "eye_phi": [3.] * 5}, - index=pd.Index(name="Time (s)", - data=[3., 4., 5., 6., 7.])), - # expected_gaze_data - pd.DataFrame({"raw_eye_area": [3., 5., 7., 9., 11.], - "raw_pupil_area": [2., 4., 6., 8., 10.], - "raw_screen_coordinates_x_cm": [3., 5., 7., 9., 11.], - "raw_screen_coordinates_y_cm": [2., 4., 6., 8., 10.], - "raw_screen_coordinates_spherical_x_deg": [3., - 5., - 7., - 9., - 11.], - "raw_screen_coordinates_spherical_y_deg": [2., - 4., - 6., - 8., - 10.], - "filtered_eye_area": [3., 5., np.nan, 9., 11.], - "filtered_pupil_area": [2., 4., np.nan, 8., 10.], - "filtered_screen_coordinates_x_cm": [3., - 5., - np.nan, - 9., - 11.], - "filtered_screen_coordinates_y_cm": [2., - 4., - np.nan, - 8., - 10.], - "filtered_screen_coordinates_spherical_x_deg": [3., - 5., - np.nan, - 9., - 11.], - "filtered_screen_coordinates_spherical_y_deg": [2., - 4., - np.nan, - 8., - 10.]}, - index=pd.Index(name="Time (s)", - data=[3., 4., 5., 6., 7.])) + "eye_tracking_frame_times, eye_dlc_tracking_data, " + "eye_gaze_data, expected_pupil_data, " + "expected_gaze_data" ), -]) -def test_add_eye_tracking_data_to_nwbfile(nwbfile, - roundtripper, - roundtrip, - eye_tracking_frame_times, - eye_dlc_tracking_data, - eye_gaze_data, - expected_pupil_data, - expected_gaze_data): - nwbfile = \ - write_nwb.add_eye_tracking_data_to_nwbfile(nwbfile, - eye_tracking_frame_times, - eye_dlc_tracking_data, - eye_gaze_data) + [ + ( + # eye_tracking_frame_times + pd.Series([3.0, 4.0, 5.0, 6.0, 7.0]), + # eye_dlc_tracking_data + { + "pupil_params": create_preload_eye_tracking_df( + np.full((5, 5), 1.0) + ), + "cr_params": create_preload_eye_tracking_df( + np.full((5, 5), 2.0) + ), + "eye_params": create_preload_eye_tracking_df( + np.full((5, 5), 3.0) + ), + }, + # eye_gaze_data + { + "raw_pupil_areas": pd.Series([2.0, 4.0, 6.0, 8.0, 10.0]), + "raw_eye_areas": pd.Series([3.0, 5.0, 7.0, 9.0, 11.0]), + "raw_screen_coordinates": pd.DataFrame( + { + "y": [2.0, 4.0, 6.0, 8.0, 10.0], + "x": [3.0, 5.0, 7.0, 9.0, 11.0], + } + ), + "raw_screen_coordinates_spherical": pd.DataFrame( + { + "y": [2.0, 4.0, 6.0, 8.0, 10.0], + "x": [3.0, 5.0, 7.0, 9.0, 11.0], + } + ), + "new_pupil_areas": pd.Series([2.0, 4.0, np.nan, 8.0, 10.0]), + "new_eye_areas": pd.Series([3.0, 5.0, np.nan, 9.0, 11.0]), + "new_screen_coordinates": pd.DataFrame( + { + "y": [2.0, 4.0, np.nan, 8.0, 10.0], + "x": [3.0, 5.0, np.nan, 9.0, 11.0], + } + ), + "new_screen_coordinates_spherical": pd.DataFrame( + { + "y": [2.0, 4.0, np.nan, 8.0, 10.0], + "x": [3.0, 5.0, np.nan, 9.0, 11.0], + } + ), + "synced_frame_timestamps": pd.Series( + [3.0, 4.0, 5.0, 6.0, 7.0] + ), + }, + # expected_pupil_data + pd.DataFrame( + { + "corneal_reflection_center_x": [2.0] * 5, + "corneal_reflection_center_y": [2.0] * 5, + "corneal_reflection_height": [4.0] * 5, + "corneal_reflection_width": [4.0] * 5, + "corneal_reflection_phi": [2.0] * 5, + "pupil_center_x": [1.0] * 5, + "pupil_center_y": [1.0] * 5, + "pupil_height": [2.0] * 5, + "pupil_width": [2.0] * 5, + "pupil_phi": [1.0] * 5, + "eye_center_x": [3.0] * 5, + "eye_center_y": [3.0] * 5, + "eye_height": [6.0] * 5, + "eye_width": [6.0] * 5, + "eye_phi": [3.0] * 5, + }, + index=pd.Index( + name="Time (s)", data=[3.0, 4.0, 5.0, 6.0, 7.0] + ), + ), + # expected_gaze_data + pd.DataFrame( + { + "raw_eye_area": [3.0, 5.0, 7.0, 9.0, 11.0], + "raw_pupil_area": [2.0, 4.0, 6.0, 8.0, 10.0], + "raw_screen_coordinates_x_cm": [3.0, 5.0, 7.0, 9.0, 11.0], + "raw_screen_coordinates_y_cm": [2.0, 4.0, 6.0, 8.0, 10.0], + "raw_screen_coordinates_spherical_x_deg": [ + 3.0, + 5.0, + 7.0, + 9.0, + 11.0, + ], + "raw_screen_coordinates_spherical_y_deg": [ + 2.0, + 4.0, + 6.0, + 8.0, + 10.0, + ], + "filtered_eye_area": [3.0, 5.0, np.nan, 9.0, 11.0], + "filtered_pupil_area": [2.0, 4.0, np.nan, 8.0, 10.0], + "filtered_screen_coordinates_x_cm": [ + 3.0, + 5.0, + np.nan, + 9.0, + 11.0, + ], + "filtered_screen_coordinates_y_cm": [ + 2.0, + 4.0, + np.nan, + 8.0, + 10.0, + ], + "filtered_screen_coordinates_spherical_x_deg": [ + 3.0, + 5.0, + np.nan, + 9.0, + 11.0, + ], + "filtered_screen_coordinates_spherical_y_deg": [ + 2.0, + 4.0, + np.nan, + 8.0, + 10.0, + ], + }, + index=pd.Index( + name="Time (s)", data=[3.0, 4.0, 5.0, 6.0, 7.0] + ), + ), + ), + ], +) +def test_add_eye_tracking_data_to_nwbfile( + nwbfile, + roundtripper, + roundtrip, + eye_tracking_frame_times, + eye_dlc_tracking_data, + eye_gaze_data, + expected_pupil_data, + expected_gaze_data, +): + nwbfile = write_nwb.add_eye_tracking_data_to_nwbfile( + nwbfile, eye_tracking_frame_times, eye_dlc_tracking_data, eye_gaze_data + ) if roundtrip: obt = roundtripper(nwbfile, EcephysNwbSessionApi) @@ -1409,7 +1676,9 @@ def test_add_eye_tracking_data_to_nwbfile(nwbfile, include_filtered_data=True ) - pd.testing.assert_frame_equal(obtained_pupil_data, - expected_pupil_data, check_like=True) - pd.testing.assert_frame_equal(obtained_screen_gaze_data, - expected_gaze_data, check_like=True) + pd.testing.assert_frame_equal( + obtained_pupil_data, expected_pupil_data, check_like=True + ) + pd.testing.assert_frame_equal( + obtained_screen_gaze_data, expected_gaze_data, check_like=True + ) diff --git a/allensdk/test/core/test_datafame_utils.py b/allensdk/test/core/test_datafame_utils.py index b7629c4d1..d330ac33d 100644 --- a/allensdk/test/core/test_datafame_utils.py +++ b/allensdk/test/core/test_datafame_utils.py @@ -1,9 +1,8 @@ -import pytest import pandas as pd - +import pytest from allensdk.core.dataframe_utils import ( patch_df_from_other, - return_one_dataframe_row_only + return_one_dataframe_row_only, ) @@ -13,10 +12,10 @@ def target_df_fixture(): Return an example target dataframe with columns 'a', 'b', 'c', 'd' """ data = [ - {'a': 1, 'b': 3.4, 'c': 'apple', 'd': None}, - {'a': 9, 'b': 4.5, 'c': 'banana', 'd': 4.6}, - {'a': 12, 'b': 7.8, 'c': 'pineapple', 'd': 'purple'}, - {'a': 17, 'b': None, 'c': 'papaya', 'd': 11} + {"a": 1, "b": 3.4, "c": "apple", "d": None}, + {"a": 9, "b": 4.5, "c": "banana", "d": 4.6}, + {"a": 12, "b": 7.8, "c": "pineapple", "d": "purple"}, + {"a": 17, "b": None, "c": "papaya", "d": 11}, ] return pd.DataFrame(data=data) @@ -27,77 +26,75 @@ def source_df_fixture(): Return an example source dataframe with columns 'a', 'b', 'd', 'e' """ data = [ - {'a': 1, 'b': 7.5, 'd': 'tree', 'e': 'frog'}, - {'a': 12, 'b': 88.9, 'd': None, 'e': 'dog'}, - {'a': 17, 'b': 35.2, 'd': 17, 'e': 'cat'} + {"a": 1, "b": 7.5, "d": "tree", "e": "frog"}, + {"a": 12, "b": 88.9, "d": None, "e": "dog"}, + {"a": 17, "b": 35.2, "d": 17, "e": "cat"}, ] return pd.DataFrame(data=data) -def test_patch_error_on_no_index( - target_df_fixture, - source_df_fixture): +def test_patch_error_on_no_index(target_df_fixture, source_df_fixture): """ Test than an exception is raise when you specify an index_column that does not exist in one or the other dataframe """ - with pytest.raises(ValueError, match='not in target_df'): + with pytest.raises(ValueError, match="not in target_df"): patch_df_from_other( source_df=source_df_fixture, target_df=target_df_fixture, - index_column='e', - columns_to_patch=['a', 'b']) + index_column="e", + columns_to_patch=["a", "b"], + ) - with pytest.raises(ValueError, match='not in source_df'): + with pytest.raises(ValueError, match="not in source_df"): patch_df_from_other( source_df=source_df_fixture, target_df=target_df_fixture, - index_column='c', - columns_to_patch=['a', 'b']) + index_column="c", + columns_to_patch=["a", "b"], + ) -def test_patch_error_on_missing_col( - source_df_fixture, - target_df_fixture): +def test_patch_error_on_missing_col(source_df_fixture, target_df_fixture): """ Test that an error is raised when you specify a column_to_patch that is not in the source_df """ - with pytest.raises(ValueError, match='not in source_df'): + with pytest.raises(ValueError, match="not in source_df"): patch_df_from_other( source_df=source_df_fixture, target_df=target_df_fixture, - index_column='a', - columns_to_patch=['c']) + index_column="a", + columns_to_patch=["c"], + ) -def test_error_on_not_unique_index( - target_df_fixture): +def test_error_on_not_unique_index(target_df_fixture): """ Test that an exception is raised when the values of index_column in source_df are not unique """ data = [ - {'a': 1, 'b': 7.5, 'e': 'frog'}, - {'a': 12, 'b': 88.9, 'e': 'dog'}, - {'a': 17, 'b': 35.2, 'e': 'cat'}, - {'a': 17, 'b': 55.7, 'e': 'mosquito'} + {"a": 1, "b": 7.5, "e": "frog"}, + {"a": 12, "b": 88.9, "e": "dog"}, + {"a": 17, "b": 35.2, "e": "cat"}, + {"a": 17, "b": 55.7, "e": "mosquito"}, ] source_df = pd.DataFrame(data=data) - with pytest.raises(ValueError, match='in source_df are not unique'): + with pytest.raises(ValueError, match="in source_df are not unique"): patch_df_from_other( source_df=source_df, target_df=target_df_fixture, - index_column='a', - columns_to_patch=['c']) + index_column="a", + columns_to_patch=["c"], + ) -@pytest.mark.parametrize('original_index', [None, 'a', 'c', 'b']) +@pytest.mark.parametrize("original_index", [None, "a", "c", "b"]) def test_patch_no_duplicates( - source_df_fixture, - target_df_fixture, - original_index): + source_df_fixture, target_df_fixture, original_index +): """ Test that we get the expected dataframe back in the case where there are no duplicate values of index_column @@ -109,10 +106,10 @@ def test_patch_no_duplicates( """ expected_data = [ - {'a': 1, 'b': 7.5, 'c': 'apple', 'd': 'tree'}, - {'a': 9, 'b': 4.5, 'c': 'banana', 'd': 4.6}, - {'a': 12, 'b': 88.9, 'c': 'pineapple', 'd': 'purple'}, - {'a': 17, 'b': 35.2, 'c': 'papaya', 'd': 17} + {"a": 1, "b": 7.5, "c": "apple", "d": "tree"}, + {"a": 9, "b": 4.5, "c": "banana", "d": 4.6}, + {"a": 12, "b": 88.9, "c": "pineapple", "d": "purple"}, + {"a": 17, "b": 35.2, "c": "papaya", "d": 17}, ] expected_df = pd.DataFrame(data=expected_data) @@ -124,19 +121,19 @@ def test_patch_no_duplicates( target_df = target_df_fixture.copy(deep=True) actual = patch_df_from_other( - source_df=source_df_fixture.copy(deep=True), - target_df=target_df, - index_column='a', - columns_to_patch=['b', 'd']) + source_df=source_df_fixture.copy(deep=True), + target_df=target_df, + index_column="a", + columns_to_patch=["b", "d"], + ) pd.testing.assert_frame_equal(actual, expected_df) -@pytest.mark.parametrize('original_index', [None, 'c', 'b']) +@pytest.mark.parametrize("original_index", [None, "c", "b"]) def test_patch_with_duplicates( - source_df_fixture, - target_df_fixture, - original_index): + source_df_fixture, target_df_fixture, original_index +): """ Test that we get the expected dataframe back in the case where there are duplicate values of index_column @@ -148,20 +145,20 @@ def test_patch_with_duplicates( """ data = [ - {'a': 1, 'b': 3.4, 'c': 'apple', 'd': None}, - {'a': 9, 'b': 4.5, 'c': 'banana', 'd': 4.6}, - {'a': 17, 'b': 11.3, 'c': 'tomato', 'd': 'blue'}, - {'a': 12, 'b': 7.8, 'c': 'pineapple', 'd': 'purple'}, - {'a': 17, 'b': None, 'c': 'papaya', 'd': 11} + {"a": 1, "b": 3.4, "c": "apple", "d": None}, + {"a": 9, "b": 4.5, "c": "banana", "d": 4.6}, + {"a": 17, "b": 11.3, "c": "tomato", "d": "blue"}, + {"a": 12, "b": 7.8, "c": "pineapple", "d": "purple"}, + {"a": 17, "b": None, "c": "papaya", "d": 11}, ] target_df = pd.DataFrame(data=data) expected_data = [ - {'a': 1, 'b': 7.5, 'c': 'apple', 'd': 'tree'}, - {'a': 9, 'b': 4.5, 'c': 'banana', 'd': 4.6}, - {'a': 17, 'b': 35.2, 'c': 'tomato', 'd': 17}, - {'a': 12, 'b': 88.9, 'c': 'pineapple', 'd': 'purple'}, - {'a': 17, 'b': 35.2, 'c': 'papaya', 'd': 17} + {"a": 1, "b": 7.5, "c": "apple", "d": "tree"}, + {"a": 9, "b": 4.5, "c": "banana", "d": 4.6}, + {"a": 17, "b": 35.2, "c": "tomato", "d": 17}, + {"a": 12, "b": 88.9, "c": "pineapple", "d": "purple"}, + {"a": 17, "b": 35.2, "c": "papaya", "d": 17}, ] expected_df = pd.DataFrame(data=expected_data) @@ -170,82 +167,90 @@ def test_patch_with_duplicates( target_df = target_df.set_index(original_index) actual = patch_df_from_other( - source_df=source_df_fixture.copy(deep=True), - target_df=target_df, - index_column='a', - columns_to_patch=['b', 'd']) + source_df=source_df_fixture.copy(deep=True), + target_df=target_df, + index_column="a", + columns_to_patch=["b", "d"], + ) pd.testing.assert_frame_equal(actual, expected_df) -def test_patch_new_column( - target_df_fixture, - source_df_fixture): +def test_patch_new_column(target_df_fixture, source_df_fixture): """ Test case where we use patch_df_from_other in the case where we are adding a column to target_df """ expected_data = [ - {'a': 1, 'b': 3.4, 'c': 'apple', 'd': None, 'e': 'frog'}, - {'a': 9, 'b': 4.5, 'c': 'banana', 'd': 4.6, 'e': None}, - {'a': 12, 'b': 7.8, 'c': 'pineapple', 'd': 'purple', 'e': 'dog'}, - {'a': 17, 'b': None, 'c': 'papaya', 'd': 11, 'e': 'cat'} + {"a": 1, "b": 3.4, "c": "apple", "d": None, "e": "frog"}, + {"a": 9, "b": 4.5, "c": "banana", "d": 4.6, "e": None}, + {"a": 12, "b": 7.8, "c": "pineapple", "d": "purple", "e": "dog"}, + {"a": 17, "b": None, "c": "papaya", "d": 11, "e": "cat"}, ] expected_df = pd.DataFrame(data=expected_data) actual = patch_df_from_other( - source_df=source_df_fixture.copy(deep=True), - target_df=target_df_fixture, - index_column='a', - columns_to_patch=['e']) + source_df=source_df_fixture.copy(deep=True), + target_df=target_df_fixture, + index_column="a", + columns_to_patch=["e"], + ) pd.testing.assert_frame_equal(actual, expected_df) def test_row_index_not_in_dataframe_error(): """Test for the correct failure row index not in dataframe""" - mock_behavior_sessions = pd.DataFrame({ - "behavior_session_id": [1, 2, 3, 4], - "ecephys_session_id": pd.Series( - [10, 11, 12, 13], - dtype='Int64'), - "mouse_id": [4, 4, 2, 1]}).set_index("behavior_session_id") + mock_behavior_sessions = pd.DataFrame( + { + "behavior_session_id": [1, 2, 3, 4], + "ecephys_session_id": pd.Series([10, 11, 12, 13], dtype="Int64"), + "mouse_id": [4, 4, 2, 1], + } + ).set_index("behavior_session_id") table_name = "behavior_session_table" index_name = "behavior_session_id" session_id = 1234 with pytest.raises( - RuntimeError, - match=f"The {table_name} should have " - "1 and only 1 entry for a given " - f"{index_name}. No indexed rows found for " - f"id={session_id}" + RuntimeError, + match=f"The {table_name} should have " + "1 and only 1 entry for a given " + f"{index_name}. No indexed rows found for " + f"id={session_id}", ): - return_one_dataframe_row_only(input_table=mock_behavior_sessions, - index_value=session_id, - table_name=table_name) + return_one_dataframe_row_only( + input_table=mock_behavior_sessions, + index_value=session_id, + table_name=table_name, + ) def test_multiple_indexes_in_dataframe(): """Test for the correct failure on multiple rows with same index dataframe""" - mock_behavior_sessions = pd.DataFrame({ - "behavior_session_id": [1, 2, 2, 3, 4], - "ecephys_session_id": pd.Series( - [10, 11, 0, 12, 13], - dtype='Int64'), - "mouse_id": [4, 4, 4, 2, 1]}).set_index("behavior_session_id") + mock_behavior_sessions = pd.DataFrame( + { + "behavior_session_id": [1, 2, 2, 3, 4], + "ecephys_session_id": pd.Series( + [10, 11, 0, 12, 13], dtype="Int64" + ), + "mouse_id": [4, 4, 4, 2, 1], + } + ).set_index("behavior_session_id") table_name = "behavior_session_table" index_name = "behavior_session_id" session_id = 2 n_rows = len(mock_behavior_sessions.loc[session_id]) with pytest.raises( - RuntimeError, - match=f"The {table_name} should have " - "1 and only 1 entry for a given " - f"{index_name}. For " - f"{session_id} " - f" there are {n_rows} entries." + RuntimeError, + match=f"The {table_name} should have " + "1 and only 1 entry for a given " + f"{index_name}. For " + f"{session_id} " + f" there are {n_rows} entries.", ): - return_one_dataframe_row_only(input_table=mock_behavior_sessions, - index_value=session_id, - table_name=table_name) \ No newline at end of file + return_one_dataframe_row_only( + input_table=mock_behavior_sessions, + index_value=session_id, + table_name=table_name, + )