<jupyter_start><jupyter_text>Calculate MI for each species and plot goodness of fit by length of analysis1. load datasets2. calculate MI<jupyter_code>import pandas as pd
import numpy as np
from parallelspaper.config.paths import DATA_DIR, FIGURE_DIR
from parallelspaper.birdsong_datasets import MI_seqs, compress_seq, BCOL_DICT
from parallelspaper import information_theory as it 
from tqdm.autonotebook import tqdm
from parallelspaper import model_fitting as mf
from parallelspaper.utils import save_fig
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>starling_seq_df = pd.read_pickle(DATA_DIR / 'song_seq_df/starling.pickle')
CAVI_CATH_seq_df = pd.read_pickle(DATA_DIR / 'song_seq_df/CAVI_CATH.pickle')
BF_seq_df = pd.read_pickle(DATA_DIR / 'song_seq_df/BF.pickle')
seq_dfs = pd.concat([starling_seq_df, CAVI_CATH_seq_df, BF_seq_df])
# sequence lengths
seq_dfs['sequence_lens'] = [len(i) for i in seq_dfs.syllables]
# recording number as integer
seq_dfs['rec_num'] = seq_dfs.rec_num.values.astype('int32')
# sort sequences
seq_dfs = seq_dfs.sort_values(by=['species', 'bird', 'rec_num'])
# get rid of unID'd birds (CAVI, CATH)
seq_dfs = seq_dfs[seq_dfs.bird != '?']
len(seq_dfs)
seq_dfs[:3]<jupyter_output><empty_output><jupyter_text>Calculate MI<jupyter_code>distances = np.arange(1,1001)
verbosity = 0; n_jobs = 20
np.unique(seq_dfs.species)
MI_DF = pd.DataFrame(columns=['species', 'type', 'MI', 'MI_shuff', 'distances',
                              'MI_var', 'MI_shuff_var'])

for species in np.unique(seq_dfs.species):
    species_df = seq_dfs[seq_dfs.species == species].sort_values(by=['bird', 'rec_num'])
    print(species)
   
    # analysis by day
    day_group = []
    for bird in np.unique(species_df.bird.values):
        bird_df = species_df[species_df.bird==bird]
        for day in np.unique(bird_df.day.values):
            day_df = bird_df[bird_df.day == day]
            day_group.append(np.concatenate(day_df.syllables.values))
    units = day_group
    (MI, var_MI), (MI_shuff, MI_shuff_var) = it.sequential_mutual_information(units, distances, n_jobs = n_jobs, verbosity = verbosity)
    MI_DF.loc[len(MI_DF)] = [species, 'session', MI, MI_shuff, distances, var_MI, MI_shuff_var]
MI_DF.to_pickle(DATA_DIR / 'MI_DF/birdsong/birdsong_MI_DF_long.pickle')
MI_DF<jupyter_output><empty_output><jupyter_text>Calculate fit<jupyter_code>from sklearn.externals.joblib import Parallel, delayed
n_jobs = 20; verbosity = 0
def get_fit(species, d, distances, sig):
    results_power, results_exp, results_pow_exp, best_fit_model = mf.fit_models(
        distances[:d], sig[:d])
    R2_exp, R2_concat, R2_power, AICc_exp, AICc_pow, AICc_concat = mf.fit_results(
        sig[:d], distances[:d],  results_exp, results_power, results_pow_exp)

    y_model = mf.get_y(mf.pow_exp_decay, results_pow_exp, distances)
    y_pow = mf.get_y(mf.powerlaw_decay, results_pow_exp, distances)
    y_exp = mf.get_y(mf.exp_decay, results_pow_exp, distances)

    R2_exp_comp = mf.r2(sig[:d] - y_pow[:d], y_exp[:d] -
                        results_pow_exp.params['intercept'].value, distances[:d], logscaled=True)
    s = sig[:d] - y_exp[:d]
    m = y_pow[:d]-results_pow_exp.params['intercept'].value
    mask = s > 0
    R2_pow_comp = mf.r2(s[mask], m[mask], distances[:d][mask], logscaled=True)
    # print(R2_pow_comp)
    #plt.plot(distances[:d], mf.residuals(s, m,distances[:d]))

    AICc_exp_comp = mf.AICc(d, len(results_exp.params), sig[:d] - y_pow[:d], y_exp[:d] -
                            results_pow_exp.params['intercept'].value, distances[:d], logscaled=True)
    AICc_pow_comp = mf.AICc(d, len(results_power.params),
                            sig[:d] - y_exp[:d], y_pow[:d]-results_pow_exp.params['intercept'].value, distances[:d], logscaled=True)
    return (species, d, R2_exp, R2_concat, R2_power, AICc_exp, AICc_pow, 
            AICc_concat, R2_pow_comp, R2_exp_comp, AICc_exp_comp, AICc_pow_comp)
# aic / r2 for individual components
fit_df = []

columns = ['species', 'd', 'R2_exp', 'R2_concat', 'R2_power', 'AICc_exp', 'AICc_pow', 
                                 'AICc_concat', 'R2_pow_comp', 'R2_exp_comp',  'AICc_exp_comp', 'AICc_pow_comp']

