plot_comparison.py

import numpy as np
import matplotlib.pyplot as plt
import math
import os
import matplotlib.patches as mpatches
import pandas as pd
import utilities

# Plot: Comparison between aCCE values, simulated from the model, and actual experimental data
# (Exp. 2 of Lohmann et al. 2019)

# 'consistent mapping' here taken as 'no drift' condition, 'variable mapping' as 'constant drift' condition

# Timesteps - rough equivalence
soa_exp_onset = "MovementOnset" # in the data
soa_exp_partway = "OneHalf"
t_model_onset = 0
t_model_partway = 50 #

exp_data_rts = "Data/experimental_Lohmannetal2019Exp2/RawVerbalRT.dat"
#exp_data_acce = "Data/experimental_Lohmannetal2019Exp2/CongruencyData.dat"
data_folder = "Exp-Data/"


def _safe_load_model_rts(testfolder, folder="Exp-Data/", alpha_lr=None, threshold=None, t_reach=None):
    """Load or compute RTs_seq for a given model testfolder using utilities.run_rts_over_time.

    Returns the RTs_seq numpy array (timesteps x 17) or None on error.
    """
    try:
        kwargs = {}
        if alpha_lr is not None:
            kwargs["alpha_lr"] = alpha_lr
        if threshold is not None:
            kwargs["threshold"] = threshold
        # Run it as long as reach for now?
        if t_reach is not None:
            kwargs["t_to_run"] = t_reach
        RTs_seq = utilities.run_rts_over_time(testfolder, folder=folder, **kwargs)
        return RTs_seq
    except Exception as e:
        print(f"Error loading RTs_seq for {testfolder}: {e}")
        return None


def model_acce_from_RTs_seq(RTs_seq, t_idx):
    """Compute aCCE values (upright and upside-down) and approximate std errors from RTs_seq at timestep t_idx.

    Returns dict: {'overall': val, 'upright': val, 'upsidedown': val, 'se_upright': se, 'se_upsidedown': se}
    """
    if RTs_seq is None:
        return None
    t_max = RTs_seq.shape[0]
    if t_idx < 0:
        t_idx = 0
    if t_idx >= t_max:
        t_idx = t_max - 1

    overall = float(RTs_seq[t_idx, 0])
    # columns: 5=mean_congruent_upright,6=mean_congruent_upsidedown,7=incongr_upright,8=incongr_upsidedown
    cong_u = float(RTs_seq[t_idx, 5])
    cong_ud = float(RTs_seq[t_idx, 6])
    incong_u = float(RTs_seq[t_idx, 7])
    incong_ud = float(RTs_seq[t_idx, 8])

    val_u = incong_u - cong_u
    val_ud = incong_ud - cong_ud

    # stds: 13=std_congruent_upright,14=std_congruent_upsidedown,15=std_incongruent_upright,16=std_incongruent_upsidedown
    try:
        s_cong_u = float(RTs_seq[t_idx, 13])
        s_cong_ud = float(RTs_seq[t_idx, 14])
        s_incong_u = float(RTs_seq[t_idx, 15])
        s_incong_ud = float(RTs_seq[t_idx, 16])
        # approximate se for difference as sqrt(s1^2 + s2^2) / sqrt(N). We don't have N here; use pooled sqrt(s1^2+s2^2)
        se_u = math.sqrt(s_cong_u ** 2 + s_incong_u ** 2)
        se_ud = math.sqrt(s_cong_ud ** 2 + s_incong_ud ** 2)
    except Exception:
        se_u = se_ud = 0.0

    return {
        'overall': overall,
        'upright': val_u,
        'upsidedown': val_ud,
        'se_upright': se_u,
        'se_upsidedown': se_ud
    }


