In [1]:
import pandas as pd
import nibabel as nib
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%load_ext autoreload
%autoreload 2
import sys
sys.path.append('..')
import sfp
import math
from scipy import stats
from scipy import optimize as opt
import torch
import glob
import os
import warnings
import itertools
import hessian
from torch.utils import data as torchdata
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
sns.set_style('whitegrid')
#df_path = "/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/derivatives/first_level_analysis/stim_class/posterior/sub-wlsubj045/ses-04/sub-wlsubj045_ses-04_task-sfprescaled_v1_e1-12_summary.csv"
df_path = "/mnt/prince_scratch/spatial_frequency_preferences/derivatives/first_level_analysis/stim_class/bayesian_posterior/sub-wlsubj045/ses-04/sub-wlsubj045_ses-04_task-sfprescaled_v1_e1-12_summary.csv"
In [2]:
sns.set_style('white', {'axes.spines.right': False, 'axes.spines.top': False})
In [5]:
# Load in data
df = pd.read_csv(df_path)
Let's look at some voxels that have good GLM $R^2$ values and pick one
In [11]:
df[(df.varea==1)&(df.GLM_R2>40)].drop_duplicates('voxel').sort_values('GLM_R2', ascending=False).head(20)[['voxel', 'GLM_R2', 'varea', 'hemi', 'angle' ,'eccen', 'precision', 'prf_vexpl']]
Out[11]:
In [12]:
# Pick a V1 voxel with a good R2
voxel_df = df[(df.voxel.isin([3269]))]#, 1421]))]
voxel_df.head()
Out[12]:
Let's examine the response of this voxel as a function of spatial frequency. In the plot below, we plot the normed amplitude estimate as a function of the local spatial frequency. We see that the response looks roughly log-Normal, but there appears to be some difference between the different stimulus classes.
In [232]:
pal = sfp.plotting.stimulus_type_palette('relative')
hue_order = sfp.plotting.stimulus_type_order('relative')
In [15]:
g = sns.FacetGrid(voxel_df[~voxel_df.stimulus_superclass.isin(['mixtures'])], hue='stimulus_superclass',palette=pal, size=8, aspect=1.5, hue_order=hue_order)
g.map(plt.scatter, 'local_sf_magnitude', 'amplitude_estimate_median', linewidth=6)
g.map(plt.plot, 'local_sf_magnitude', 'amplitude_estimate_median', linewidth=6)
# g.ax.set_xscale('log', basex=2)
g.add_legend()
g.ax.tick_params(size=0)
g.ax.set_xlim((0, 3))
g.set_xlabels('Local spatial frequency (cpd)')
g.set_ylabels('Amplitude Estimate')
Out[15]:
These classes differ in their local orientation, so we can look at a plot of the response as a function of the local spatial frequency with respect to x and y (size represents the response). This plot is difficult to parse, but the main point is that these different stimulus classes are not arbitrary and discrete: they lie on a continuum, related by the stimulus orientation, and so we can fit the response of the voxel as a 2d tuning curve.
In [16]:
def scatter_sizes(x, y, s, plot_color=False, cmap=None, size_scale=1, **kwargs):
if plot_color:
kwargs.pop('color')
if cmap is None:
cmap = 'Blues'
plt.scatter(x, y, s=s*80*size_scale, c=s, cmap=cmap, **kwargs)
else:
plt.scatter(x, y, s=s*80*size_scale, **kwargs)
with sns.axes_style('whitegrid'):
voxel_df['normalized_resp'] = voxel_df['amplitude_estimate_median'].copy()
voxel_df['normalized_resp'] = (voxel_df['normalized_resp'] - voxel_df['normalized_resp'].min()) / (voxel_df['normalized_resp'].max() - voxel_df['normalized_resp'].min())
g=sns.FacetGrid(voxel_df, size=8, aspect=1, hue='stimulus_superclass', palette=pal, hue_order=hue_order)
g.map(scatter_sizes, 'local_w_x', 'local_w_y', 'normalized_resp', plot_color=False, size_scale=2)
g.add_legend()
scatter_ax = plt.gca()
scatter_ax.set_aspect('equal')
g.ax.tick_params(size=0)
This is just that same data, but rotated and plotted on a semi-log plots.
In [17]:
with sns.axes_style('whitegrid'), sns.plotting_context('notebook'):
g=sns.FacetGrid(voxel_df, hue='stimulus_superclass', size=5, aspect=1, palette=pal, hue_order=hue_order)
g.map(scatter_sizes, 'local_w_r', 'local_w_a', 'normalized_resp')
g.add_legend()
scatter_ax = plt.gca()
scatter_ax.set_xscale('symlog', basex=2, linthreshx=2**(-4))
scatter_ax.set_yscale('symlog', basey=2, linthreshy=2**(-4))
But then the question is: how does the tuning change with orientation? Two possibilities are:
Then there's the question of how either the mode or the amplitude changes with orientation. Let's assume it changes smoothly and periodically, symmetrically about 180 degrees (because 2d orientation is runs from 0 to 180 degrees). We'll examine three possibilities:
I'll build up the complexity of the model, adding terms (and plots) to show additional possibilities that we're modeling. However, note that there are several things we are not modeling:
There are lots of ways to show this. For our purposes, we're going to show three types of plots:
We'll build up the complexity of these examples, starting with the simplest and getting more complex as we go.
As we mentioned, we noticed the response looked roughly log-normal, so we start with that:
$\hat{\beta} = A \exp\left(\frac{-(\log_2\omega+\log_2p)^2}{2\sigma^2} \right)$
Here:
GLMdenoise) of a given voxel to a given stimulusThis here is what we were fitting in our first pass analysis, as seen in notebook 04-MRI-first-level; for each eccentricity band and stimulus class, we fit all three parameters ($p$, $\sigma$, and $A$) separately, and then plotted them out as a function of eccentricity to see their relationship.
Now, we're interested in parameterizing the relationship between these parameters and the voxel's eccentricity, retinotopic angle, and the stimulus's local orientation information. The simplest way we'll do this, corresponding to the largest effect we saw, is to say that retinotopic angle and stimulus local orientation don't matter and that the preferred period is some linear function of eccentricity, $e$:
$p=ae+b$
where we fit $a$ and $b$ across all voxels and stimuli. Note that we now fit the same $\sigma$ to all voxels (but since it's in octaves, there's an implicit scaling here). We don't actually fit $A$ here, since we normalize each voxel's response in our loss function, so any constant value for $A$ will be equivalent.
In [35]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=0)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
For this simple model, there's on effect of retinotopic angle or stimulus orientation, so we won't look at any other plots.
However, we can look at the effect of changing these two parameters, $a$ and $b$, which do what you expect:
In [34]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.5, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=0)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Increased $a$', ylim=(-.1, 2.5))
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.05, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=0)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Decreased $a$', ylim=(-.1, 2.5))
In [36]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.4, orientation_type='full', abs_mode_cardinals=0)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Increased $b$', ylim=(-.1, 2.5))
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=0, orientation_type='full', abs_mode_cardinals=0)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Decreased $b$', ylim=(-.1, 2.5))
Alright, now we're going to make things more complicated by parameterizing an effect of stimulus orientation and voxel retinotopic angle on the preferred period:
$\hat{\beta} = A \exp\left(\frac{-(\log_2\omega+\log_2p)^2}{2\sigma^2} \right)$
$p=(ae+b)(1+p_1\cos 2\theta + p_2\cos 4\theta + p_3\cos 2(\theta-\phi) + p_4\cos 4(\theta-\phi))$
where all parameters are as before and we now have two new variables:
and the new $p_i$ parameters control the size of the effects on the preferred period:
Note that we're saying the effect of eccentricity and retinotopic angle / orientation are multiplied together, that they have separable effects on the preferred period.
Let's first consider a model that has $p_3=p_4=0$ and $p_1>p_2>0$, that is, where there's only an effect of absolute orientation and the effect of vertical vs. horizontal is larger than that of cardinals vs. obliques.
In [562]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=.2, abs_mode_obliques=.1)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
Now we'll look at these plots again, to see what happens when we increase $p_1$
In [69]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=.4, abs_mode_obliques=.1)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
... when we make $p_1=p_2$
In [70]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=.2, abs_mode_obliques=.2)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
... and we make $p_1<p_2$
In [71]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', abs_mode_cardinals=.2, abs_mode_obliques=.4)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
And now let's set $p_1=p_2=0$ and look at the effects of $p_3$ and $p_4$, with $p_3 > p_4$
In [72]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', rel_mode_cardinals=.2, rel_mode_obliques=.1)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
... when we increase $p_3$
In [73]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', rel_mode_cardinals=.4, rel_mode_obliques=.1)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
... when $p_3 = p_4$
In [74]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', rel_mode_cardinals=.2, rel_mode_obliques=.2)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
... when $p_3 < p_4$
In [75]:
m = sfp.model.LogGaussianDonut(vary_amplitude=False, sf_ecc_slope=.2, sf_ecc_intercept=.2, orientation_type='full', rel_mode_cardinals=.2, rel_mode_obliques=.4)
tmp = sfp.analyze_model.create_preferred_period_df(m, reference_frame='relative',)
g = sfp.plotting.feature_df_plot(tmp, col=None, title='Preferred period averaged across all retinotopic angles', ylim=(-.1, 2.5))
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='relative', orientation=np.linspace(0, np.pi, 2, endpoint=False),)
g = sfp.plotting.feature_df_polar_plot(tmp, )
tmp = sfp.analyze_model.create_preferred_period_contour_df(m, reference_frame='absolute', orientation=np.linspace(0, np.pi, 2, endpoint=False))
g = sfp.plotting.feature_df_polar_plot(tmp)
And we do the same thing with amplitude, parameterizing it as:
$\hat{\beta} = A \exp\left(\frac{-(\log_2\omega+\log_2p)^2}{2\sigma^2} \right)$
$p=(ae+b)(1+p_1\cos 2\theta + p_2\cos 4\theta + p_3\cos 2(\theta-\phi) + p_4\cos 4(\theta-\phi))$
$A=(1+A_1\cos 2\theta + A_2\cos 4\theta + A_3\cos 2(\theta-\phi) + A_4\cos 4(\theta-\phi))$
In [2]:
model_names = ['constant iso', 'scaling iso', 'full iso', 'full absolute', 'full relative', 'full full', 'full absolute amps', 'full relative amps', 'full full amps']
parameters = [r'$\sigma$', r'$a$', r'$b$', r'$p_1$', r'$p_2$', r'$p_3$', r'$p_4$', r'$A_1$', r'$A_2$', r'$A_3$', r'$A_4$']
model_variants = np.zeros((len(model_names), len(parameters))).astype(bool)
# everyone fits sigma
model_variants[:, 0] = True
model_variants[1:, 1] = True
model_variants[0, 2] = True
model_variants[2:, 2] = True
model_variants[3, [3, 4]] = True
model_variants[4, [5, 6]] = True
model_variants[5, [3, 4, 5, 6]] = True
model_variants[6, [3, 4, 7, 8]] = True
model_variants[7, [5, 6, 9, 10]] = True
model_variants[8, 3:] = True
In [5]:
model_variants = pd.DataFrame(model_variants, model_names, parameters)
green, red = sns.color_palette('deep', 4)[2:]
pal = sns.blend_palette([red, green])
In [6]:
with sns.plotting_context('poster', font_scale=1.5):
fig = plt.figure(figsize=(12,10))
ax = sns.heatmap(model_variants, cmap=pal, cbar=False, )
ax.set_yticklabels(model_names, rotation=0)
One of the things we care about is the noisiness of our voxels. We'll use this in our loss function, so we want to know whether we can just use the noise we've calculated on a voxel-by-voxel basis or whether we have to worry about differing levels of noise for different stimulus classes.
In [7]:
df = pd.read_csv(df_path)
What's unclear from the following plot is whether the apparent relationship between the median amplitude estimate and the error is because of differences between voxels or within them. If it's between voxels, we're fine, but it's bad if it's within them (this would be like for a Poisson-process neuron, whose variance increases with the mean firing rate)
In [8]:
sns.relplot(x='amplitude_estimate_median', y='amplitude_estimate_std_error', data=df)
Out[8]:
So to answer this question, we do a simple linear regression between the median amplitude estimate and the standard error.
In [9]:
def linear_regression(data, x='amplitude_estimate_median', y='amplitude_estimate_std_error'):
s,i,r,p,e = stats.linregress(data[x], data[y])
return pd.Series({'slope': s, 'intercept': i, 'r_squared': r**2, 'p_value': p, 'stderr': e})
linregress = df.groupby('voxel').apply(linear_regression)
And what we see here is that there doesn't appear to be a consistent relationship between the median amplitude estimate and the error of those estimates: the slope looks like it's distributed around zero, as does the relationship between slope and intercept. This is good, so we can treat precision as a voxel-by-voxel number
In [10]:
sns.pairplot(linregress, vars=['slope', 'intercept'], height=5)
Out[10]:
Since we're happy with that, we'll use the following loss function, evaluated for each voxel $v$ (NOTE: probably should add $_v$ to all these values):
$L(\beta,\hat{\beta})=\frac{1}{\sigma_v^2}\sum_{i=1}^{n}\frac{1}{n}\left(\frac{\beta_i}{||\beta||_2} - \frac{\hat{\beta}_i}{||\hat{\beta}||_2} \right)^2$
where $i$ indexes the $n$ different stimulus classes, $\beta_i$ is the response (estimated using GLMdenoise) to stimulus class $i$, $\hat{\beta}_i$ is the response to stimulus class $i$ predicted by our model, $||\beta||_2$ is the L2-norm of $\beta$ (for this voxel, across all stimulus classes), and $\sigma_v^2$ is the variance of the voxel's response (that is, $\sigma_v^2 = \frac{1}{n}\sum_{i=1}^n\sigma^2_{vi}$, where $\sigma_{vi}$ is half of the 68 percentile range of the response of voxel $v$ to stimulus class $i$, as estimated by GLMdenoise).
If we didn't have that $\frac{1}{\sigma_v}$ in the beginning, this would be $\frac{2}{n}(1-\cos(\beta,\hat{\beta}))$, where $\cos(A,B)$ is the cosine similarity between two vectors $A$ and $B$. It's similar to the Pearson correlation, except that it's not de-meaned (see this StackOverflow post for a bit more information)
To see how:
$\sum_i\frac{1}{n}\left(\frac{\beta_i}{||\beta||} - \frac{\hat{\beta}_i}{||\hat{\beta}||} \right)^2$
$\frac{1}{n}\sum_i\left(\frac{\beta_i^2}{||\beta||^2} - \frac{2\beta_i\hat{\beta}_i}{||\hat{\beta}|| ||\beta||} + \frac{\hat{\beta}_i^2}{||\hat{\beta}||^2}\right)$
$\frac{1}{n}\left(\sum_i\frac{\beta_i^2}{||\beta||^2} + \sum_i\frac{\hat{\beta}_i^2}{||\hat{\beta}||^2} - \sum_i\frac{2\beta_i\hat{\beta}_i}{||\hat{\beta}|| ||\beta||} \right)$
$\frac{1}{n}\left(\frac{\sum_i\beta_i^2}{||\beta||^2} + \frac{\sum_i\hat{\beta}_i^2}{||\hat{\beta}||^2} - \frac{2\sum_i\beta_i\hat{\beta}_i}{||\hat{\beta}|| ||\beta||} \right)$
$\frac{1}{n}\left(\frac{||\beta||^2}{||\beta||^2} + \frac{ ||\hat{\beta}||^2}{||\hat{\beta}||^2} - \frac{2\sum_i\beta_i\hat{\beta}_i}{||\hat{\beta}|| ||\beta||} \right)$
$\frac{1}{n}\left(2 - 2\frac{\sum_i\beta_i\hat{\beta}_i}{||\hat{\beta}|| ||\beta||} \right)$
$\frac{2}{n}\left(1 - \cos(\beta,\hat{\beta}) \right)$
For cross-validated loss, we fit model to 44 of the 48 classes using standard procedure (with above loss), then get predictions for the 4 held out classes. We do this for each of the 12 subsets, which gets us a complete $\hat{\beta}_v$ that we can compare against the actual $\beta_v$ using the above loss. Therefore, interpretation should be the same.
In the 05a-Simulations notebook, we look at some simulated data to show how we picked the learning hyperparameters (learning rate and batch size), as well as some simulations that reassure us about model recovery (finding the correct model and the correct parameter values). Thus reassured, we examine the actual data.
Now we look at all the cross-validated fits to subjects and sessions, using the pre-shrunk data computed by the snakemake rule summarize_model_cv (because the files are too big to hold all in memory).
In [151]:
base_str = '/users/broderick/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/derivatives/tuning_2d_model/stim_class/bayesian_posterior/initial_cv/v1_e1-12_summary_b10_r0.001_g0_s0_all_'
models = pd.read_csv(base_str + 'models.csv')
timing_df = pd.read_csv(base_str + 'timing.csv')
avg_loss = pd.read_csv(base_str + 'loss.csv')
cv_loss = pd.read_csv(base_str + 'cv_loss.csv')
In [111]:
model_order = ['constant_donut_iso_amps-constant', 'scaling_donut_iso_amps-constant', 'full_donut_iso_amps-constant',
'full_donut_absolute_amps-constant', 'full_donut_relative_amps-constant', 'full_donut_full_amps-constant',
'full_donut_absolute_amps-vary', 'full_donut_relative_amps-vary', 'full_donut_full_amps-vary']
In [112]:
avg_loss['indicator'] = avg_loss.apply(lambda x: str((x.subject, x.session)), 1)
models['indicator'] = models.apply(lambda x: str((x.subject, x.session)), 1)
cv_loss['indicator'] = cv_loss.apply(lambda x: str((x.subject, x.session)), 1)
timing_df['indicator'] = timing_df.apply(lambda x: str((x.subject, x.session)), 1)
indicator_order = cv_loss.indicator.unique()
In [113]:
pd.pivot_table(cv_loss, 'cv_loss', ['fit_model_type'], 'indicator')
Out[113]:
First thing we note is that the cross-validated loss varies a lot across subjects / scanning sessions and is generally pretty tightly packed within them. In order to examine which model fit the data best, we will demean them -- average the cross-validated loss for a given subject / scanning session across model types, and subtract that off.
In [116]:
g = sns.catplot('indicator', 'cv_loss', 'fit_model_type', data=cv_loss, hue_order=model_order,
legend='full', height=8, aspect=1.5, order=indicator_order, palette=sns.color_palette('deep', 9), )
for ax in g.axes.flatten():
labels = ax.get_xticklabels()
if labels:
ax.set_xticklabels(labels, rotation=25)
#g.set(ylim=(0, .15))
g.fig.subplots_adjust(bottom=.25)
Once we do that, we can start to see some more obvious patterns, but we want to view these as a function of model type, not of subjects / scanning sessions.
In [117]:
demeaned_loss = cv_loss.set_index(['indicator'])
demeaned_loss['loss_mean'] = cv_loss.groupby(['indicator']).cv_loss.mean()
demeaned_loss['demeaned_loss'] = (demeaned_loss.cv_loss - demeaned_loss.loss_mean)# / demeaned_all_loss.loss_mean
demeaned_loss = demeaned_loss.reset_index()
g = sns.catplot('indicator', 'demeaned_loss', 'fit_model_type', data=demeaned_loss, hue_order=model_order,
legend='full', height=6, aspect=1.5, order=indicator_order, palette=sns.color_palette('deep', 9),)
for ax in g.axes.flatten():
labels = ax.get_xticklabels()
if labels:
ax.set_xticklabels(labels, rotation=25)
ax.axhline(color='grey', linestyle='dashed')
#g.set(ylim=(-.005, .01))
g.fig.subplots_adjust(bottom=.25)
Here the patterns become more clear...
In [118]:
with sns.plotting_context('poster'), sns.color_palette('husl', demeaned_loss.indicator.nunique()):
g = sns.catplot('fit_model_type', 'demeaned_loss', 'indicator', data=demeaned_loss, s=15,
legend='full', height=17, aspect=.75, kind='strip', ci=95, estimator=np.median, order=model_order)
for ax in g.axes.flatten():
labels = ax.get_xticklabels()
if labels:
ax.set_xticklabels(labels, rotation=25)
ax.axhline(color='grey', linestyle='dashed')
# g.set(ylim=(-.002, .006))
g.fig.subplots_adjust(bottom=.25)
g.set(xlabel='Model type', ylabel='Cross-validated Loss (demeaned by subject)')
#g.fig.savefig('cv_loss_by_type.svg')
And by plotting the median with 95% bootstrapped confidence intervals, we can get a clearer picture: we get a very big gain by fitting both parameters $a$ and $b$ (instead of one or the other), which we expected from the 1d results seen in the 04-MRI-first-level notebook, and fitting all parameters does seem to be the best, though it appears that the absolute amplitudes are doing the heavy lifting.
In [119]:
with sns.plotting_context('poster'):
g = sns.catplot('fit_model_type', 'demeaned_loss', 'fit_model_type', data=demeaned_loss, s=15,
legend='full', height=17, aspect=.75, kind='point', ci=95, estimator=np.median, order=model_order, hue_order=model_order)
for ax in g.axes.flatten():
labels = ax.get_xticklabels()
if labels:
ax.set_xticklabels(labels, rotation=25)
ax.axhline(color='grey', linestyle='dashed')
# g.set(ylim=(-.002, .006))
g.fig.subplots_adjust(bottom=.25)
g.set(xlabel='Model type', ylabel='Cross-validated Loss (demeaned by subject)')
#g.fig.savefig('cv_loss_by_type.svg')
In [5]:
models = pd.read_csv('/users/broderick/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/derivatives/tuning_2d_model/stim_class/bayesian_posterior/initial/v1_e1-12_summary_b10_r0.001_g0_full_full_vary_all_models.csv')
bootstrap_models = pd.read_csv('/users/broderick/mnt/winawerlab/Projects/spatial_frequency_preferences/BIDS/derivatives/tuning_2d_model/stim_class/bayesian_posterior/bootstrap/v1_e1-12_full_b10_r0.001_g0_full_full_vary_all_models.csv')
In [6]:
indicator_order = models.indicator.unique()
In [7]:
parameters = [r'$\sigma$', r'$a$', r'$b$', r'$p_1$', r'$p_2$', r'$p_3$', r'$p_4$', r'$A_1$', r'$A_2$', r'$A_3$', r'$A_4$']
original_param_names = ['sigma', 'sf_ecc_slope', 'sf_ecc_intercept', 'abs_mode_cardinals', 'abs_mode_obliques', 'rel_mode_cardinals', 'rel_mode_obliques',
'abs_amplitude_cardinals', 'abs_amplitude_obliques', 'rel_amplitude_cardinals', 'rel_amplitude_obliques']
rename_params = dict((k, v) for k, v in zip(original_param_names, parameters))
models_plot = models.copy()#.query('fit_model_type==@fit_model_type')
models_plot = models_plot.set_index('model_parameter')
models_plot.loc['sigma', 'param_category'] = 'sigma'
models_plot.loc[['sf_ecc_slope', 'sf_ecc_intercept'], 'param_category'] = 'eccen'
models_plot.loc[['abs_mode_cardinals', 'abs_mode_obliques', 'rel_mode_cardinals', 'rel_mode_obliques',
'abs_amplitude_cardinals', 'abs_amplitude_obliques', 'rel_amplitude_cardinals', 'rel_amplitude_obliques'], 'param_category'] = 'orientation'
models_plot = models_plot.reset_index()
models_plot['model_parameter'] = models_plot.model_parameter.map(rename_params)
pal = dict(zip(models.indicator.unique(), sns.color_palette('husl', models.indicator.nunique())))
Model fit to the median responses (across bootstraps) for each subject/session, then this data shows the median across subject/sessions, with confidence intervals giving the 95% bootstrapped (across subject/sessions) confidence intervals
In [164]:
with sns.plotting_context('poster', font_scale=1):
fig, axes = plt.subplots(1, 3, figsize=(20, 10), gridspec_kw={'width_ratios': [.15, .3, .6]})
for i, ax in enumerate(axes):
cat = ['sigma', 'eccen', 'orientation'][i]
models_plot_order = [i for i in parameters if i in models_plot.query("param_category==@cat").model_parameter.unique()]
sns.pointplot('model_parameter', 'fit_value', 'model_parameter', data=models_plot.query("param_category==@cat"), estimator=np.median, s=15, ax=ax, order=models_plot_order)
if ax.legend_:
ax.legend_.remove()
#if i==2:
# ax.set(ylim=(-.25, .25))
ax.axhline(color='grey', linestyle='dashed')
ax.set(ylabel='Fit value', xlabel='Parameter')
#fig.savefig('param_vals.svg')
Same as above, except now we're showing each subject/session as a different point, also indicated by hue
In [124]:
# models_plot = models_plot.query("session in ['ses-03', 'ses-04']")
with sns.plotting_context('poster', font_scale=1):
fig, axes = plt.subplots(1, 3, figsize=(20, 10), gridspec_kw={'width_ratios': [.15, .3, .6]})
for i, ax in enumerate(axes):
cat = ['sigma', 'eccen', 'orientation'][i]
models_plot_order = [i for i in parameters if i in models_plot.query("param_category==@cat").model_parameter.unique()]
sns.stripplot('model_parameter', 'fit_value', 'indicator', data=models_plot.query("param_category==@cat"), s=15, ax=ax, order=models_plot_order, palette=pal)
if ax.legend_:
ax.legend_.remove()
#if i==2:
# ax.set(ylim=(-.25, .25))
ax.axhline(color='grey', linestyle='dashed')
ax.set(ylabel='Fit value', xlabel='Parameter')
#fig.savefig('param_vals.svg')
In [18]:
parameters = [r'$\sigma$', r'$a$', r'$b$', r'$p_1$', r'$p_2$', r'$p_3$', r'$p_4$', r'$A_1$', r'$A_2$', r'$A_3$', r'$A_4$']
original_param_names = ['sigma', 'sf_ecc_slope', 'sf_ecc_intercept', 'abs_mode_cardinals', 'abs_mode_obliques', 'rel_mode_cardinals', 'rel_mode_obliques',
'abs_amplitude_cardinals', 'abs_amplitude_obliques', 'rel_amplitude_cardinals', 'rel_amplitude_obliques']
rename_params = dict((k, v) for k, v in zip(original_param_names, parameters))
models_plot = bootstrap_models.copy()#.query('fit_model_type==@fit_model_type')
models_plot = models_plot.set_index('model_parameter')
models_plot.loc['sigma', 'param_category'] = 'sigma'
models_plot.loc[['sf_ecc_slope', 'sf_ecc_intercept'], 'param_category'] = 'eccen'
models_plot.loc[['abs_mode_cardinals', 'abs_mode_obliques', 'rel_mode_cardinals', 'rel_mode_obliques',
'abs_amplitude_cardinals', 'abs_amplitude_obliques', 'rel_amplitude_cardinals', 'rel_amplitude_obliques'], 'param_category'] = 'orientation'
models_plot = models_plot.reset_index()
models_plot['model_parameter'] = models_plot.model_parameter.map(rename_params)
pal = dict(zip(models.indicator.unique(), sns.color_palette('husl', models.indicator.nunique())))
#for k, v in zip(models.indicator.unique(), sns.color_palette('husl', models.indicator.nunique())):
# pal[k] = np.expand_dims(list(v)+[1], 0)
Same as above, except now we're showing the values from when we fit the model to each bootstrap separately (above was just fit to the median across bootstraps), with error bars showing the 95% confidence intervals across bootstraps.
You might look at this, compare it to the previous plot, and think -- these two don't look the same
In [126]:
with sns.plotting_context('poster', font_scale=1):
fig, axes = plt.subplots(1, 3, figsize=(20, 10), gridspec_kw={'width_ratios': [.15, .3, .6]})
for i, ax in enumerate(axes):
cat = ['sigma', 'eccen', 'orientation'][i]
models_plot_order = [i for i in parameters if i in models_plot.query("param_category==@cat").model_parameter.unique()]
for n, g in models_plot.query("param_category==@cat & session not in ['ses-pilot00', 'ses-pilot01'] & task!='task-sfpconstant'").groupby('indicator'):
sfp.plotting.scatter_ci_dist('model_parameter', 'fit_value', data=g, label=n, ax=ax, x_jitter=.2, color=pal[n], x_order=models_plot_order)
if ax.legend_:
ax.legend_.remove()
#if i==2:
# ax.set(ylim=(-.25, .25))
ax.axhline(color='grey', linestyle='dashed')
ax.set(ylabel='Fit value', xlabel='Parameter')
ax.set_xlim((-.5, models_plot.query("param_category==@cat").model_parameter.nunique()-.5))
#fig.savefig('param_vals.svg')
In [79]:
compare_df = models[['model_parameter', 'indicator', 'fit_value', 'session', 'subject', 'task']]
tmp = bootstrap_models[['model_parameter', 'indicator', 'fit_value', 'session', 'subject', 'task']].rename(columns={'fit_value': 'fit_value_bs'})
compare_df = pd.merge(tmp, compare_df, on=['model_parameter', 'indicator', 'session', 'subject', 'task'])
And you'd be correct, the median parameter value across bootstraps and the parameter value when fit to the median across bootstraps are not always the same. The following has each parameter on a separate subplot, with different subject/sessions on the x-axis, and there are two points for each: the one with error bars is the median parameter value across bootstraps (with 90% CI), while the one without is the value fit to the median across bootstraps. Hue shows the different tasks: regular in blue, constant stimuli in orange, and rescaled stimuli in green.
It looks like a handful of subjects (notably sub-wlsubj042) that have the most issues.
In [128]:
g = sns.FacetGrid(compare_df.query("session not in ['ses-pilot00', 'ses-pilot01']"), col='model_parameter', hue='task', col_wrap=4, sharey=False, aspect=2.5, height=3, col_order=original_param_names)
g.map_dataframe(sfp.plotting.scatter_ci_dist, 'indicator', 'fit_value_bs', ci_vals=[5, 95])
for ax in g.axes:
ax.set_xticklabels(ax.get_xticklabels(), rotation=90)
g.map_dataframe(plt.scatter, 'indicator', 'fit_value')
Out[128]:
In the above, it's easier to see which subject/sessions differ the most between these two, the following is an attempt to see for which parameters it differs the most.
In [130]:
id_vars = ['model_parameter', 'indicator', 'session', 'subject', 'task']
tmp = pd.melt(compare_df, id_vars,)
In [131]:
g = sns.FacetGrid(tmp.query("session not in ['ses-pilot00', 'ses-pilot01']"), row='subject', col='session', hue='model_parameter', sharey=False, aspect=2.5, height=3, )
g.map_dataframe(sfp.plotting.scatter_ci_dist, 'variable', 'value', ci_vals=[5, 95], x_jitter=True, join=True)
g.add_legend()
g.set(xlim=(-1, 2))
Out[131]:
In [132]:
g = sns.FacetGrid(tmp.query("session not in ['ses-pilot00', 'ses-pilot01']"), col='model_parameter', hue='indicator', col_wrap=4, sharey=False, aspect=2.5, height=3, col_order=original_param_names)
g.map_dataframe(sfp.plotting.scatter_ci_dist, 'variable', 'value', ci_vals=[5, 95], x_jitter=True, join=True)
g.add_legend()
g.set(xlim=(-1, 2))
Out[132]:
Looking at parameter values, it's hard to interpret what they mean. So we'll make some plots showing the variability in preferred period (as function of eccentricity, angle, and orientation) and amplitude (as function of angle and orientation), across and within subjects.
First, bootstrapping across subjects (using their median)
In [248]:
abs_period = sfp.analyze_model.create_feature_df(models.query("task=='task-sfprescaled'"), )
rel_period = sfp.analyze_model.create_feature_df(models.query("task=='task-sfprescaled'"), reference_frame='relative')
In [250]:
sfp.plotting.feature_df_plot(rel_period)
Out[250]:
In [252]:
sfp.plotting.feature_df_plot(rel_period, col='subject', pre_boot_gb_func='mean', col_wrap=8, )
Out[252]:
In [251]:
sfp.plotting.feature_df_plot(rel_period, col=None, pre_boot_gb_func=np.mean)
Out[251]:
In [249]:
sfp.plotting.feature_df_plot(abs_period)
Out[249]:
In [253]:
sfp.plotting.feature_df_plot(abs_period, col='subject', pre_boot_gb_func='mean', col_wrap=8, )
Out[253]:
In [246]:
sfp.plotting.feature_df_plot(abs_period, col=None, row=None, pre_boot_gb_func='mean', )
Out[246]:
In [320]:
rel_period = sfp.analyze_model.create_feature_df(bootstrap_models.query("task=='task-sfprescaled'"), reference_frame='relative', )
In [366]:
g = sfp.plotting.feature_df_plot(rel_period, col='subject', pre_boot_gb_func='mean', plot_func=sfp.plotting.scatter_ci_dist, draw_ctr_pts=False, ci_mode='fill', join=True)
In [384]:
period = sfp.analyze_model.create_feature_df(bootstrap_models.query("task=='task-sfprescaled'"), eccentricity=[2, 4, 10], retinotopic_angle=np.linspace(0, 2*np.pi, 49),)
In [385]:
period.head()
Out[385]:
In [435]:
g = sfp.plotting.feature_df_polar_plot(period, row='Eccentricity (deg)' , col='subject', plot_func=sfp.plotting.scatter_ci_dist, draw_ctr_pts=False, ci_mode='fill', join=True, r='Preferred period (deg)',
top=.9, hspace=.3, wspace=.1)
In [400]:
period.columns
Out[400]:
In [408]:
a = {'a': 1}
In [409]:
a.update({'a': 2})
In [433]:
g = sfp.plotting.feature_df_polar_plot(period[period['Eccentricity (deg)']==2], row='Eccentricity (deg)' , col='subject', plot_func=sfp.plotting.scatter_ci_dist, draw_ctr_pts=False, ci_mode='fill', join=True, r='Preferred period (deg)',
all_tick_labels=['r', 'theta'], top=.9)
In [390]:
abs_period_contour = sfp.analyze_model.create_feature_df(models, 'preferred_period_contour', )
rel_period_contour = sfp.analyze_model.create_feature_df(models, 'preferred_period_contour', reference_frame='relative')
In [5]:
sfp.plotting.feature_df_polar_plot(abs_period_contour)
Out[5]:
In [6]:
sfp.plotting.feature_df_polar_plot(abs_period_contour.query("`Stimulus type` in ['vertical', 'horizontal']"))
Out[6]:
In [7]:
sfp.plotting.feature_df_polar_plot(rel_period_contour)
Out[7]:
In [8]:
sfp.plotting.feature_df_polar_plot(rel_period_contour.query("`Stimulus type` in ['angular', 'radial']"))
Out[8]:
In [17]:
abs_ampl = sfp.analyze_model.create_feature_df(models, 'max_amplitude', )
rel_ampl = sfp.analyze_model.create_feature_df(models, 'max_amplitude', reference_frame='relative')
In [18]:
sfp.plotting.feature_df_polar_plot(abs_ampl, col=None, r='Max amplitude')
Out[18]:
In [19]:
sfp.plotting.feature_df_polar_plot(abs_ampl.query("`Stimulus type` in ['vertical', 'horizontal']"), col=None, r='Max amplitude')
Out[19]:
In [20]:
sfp.plotting.feature_df_polar_plot(rel_ampl, col=None, r='Max amplitude')
Out[20]:
In [21]:
sfp.plotting.feature_df_polar_plot(rel_ampl.query("`Stimulus type` in ['angular', 'radial']"), col=None, r='Max amplitude')
Out[21]:
In [96]:
b = sfp.analyze_model.bootstrap_features(feat, 100)
In [391]:
fit_model_type = 'full_donut_full_amps-vary'
angles = np.linspace(0, 2*np.pi, 49)
oris = np.linspace(0, np.pi, 4, endpoint=False)
eccens = np.linspace(0, 11, 11)
period_target = [.5, 1, 1.5]
ori_labels = ['radial', 'forward spiral', 'angular', 'reverse spiral']
abs_ori_labels = ['vertical', 'forward diagonal', 'horizontal', 'reverse diagonal']
angle_labels = ['0', r'$\frac{\pi}{4}$', r'$\frac{\pi}{2}$', r'$\frac{3\pi}{4}$']
pal = sns.color_palette('deep', len(ori_labels))
abs_pal = sns.color_palette('cubehelix', len(abs_ori_labels))
indicators = models.indicator.unique()
all_contours = np.empty((len(indicators), len(period_target), len(angles), len(oris)))
all_rel_contours = np.empty_like(all_contours)
all_periods = np.empty((len(indicators), len(angles), len(oris), len(eccens)))
all_rel_periods = np.empty_like(all_periods)
all_amps = np.empty((len(indicators), len(angles), len(oris)))
all_rel_amps = np.empty_like(all_amps)
for j, ind in enumerate(indicators):
single_model_df = models.query('indicator==@ind & fit_model_type==@fit_model_type')
model = sfp.model.LogGaussianDonut.init_from_df(single_model_df)
for k, p in enumerate(period_target):
all_contours[j, k, :, :] = model.preferred_period_contour(p, angles, oris).detach().numpy()
all_rel_contours[j, k, :, :] = model.preferred_period_contour(p, angles, rel_sf_angle=oris).detach().numpy()
all_amps[j, :, :] = model.max_amplitude(angles, oris).detach().numpy()
all_rel_amps[j, :, :] = model.max_amplitude(angles, rel_sf_angle=oris).detach().numpy()
for k, a in enumerate(oris):
all_rel_periods[j, :, k, :] = model.preferred_period(eccens, angles, rel_sf_angle=a).detach().numpy()
all_periods[j, :, k, :] = model.preferred_period(eccens, angles, a).detach().numpy()
In [392]:
n_bootstraps = 10000
bootstraps = np.random.randint(0, len(indicators), size=(n_bootstraps, len(indicators)))
contours = np.empty((n_bootstraps, len(period_target), len(angles), len(oris)))
rel_contours = np.empty_like(contours)
periods = np.empty((n_bootstraps, len(angles), len(oris), len(eccens)))
rel_periods = np.empty_like(periods)
amps = np.empty((n_bootstraps, len(angles), len(oris)))
rel_amps = np.empty_like(amps)
for i, b in enumerate(bootstraps):
periods[i, :, :, :] = np.mean(all_periods[b], 0)
rel_periods[i, :, :, :] = np.mean(all_rel_periods[b], 0)
contours[i, :, :, :] = np.mean(all_contours[b], 0)
rel_contours[i, :, :, :] = np.mean(all_rel_contours[b], 0)
amps[i, :, :] = np.mean(all_amps[b], 0)
rel_amps[i, :, :] = np.mean(all_rel_amps[b], 0)
In [208]:
periods_ci = np.percentile(periods, [16, 84], 0)
periods_mn = np.median(periods, 0)
rel_periods_ci = np.percentile(rel_periods, [16, 84], 0)
rel_periods_mn = np.median(rel_periods, 0)
angle_labels_tmp = angle_labels[::2]
fig, axes = plt.subplots(2, 2, figsize=(15, 10), sharey=True, sharex=True, subplot_kw={'ylim': (0, 2.25)})
for j_idx, j in enumerate(np.where(np.in1d(angles, [0, np.pi/2]))[0]):
for i in range(4):
axes[0, j_idx].plot(eccens, periods_mn[j,i,:], label=abs_ori_labels[i], c=abs_pal[i])
axes[0, j_idx].fill_between(eccens, *periods_ci[:,j,i,:], alpha=.2, color=abs_pal[i])
axes[1, j_idx].plot(eccens, rel_periods_mn[j,i,:], label=ori_labels[i], c=pal[i])
axes[1, j_idx].fill_between(eccens, *rel_periods_ci[:,j,i,:], alpha=.2, color=pal[i])
axes[0,j_idx].set_title('retinotopic angle %s' % angle_labels_tmp[j_idx])
axes[0,j_idx].legend()
axes[1,j_idx].legend()
axes[1, j_idx].set_xlabel("Eccentricity (deg)")
axes[0, 0].set_ylabel("Preferred Period (deg)")
axes[1, 0].set_ylabel("Preferred Period (deg)")
plt.tight_layout()
fig.savefig('all_pref_period.svg')
In [105]:
periods_ci = np.percentile(periods.mean(1), [16, 84], 0)
periods_mn = np.median(periods.mean(1), 0)
rel_periods_ci = np.percentile(rel_periods.mean(1), [16, 84], 0)
rel_periods_mn = np.median(rel_periods.mean(1), 0)
fig, axes = plt.subplots(1, 2, figsize=(15, 5), sharex=True, sharey=True, subplot_kw={'ylim': (0, 2.25)})
for i in range(4):
axes[0].plot(eccens, periods_mn[i,:], label=abs_ori_labels[i], c=abs_pal[i])
axes[0].fill_between(eccens, *periods_ci[:,i,:], alpha=.2, color=abs_pal[i])
axes[1].plot(eccens, rel_periods_mn[i,:], label=ori_labels[i], c=pal[i])
axes[1].fill_between(eccens, *rel_periods_ci[:,i,:], alpha=.2, color=pal[i])
axes[0].set_title('Averaged across retinotopic angles')
axes[0].legend()
axes[1].set_title('Averaged across retinotopic angles')
axes[1].legend()
axes[0].set_ylabel("Preferred Period (deg)")
axes[0].set_xlabel("Eccentricity (deg)")
axes[1].set_xlabel("Eccentricity (deg)")
plt.tight_layout()
fig.savefig('pref_period.svg')
In [359]:
periods_ci = np.percentile(periods, [16, 84], 0)
periods_mn = np.median(periods, 0)
rel_periods_ci = np.percentile(rel_periods, [16, 84], 0)
rel_periods_mn = np.median(rel_periods, 0)
angle_labels_tmp = angle_labels[::2]
with sns.plotting_context('poster'):
fig, axes = plt.subplots(1, 2, figsize=(14, 6), sharey=True, sharex=True, )#subplot_kw={'ylim': (0, 2.25)})
for j_idx, j in enumerate(np.where(np.in1d(angles, [0, np.pi/2]))[0]):
for i in range(4):
axes[j_idx].plot(eccens, rel_periods_mn[j,i,:], label=ori_labels[i], c=pal[i])
axes[j_idx].fill_between(eccens, *rel_periods_ci[:,j,i,:], alpha=.2, color=pal[i])
axes[j_idx].set_title('retinotopic angle %s' % angle_labels_tmp[j_idx])
axes[j_idx].set_xlabel("Eccentricity (deg)")
axes[j_idx].axhline(color='k', linestyle='--')
axes[j_idx].axvline(color='k', linestyle='--')
axes[0].set_ylabel("Preferred Period (deg)")
plt.tight_layout()
fig.savefig('all_pref_period.svg')
In [365]:
periods_ci = np.percentile(periods.mean(1), [16, 84], 0)
periods_mn = np.median(periods.mean(1), 0)
rel_periods_ci = np.percentile(rel_periods.mean(1), [16, 84], 0)
rel_periods_mn = np.median(rel_periods.mean(1), 0)
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7, 6), sharex=True, sharey=True,)#subplot_kw={'ylim': (0, 2.5)})
for i in range(4):
ax.plot(eccens, rel_periods_mn[i,:], label=ori_labels[i], c=pal[i])
ax.fill_between(eccens, *rel_periods_ci[:,i,:], alpha=.2, color=pal[i])
ax.set_title('Averaged across retinotopic angles')
ax.legend()
ax.axhline(color='k', linestyle='--')
ax.axvline(color='k', linestyle='--')
ax.set_ylabel("Preferred Period (deg)")
ax.set_xlabel("Eccentricity (deg)")
plt.tight_layout()
fig.savefig('pref_period.svg')
In [366]:
contours_ci = np.percentile(contours, [16, 84], 0)
contours_mn = contours.mean(0)
# with sns.plotting_context('poster'):
# fig, axes = plt.subplots(1, 3, figsize=(24, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
# for j, ax in enumerate(axes.flatten()):
# for i in [0, 2, ]:
# ax.plot(angles, contours_mn[j,:,i], label=abs_ori_labels[i], color=abs_pal[i])
# ax.fill_between(angles, contours_ci[0,j,:,i], contours_ci[1,j,:,i], alpha=.2, color=abs_pal[i])
# ax.set_title('Preferred period %s' % period_target[j])
# ax.legend()
# fig.savefig('abs_contour.svg')
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
for i in [0, 2, ]:
ax.plot(angles, contours_mn[1,:,i], label=abs_ori_labels[i], color=abs_pal[i])
ax.fill_between(angles, *contours_ci[:,1,:,i], alpha=.2, color=abs_pal[i])
ax.set_title('Preferred period=%s contour' % period_target[1])
ax.legend()
fig.savefig('abs_contour.svg')
In [130]:
contours_ci = np.percentile(rel_contours, [16, 84], 0)
contours_mn = rel_contours.mean(0)
# with sns.plotting_context('poster'):
# fig, axes = plt.subplots(1, 3, figsize=(24, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
# for j, ax in enumerate(axes.flatten()):
# for i in [0, 2]:
# ax.plot(angles, contours_mn[j,:,i], label=ori_labels[i], color=pal[i])
# ax.fill_between(angles, contours_ci[0,j,:,i], contours_ci[1,j,:,i], alpha=.2, color=pal[i])
# ax.set_title('Preferred period %s' % period_target[j])
# ax.legend()
# fig.savefig('rel_contour.svg')
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
for i in [0, 2]:
ax.plot(angles, contours_mn[1,:,i], label=ori_labels[i], color=pal[i])
ax.fill_between(angles, *contours_ci[:,1,:,i], alpha=.2, color=pal[i])
ax.set_title('Preferred period %s' % period_target[1])
ax.legend()
fig.savefig('rel_contour.svg')
In [219]:
contours_ci = np.percentile(rel_periods, [16, 84], 0)
contours_mn = rel_periods.mean(0)
# with sns.plotting_context('poster'):
# fig, axes = plt.subplots(1, 3, figsize=(24, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
# for j, ax in enumerate(axes.flatten()):
# for i in [0, 2]:
# ax.plot(angles, contours_mn[j,:,i], label=ori_labels[i], color=pal[i])
# ax.fill_between(angles, contours_ci[0,j,:,i], contours_ci[1,j,:,i], alpha=.2, color=pal[i])
# ax.set_title('Preferred period %s' % period_target[j])
# ax.legend()
# fig.savefig('rel_contour.svg')
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
for i in [0, 2]:
ax.plot(angles, contours_mn[:,i,4], label=ori_labels[i], color=pal[i])
ax.fill_between(angles, *contours_ci[:,:,i,4], alpha=.2, color=pal[i])
ax.set_title('Preferred period at eccentricity %.01f' % eccens[4])
ax.legend()
fig.savefig('rel_pref_period.svg')
In [225]:
contours_ci = np.percentile(periods, [16, 84], 0)
contours_mn = periods.mean(0)
# with sns.plotting_context('poster'):
# fig, axes = plt.subplots(1, 3, figsize=(24, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
# for j, ax in enumerate(axes.flatten()):
# for i in [0, 2]:
# ax.plot(angles, contours_mn[j,:,i], label=ori_labels[i], color=pal[i])
# ax.fill_between(angles, contours_ci[0,j,:,i], contours_ci[1,j,:,i], alpha=.2, color=pal[i])
# ax.set_title('Preferred period %s' % period_target[j])
# ax.legend()
# fig.savefig('rel_contour.svg')
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': 'polar'}, sharex=True, sharey=True)
for i in [0, 2]:
ax.plot(angles, contours_mn[:,i,4], label=abs_ori_labels[i], color=abs_pal[i])
ax.fill_between(angles, *contours_ci[:,:,i,4], alpha=.2, color=abs_pal[i])
ax.set_title('Preferred period at eccentricity %.01f' % eccens[4])
ax.legend()
fig.savefig('abs_pref_period.svg')
In [449]:
amps_ci = np.percentile(amps, [16, 84], 0)
amps_mn = amps.mean(0)
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7, 7), subplot_kw={'projection': 'polar'})
for i in [0, 2,1,3]:
ax.plot(angles, amps_mn[:,i], label=abs_ori_labels[i], color=abs_pal[i])
ax.fill_between(angles, *amps_ci[:,:,i], alpha=.2, color=abs_pal[i])
ax.set_title('Max amplitude')
ax.legend()
ax.set_ylim((0, 1.2))
ax.set(yticklabels=['', .4, '', '', 1])
fig.savefig('abs_amps.svg')
In [450]:
amps_ci = np.percentile(rel_amps, [16, 84], 0)
amps_mn = rel_amps.mean(0)
with sns.plotting_context('poster'):
fig, ax = plt.subplots(1, 1, figsize=(7,7), subplot_kw={'projection': 'polar'})
for i in [0, 2,1,3]:
ax.plot(angles, amps_mn[:,i], label=ori_labels[i], color=pal[i])
ax.fill_between(angles, *amps_ci[:,:,i], alpha=.2, color=pal[i])
ax.set_title('Max amplitude')
ax.legend()
ax.set_ylim((0, 1.2))
ax.set(yticklabels=['', .4, '', '', 1])
fig.savefig('rel_amps.svg')
In [245]:
with sns.plotting_context('poster'):
fig, axes = plt.subplots(1, 7, figsize=(49, 7), subplot_kw={'projection': 'polar'})
for j, ax in enumerate(axes.flatten()):
for i in [0, 2]:
ax.plot(angles, all_amps[j,:,i], label=abs_ori_labels[i], color=abs_pal[i])
ax.set_title('Max amplitude')
ax.legend()
ax.set_ylim((0, 1.2))
ax.set(yticklabels=['', .4, '', '', 1])
In [394]:
with sns.plotting_context('poster'):
fig, axes = plt.subplots(1, 7, figsize=(49, 7), subplot_kw={'projection': 'polar'})
for j, ax in enumerate(axes.flatten()):
for i in [0, 2]:
ax.plot(angles, all_rel_amps[j,:,i], label=ori_labels[i], color=pal[i])
ax.set_title('Max amplitude')
ax.legend()
ax.set_ylim((0, 1.2))
ax.set(yticklabels=['', .4, '', '', 1])
In [246]:
with sns.plotting_context('poster'):
fig, axes = plt.subplots(1, 7, figsize=(49, 7), subplot_kw={'projection': 'polar'})
for j, ax in enumerate(axes.flatten()):
for i in [0, 2]:
ax.plot(angles, all_contours[j,1,:,i], label=abs_ori_labels[i], color=abs_pal[i])
In [393]:
with sns.plotting_context('poster'):
fig, axes = plt.subplots(1, 7, figsize=(49, 7), subplot_kw={'projection': 'polar'})
for j, ax in enumerate(axes.flatten()):
for i in [0, 2]:
ax.plot(angles, all_rel_contours[j,1,:,i], label=ori_labels[i], color=pal[i])
In [383]:
original_param_names = ['sigma', 'sf_ecc_slope', 'sf_ecc_intercept', 'abs_mode_cardinals', 'abs_mode_obliques', 'rel_mode_cardinals', 'rel_mode_obliques',
'abs_amplitude_cardinals', 'abs_amplitude_obliques', 'rel_amplitude_cardinals', 'rel_amplitude_obliques']
rename_params = dict((k, v) for k, v in zip(original_param_names, parameters))
tmp = models.query('fit_model_type==@fit_model_type')
tmp = tmp.set_index('model_parameter')
tmp.loc['sigma', 'param_category'] = 'sigma'
tmp.loc[['sf_ecc_slope', 'sf_ecc_intercept'], 'param_category'] = 'eccen'
tmp.loc[['abs_mode_cardinals', 'abs_mode_obliques', 'rel_mode_cardinals', 'rel_mode_obliques',
'abs_amplitude_cardinals', 'abs_amplitude_obliques', 'rel_amplitude_cardinals', 'rel_amplitude_obliques'], 'param_category'] = 'orientation'
tmp = tmp.reset_index()
tmp['model_parameter'] = tmp.model_parameter.map(rename_params)
In [448]:
with sns.plotting_context('poster', font_scale=1):
fig, axes = plt.subplots(1, 3, figsize=(20, 10), gridspec_kw={'width_ratios': [.15, .3, .6]})
for i, ax in enumerate(axes):
cat = ['sigma', 'eccen', 'orientation'][i]
tmp_order = [i for i in parameters if i in tmp.query("param_category==@cat").model_parameter.unique()]
sns.stripplot('model_parameter', 'fit_value', 'indicator', data=tmp.query("param_category==@cat"), s=15, ax=ax, order=tmp_order)
if ax.legend_:
ax.legend_.remove()
if i==2:
ax.set(ylim=(-.25, .25))
ax.axhline(color='grey', linestyle='dashed')
ax.set(ylabel='Fit value', xlabel='Parameter')
#fig.savefig('param_vals.svg')