for axi, (idx, row) in enumerate(MI_DF.sort_values(by=['species']).iterrows()):
    species = row.species
    sig = row.MI-row.MI_shuff
    with Parallel(n_jobs=n_jobs, verbose=verbosity) as parallel:
        x = parallel(
            delayed(get_fit)(species, d, row.distances, sig)
                 for d in tqdm(np.unique(np.linspace(16,1000, 200).astype(int))))
    
    fit_df_lang = pd.DataFrame(x, columns = columns)
    fit_df.append(fit_df_lang)
fit_df = pd.concat(fit_df)
fit_df.to_pickle(DATA_DIR / 'MI_DF/birdsong/birdsong_fit_df_long.pickle')
fit_df[:3]<jupyter_output><empty_output><jupyter_text>Plot results<jupyter_code>import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>R2 full concatenative model<jupyter_code>ncol = len(np.unique(fit_df.species))
zoom = 4
fig, axs = plt.subplots(ncols=ncol, figsize= (ncol*zoom, zoom))
for sax, species in enumerate(np.unique(fit_df.species)):
    color = BCOL_DICT[species]
    ax = axs.flatten()[sax]
    spec_fit_df = fit_df[fit_df.species == species]
    
    ax.plot(spec_fit_df.d, spec_fit_df.R2_concat.values, lw=4, color=color)
    ax.set_title(species)
    ax.set_ylim([0.5,1.01])
    ax.set_xlim([np.min(spec_fit_df.d), np.max(spec_fit_df.d)])
    spec_fit_df = spec_fit_df[spec_fit_df.d > 100]
    d = spec_fit_df.d.values[np.where(spec_fit_df.R2_concat.values > (spec_fit_df.R2_concat.values[0]*.999))[0][-1]]
    print(d)
    #ax.axvline(d,color='k', ls='dashed', alpha=0.5)
    #ax.set_ylim(np.exp(sig_lims))
    ax.tick_params(which='both', direction='in')
    ax.tick_params(which='major', length=10, width =3)
    ax.tick_params(which='minor', length=5, width =2)
    ax.set_xlabel('Distance (syllables)', fontsize=18)
    ax.set_xscale( "log" , basex=10)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(3)
        ax.spines[axis].set_color('k')

axs[0].set_ylabel('$r^2$ concat. model', fontsize=18)
plt.tight_layout()

save_fig(FIGURE_DIR/'R2_song')<jupyter_output><empty_output><jupyter_text>R2 of power-law compenent<jupyter_code>ncol = len(np.unique(fit_df.species))
zoom = 4
fig, axs = plt.subplots(ncols=ncol, figsize= (ncol*zoom, zoom))
for sax, species in enumerate(np.unique(fit_df.species)):
    color = BCOL_DICT[species]
    ax = axs.flatten()[sax]
    spec_fit_df = fit_df[fit_df.species == species]
    ax.plot(spec_fit_df.d, spec_fit_df.R2_pow_comp.values, lw=4, color=color)
    #ax.plot(spec_fit_df.d, spec_fit_df.R2_pow_comp.values, lw=4, color=color)
    ax.set_title(species)
    ax.set_ylim([0.5,1.01])
    ax.set_xlim([np.min(spec_fit_df.d), np.max(spec_fit_df.d)])
    
    #ax.axvline(d,color='k', ls='dashed', alpha=0.5)
    #ax.set_ylim(np.exp(sig_lims))
    ax.tick_params(which='both', direction='in')
    ax.tick_params(which='major', length=10, width =3)
    ax.tick_params(which='minor', length=5, width =2)
    ax.set_xlabel('Distance (syllables)', fontsize=18)
    ax.set_xscale( "log" , basex=10)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(3)
        ax.spines[axis].set_color('k')

axs[0].set_ylabel('$r^2$ power-law component', fontsize=18)
plt.tight_layout()


save_fig(FIGURE_DIR/'r2_powerlaw_song')<jupyter_output><empty_output><jupyter_text>$\Delta$AICc concat vs expon<jupyter_code>ncol = len(np.unique(fit_df.species))
zoom = 4
fig, axs = plt.subplots(ncols=ncol, figsize= (ncol*zoom, zoom))
for sax, species in enumerate(np.unique(fit_df.species)):
    color = BCOL_DICT[species]
    ax = axs.flatten()[sax]
    spec_fit_df = fit_df[fit_df.species == species]
    ax.plot(spec_fit_df.d, spec_fit_df.AICc_concat.values - spec_fit_df.AICc_exp.values, lw=4, color=color)
    ax.set_title(species)
    ax.tick_params(which='both', direction='in')
    ax.tick_params(which='major', length=10, width =3)
    ax.tick_params(which='minor', length=5, width =2)
    ax.set_xlabel('Distance (syllables)', fontsize=18)
    ax.set_xscale( "log" , basex=10)
    for axis in ['top','bottom','left','right']:
        ax.spines[axis].set_linewidth(3)
        ax.spines[axis].set_color('k')
    ax.set_xlim([np.min(spec_fit_df.d), np.max(spec_fit_df.d)])

axs[0].set_ylabel('$\Delta$AICc (concat.- exp.)', fontsize=18)
plt.tight_layout()

save_fig(FIGURE_DIR/'delta_AIC_song')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14952614952614954, 'fraction_numerical': 0.014742014742014743, 'mean_word_length': 4.74076561450638, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '606872', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3564, 'n_tokens_neox': 3435, 'n_words': 528, 'path': 'notebooks/birdsong/07.0-birdsong-plot-goodness-of-fit-by-analysis-length.ipynb', 'repo_name': 'timsainb/ParallelsBirdsongLanguagePaper', 'size': '8547'}