Port some of the schrodinger's plots

[msultan]

Someone on mdtraj linked to pdf http://content.schrodinger.com/Resources/SID/3p6h_desmond_npt_10ns.pdf . There are some nice plots in there that we might want to consider bringing in. I like the dihedral ones in particular.

[cxhernandez]

Whaddya think?

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import seaborn.apionly as sns
import numpy as np

N = 25

gs = gridspec.GridSpec(2, 6)
ax1 = plt.subplot(gs[:1, :2], polar=True)
ax2 = plt.subplot(gs[:1, 2:])

sns.distplot(x, bins=N, ax=ax2)
radii, theta = np.histogram(x, bins=N, normed=True)
ax1.set_yticklabels([])
ax2.set_yticks([])
_=plt.xlim([0, 360])
width = (2*np.pi) / N

bars = ax1.bar(np.pi * theta[1:] / 180., radii, width=width,)

plt.tight_layout()

[msultan]

that looks sweet!

[msultan]

One small nitpicky thing, maybe have the limits be -pi to pi?

[jeiros]

Hi @cxhernandez That's a great plot! I've implemented it in a function that allows for wrapping between -180 to 180, or constraint from 0 to 360.

However, I'm running into some odd results that I can't quite figure out. Once I can clear the confusion I think it could be a great addition to the package.

This is how I've done it. I first define a couple of functions that wrap angles back to the desired range, even if the angles are several phases larger than the initial range.

def wrapAngle(x):
    """Wraps an angle between -180 and 180 degrees"""
    x = (x+180) % 360
    if x < 0:
        x += 360
    return x - 180

def constrainAngle(x):
    """Constrains an angle between 0 and 360 degrees"""
    x = x % 360
    if x < 0:
        x += 360
    return x

def cleanup_top_right_axes(ax):
    """
    Removes the top and right axis lines of a matplotlib axis
    """
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    return ax

Then this is your function with some additions:

def plot_angle(data, N=50, title=None, ax1=None, ax2=None, color=None, wrap=True):
    if ax1 is None or ax2 is None:
        gs = gridspec.GridSpec(2, 6)
        ax1 = plt.subplot(gs[:1, :2], polar=True)
        ax2 = plt.subplot(gs[:1, 2:])

    if wrap:
        vf = np.vectorize(wrapAngle)        
    else:
        vf = np.vectorize(constrainAngle)
    x = vf(data)

    sns.distplot(x, bins=N, ax=ax2, color=color, kde=True)
    radii, theta = np.histogram(x, bins=N, normed=True)
    ax1.set_yticklabels([])

    if wrap:
        ax2.set_xlim(-180, 180)
        ax1.set_xticklabels(['{}°'.format(x) for x in [0, 45, 90, 135, 180, -135, -90, -45]])
        ax2.set_xticklabels(['{}°'.format(x) for x in [-180, -135, -90, -45, 0, 45, 90, 135, 180]])

    else:
        ax2.set_xlim(0, 360)
        ax2.set_xticklabels(['{}°'.format(x) for x in [0, 45, 90, 135, 180, 225, 270, 315, 360]])

    ax2.set_yticks([])
    ax2.set(xlabel='Angle', ylabel='Density')
    cleanup_top_right_axes(ax2)

    width = (2 * np.pi) / N

    bars = ax1.bar(np.deg2rad(theta[1:]), radii, width=width, color=color, alpha=.5)

    if title is not None:
        plt.suptitle(title)

    plt.tight_layout()

    f = plt.gcf()
    return f, (ax1, ax2)

So let's test it:

# Generate a bimodal angle
a = np.random.normal(135, size=1000)
b = np.random.normal(270, size=1000)
c = np.concatenate([a, b])

# plot with no wrap, from 0 to 360
plot_angle(c, wrap=False)

# wrap from -180 to 180
plot_angle(c, wrap=True)

There's a couple of weird things going on. The polar plots recover the actual values for the angles, however the distribution plots from seaborn seem off?

In the case of constraining from 0 to 360, since the values I'm generating are already in that bound, I've made sure that the numbers match:

np.allclose(c, np.vectorize(constrainAngle)(c))
True

And then if I just call sns.distplot on it's own, it correctly places the distribution horizontally, contrarily to the one in the function:

sns.distplot(np.vectorize(constrainAngle)(c))
sns.distplot(c)

Also a minor thing that I can't figure out is even when I specify the xlimits with ax2.set_xlim this seems to be ignored.

Sorry for the long post, as you can see I'm a bit lost as to what it going on and could use some help 🤣

[cxhernandez]

Nice plots! I think the issue here is just that you need to set the xticks to make sure that there are as many as the xticklabels. If you comment out the set_xticklabels lines, I think it'll work just fine.

[jeiros]

Oh my god I was going insane trying to figure it out! Yes, it was the tick labels being in the wrong place so it was confusing. Thanks!!

I'll clean the function a bit and then put together a PR to add it if that is OK!

[cxhernandez]

done in #119