def is_congruent_exp_row(orient, light, tactile):
    """Return True if the trial (experimental row) is congruent according to the mapping described.

    Accepts orient: can be 'upright','upside-down', 1, -1, etc.; light: 'left'/'right' or numeric; tactile: 'thumb'/'index' or numeric.
    """
    # normalize inputs to simple strings when possible
    def norm(x):
        if x is None:
            return None
        if isinstance(x, str):
            return x.strip().lower()
        try:
            return int(x)
        except Exception:
            return x

    o = norm(orient)
    l = norm(light)
    t = norm(tactile)

    # interpret orientation
    upr = None
    if isinstance(o, int):
        upr = (o == 1)
    elif isinstance(o, str):
        upr = ('upr' in o) or (o == 'up')

    # interpret light and tactile
    left = None
    if isinstance(l, int):
        left = (l == 1)
    elif isinstance(l, str):
        left = ('left' in l)

    is_thumb = None
    if isinstance(t, int):
        is_thumb = (t == 1)
    elif isinstance(t, str):
        is_thumb = ('thumb' in t) or ('thal' in t)

    # apply congruency rules
    if upr is None or left is None or is_thumb is None:
        return False
    if upr:
        # congruent: left light + thumb OR right light + index
        if left and is_thumb:
            return True
        if (not left) and (not is_thumb):
            return True
        return False
    else:
        # upside-down: congruent is reversed
        if left and (not is_thumb):
            return True
        if (not left) and is_thumb:
            return True
        return False


def load_experiment_acce(path, soa_values=(soa_exp_onset, soa_exp_partway)):
    """Load experimental congruency table and compute aCCE (incongr - congr) for each SOA and bottle orientation.

    Returns dict keyed by soa -> {'upright': (mean, sem, n), 'upsidedown': (mean, sem, n)}
    """
    try:
        df = pd.read_csv(path, sep=None, engine='python')
    except Exception as e:
        print(f"Could not read experimental data from {path}: {e}")
        return None

    # ensure expected columns
    cols = [c.lower() for c in df.columns]
    # try to find column names in a case-insensitive way
    def find_col(prefixes):
        for p in prefixes:
            for c in df.columns:
                if p.lower() == c.lower():
                    return c
        return None

    col_rt = find_col(['RT', 'Rt', 'rt', 'ReactionTime'])
    col_orient = find_col(['BottleOrientation', 'Bottle', 'Orientation', 'bottleorientation'])
    col_light = find_col(['LightStimulation', 'Light', 'lightstimulation'])
    col_tac = find_col(['TactileStimulation', 'Tactile', 'tactilestimulation'])
    col_soa = find_col(['SOACondition', 'SOA', 'SoaCondition'])
    col_pid = find_col(['ParticipantID', 'Participant', 'PID'])
    col_map = find_col(['Mapping', 'mapping', 'Map', 'map', 'MappingCondition'])

    if col_rt is None or col_orient is None or col_light is None or col_tac is None or col_soa is None:
        print("Experimental file missing one of required columns: RT, BottleOrientation, LightStimulation, TactileStimulation, SOACondition")
        return None

    results = {}
    for soa in soa_values:
        df_soa = df[df[col_soa].astype(str) == str(soa)]
        if df_soa.shape[0] == 0:
            # try matching by substring
            df_soa = df[df[col_soa].astype(str).str.contains(str(soa), case=False, na=False)]

        # compute congruency per row
        df_soa = df_soa.copy()
        df_soa['is_congr'] = df_soa.apply(lambda r: is_congruent_exp_row(r[col_orient], r[col_light], r[col_tac]), axis=1)

        # determine mapping groups
        if col_map is not None:
            df_soa['_mapping'] = df_soa[col_map].astype(str).str.strip().str.lower()
            mapping_vals = sorted(df_soa['_mapping'].dropna().unique())
        else:
            df_soa['_mapping'] = 'all'
            mapping_vals = ['all']

        out_map = {}
        for map_val in mapping_vals:
            d_map = df_soa[df_soa['_mapping'] == map_val]
            out = {}
            for orient_val in d_map[col_orient].unique():
                d_or = d_map[d_map[col_orient] == orient_val]
                cong = d_or[d_or['is_congr'] == True]
                incong = d_or[d_or['is_congr'] == False]
                n_cong = 0
                n_incong = 0
                try:
                    mean_cong = cong[col_rt].astype(float).mean()
                    mean_incong = incong[col_rt].astype(float).mean()
                    n_cong = cong.shape[0]
                    n_incong = incong.shape[0]
                    sem_cong = cong[col_rt].astype(float).std(ddof=1) / math.sqrt(n_cong) if n_cong > 1 else 0.0
                    sem_incong = incong[col_rt].astype(float).std(ddof=1) / math.sqrt(n_incong) if n_incong > 1 else 0.0
                    acce_val = float(mean_incong - mean_cong) if (not math.isnan(mean_incong)) and (not math.isnan(mean_cong)) else float('nan')
                    sem_diff = math.sqrt((sem_cong ** 2) + (sem_incong ** 2))
                except Exception:
                    acce_val = float('nan')
                    sem_diff = 0.0
                key = 'upright' if ('upr' in str(orient_val).lower()) or (str(orient_val) in ['1', '1.0']) else 'upsidedown'
                out[key] = (acce_val, sem_diff, int(n_cong + n_incong))
            out_map[map_val] = out

        results[soa] = out_map

    return results


