diff --git a/docs/MIGRATION.md b/docs/MIGRATION.md index 473df79..bd61491 100644 --- a/docs/MIGRATION.md +++ b/docs/MIGRATION.md @@ -47,7 +47,7 @@ and **Plan Analysis** v0.5. | `PsychometricFitConfig` | Lookup | One versioned psychometric eligibility configuration | none | minimum choices/stimulus values, analysis version | | `PsychometricSessionFit` | Computed | One fit for one upstream TrialSet and config | `DecisionTask.TrialSet` + config | status, fit sample size, curve/parameters/diagnostics | | `PsychometricSubjectFit` | Computed | One pooled fit for one analysis set, subject, condition, and config | analysis set + `Subject` + config | status, fit sample size, curve/parameters/diagnostics | -| `PsychophysicalKernelFitConfig` | Lookup | One versioned pooled-kernel configuration | none | binning method, observation window, evidence encoding, CV, regularization, version | +| `PsychophysicalKernelFitConfig` | Lookup | One integer-keyed pooled-kernel configuration | none | binning method, observation window, evidence model, CV, regularization | | `PsychophysicalKernel` | Computed | One pooled kernel for one analysis set, subject, condition, and config | analysis set + `Subject` + kernel config | status, timing source, fit sample size, weights, held-out scores, bias | Keep: @@ -133,9 +133,12 @@ subpage. ## Kernel timing-source extension -The canonical pooled kernel now supports fixed 100 ms windows in addition to -the migrated variable-width `v1_10bin_10fold` fit. For fixed-window configs, -each selected session: +The canonical pooled kernel uses fixed 100 ms bins and the Odoemene et al. +(2018) Equation 5 evidence model: each bin's flash count is centered on that +trial's generative stimulus rate, which is also included as a separate nuisance +regressor. Config IDs are opaque integers: config 0 uses the center-exit window +and config 1 uses the response window. The lookup parameters, not the ID, +define each analysis. For each selected session: 1. uses NIDAQ/OneBox visual-flash and port-event times when a `visual_stim` `EventMapping` exists and all required mappings validate; @@ -144,6 +147,10 @@ each selected session: 3. raises on incomplete or invalid NIDAQ mappings rather than silently falling back to Bpod. +Within a NIDAQ-timed session, mapped flashes and port edges must agree with the +Bpod trial alignment within 100 ms. Trials with missing or misaligned hardware +events are excluded rather than replaced with Bpod timestamps. + `PsychophysicalKernel.timing_source` records `nidaq`, `bpod`, or `mixed` for the pooled result. Detailed session provenance remains derivable from `BehaviorAnalysisSet.TrialSet` and the ephys `EventMapping`; no duplicate @@ -156,7 +163,6 @@ uv run python scripts/analyses/migrate_kernel_timing_source.py uv run python scripts/analyses/migrate_kernel_timing_source.py --apply ``` -The first command is read-only and prints the missing `ALTER TABLE` statements. -The second adds the new nullable/defaulted columns and inserts the two v2 -fixed-window configs while preserving all existing v1 rows as Bpod-timed -legacy fits. +The first command is read-only and reports the rows that will be removed. The +second recreates the disposable kernel result and config tables with the +canonical schema. Existing kernel rows are not migrated. diff --git a/labdata_plugin/analysisschema.py b/labdata_plugin/analysisschema.py index e913984..ccf96e8 100644 --- a/labdata_plugin/analysisschema.py +++ b/labdata_plugin/analysisschema.py @@ -148,65 +148,45 @@ def make(self, key): @rojasbowe_schema class PsychophysicalKernelFitConfig(dj.Lookup): - """Versioned settings for pooled psychophysical-kernel fits.""" + """Settings for pooled psychophysical-kernel fits.""" definition = """ - kernel_fit_config_id : varchar(48) + kernel_fit_config_id : int --- timebins : int + binning_method : enum('fixed_width') + bin_width_s : float + observation_window : enum('center_exit', 'response') + evidence_model : enum('trial_rate_residual') # Odoemene Eq. 5 + min_trials_per_bin : int cv_splits : int random_state : int - max_rate_hz : float # calibration rate regularization_c : float - kernel_method = 'legacy_variable' : enum('legacy_variable', 'fixed_window') - bin_width_s = NULL : float - evidence_encoding = NULL : enum('max_rate', 'trial_rate') - min_trials_per_bin = NULL : int - observation_window = NULL : enum('center_exit', 'response') - analysis_version : varchar(32) """ contents = [ # noqa: RUF012 ( - "v1_10bin_10fold", + 0, 10, + "fixed_width", + 0.1, + "center_exit", + "trial_rate_residual", + 50, 10, 0, - 20.0, 1.0, - "legacy_variable", - None, - None, - None, - None, - "v1", ), ( - "v2_100ms_10bin_center_rate", - 10, + 1, 10, - 0, - 20.0, - 1.0, - "fixed_window", + "fixed_width", 0.1, - "trial_rate", + "response", + "trial_rate_residual", 50, - "center_exit", - "v2", - ), - ( - "v2_100ms_10bin_response_rate", - 10, 10, 0, - 20.0, 1.0, - "fixed_window", - 0.1, - "trial_rate", - 50, - "response", - "v2", ), ] @@ -251,10 +231,7 @@ def key_source(self): def make(self, key): config = (PsychophysicalKernelFitConfig() & key).fetch1() trialset_keys = _selected_trialset_keys(key) - if config["kernel_method"] == "legacy_variable": - payload = _legacy_kernel_payload(key, config, trialset_keys) - else: - payload = _fixed_kernel_payload(key, config, trialset_keys) + payload = _kernel_payload(key, config, trialset_keys) self.insert1( { **key, @@ -284,73 +261,11 @@ def _fetch_trialset_rows_for_subject(key): return list((DecisionTask.TrialSet() & trialset_keys).fetch(as_dict=True)) -def _legacy_kernel_payload(key, config, trialset_keys): - from behavior_analyses.io import get_chipmunk_table - from behavior_analyses.kernels import ( - fit_psychophysical_kernel, - interpret_kernel_weights, - ) - - Chipmunk = get_chipmunk_table() - relation = ( - Chipmunk() * Chipmunk.Trial() * Chipmunk.TrialParameters() - & trialset_keys - & {"rewarded_modality": key["trialset_description"]} - ) - stim_events, response_values = relation.fetch("stim_events", "response") - result = fit_psychophysical_kernel( - stim_events, - response_values, - timebins=int(config["timebins"]), - cv_splits=int(config["cv_splits"]), - random_state=int(config["random_state"]), - max_rate_hz=float(config["max_rate_hz"]), - regularization_c=float(config["regularization_c"]), - ) - n_trials_fit = int(result["choice_right"].size) - if result["weights"].size == 0: - return { - "fit_status": "skipped", - "fit_message": "insufficient trials or response classes for CV", - "timing_source": "bpod", - "n_trials_fit": n_trials_fit, - "n_bins_fit": 0, - } - - n_observed = np.full(int(config["timebins"]), n_trials_fit, dtype=int) - score_mean = float(np.mean(result["scores"])) - majority_accuracy = float( - max(np.mean(result["choice_right"]), 1.0 - np.mean(result["choice_right"])) - ) - weights_mean = np.mean(result["weights"], axis=0) - return { - "fit_status": "fit", - "timing_source": "bpod", - "n_trials_fit": n_trials_fit, - "n_bins_fit": int(config["timebins"]), - "n_observed_per_bin": n_observed, - "weights": result["weights"], - "weights_mean": weights_mean, - "weights_error": np.mean(result["error"], axis=0), - "scores": result["scores"], - "score_mean": score_mean, - "majority_accuracy": majority_accuracy, - "score_above_majority": score_mean - majority_accuracy, - "bias": result["bias"], - "bias_mean": float(np.mean(result["bias"])), - "interpretation": interpret_kernel_weights( - weights_mean, - n_observed, - min_trials_per_bin=int(config["cv_splits"]), - ), - } - - -def _fixed_kernel_payload(key, config, trialset_keys): +def _kernel_payload(key, config, trialset_keys): from behavior_analyses.kernel_timing import fetch_pooled_kernel_inputs from behavior_analyses.kernels import ( - build_fixed_residual_rate_matrix, - fit_fixed_psychophysical_kernel, + build_residual_rate_matrix, + fit_psychophysical_kernel, interpret_kernel_weights, ) @@ -359,22 +274,18 @@ def _fixed_kernel_payload(key, config, trialset_keys): key["trialset_description"], observation_window=str(config["observation_window"]), ) - trial_rate_hz = ( - inputs["trial_rate_hz"] if config["evidence_encoding"] == "trial_rate" else None - ) residual, choices, n_observed, bin_centers, expected_counts = ( - build_fixed_residual_rate_matrix( + build_residual_rate_matrix( inputs["stim_times_per_trial"], inputs["first_stim_times"], inputs["observation_end_times"], inputs["response_values"], timebins=int(config["timebins"]), bin_width_s=float(config["bin_width_s"]), - max_rate_hz=float(config["max_rate_hz"]), - trial_rate_hz=trial_rate_hz, + trial_rate_hz=inputs["trial_rate_hz"], ) ) - result = fit_fixed_psychophysical_kernel( + result = fit_psychophysical_kernel( residual, choices, n_observed_per_bin=n_observed, @@ -382,9 +293,7 @@ def _fixed_kernel_payload(key, config, trialset_keys): random_state=int(config["random_state"]), min_trials_per_bin=int(config["min_trials_per_bin"]), regularization_c=float(config["regularization_c"]), - expected_counts=( - expected_counts if config["evidence_encoding"] == "trial_rate" else None - ), + expected_counts=expected_counts, ) base = { "timing_source": inputs["timing_source"], diff --git a/scripts/analyses/migrate_behavior_analysis_schema.py b/scripts/analyses/migrate_behavior_analysis_schema.py index 2ae40c6..5408340 100644 --- a/scripts/analyses/migrate_behavior_analysis_schema.py +++ b/scripts/analyses/migrate_behavior_analysis_schema.py @@ -30,7 +30,6 @@ "__psychometric_subject_fit", "__psychophysical_kernel", } -EXPECTED_KERNEL_CONFIG = (10, 10, 0) TRIALSET_KEY_FIELDS = ( "subject_name", "session_name", @@ -64,11 +63,9 @@ def main() -> None: database = f"{labdata_schema.dbase_name}_user" if args.resume: _validate_resume_state(connection, database) - _validate_kernel_configs(connection, database, archived=True) _print_archive_counts(connection, database) else: _validate_table_state(connection, database) - _validate_kernel_configs(connection, database) _print_source_counts(connection, database) if not args.apply: print("Dry run only. Re-run with --apply after exact live-write approval.") @@ -86,14 +83,9 @@ def main() -> None: PsychometricSessionFit, PsychometricSubjectFit, PsychophysicalKernel, - PsychophysicalKernelFitConfig, ) PsychometricFitConfig.insert1(("v1", 100, 6, "v1"), skip_duplicates=True) - PsychophysicalKernelFitConfig.insert1( - ("v1_10bin_10fold", 10, 10, 0, 20.0, 1.0, "v1"), - skip_duplicates=True, - ) old_master = _archive_table(connection, database, "behavior_session_set") BehaviorAnalysisSet.insert( @@ -166,35 +158,6 @@ def main() -> None: allow_direct_insert=True, ) - old_kernels = _archive_table(connection, database, "__psychophysical_kernel") - PsychophysicalKernel.insert( - [ - { - "analysis_set_id": row["session_set_id"], - "subject_name": row["subject_name"], - "trialset_description": row["trialset_description"], - "kernel_fit_config_id": "v1_10bin_10fold", - "fit_status": "fit", - "n_trials_fit": row["n_trials"], - **{ - field: row[field] - for field in ( - "weights", - "weights_mean", - "weights_error", - "scores", - "score_mean", - "bias", - "bias_mean", - ) - }, - } - for row in old_kernels.fetch(as_dict=True) - ], - skip_duplicates=True, - allow_direct_insert=True, - ) - print(f"BehaviorAnalysisSet: {len(BehaviorAnalysisSet())}") print(f"BehaviorAnalysisSet.TrialSet: {len(BehaviorAnalysisSet.TrialSet())}") print(f"PsychometricSessionFit copied: {len(PsychometricSessionFit())}") @@ -241,35 +204,6 @@ def _validate_resume_state(connection, database): ) -def _validate_kernel_configs(connection, database, *, archived=False): - master = ( - ARCHIVE_TABLES["behavior_session_set"] if archived else "behavior_session_set" - ) - kernel = ( - ARCHIVE_TABLES["__psychophysical_kernel"] - if archived - else "__psychophysical_kernel" - ) - queries = { - master: ( - "SELECT DISTINCT kernel_timebins, kernel_cv_splits, " - f"kernel_random_state FROM `{database}`.`{master}`" - ), - kernel: ( - "SELECT DISTINCT timebins, cv_splits, random_state " - f"FROM `{database}`.`{kernel}`" - ), - } - incompatible = {} - for table, query in queries.items(): - configs = {tuple(map(int, row)) for row in connection.query(query).fetchall()} - unexpected = configs - {EXPECTED_KERNEL_CONFIG} - if unexpected: - incompatible[table] = sorted(unexpected) - if incompatible: - raise RuntimeError(f"Incompatible legacy kernel configs: {incompatible}") - - def _print_source_counts(connection, database): for source, archive in ARCHIVE_TABLES.items(): count = connection.query( diff --git a/scripts/analyses/migrate_kernel_timing_source.py b/scripts/analyses/migrate_kernel_timing_source.py index 8e9da67..5c1fe14 100644 --- a/scripts/analyses/migrate_kernel_timing_source.py +++ b/scripts/analyses/migrate_kernel_timing_source.py @@ -7,38 +7,6 @@ CONFIG_TABLE = "#psychophysical_kernel_fit_config" KERNEL_TABLE = "__psychophysical_kernel" -CONFIG_COLUMNS = ( - ( - "kernel_method", - "enum('legacy_variable','fixed_window') NOT NULL DEFAULT 'legacy_variable'", - "regularization_c", - ), - ("bin_width_s", "float DEFAULT NULL", "kernel_method"), - ( - "evidence_encoding", - "enum('max_rate','trial_rate') DEFAULT NULL", - "bin_width_s", - ), - ("min_trials_per_bin", "int DEFAULT NULL", "evidence_encoding"), - ( - "observation_window", - "enum('center_exit','response') DEFAULT NULL", - "min_trials_per_bin", - ), -) -KERNEL_COLUMNS = ( - ( - "timing_source", - "enum('nidaq','bpod','mixed') NOT NULL DEFAULT 'bpod'", - "fit_message", - ), - ("n_bins_fit", "int DEFAULT NULL", "n_trials_fit"), - ("n_observed_per_bin", "longblob DEFAULT NULL", "n_bins_fit"), - ("bin_centers_s", "longblob DEFAULT NULL", "n_observed_per_bin"), - ("majority_accuracy", "float DEFAULT NULL", "score_mean"), - ("score_above_majority", "float DEFAULT NULL", "majority_accuracy"), - ("interpretation", "varchar(32) DEFAULT NULL", "bias_mean"), -) def parse_args() -> argparse.Namespace: @@ -46,7 +14,7 @@ def parse_args() -> argparse.Namespace: parser.add_argument( "--apply", action="store_true", - help="Alter the live kernel tables and insert fixed-window configs.", + help="Recreate the disposable kernel tables with the canonical schema.", ) return parser.parse_args() @@ -58,26 +26,16 @@ def main() -> None: connection = dj.conn() database = f"{labdata_schema.dbase_name}_user" - table_columns = { - CONFIG_TABLE: CONFIG_COLUMNS, - KERNEL_TABLE: KERNEL_COLUMNS, - } - existing = { - table: _existing_columns(connection, database, table) for table in table_columns - } - statements = [ - statement - for table, columns in table_columns.items() - for statement in _missing_column_statements( - database, table, columns, existing[table] - ) - ] + statements = _reset_statements(database) if not args.apply: + for table in (KERNEL_TABLE, CONFIG_TABLE): + count = connection.query( + f"SELECT COUNT(*) FROM `{database}`.`{table}`" + ).fetchone()[0] + print(f"{table}: {count} rows will be deleted") for statement in statements: print(statement) - if not statements: - print("Kernel timing-source columns are already present.") print("Dry run only. Re-run with --apply after exact live-write approval.") return @@ -90,8 +48,15 @@ def main() -> None: PsychophysicalKernelFitConfig.contents, skip_duplicates=True, ) - _validate_migration(connection, database) - print("Kernel timing-source schema migration complete.") + _validate_schema(connection, database) + print("Kernel schema reset complete.") + + +def _reset_statements(database: str) -> list[str]: + return [ + f"DROP TABLE IF EXISTS `{database}`.`{KERNEL_TABLE}`", + f"DROP TABLE IF EXISTS `{database}`.`{CONFIG_TABLE}`", + ] def _existing_columns(connection, database: str, table: str) -> set[str]: @@ -104,40 +69,55 @@ def _existing_columns(connection, database: str, table: str) -> set[str]: } -def _missing_column_statements( - database: str, - table: str, - columns: tuple[tuple[str, str, str], ...], - existing: set[str], -) -> list[str]: - return [ - f"ALTER TABLE `{database}`.`{table}` ADD COLUMN `{name}` " - f"{definition} AFTER `{after}`" - for name, definition, after in columns - if name not in existing - ] - - -def _validate_migration(connection, database: str) -> None: - expected = { - CONFIG_TABLE: {name for name, _, _ in CONFIG_COLUMNS}, - KERNEL_TABLE: {name for name, _, _ in KERNEL_COLUMNS}, +def _validate_schema(connection, database: str) -> None: + config_columns = _existing_columns(connection, database, CONFIG_TABLE) + required_config_columns = { + "kernel_fit_config_id", + "timebins", + "binning_method", + "bin_width_s", + "observation_window", + "evidence_model", + "min_trials_per_bin", + "cv_splits", + "random_state", + "regularization_c", } - missing = { - table: sorted(columns - _existing_columns(connection, database, table)) - for table, columns in expected.items() + if config_columns != required_config_columns: + raise RuntimeError( + "Kernel config schema mismatch: " + f"expected={sorted(required_config_columns)}, " + f"found={sorted(config_columns)}" + ) + + kernel_columns = _existing_columns(connection, database, KERNEL_TABLE) + required_kernel_columns = { + "timing_source", + "n_bins_fit", + "n_observed_per_bin", + "bin_centers_s", + "majority_accuracy", + "score_above_majority", + "interpretation", } - missing = {table: columns for table, columns in missing.items() if columns} - if missing: - raise RuntimeError(f"Kernel timing-source migration incomplete: {missing}") - - invalid_sources = connection.query( - f"SELECT COUNT(*) FROM `{database}`.`{KERNEL_TABLE}` " - "WHERE timing_source NOT IN ('nidaq','bpod','mixed')" - ).fetchone()[0] - if invalid_sources: + missing_kernel_columns = required_kernel_columns - kernel_columns + if missing_kernel_columns: + raise RuntimeError( + f"Kernel result schema is missing: {sorted(missing_kernel_columns)}" + ) + + rows = connection.query( + "SELECT `kernel_fit_config_id`, `binning_method`, " + "`observation_window`, `evidence_model` " + f"FROM `{database}`.`{CONFIG_TABLE}` ORDER BY `kernel_fit_config_id`" + ).fetchall() + expected = [ + (0, "fixed_width", "center_exit", "trial_rate_residual"), + (1, "fixed_width", "response", "trial_rate_residual"), + ] + if list(rows) != expected: raise RuntimeError( - f"Kernel timing-source migration found {invalid_sources} invalid rows" + f"Kernel config rows mismatch: expected={expected}, found={rows}" ) diff --git a/src/behavior_analyses/kernel_timing.py b/src/behavior_analyses/kernel_timing.py index 7146d31..1bf1087 100644 --- a/src/behavior_analyses/kernel_timing.py +++ b/src/behavior_analyses/kernel_timing.py @@ -13,6 +13,7 @@ "center_port", "right_port", ) +MAX_NIDAQ_ALIGNMENT_ERROR_S = 0.1 def fetch_pooled_kernel_inputs( @@ -87,15 +88,21 @@ def extract_bpod_kernel_inputs( continue response = row["response"] observation_end = row.get(end_field) + trial_sync = row.get("t_sync") stims = np.asarray(row.get("stim_events", []), dtype=float) stims = stims[np.isfinite(stims)] if ( response not in (-1, 1) or observation_end is None or not np.isfinite(observation_end) + or trial_sync is None + or not np.isfinite(trial_sync) or stims.size == 0 ): continue + # Chipmunk stimulus events are relative to the Bpod sync pulse, while + # state-transition timestamps are absolute within the session. + observation_end = float(observation_end) - float(trial_sync) stims = stims[stims < float(observation_end)] if stims.size == 0 or observation_end <= stims[0]: continue @@ -136,7 +143,6 @@ def extract_nidaq_kernel_inputs( nidaq_sync = nidaq_sync[order] stims = np.asarray(aligned_events["visual_stim"], dtype=float) - center_entries = np.asarray(aligned_events["center_port"], dtype=float) center_exits = np.asarray(aligned_events["center_port_exit"], dtype=float) left_entries = np.asarray(aligned_events["left_port"], dtype=float) right_entries = np.asarray(aligned_events["right_port"], dtype=float) @@ -148,54 +154,71 @@ def extract_nidaq_kernel_inputs( or row["response"] not in (-1, 1) or row.get("t_react") is None or not np.isfinite(row["t_react"]) + or row.get("t_sync") is None + or not np.isfinite(row["t_sync"]) ): continue + bpod_stims = np.asarray(row.get("stim_events", []), dtype=float) + bpod_stims = bpod_stims[np.isfinite(bpod_stims)] + if bpod_stims.size == 0: + continue trial_number = int(row["trial_num"]) if trial_number >= trial_starts.size: - raise ValueError(f"NIDAQ trial_start is missing trial_num={trial_number}") + continue trial_start = trial_starts[trial_number] trial_end = ( trial_starts[trial_number + 1] if trial_number + 1 < trial_starts.size else np.inf ) + interpolated_first_stim = float( + np.interp( + float(row["t_sync"]) + float(bpod_stims[0]), + bpod_sync, + nidaq_sync, + ) + ) interpolated_exit = float( np.interp(float(row["t_react"]), bpod_sync, nidaq_sync) ) trial_exits = center_exits[ (center_exits > trial_start) & (center_exits < trial_end) ] - if trial_exits.size: - center_exit = float( - trial_exits[np.argmin(np.abs(trial_exits - interpolated_exit))] - ) - else: - center_exit = interpolated_exit - - center_mask = ( - (center_entries > trial_start) - & (center_entries < center_exit) - & (center_entries < trial_end) - ) - if not center_mask.any(): - raise ValueError(f"No NIDAQ center-port entry for trial_num={trial_number}") - center_entry = float(center_entries[center_mask][-1]) + center_exit = _nearest_aligned_event(trial_exits, interpolated_exit) + if center_exit is None: + continue - response_entries = right_entries if int(row["response"]) == 1 else left_entries - response_mask = (response_entries > center_exit) & ( - response_entries < trial_end - ) - if not response_mask.any(): - raise ValueError( - f"No NIDAQ response-port entry for trial_num={trial_number}" + observation_end = center_exit + if observation_window == "response": + response_time = row.get("t_response") + if response_time is None or not np.isfinite(response_time): + continue + interpolated_response = float( + np.interp(float(response_time), bpod_sync, nidaq_sync) ) - response_entry = float(response_entries[response_mask][0]) - observation_end = ( - center_exit if observation_window == "center_exit" else response_entry - ) - trial_stims = stims[(stims >= center_entry) & (stims < observation_end)] + response_entries = ( + right_entries if int(row["response"]) == 1 else left_entries + ) + trial_responses = response_entries[ + (response_entries > center_exit) & (response_entries < trial_end) + ] + response_entry = _nearest_aligned_event( + trial_responses, interpolated_response + ) + if response_entry is None: + continue + observation_end = response_entry + + trial_stims = stims[ + (stims >= trial_start) & (stims < observation_end) & (stims < trial_end) + ] if trial_stims.size == 0: continue + # Use the Bpod schedule only to identify the first hardware-timed flash. + first_stim = _nearest_aligned_event(trial_stims, interpolated_first_stim) + if first_stim is None: + continue + trial_stims = trial_stims[trial_stims >= first_stim] _append_kernel_trial(result, trial_stims, observation_end, row) return result @@ -375,6 +398,15 @@ def _append_kernel_trial( result["trial_rate_hz"].append(float(rate)) +def _nearest_aligned_event(events: np.ndarray, target: float) -> float | None: + if events.size == 0: + return None + event = float(events[np.argmin(np.abs(events - target))]) + if abs(event - target) > MAX_NIDAQ_ALIGNMENT_ERROR_S: + return None + return event + + def _digital_onsets(timestamps: np.ndarray, values: np.ndarray | None) -> np.ndarray: return timestamps[::2] if values is None else timestamps[values == 1] diff --git a/src/behavior_analyses/kernels.py b/src/behavior_analyses/kernels.py index 979fbd1..a215c03 100644 --- a/src/behavior_analyses/kernels.py +++ b/src/behavior_analyses/kernels.py @@ -9,100 +9,16 @@ def build_residual_rate_matrix( - stim_events, - response_values, - *, - timebins: int = 10, - max_rate_hz: float = 20.0, -) -> tuple[np.ndarray, np.ndarray]: - rows = [] - choices = [] - for events, response in zip(stim_events, response_values): - if response not in (-1, 1): - continue - events = np.asarray(events, dtype=float) - events = events[np.isfinite(events)] - if events.size < 2: - continue - bins = np.linspace(events[0], events[-1], num=timebins + 1) - specific_rate = max_rate_hz / len(bins) - instantaneous_rate, _ = np.histogram(events, bins=bins) - rows.append(instantaneous_rate - specific_rate) - choices.append(response == 1) - if not rows: - return np.empty((0, timebins)), np.empty((0,), dtype=int) - return np.asarray(rows, dtype=float), np.asarray(choices, dtype=int) - - -def fit_psychophysical_kernel( - stim_events, - response_values, - *, - timebins: int = 10, - cv_splits: int = 10, - random_state: int = 0, - max_rate_hz: float = 20.0, - regularization_c: float = 1.0, -) -> dict: - x, y = build_residual_rate_matrix( - stim_events, response_values, timebins=timebins, max_rate_hz=max_rate_hz - ) - if x.shape[0] < cv_splits or np.unique(y).size < 2: - return { - "design_matrix": x, - "choice_right": y, - "weights": np.empty((0, timebins)), - "scores": np.empty((0,)), - "bias": np.empty((0,)), - "error": np.empty((0, timebins)), - } - - splitter = StratifiedKFold( - n_splits=cv_splits, shuffle=True, random_state=random_state - ) - weights = [] - scores = [] - biases = [] - errors = [] - for train_index, test_index in splitter.split(x, y): - x_train, x_test = x[train_index], x[test_index] - y_train, y_test = y[train_index], y[test_index] - model = LogisticRegression( - penalty="l2", - solver="liblinear", - C=regularization_c, - fit_intercept=True, - ).fit(x_train, y_train) - predict_prob = model.predict_proba(x_train) - variance = np.prod(predict_prob, axis=1) - covariance = np.linalg.pinv(np.dot(x_train.T * variance, x_train)) - errors.append(np.sqrt(np.diag(covariance))) - weights.append(model.coef_[0]) - scores.append(model.score(x_test, y_test)) - biases.append(model.intercept_[0]) - - return { - "design_matrix": x, - "choice_right": y, - "weights": np.asarray(weights, dtype=float), - "scores": np.asarray(scores, dtype=float), - "bias": np.asarray(biases, dtype=float), - "error": np.asarray(errors, dtype=float), - } - - -def build_fixed_residual_rate_matrix( stim_times_per_trial: Sequence[np.ndarray], first_stim_times: Sequence[float], observation_end_times: Sequence[float], response_values: Sequence[Any], + trial_rate_hz: Sequence[float], *, timebins: int = 10, bin_width_s: float = 0.1, - max_rate_hz: float = 20.0, - trial_rate_hz: Sequence[float] | None = None, ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: - """Build fixed-width residual flash counts for a timing-aware kernel.""" + """Build Odoemene Eq. 5 residual counts in fixed-width time bins.""" if timebins < 1: raise ValueError("timebins must be >= 1") if bin_width_s <= 0: @@ -119,14 +35,11 @@ def build_fixed_residual_rate_matrix( first_stim = np.asarray(first_stim_times, dtype=float)[choice_mask] observation_end = np.asarray(observation_end_times, dtype=float)[choice_mask] - if trial_rate_hz is None: - expected_rates = np.full(choice_mask.shape, max_rate_hz, dtype=float) - else: - expected_rates = np.asarray(trial_rate_hz, dtype=float) - if expected_rates.shape != choice_mask.shape: - raise ValueError("trial_rate_hz must match response_values") - if np.any(~np.isfinite(expected_rates)) or np.any(expected_rates < 0): - raise ValueError("trial_rate_hz must contain finite nonnegative rates") + expected_rates = np.asarray(trial_rate_hz, dtype=float) + if expected_rates.shape != choice_mask.shape: + raise ValueError("trial_rate_hz must match response_values") + if np.any(~np.isfinite(expected_rates)) or np.any(expected_rates < 0): + raise ValueError("trial_rate_hz must contain finite nonnegative rates") expected_rates = expected_rates[choice_mask] bin_edges = np.arange(timebins + 1, dtype=float) * bin_width_s @@ -166,18 +79,18 @@ def build_fixed_residual_rate_matrix( return residual, choices, n_observed_per_bin, bin_centers_s, expected_counts -def fit_fixed_psychophysical_kernel( +def fit_psychophysical_kernel( residual: np.ndarray, choice_right: np.ndarray, *, + expected_counts: np.ndarray, n_observed_per_bin: np.ndarray | None = None, cv_splits: int = 10, random_state: int = 0, min_trials_per_bin: int = 50, regularization_c: float = 1.0, - expected_counts: np.ndarray | None = None, ) -> dict[str, Any]: - """Fit an L2 logistic kernel using the longest complete-case bin prefix.""" + """Fit the Odoemene Eq. 5 kernel over the longest complete-case prefix.""" residual = np.asarray(residual, dtype=float) choice_right = np.asarray(choice_right, dtype=int) n_bins = residual.shape[1] if residual.ndim == 2 else 0 @@ -222,7 +135,7 @@ def fit_fixed_psychophysical_kernel( if n_splits < 2: return empty - design, coefficient_to_weights = _fixed_kernel_design( + design, coefficient_to_weights = _kernel_design( residual, x, complete, @@ -335,16 +248,13 @@ def _complete_case_prefix( return prefix -def _fixed_kernel_design( +def _kernel_design( residual: np.ndarray, complete_residual: np.ndarray, complete: np.ndarray, n_bins_fit: int, - expected_counts: np.ndarray | None, + expected_counts: np.ndarray, ) -> tuple[np.ndarray, np.ndarray]: - if expected_counts is None: - return complete_residual, np.eye(n_bins_fit) - expected_counts = np.asarray(expected_counts, dtype=float) if expected_counts.shape != residual.shape: raise ValueError("expected_counts must match residual") diff --git a/tests/test_analysis_functions.py b/tests/test_analysis_functions.py index c1bd951..47f54a1 100644 --- a/tests/test_analysis_functions.py +++ b/tests/test_analysis_functions.py @@ -86,35 +86,17 @@ def test_summarize_trialset_counts_choice_and_means(self): class KernelTests(unittest.TestCase): - def test_kernel_design_matrix_skips_no_choice_and_short_stim_trials(self): + def test_kernel_keeps_unobserved_late_bins_as_nan(self): from behavior_analyses.kernels import build_residual_rate_matrix - stim_events = [ - np.array([0.0, 0.1, 0.2]), - np.array([0.0]), - np.array([0.0, 0.2, 0.4]), - ] - responses = np.array([1, -1, 0]) - - x, y = build_residual_rate_matrix(stim_events, responses, timebins=2) - - self.assertEqual(x.shape, (1, 2)) - np.testing.assert_allclose(x[0], [1 - 20 / 3, 2 - 20 / 3]) - np.testing.assert_array_equal(y, [1]) - - def test_fixed_kernel_keeps_unobserved_late_bins_as_nan(self): - from behavior_analyses.kernels import build_fixed_residual_rate_matrix - - residual, choices, n_observed, centers, expected = ( - build_fixed_residual_rate_matrix( - [np.array([1.0, 1.05, 1.12])], - [1.0], - [1.15], - [1], - timebins=3, - bin_width_s=0.1, - trial_rate_hz=[20.0], - ) + residual, choices, n_observed, centers, expected = build_residual_rate_matrix( + [np.array([1.0, 1.05, 1.12])], + [1.0], + [1.15], + [1], + [20.0], + timebins=3, + bin_width_s=0.1, ) np.testing.assert_array_equal(choices, [1]) @@ -123,26 +105,29 @@ def test_fixed_kernel_keeps_unobserved_late_bins_as_nan(self): self.assertTrue(np.isnan(residual[0, 2])) self.assertTrue(np.isnan(expected[0, 2])) - def test_fixed_kernel_fit_is_deterministic(self): - from behavior_analyses.kernels import fit_fixed_psychophysical_kernel + def test_kernel_fit_is_deterministic(self): + from behavior_analyses.kernels import fit_psychophysical_kernel rng = np.random.default_rng(4) residual = rng.normal(size=(120, 4)) choices = (residual[:, 0] > 0).astype(int) residual[60:, 3] = np.nan n_observed = np.sum(np.isfinite(residual), axis=0) + expected = np.where(np.isfinite(residual), 2.0, np.nan) - first = fit_fixed_psychophysical_kernel( + first = fit_psychophysical_kernel( residual, choices, + expected_counts=expected, n_observed_per_bin=n_observed, cv_splits=5, random_state=7, min_trials_per_bin=50, ) - second = fit_fixed_psychophysical_kernel( + second = fit_psychophysical_kernel( residual, choices, + expected_counts=expected, n_observed_per_bin=n_observed, cv_splits=5, random_state=7, diff --git a/tests/test_cli_and_migration_contracts.py b/tests/test_cli_and_migration_contracts.py index c34178c..52f2a4d 100644 --- a/tests/test_cli_and_migration_contracts.py +++ b/tests/test_cli_and_migration_contracts.py @@ -51,19 +51,17 @@ def populate(*_args, **_kwargs): class CliContractTests(unittest.TestCase): - def test_kernel_timing_migration_only_adds_missing_columns(self): + def test_kernel_timing_migration_drops_result_before_config(self): module = runpy.run_path(str(SCRIPTS / "migrate_kernel_timing_source.py")) - statements = module["_missing_column_statements"]( - "labdata_user", - module["KERNEL_TABLE"], - module["KERNEL_COLUMNS"], - {"timing_source", "n_trials_fit", "fit_message"}, - ) + statements = module["_reset_statements"]("labdata_user") - self.assertEqual(len(statements), len(module["KERNEL_COLUMNS"]) - 1) - self.assertTrue(all("timing_source" not in sql for sql in statements)) - self.assertTrue(all(sql.startswith("ALTER TABLE") for sql in statements)) + self.assertEqual(len(statements), 2) + self.assertIn(module["KERNEL_TABLE"], statements[0]) + self.assertIn(module["CONFIG_TABLE"], statements[1]) + self.assertTrue( + all(sql.startswith("DROP TABLE IF EXISTS") for sql in statements) + ) def test_schema_migration_archive_names_leave_room_for_foreign_keys(self): module = runpy.run_path(str(SCRIPTS / "migrate_behavior_analysis_schema.py")) @@ -86,14 +84,6 @@ def test_schema_migration_rejects_occupied_target_tables(self): with self.assertRaisesRegex(RuntimeError, "occupied_targets"): module["_validate_table_state"](connection, "labdata_user") - expected = MagicMock() - expected.fetchall.return_value = [(10, 10, 0)] - incompatible = MagicMock() - incompatible.fetchall.return_value = [(8, 10, 0)] - connection.query.side_effect = [expected, incompatible] - with self.assertRaisesRegex(RuntimeError, "Incompatible legacy kernel"): - module["_validate_kernel_configs"](connection, "labdata_user") - def test_schema_migration_accepts_expected_resume_state(self): module = runpy.run_path(str(SCRIPTS / "migrate_behavior_analysis_schema.py")) connection = MagicMock() diff --git a/tests/test_kernel_timing.py b/tests/test_kernel_timing.py index bf32979..095b6f1 100644 --- a/tests/test_kernel_timing.py +++ b/tests/test_kernel_timing.py @@ -20,17 +20,17 @@ def setUp(self): "stim_events": np.array([0.2, 0.35, 0.7]), "stim_rate_vision": 12.0, "response": 1, - "t_sync": 0.0, - "t_react": 0.5, - "t_response": 0.9, + "t_sync": 10.0, + "t_react": 10.5, + "t_response": 10.9, }, { "trial_num": 1, "rewarded_modality": "visual", - "stim_events": np.array([2.2]), + "stim_events": np.array([0.2]), "stim_rate_vision": 8.0, "response": 0, - "t_sync": 2.0, + "t_sync": 12.0, "t_react": None, "t_response": None, }, @@ -101,6 +101,37 @@ def test_nidaq_window_uses_mapped_flash_and_port_times(self): np.testing.assert_allclose(center["observation_end_times"], [0.5]) np.testing.assert_allclose(response["observation_end_times"], [0.9]) + def test_bpod_and_nidaq_match_for_aligned_trial(self): + from behavior_analyses.kernel_timing import ( + extract_bpod_kernel_inputs, + extract_nidaq_kernel_inputs, + resolve_nidaq_event_arrays, + ) + + mapping_rows, event_rows = self._mapped_event_fixture() + aligned = resolve_nidaq_event_arrays( + event_rows, mapping_rows, "GRB006", "session" + ) + + for window in ("center_exit", "response"): + bpod = extract_bpod_kernel_inputs( + self.trial_rows, "visual", observation_window=window + ) + nidaq = extract_nidaq_kernel_inputs( + aligned, + self.trial_rows, + "visual", + observation_window=window, + ) + + self.assertEqual(bpod["response_values"], nidaq["response_values"]) + np.testing.assert_allclose( + bpod["observation_end_times"], nidaq["observation_end_times"] + ) + np.testing.assert_allclose( + bpod["stim_times_per_trial"][0], nidaq["stim_times_per_trial"][0] + ) + def test_combined_provenance_is_mixed(self): from behavior_analyses.kernel_timing import ( combine_kernel_inputs, @@ -124,15 +155,15 @@ def test_combined_provenance_is_mixed(self): @staticmethod def _mapped_event_fixture(): sources = { - "visual_stim": ("ai0", [0.2, 0.35, 0.7, 2.2], None), + "visual_stim": ("ai0", [0.15, 0.2, 0.35, 0.7, 2.2], None), "trial_start": ("line0", [0.0, 2.0], [1, 1]), "left_port": ("line1", [2.3, 2.4], [1, 0]), "center_port": ( "line2", - [0.1, 0.5, 2.1, 2.5], - [1, 0, 1, 0], + [0.1, 0.45, 0.499, 0.5, 2.1, 2.5], + [1, 0, 1, 0, 1, 0], ), - "right_port": ("line3", [0.9, 1.0], [1, 0]), + "right_port": ("line3", [0.55, 0.56, 0.9, 1.0], [1, 0, 1, 0]), } mapping_rows = [] event_rows = [] diff --git a/tests/test_schema_imports.py b/tests/test_schema_imports.py index aed5355..339591e 100644 --- a/tests/test_schema_imports.py +++ b/tests/test_schema_imports.py @@ -71,6 +71,26 @@ def test_analysis_schema_imports_with_fake_labdata(self): self.assertTrue(hasattr(module, "PsychometricFitConfig")) self.assertTrue(hasattr(module, "PsychophysicalKernelFitConfig")) self.assertTrue(hasattr(module, "PsychophysicalKernel")) + self.assertIn( + "kernel_fit_config_id : int", + module.PsychophysicalKernelFitConfig.definition, + ) + self.assertNotIn( + "analysis_version", + module.PsychophysicalKernelFitConfig.definition, + ) + self.assertNotIn( + "kernel_method", + module.PsychophysicalKernelFitConfig.definition, + ) + self.assertNotIn( + "evidence_encoding", + module.PsychophysicalKernelFitConfig.definition, + ) + self.assertEqual( + [row[0] for row in module.PsychophysicalKernelFitConfig.contents], + [0, 1], + ) self.assertFalse(hasattr(module, "LearningSessionMetrics"))