Spatial Map¶

epoch_df = loading.load_epoch(basepath)
beh_epochs = nel.EpochArray(np.array([epoch_df.startTime, epoch_df.stopTime]).T)

# you can change this based on which session you want to look at
behavior_idx = 11

# print out data frame
epoch_df

spatial_maps = maps.SpatialMap(
    pos[beh_epochs[behavior_idx]],
    st[beh_epochs[behavior_idx]],
    s_binsize=3,
    speed=speed,
    tuning_curve_sigma=3,
    place_field_min_size=15,
    place_field_max_size=1000,
    place_field_sigma=3,
)
spatial_maps.find_fields()

spatial_maps

[k for k in dir(spatial_maps.tc) if not k.startswith('_')]

# pull out velocity restricted position
current_pos = pos[beh_epochs[behavior_idx]][spatial_maps.run_epochs]

# iterate over the tuning curves
for ratemap_i, ratemap__ in enumerate(spatial_maps.tc.ratemap):
    ratemap_ = ratemap__.copy()
    ratemap_[spatial_maps.tc.occupancy < 0.01] = np.nan
    fig, ax = plt.subplots(1, 2, figsize=(10, 5))
    sns.heatmap(ratemap_.T, ax=ax[0])

    field_ids = np.unique(spatial_maps.tc.field_mask[ratemap_i])

    # plot field boundaries over the heat map
    if len(field_ids) > 1:
        for field_i in range(len(field_ids) - 1):
            bc = measure.find_contours(
                (spatial_maps.tc.field_mask[ratemap_i] == field_i + 1).T,
                0,
                fully_connected="low",
                positive_orientation="low",
            )
            for c in bc:
                ax[0].plot(c[:, 1], c[:, 0], linewidth=4)
    ax[0].invert_yaxis()

    # plot xy coords for animal position
    npl.plot2d(current_pos, lw=2, c="0.8", ax=ax[1])

    # plot xy coords for each spike
    x_time, pos_at_spikes = pos.asarray(
        at=st[beh_epochs[behavior_idx]][spatial_maps.run_epochs].data[ratemap_i]
    )
    ax[1].plot(pos_at_spikes[0, :], pos_at_spikes[1, :], ".", color="k")
    ax[1].set_aspect("equal")
    ax[0].set_aspect("equal")
    ax[1].set_title(f"{cm.iloc[ratemap_i].brainRegion} {cm.iloc[ratemap_i].UID}")
    ax[0].set_title(
        f"field width, {spatial_maps.tc.field_width[ratemap_i]:.2f}. peak rate {spatial_maps.tc.field_peak_rate[ratemap_i]:.2f}."
    )
    plt.show()

# change to the first linear track session
behavior_idx = 4

# print out data frame
epoch_df

pos = nel.AnalogSignalArray(
    data=position_df["linearized"].values.T,
    timestamps=position_df.timestamps.values,
)
pos = pos[beh_epochs[behavior_idx]]
speed = nel.utils.ddt_asa(pos, smooth=True, sigma=0.250, norm=True)
pos

plt.figure(figsize=(20, 4))
plt.plot(pos.data.T)

# get outbound and inbound epochs
(outbound_epochs, inbound_epochs) = get_linear_track_lap_epochs(
    pos.abscissa_vals, pos.data[0], newLapThreshold=20
)


outbound_epochs, inbound_epochs

plt.figure(figsize=(20, 4))
plt.scatter(
    pos[outbound_epochs].abscissa_vals,
    pos[outbound_epochs].data,
    color="r",
    s=1,
    label="outbound_epochs",
)
plt.scatter(
    pos[inbound_epochs].abscissa_vals,
    pos[inbound_epochs].data,
    color="b",
    s=1,
    label="inbound_epochs",
)
plt.legend()

tc = spatial_maps.tc

w, h = set_size("thesis", fraction=0.3, subplots=(3, 1))

with npl.FigureManager(show=True, figsize=(w, h)) as (fig_, ax):
    npl.utils.skip_if_no_output(fig_)
    ax1 = npl.plot_tuning_curves1D(
        tc.reorder_units(), normalize=True, pad=1, fill=True, alpha=0.5
    )

    leg_lines = ax1.get_lines()
    plt.setp(leg_lines, linewidth=0.5)
    ax.set_xlabel("position (cm)")
    ax.set_ylabel("n pyr cells")
    ax.set_xticks([0, tc.bins.max()])