def plot_compare(model_no_drift_folder, model_drift_folder, exp_table_path, alpha_lr=None, threshold=None, t_reach=None,
                 timesteps=(t_model_onset, t_model_partway), title_part=None, use_two_axes=True,
                 model_time_scale=None, use_e_bars=True):
    """Compute aCCE for models and experiment and produce comparison plots for the two timesteps/SOAs.
    """
    # load experimental
    exp_acce = load_experiment_acce(exp_table_path, soa_values=(soa_exp_onset, soa_exp_partway))
    if exp_acce is None:
        print("Failed to load experimental data - aborting plot")
        return

    # load models
    RTs_no = _safe_load_model_rts(model_no_drift_folder, folder=data_folder, alpha_lr=alpha_lr, threshold=threshold, t_reach=t_reach)
    RTs_c = _safe_load_model_rts(model_drift_folder, folder=data_folder, alpha_lr=alpha_lr, threshold=threshold, t_reach=t_reach)

    # compute global y-limits for experiment and model so both panels use same scales
    # collect both values and their errors (sem) so axis limits can include error bars when requested
    exp_vals = []
    exp_errs = []
    model_vals = []
    model_errs = []
    # collect experiment values across both SOAs
    for soa_k in (soa_exp_onset, soa_exp_partway):
        exp_for_soa = exp_acce.get(soa_k, {})
        if isinstance(exp_for_soa, dict) and exp_for_soa:
            # mapping case
            first_val = next(iter(exp_for_soa.values()))
            if isinstance(first_val, dict):
                for mkey, mp in exp_for_soa.items():
                    for orient in ('upright', 'upsidedown'):
                        t = mp.get(orient, (np.nan, 0.0, 0))
                        try:
                            v = float(t[0])
                            se = float(t[1]) if len(t) > 1 else 0.0
                            if not (isinstance(v, float) and math.isnan(v)):
                                exp_vals.append(v)
                                exp_errs.append(se)
                        except Exception:
                            pass
            else:
                for orient in ('upright', 'upsidedown'):
                    t = exp_for_soa.get(orient, (np.nan, 0.0, 0))
                    try:
                        v = float(t[0])
                        se = float(t[1]) if len(t) > 1 else 0.0
                        if not (isinstance(v, float) and math.isnan(v)):
                            exp_vals.append(v)
                            exp_errs.append(se)
                    except Exception:
                        pass
    # collect model values across both timesteps
    for t_idx in timesteps:
        m_no_tmp = model_acce_from_RTs_seq(RTs_no, int(t_idx))
        m_c_tmp = model_acce_from_RTs_seq(RTs_c, int(t_idx))
        for mtmp in (m_no_tmp, m_c_tmp):
            if mtmp:
                try:
                    v_u = float(mtmp.get('upright', np.nan))
                    se_u = float(mtmp.get('se_upright', 0.0))
                    if model_time_scale is not None:
                        try:
                            v_u = v_u * float(model_time_scale)
                            se_u = se_u * float(model_time_scale)
                        except Exception:
                            pass
                    if not (isinstance(v_u, float) and math.isnan(v_u)):
                        model_vals.append(v_u)
                        model_errs.append(se_u)
                except Exception:
                    pass
                try:
                    v_ud = float(mtmp.get('upsidedown', np.nan))
                    se_ud = float(mtmp.get('se_upsidedown', 0.0))
                    if model_time_scale is not None:
                        try:
                            v_ud = v_ud * float(model_time_scale)
                            se_ud = se_ud * float(model_time_scale)
                        except Exception:
                            pass
                    if not (isinstance(v_ud, float) and math.isnan(v_ud)):
                        model_vals.append(v_ud)
                        model_errs.append(se_ud)
                except Exception:
                    pass

    def compute_limits(vals_list, errs_list=None, pad_fraction=0.1):
        # vals_list: list of floats; errs_list: optional list of sems of same length
        clean_vals = []
        clean_errs = []
        for i, v in enumerate(vals_list):
            if v is None:
                continue
            try:
                if isinstance(v, float) and math.isnan(v):
                    continue
            except Exception:
                pass
            clean_vals.append(float(v))
            if errs_list is not None and i < len(errs_list):
                try:
                    clean_errs.append(float(errs_list[i]))
                except Exception:
                    clean_errs.append(0.0)
            else:
                clean_errs.append(0.0)

        if len(clean_vals) == 0:
            return (-1.0, 1.0)

        # if errors provided and positive, expand limits to include error bars
        if errs_list is not None and any([e for e in clean_errs if (not (isinstance(e, float) and math.isnan(e))) and e != 0]):
            lowers = [v - e for v, e in zip(clean_vals, clean_errs)]
            uppers = [v + e for v, e in zip(clean_vals, clean_errs)]
            vmin = float(np.nanmin(lowers))
            vmax = float(np.nanmax(uppers))
        else:
            vmin = float(np.nanmin(clean_vals))
            vmax = float(np.nanmax(clean_vals))

        if vmin == vmax:
            v = abs(vmin) if vmin != 0 else 1.0
            return (-v, v)
        rng = vmax - vmin
        return (vmin - pad_fraction * rng, vmax + pad_fraction * rng)

    exp_ylim = compute_limits(exp_vals, exp_errs if use_e_bars else None)
    model_ylim = compute_limits(model_vals, model_errs if use_e_bars else None)

    # Ensure zero is included in both ranges (so we can align the zero point)
    def include_zero(lim):
        lo, hi = float(lim[0]), float(lim[1])
        if lo <= 0.0 <= hi:
            return (lo, hi)
        return (min(lo, 0.0), max(hi, 0.0))

    exp_ylim = include_zero(exp_ylim)
    model_ylim = include_zero(model_ylim)

    # Align zero point: compute fraction f of zero in exp axis and expand model axis so zero has same fraction
    exp_lo, exp_hi = float(exp_ylim[0]), float(exp_ylim[1])
    model_lo, model_hi = float(model_ylim[0]), float(model_ylim[1])
    exp_span = exp_hi - exp_lo if (exp_hi - exp_lo) != 0 else 1.0
    model_span = model_hi - model_lo if (model_hi - model_lo) != 0 else 1.0

    # fraction of height where zero lies on experiment axis
    f = (0.0 - exp_lo) / exp_span
    # clamp f to avoid exact 0 or 1 which would cause division issues
    eps = 1e-6
    if f < eps:
        f = eps
    if f > 1.0 - eps:
        f = 1.0 - eps

    # Determine minimal model span needed to keep existing model_lo/model_hi while making zero at fraction f
    # Constraints: new_model_lo = -f * new_span <= model_lo  => new_span >= (-model_lo) / f
    #              new_model_hi = new_model_lo + new_span >= model_hi => new_span >= model_hi / (1-f)
    needed_span1 = (-model_lo) / f if f > 0 else float('inf')
    needed_span2 = model_hi / (1.0 - f) if (1.0 - f) > 0 else float('inf')
    needed_span = max(model_span, needed_span1, needed_span2)

    new_model_lo = -f * needed_span
    new_model_hi = new_model_lo + needed_span
    model_ylim = (new_model_lo, new_model_hi)

    # If model_time_scale is provided we likely want both axes to show the same numeric scale
    # so that bars look identical; in that case, override model_ylim to equal exp_ylim when using two axes
    if (model_time_scale is not None) and use_two_axes:
        model_ylim = exp_ylim

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    for ax_idx, (timestep, soa) in enumerate(zip(timesteps, (soa_exp_onset, soa_exp_partway))):
        ax = axes[ax_idx]
        # secondary axis for model values (or shared axis if requested)
        if use_two_axes:
            ax2 = ax.twinx()
        else:
            ax2 = ax
        # model values
        m_no = model_acce_from_RTs_seq(RTs_no, int(timestep))
        m_c = model_acce_from_RTs_seq(RTs_c, int(timestep))

        # Clusters: 4 clusters (Exp consistent, Exp variable, Model (no drift), Model (drift))
        clusters = ['Exp consistent', 'Exp variable', 'Model (no drift)', 'Model (drift)']
        n_clusters = len(clusters)
        cluster_centers = np.arange(n_clusters)

        # orientations bars per cluster: upright, upsidedown, and mean (mean will be computed later)
        orient_base = ['upright', 'upsidedown']
        orientations = orient_base + ['mean']
        n_orient = len(orientations)

        exp_for_soa = exp_acce.get(soa, {})
        has_mapping = False
        mapping_dict = {}
        if isinstance(exp_for_soa, dict) and exp_for_soa:
            first_val = next(iter(exp_for_soa.values()))
            if isinstance(first_val, dict):
                has_mapping = True
                mapping_dict = exp_for_soa

        def pick_mapping_key(mapping_keys, want='consistent'):
            want = want.lower()
            for k in mapping_keys:
                if str(k).strip().lower() == want:
                    return k
            for k in mapping_keys:
                kl = str(k).lower()
                if want == 'consistent' and ('cons' in kl or 'fix' in kl or 'stable' in kl):
                    return k
                if want == 'variable' and ('var' in kl or 'rand' in kl or 'instab' in kl or 'vari' in kl):
                    return k
            if want == 'consistent':
                return mapping_keys[0]
            return mapping_keys[1] if len(mapping_keys) > 1 else mapping_keys[0]

        # build values arrays per cluster and orientation
        vals = np.full((n_clusters, n_orient), np.nan, dtype=float)
        sems = np.zeros((n_clusters, n_orient), dtype=float)

        # fill experimental clusters (only upright and upsidedown; mean computed later)
        if has_mapping:
            keys = sorted(mapping_dict.keys())
            k_cons = pick_mapping_key(keys, 'consistent') if keys else None
            k_var = pick_mapping_key(keys, 'variable') if keys else None
            for oi, orient in enumerate(orient_base):
                if k_cons in mapping_dict:
                    t = mapping_dict[k_cons].get(orient, (np.nan, 0.0, 0))
                    vals[0, oi] = t[0]; sems[0, oi] = t[1]
                if k_var in mapping_dict:
                    t = mapping_dict[k_var].get(orient, (np.nan, 0.0, 0))
                    vals[1, oi] = t[0]; sems[1, oi] = t[1]
        else:
            # no mapping: put available experiment values into first cluster (consistent)
            for oi, orient in enumerate(orient_base):
                t = exp_for_soa.get(orient, (np.nan, 0.0, 0))
                vals[0, oi] = t[0]; sems[0, oi] = t[1]
                vals[1, oi] = np.nan; sems[1, oi] = 0.0

        # fill model clusters: cluster 2 -> model no drift (m_no), cluster 3 -> model drift (m_c)
        for oi, orient in enumerate(orient_base):
            # model no drift -> cluster index 2
            if m_no:
                tval = m_no['upright'] if orient == 'upright' else m_no['upsidedown']
                tse = m_no['se_upright'] if orient == 'upright' else m_no['se_upsidedown']
                try:
                    if model_time_scale is not None:
                        tval = float(tval) * float(model_time_scale)
                        tse = float(tse) * float(model_time_scale)
                except Exception:
                    pass
                vals[2, oi] = tval; sems[2, oi] = tse
            else:
                vals[2, oi] = np.nan; sems[2, oi] = 0.0
            # model drift -> cluster index 3
            if m_c:
                tval = m_c['upright'] if orient == 'upright' else m_c['upsidedown']
                tse = m_c['se_upright'] if orient == 'upright' else m_c['se_upsidedown']
                try:
                    if model_time_scale is not None:
                        tval = float(tval) * float(model_time_scale)
                        tse = float(tse) * float(model_time_scale)
                except Exception:
                    pass
                vals[3, oi] = tval; sems[3, oi] = tse
            else:
                vals[3, oi] = np.nan; sems[3, oi] = 0.0

        # plotting: for each cluster, draw bars for upright, upside-down and mean
        total_width = 0.8
        bar_width = total_width / n_orient / 1.1  # spacing
        # centered offsets for n_orient bars
        orientation_offsets = np.linspace(-total_width/2 + bar_width/2, total_width/2 - bar_width/2, n_orient)
        cluster_spacing = 1.0

        # order of colors matches clusters above: last two are model no drift, model drift
        colors = ['#777777', '#bbbbbb', 'gold', 'orange']
        labels_cluster = ['Exp. (consistent)', 'Exp. (variable)', 'Model (no drift)', 'Model (drift)']

        # compute mean values and sem for each cluster (mean of upright and upside-down)
        for ci in range(n_clusters):
            v_u = vals[ci, 0]
            v_ud = vals[ci, 1]
            se_u = sems[ci, 0]
            se_ud = sems[ci, 1]
            if not (math.isnan(v_u) and math.isnan(v_ud)):
                # if one is nan, take the other as mean
                if math.isnan(v_u):
                    vals[ci, 2] = v_ud
                    sems[ci, 2] = se_ud
                elif math.isnan(v_ud):
                    vals[ci, 2] = v_u
                    sems[ci, 2] = se_u
                else:
                    vals[ci, 2] = (v_u + v_ud) / 2.0
                    # approximate sem of mean: sqrt(se_u^2 + se_ud^2)/2
                    sems[ci, 2] = math.sqrt(se_u ** 2 + se_ud ** 2) / 2.0

        for ci, center in enumerate(cluster_centers):
            for oi, orient in enumerate(orientations):
                x_pos = center + orientation_offsets[oi]
                v = vals[ci, oi]
                se = sems[ci, oi]
                # label the cluster only once (first orientation)
                label = labels_cluster[ci] if oi == 0 else None
                # slightly different alpha for orientation types: upright (0.95), upside-down (0.6), mean (0.85)
                alpha = 0.95 if orient == 'upright' else (0.6 if orient == 'upsidedown' else 0.85)
                # use edge for mean to make it visually distinct
                edge = 'k' if orient == 'mean' else None
                # choose axis: experiment clusters (0,1) -> left axis; model clusters (2,3) -> right axis
                target_ax = ax if ci < 2 else ax2
                # decide whether to show error bars
                yerr = [se] if (use_e_bars and (se is not None) and (not (isinstance(se, float) and math.isnan(se)) ) and se != 0) else None
                target_ax.bar(x_pos, [v], width=bar_width, color=colors[ci], alpha=alpha,
                              yerr=yerr, capsize=4, label=label, edgecolor=edge)

        # set consistent y-limits across panels
        if use_two_axes:
            # experiment on left axis, model on right axis
            ax.set_ylim(exp_ylim)
            ax2.set_ylim(model_ylim)
            ax.set_ylabel('aCCE (experiment RTs)')
            if model_time_scale is not None:
                ax2.set_ylabel(f'aCCE (model units scaled, 1 timestep = {model_time_scale} ms)')
            else:
                ax2.set_ylabel('aCCE (model units)')
        else:
            # single shared axis: combine limits so both experiment and model values fit
            combined_vals = list(exp_vals) + list(model_vals)
            combined_errs = None
            if use_e_bars:
                combined_errs = list(exp_errs) + list(model_errs)
            combined_ylim = compute_limits(combined_vals, combined_errs)
            combined_ylim = include_zero(combined_ylim)
            ax.set_ylim(combined_ylim)
            # ensure ax2 (same as ax) also has same limits
            ax2.set_ylim(combined_ylim)
            ax.set_ylabel('aCCE (shared units)')

        ax.set_xticks(cluster_centers)
        ax.set_xticklabels(clusters)
        ax.set_title(f'SOA / Timestep: {soa} / t={timestep}')
        # legend for clusters (conditions)
        handles_cluster = [mpatches.Patch(color=colors[i], label=labels_cluster[i]) for i in range(n_clusters)]
        leg1 = ax.legend(handles=handles_cluster, loc='upper left')
        ax.add_artist(leg1)
        # add orientation legend (alpha/hatch indicates upright vs upside-down and mean)
        patch_u = mpatches.Patch(color='gray', alpha=0.95, label='Upright')
        patch_ud = mpatches.Patch(color='gray', alpha=0.6, label='Upside-down')
        patch_mean = mpatches.Patch(color='gray', alpha=0.85, hatch='//', label='Mean (upr+ud)/2')
        ax.legend(handles=[patch_u, patch_ud, patch_mean], loc='upper right')

    plt.tight_layout()
    # save plot
    savedir = "Exp-Results-search/Plots/"
    plotname = "compare_plot"
    if title_part is not None:
        plotname = plotname + "_" + title_part
    # If we already have a plot saved under this name, 'count up', with compare_plot_1.png, compare_plot_2.png etc.
    if os.path.exists(os.path.join(savedir, plotname + ".png")):
        i = 1
        while os.path.exists(os.path.join(savedir, f"{plotname}_{i}.png")):
            i += 1
        #plt.savefig(os.path.join(savedir, f"{plotname}_{i}.png"))
        plotname = f"{plotname}_{i}"
    plt.savefig(savedir + plotname + ".png", bbox_inches='tight', dpi=300)
    plt.show()


if __name__ == '__main__':

    # experimental file - try the path in this script first, otherwise ask user to place file there
    exp_file = exp_data_rts
    model_no = "basic_sampling_dlr0.03_prs0.01_div_0_r_1_decay0.1_b3_100bins_t_lr_0.01_sfn_0.01_nb_0.005_t_reach_100_100_t_bottle_0_stim_0"
    model_with = "basic_sampling_dlr0.03_prs0.1_div_0.1_r_1_decay0.1_b3_100bins_t_lr_0.01_sfn_0.01_nb_0.005_t_reach_100_100_t_bottle_0_stim_0"
    id_f = "base_t-100-100-100_t0_t50"
    alpha_lr = 0.05
    threshold = 0.9
    t_reach = 100
    print('Using model folders:', model_no, model_with)
    print('Experimental data file:', exp_file)
    # Toggle: use two y-axes (experiment left, model right) or a single shared y-axis for both
    use_two_axes =  False #True  # set to False to plot both experiment and model on the same y-axis
    # Scale the model automatically? Set to a number (ms per model timestep) to scale model -> experiment units,
    # or None to keep model units as-is. If provided and use_two_axes=True, the model axis will use the same limits
    # as the experiment axis so bars look identical (right axis then for illustration).
    model_t_scale = 3.0
    if model_t_scale is not None:
        id_f += f"_scale{model_t_scale}"
    # Toggle: Using error bars or not
    use_e_bars = False
    plot_compare(model_no, model_with, exp_file, timesteps=(t_model_onset, t_model_partway), title_part=id_f,
                 threshold=threshold, alpha_lr=alpha_lr, t_reach=t_reach, use_two_axes=use_two_axes,
                 model_time_scale=model_t_scale, use_e_bars=use_e_bars)