for ratemap_i, ratemap__ in enumerate(spatial_maps.tc.ratemap):
    ratemap_ = ratemap__.copy()
    x = np.arange(len(ratemap_)) * spatial_maps.s_binsize
    plt.figure(figsize=(5, 2))
    plt.plot(x, ratemap_, color="k")
    # field_ids = np.unique(spatial_maps.tc.field_mask[ratemap_i])
    plt.plot(x, spatial_maps.tc.field_mask[ratemap_i], color="r")
    plt.title(
        f"field width, {spatial_maps.tc.field_width[ratemap_i]:.2f}. peak rate {spatial_maps.tc.field_peak_rate[ratemap_i]:.2f}."
    )
    plt.show()

# pull out velocity restricted position
current_pos = pos[beh_epochs[behavior_idx]][spatial_maps.run_epochs]

# iterate over the tuning curves
for ratemap_i, ratemap__ in enumerate(spatial_maps.tc.ratemap):
    ratemap_ = ratemap__.copy()
    ratemap_[spatial_maps.tc.occupancy < 0.01] = np.nan
    fig, ax = plt.subplots(
        2, 1, figsize=(5, 3), sharex=True, gridspec_kw={"height_ratios": [1, 3]}
    )

    # plot tuning curve
    x = (spatial_maps.x_edges + np.abs(np.diff(spatial_maps.x_edges)[0] / 2))[:-1]
    ax[0].plot(x, ratemap_, color="k")

    # plot xy coords for animal position
    ax[1].plot(current_pos.data.T, current_pos.abscissa_vals, "grey")

    # plot xy coords for each spike
    x_time, pos_at_spikes = pos.asarray(
        at=st[beh_epochs[behavior_idx]][spatial_maps.run_epochs].data[ratemap_i]
    )
    ax[1].plot(pos_at_spikes.flatten(), x_time, ".", color="k")

    ax[1].set_xlim(x.min(), x.max())

    # add figure labels
    ax[0].set_ylabel("rate (Hz)")
    ax[1].set_ylabel("time (sec)")
    ax[1].set_xlabel("position (cm)")

    ax[0].set_title(
        f"field width, {spatial_maps.tc.field_width[ratemap_i]:.2f}. peak rate {spatial_maps.tc.field_peak_rate[ratemap_i]:.2f} {cm.iloc[ratemap_i].brainRegion} {cm.iloc[ratemap_i].UID}"
    )
    sns.despine()
    plt.show()

basepath = r"Z:\Data\HMC1\day8"

# load position
position_df = loading.load_animal_behavior(basepath)

# put position into a nelpy position array for ease of use
pos = nel.AnalogSignalArray(
    data=position_df["linearized"].values.T,
    timestamps=position_df.timestamps.values,
)

behavior_idx = 1

# load epoch data
epoch_df = loading.load_epoch(basepath)
beh_epochs = nel.EpochArray(np.array([epoch_df.startTime, epoch_df.stopTime]).T)

pos = pos[beh_epochs[behavior_idx]]

# get outbound and inbound epochs
outbound_epochs, inbound_epochs = get_linear_track_lap_epochs(
    pos.abscissa_vals, pos.data[0], newLapThreshold=20
)

# calculate speed
speed = nel.utils.ddt_asa(pos, smooth=True, sigma=0.250, norm=True)
run_epochs = nel.utils.get_run_epochs(speed, v1=4, v2=4).merge()


assembly_react = AssemblyReact(
    basepath=basepath,
    brainRegion="CA1",
    putativeCellType="Pyr",
    z_mat_dt=0.01,
)

# load need data (spikes, ripples, epochs)
assembly_react.load_data()

# load theta epochs
state_dict = loading.load_SleepState_states(basepath)
theta_epochs = nel.EpochArray(
    state_dict["THETA"],
)

assembly_react.get_weights(epoch=beh_epochs[behavior_idx] & theta_epochs)
assembly_react

assembly_react.plot()
plt.show()

assembly_act = assembly_react.get_assembly_act(epoch=beh_epochs[behavior_idx])
assembly_act

for tc in [spatial_map_outbound_assembly.tc, spatial_map_inbound_assembly.tc]:
    w, h = set_size("thesis", fraction=0.3, subplots=(3, 1))

    with npl.FigureManager(show=True, figsize=(w, h)) as (fig_, ax):
        npl.utils.skip_if_no_output(fig_)
        ax1 = npl.plot_tuning_curves1D(
            tc.reorder_units(), normalize=True, pad=1, fill=True, alpha=0.5
        )

        leg_lines = ax1.get_lines()
        plt.setp(leg_lines, linewidth=0.5)
        ax.set_xlabel("position (cm)")
        ax.set_ylabel("n assemblies")
        ax.set_xticks([0, tc.bins.max()])