/usr/share/psi/python/qcdb/mpl.py

"""Module with matplotlib plotting routines. These are not hooked up to
any particular qcdb data structures but can be called with basic
arguments.

"""
import os
#import matplotlib
#matplotlib.use('Agg')


def expand_saveas(saveas, def_filename, def_path=os.path.abspath(os.curdir), def_prefix='', relpath=False):
    """Analyzes string *saveas* to see if it contains information on
    path to save file, name to save file, both or neither (*saveas*
    ends in '/' to indicate directory only) (able to expand '.'). A full
    absolute filename is returned, lacking only file extension. Based on
    analysis of missing parts of *saveas*, path information from *def_path*
    and/or filename information from *def_prefix* + *def_filename* is
    inserted. *def_prefix* is intended to be something like 'mplthread_'
    to identify the type of figure.

    """
    defname = def_prefix + def_filename.replace(' ', '_')
    if saveas is None:
        pth = def_path
        fil = defname
    else:
        pth, fil = os.path.split(saveas)
        pth = pth if pth != '' else def_path
        fil = fil if fil != '' else defname

    abspathfile = os.path.join(os.path.abspath(pth), fil)
    if relpath:
        return abspathfile.split(os.path.commonprefix([abspathfile, os.getcwd()]) + os.sep)[1]
    else:
        return abspathfile


def segment_color(argcolor, saptcolor):
    """Find appropriate color expression between overall color directive
    *argcolor* and particular color availibility *rxncolor*.

    """
    import matplotlib

    # validate any sapt color
    if saptcolor is not None:
        if saptcolor < 0.0 or saptcolor > 1.0:
            saptcolor = None

    if argcolor is None:
        # no color argument, so take from rxn
        if rxncolor is None:
            clr = 'grey'
        elif saptcolor is not None:
            clr = matplotlib.cm.jet(saptcolor)
        else:
            clr = rxncolor
    elif argcolor == 'sapt':
        # sapt color from rxn if available
        if saptcolor is not None:
            clr = matplotlib.cm.jet(saptcolor)
        else:
            clr = 'grey'
    elif argcolor == 'rgb':
        # HB/MX/DD sapt color from rxn if available
        if saptcolor is not None:
            if saptcolor < 0.333:
                clr = 'blue'
            elif saptcolor < 0.667:
                clr = 'green'
            else:
                clr = 'red'
        else:
            clr = 'grey'
    else:
        # color argument is name of mpl color
        clr = argcolor

    return clr


def bars(data, title='', saveas=None, relpath=False, graphicsformat=['pdf']):
    """Generates a 'gray-bars' diagram between model chemistries with error
    statistics in list *data*, which is supplied as part of the dictionary
    for each participating bar/modelchem, along with *mc* keys in argument
    *data*. The plot is labeled with *title* and each bar with *mc* key and
    plotted at a fixed scale to facilitate comparison across projects.

    """
    import hashlib
    import matplotlib.pyplot as plt

    # initialize plot, fix dimensions for consistent Illustrator import
    fig, ax = plt.subplots(figsize=(12, 7))
    plt.ylim([0, 4.86])
    plt.xlim([0, 6])
    plt.xticks([])

    # label plot and tiers
    ax.text(0.4, 4.6, title,
        verticalalignment='bottom', horizontalalignment='left',
        family='Times New Roman', weight='bold', fontsize=12)

    widths = [0.15, 0.02, 0.02, 0.02]  # TT, HB, MX, DD
    xval = 0.1  # starting posn along x-axis

    # plot bar sets
    for bar in data:
        if bar is not None:
            lefts = [xval, xval + 0.025, xval + 0.065, xval + 0.105]

            rect = ax.bar(lefts, bar['data'], widths, linewidth=0)
            rect[0].set_color('grey')
            rect[1].set_color('red')
            rect[2].set_color('green')
            rect[3].set_color('blue')

            ax.text(xval + .08, 4.3, bar['mc'],
                verticalalignment='center', horizontalalignment='right', rotation='vertical',
                family='Times New Roman', fontsize=8)
        xval += 0.20

    # save and show
    pltuid = title + '_' + hashlib.sha1(title + repr([bar['mc'] for bar in data if bar is not None])).hexdigest()
    pltfile = expand_saveas(saveas, pltuid, def_prefix='bar_', relpath=relpath)
    files_saved = {}
    for ext in graphicsformat:
        savefile = pltfile + '.' + ext.lower()
        plt.savefig(savefile, transparent=True, format=ext, bbox_inches='tight')
        files_saved[ext.lower()] = savefile
    plt.show()
    return files_saved


def flat(data, color=None, title='', xlimit=4.0, mae=None, mape=None, view=True,
    saveas=None, relpath=False, graphicsformat=['pdf']):
    """Generates a slat diagram between model chemistries with errors in
    single-item list *data*, which is supplied as part of the dictionary
    for each participating reaction, along with *dbse* and *rxn* keys in
    argument *data*. Limits of plot are *xlimit* from the zero-line. If
    *color* is None, slats are black, if 'sapt', colors are taken from
    sapt_colors module. Summary statistic *mae* is plotted on the
    overbound side and relative statistic *mape* on the underbound side.
    Saves a file with name *title* and plots to screen if *view*.

    """
    import matplotlib.pyplot as plt

    Nweft = 1
    positions = range(-1, -1 * Nweft - 1, -1)

    # initialize plot
    fig, ax = plt.subplots(figsize=(12, 0.33))
    plt.xlim([-xlimit, xlimit])
    plt.ylim([-1 * Nweft - 1, 0])
    plt.yticks([])
    plt.xticks([])
    fig.patch.set_visible(False)
    ax.patch.set_visible(False)
    ax.axis('off')

    plt.axvline(-1.0, color='grey', linewidth=4)
    plt.axvline(-0.3, color='grey', linewidth=4)
    plt.axvline(0.0, color='grey', linewidth=4)
    plt.axvline(0.3, color='grey', linewidth=4)
    plt.axvline(1.0, color='grey', linewidth=4)

    # plot reaction errors and threads
    for rxn in data:
        xvals = rxn['data']
        clr = segment_color(color, rxn['color'] if 'color' in rxn else None)

        ax.plot(xvals, positions, '|', color=clr, markersize=13.0)

    # plot trimmings
    if mae is not None:
        plt.axvline(-1 * mae, color='black', linewidth=12)
    if mape is not None:  # equivalent to MAE for a 10 kcal/mol interaction energy
        ax.plot(0.025 * mape, positions, 'o', color='black', markersize=15.0)

    # save and show
    pltuid = title  # simple (not really unique) filename for LaTeX integration
    pltfile = expand_saveas(saveas, pltuid, def_prefix='flat_', relpath=relpath)
    files_saved = {}
    for ext in graphicsformat:
        savefile = pltfile + '.' + ext.lower()
        plt.savefig(savefile, transparent=True, format=ext)  # , bbox_inches='tight')
        files_saved[ext.lower()] = savefile
    plt.show()
    if not view:
        plt.close()
    return files_saved


#def mpl_distslat_multiplot_files(pltfile, dbid, dbname, xmin, xmax, mcdats, labels, titles):
#    """Saves a plot with basename *pltfile* with a slat representation
#    of the modelchems errors in *mcdat*. Plot is in PNG, PDF, & EPS
#    and suitable for download, no mouseover properties. Both labeled
#    and labelless (for pub) figures are constructed.
#
#    """
#    import matplotlib as mpl
#    from matplotlib.axes import Subplot
#    import sapt_colors
#    from matplotlib.figure import Figure
#
#    nplots = len(mcdats)
#    fht = nplots * 0.8
#    fig, axt = plt.subplots(figsize=(12.0, fht))
#    plt.subplots_adjust(left=0.01, right=0.99, hspace=0.3)
#
#    axt.set_xticks([])
#    axt.set_yticks([])
#    plt.axis('off')
#
#    for item in range(nplots):
#        mcdat = mcdats[item]
#        label = labels[item]
#        title = titles[item]
#
#        erdat = np.array(mcdat)
#        yvals = np.ones(len(mcdat))
#        y = np.array([sapt_colors.sapt_colors[dbname][i] for i in label])
#
#        ax = Subplot(fig, nplots, 1, item + 1)
#        fig.add_subplot(ax)
#        sc = ax.scatter(erdat, yvals, c=y, s=3000, marker="|", cmap=mpl.cm.jet, vmin=0, vmax=1)
#
#        ax.set_yticks([])
#        ax.set_xticks([])
#        ax.set_frame_on(False)
#        ax.set_xlim([xmin, xmax])
#
#    # Write files with only slats
#    plt.savefig('scratch/' + pltfile + '_plain' + '.png', transparent=True, format='PNG')
#    plt.savefig('scratch/' + pltfile + '_plain' + '.pdf', transparent=True, format='PDF')
#    plt.savefig('scratch/' + pltfile + '_plain' + '.eps', transparent=True, format='EPS')
#
#    # Rewrite files with guides and labels
#    for item in range(nplots):
#        ax_again = fig.add_subplot(nplots, 1, item + 1)
#        ax_again.set_title(titles[item], fontsize=8)
#        ax_again.text(xmin + 0.3, 1.0, stats(np.array(mcdats[item])), fontsize=7, family='monospace', verticalalignment='center')
#        ax_again.plot([0, 0], [0.9, 1.1], color='#cccc00', lw=2)
#        ax_again.set_frame_on(False)
#        ax_again.set_yticks([])
#        ax_again.set_xticks([-12.0, -8.0, -4.0, -2.0, -1.0, 0.0, 1.0, 2.0, 4.0, 8.0, 12.0])
#        ax_again.tick_params(axis='both', which='major', labelbottom='off', bottom='off')
#    ax_again.set_xticks([-12.0, -8.0, -4.0, -2.0, -1.0, 0.0, 1.0, 2.0, 4.0, 8.0, 12.0])
#    ax_again.tick_params(axis='both', which='major', labelbottom='on', bottom='off')
#
#    plt.savefig('scratch/' + pltfile + '_trimd' + '.png', transparent=True, format='PNG')
#    plt.savefig('scratch/' + pltfile + '_trimd' + '.pdf', transparent=True, format='PDF')
#    plt.savefig('scratch/' + pltfile + '_trimd' + '.eps', transparent=True, format='EPS')


def disthist(data, title='', xtitle='', xmin=None, xmax=None,
    me=None, stde=None, saveas=None, relpath=False, graphicsformat=['pdf']):
    """Saves a plot with name *saveas* with a histogram representation
    of the reaction errors in *data*. Also plots a gaussian distribution
    with mean *me* and standard deviation *stde*. Plot has x-range
    *xmin* to *xmax*, x-axis label *xtitle* and overall title *title*.

    """
    import hashlib
    import numpy as np
    import matplotlib.pyplot as plt

    def gaussianpdf(u, v, x):
        """*u* is mean, *v* is variance, *x* is value, returns probability"""
        return 1.0 / np.sqrt(2.0 * np.pi * v) * np.exp(-pow(x - u, 2) / 2.0 / v)

    me = me if me is not None else np.mean(data)
    stde = stde if stde is not None else np.std(data, ddof=1)
    xmin = xmin if xmin is not None else me - 4.0 * stde
    xmax = xmax if xmax is not None else me + 4.0 * stde

    dx = (xmax - xmin) / 40.
    nx = int(round((xmax - xmin) / dx)) + 1
    pdfx = []
    pdfy = []
    for i in xrange(nx):
        ix = xmin + i * dx
        pdfx.append(ix)
        pdfy.append(gaussianpdf(me, pow(stde, 2), ix))

    fig, ax = plt.subplots(figsize=(8, 6))
    plt.axvline(0.0, color='#cccc00')
    ax1 = fig.add_subplot(111)
    ax1.set_xlim(xmin, xmax)
    ax1.hist(data, bins=30, range=(xmin, xmax), color='#224477', alpha=0.7)
    ax1.set_xlabel(xtitle)
    ax1.set_ylabel('Count')

    ax2 = ax1.twinx()
    ax2.set_xlim(xmin, xmax)
    ax2.fill(pdfx, pdfy, color='k', alpha=0.2)
    ax2.set_ylabel('Probability Density')

    plt.title(title)

    # save and show
    pltuid = title + '_' + hashlib.sha1(title + str(me) + str(stde) + str(xmin) + str(xmax)).hexdigest()
    pltfile = expand_saveas(saveas, pltuid, def_prefix='disthist_', relpath=relpath)
    files_saved = {}
    for ext in graphicsformat:
        savefile = pltfile + '.' + ext.lower()
        plt.savefig(savefile, transparent=True, format=ext, bbox_inches='tight')
        files_saved[ext.lower()] = savefile
    plt.show()
    return files_saved


#def thread(data, labels, color=None, title='', xlimit=4.0, mae=None, mape=None):
#    """Generates a tiered slat diagram between model chemistries with
#    errors (or simply values) in list *data*, which is supplied as part of the
#    dictionary for each participating reaction, along with *dbse* and *rxn* keys
#    in argument *data*. The plot is labeled with *title* and each tier with
#    an element of *labels* and plotted at *xlimit* from the zero-line. If
#    *color* is None, slats are black, if 'sapt', colors are taken from *color*
#    key in *data* [0, 1]. Summary statistics *mae* are plotted on the
#    overbound side and relative statistics *mape* on the underbound side.
#
#    """
#    from random import random
#    import matplotlib.pyplot as plt
#
#    # initialize tiers/wefts
#    Nweft = len(labels)
#    lenS = 0.2
#    gapT = 0.04
#    positions = range(-1, -1 * Nweft - 1, -1)
#    posnS = []
#    for weft in range(Nweft):
#        posnS.extend([positions[weft] + lenS, positions[weft] - lenS, None])
#    posnT = []
#    for weft in range(Nweft - 1):
#        posnT.extend([positions[weft] - lenS - gapT, positions[weft + 1] + lenS + gapT, None])
#
#    # initialize plot
#    fht = Nweft * 0.8
#    fig, ax = plt.subplots(figsize=(12, fht))
#    plt.subplots_adjust(left=0.01, right=0.99, hspace=0.3)
#    plt.xlim([-xlimit, xlimit])
#    plt.ylim([-1 * Nweft - 1, 0])
#    plt.yticks([])
#
#    # label plot and tiers
#    ax.text(-0.9 * xlimit, -0.25, title,
#        verticalalignment='bottom', horizontalalignment='left',
#        family='Times New Roman', weight='bold', fontsize=12)
#    for weft in labels:
#        ax.text(-0.9 * xlimit, -(1.2 + labels.index(weft)), weft,
#            verticalalignment='bottom', horizontalalignment='left',
#            family='Times New Roman', weight='bold', fontsize=18)
#
#    # plot reaction errors and threads
#    for rxn in data:
#
#        # preparation
#        xvals = rxn['data']
#        clr = segment_color(color, rxn['color'] if 'color' in rxn else None)
#        slat = []
#        for weft in range(Nweft):
#            slat.extend([xvals[weft], xvals[weft], None])
#        thread = []
#        for weft in range(Nweft - 1):
#            thread.extend([xvals[weft], xvals[weft + 1], None])
#
#        # plotting
#        ax.plot(slat, posnS, color=clr, linewidth=1.0, solid_capstyle='round')
#        ax.plot(thread, posnT, color=clr, linewidth=0.5, solid_capstyle='round',
#            alpha=0.3)
#
#        # labeling
#        try:
#            toplblposn = next(item for item in xvals if item is not None)
#            botlblposn = next(item for item in reversed(xvals) if item is not None)
#        except StopIteration:
#            pass
#        else:
#            ax.text(toplblposn, -0.75 + 0.6 * random(), rxn['sys'],
#                verticalalignment='bottom', horizontalalignment='center',
#                family='Times New Roman', fontsize=8)
#            ax.text(botlblposn, -1 * Nweft - 0.75 + 0.6 * random(), rxn['sys'],
#                verticalalignment='bottom', horizontalalignment='center',
#                family='Times New Roman', fontsize=8)
#
#    # plot trimmings
#    if mae is not None:
#        ax.plot([-x for x in mae], positions, 's', color='black')
#    if mape is not None:  # equivalent to MAE for a 10 kcal/mol IE
#        ax.plot([0.025 * x for x in mape], positions, 'o', color='black')
#
#    plt.axvline(0, color='black')
#    plt.show()


def threads(data, labels, color=None, title='', xlimit=4.0, mae=None, mape=None,
    mousetext=None, mouselink=None, mouseimag=None, mousetitle=None,
    saveas=None, relpath=False, graphicsformat=['pdf']):
    """Generates a tiered slat diagram between model chemistries with
    errors (or simply values) in list *data*, which is supplied as part of the
    dictionary for each participating reaction, along with *dbse* and *rxn* keys
    in argument *data*. The plot is labeled with *title* and each tier with
    an element of *labels* and plotted at *xlimit* from the zero-line. If
    *color* is None, slats are black, if 'sapt', colors are taken from *color*
    key in *data* [0, 1]. Summary statistics *mae* are plotted on the
    overbound side and relative statistics *mape* on the underbound side.
    HTML code for mouseover if mousetext or mouselink or mouseimag specified
    based on recipe of Andrew Dalke from
    http://www.dalkescientific.com/writings/diary/archive/2005/04/24/interactive_html.html

    """
    import random
    import hashlib
    import matplotlib.pyplot as plt

    # initialize tiers/wefts
    Nweft = len(labels)
    lenS = 0.2
    gapT = 0.04
    positions = range(-1, -1 * Nweft - 1, -1)
    posnS = []
    for weft in range(Nweft):
        posnS.extend([positions[weft] + lenS, positions[weft] - lenS, None])
    posnT = []
    for weft in range(Nweft - 1):
        posnT.extend([positions[weft] - lenS - gapT, positions[weft + 1] + lenS + gapT, None])
    posnM = []

    # initialize plot
    fht = Nweft * 0.8
    fig, ax = plt.subplots(figsize=(12, fht))
    plt.subplots_adjust(left=0.01, right=0.99, hspace=0.3)
    plt.xlim([-xlimit, xlimit])
    plt.ylim([-1 * Nweft - 1, 0])
    plt.yticks([])

    # label plot and tiers
    ax.text(-0.9 * xlimit, -0.25, title,
        verticalalignment='bottom', horizontalalignment='left',
        family='Times New Roman', weight='bold', fontsize=12)
    for weft in labels:
        ax.text(-0.9 * xlimit, -(1.2 + labels.index(weft)), weft,
            verticalalignment='bottom', horizontalalignment='left',
            family='Times New Roman', weight='bold', fontsize=18)

    # plot reaction errors and threads
    for rxn in data:

        # preparation
        xvals = rxn['data']
        clr = segment_color(color, rxn['color'] if 'color' in rxn else None)
        slat = []
        for weft in range(Nweft):
            slat.extend([xvals[weft], xvals[weft], None])
        thread = []
        for weft in range(Nweft - 1):
            thread.extend([xvals[weft], xvals[weft + 1], None])

        # plotting
        ax.plot(slat, posnS, color=clr, linewidth=1.0, solid_capstyle='round')
        ax.plot(thread, posnT, color=clr, linewidth=0.5, solid_capstyle='round', alpha=0.3)

        # converting into screen coordinates for image map
        xyscreen = ax.transData.transform(zip(xvals, positions))
        xscreen, yscreen = zip(*xyscreen)
        posnM.extend(zip([rxn['db']] * Nweft, [rxn['sys']] * Nweft,
            xvals, xscreen, yscreen))

        # labeling
        if not(mousetext or mouselink or mouseimag):
            try:
                toplblposn = next(item for item in xvals if item is not None)
                botlblposn = next(item for item in reversed(xvals) if item is not None)
            except StopIteration:
                pass
            else:
                ax.text(toplblposn, -0.75 + 0.6 * random.random(), rxn['sys'],
                    verticalalignment='bottom', horizontalalignment='center',
                    family='Times New Roman', fontsize=8)
                ax.text(botlblposn, -1 * Nweft - 0.75 + 0.6 * random.random(), rxn['sys'],
                    verticalalignment='bottom', horizontalalignment='center',
                    family='Times New Roman', fontsize=8)

    # plot trimmings
    if mae is not None:
        ax.plot([-x for x in mae], positions, 's', color='black')
    if mape is not None:  # equivalent to MAE for a 10 kcal/mol IE
        ax.plot([0.025 * x for x in mape], positions, 'o', color='black')
    plt.axvline(0, color='black')

    # save and show
    pltuid = title + '_' + hashlib.sha1(title + repr(labels) + repr(xlimit)).hexdigest()
    pltfile = expand_saveas(saveas, pltuid, def_prefix='thread_', relpath=relpath)
    files_saved = {}
    for ext in graphicsformat:
        savefile = pltfile + '.' + ext.lower()
        plt.savefig(savefile, transparent=True, format=ext, bbox_inches='tight')
        files_saved[ext.lower()] = savefile
    plt.show()

    if not (mousetext or mouselink or mouseimag):
        return files_saved, None
    else:
        dpi = 80
        img_width = fig.get_figwidth() * dpi
        img_height = fig.get_figheight() * dpi

        htmlcode = """<SCRIPT>\n"""
        htmlcode += """function mouseshow(db, rxn, val) {\n"""
        if mousetext or mouselink:
            htmlcode += """   var cid = document.getElementById("cid");\n"""
            if mousetext:
                htmlcode += """   cid.innerHTML = %s;\n""" % (mousetext)
            if mouselink:
                htmlcode += """   cid.href = %s;\n""" % (mouselink)
        if mouseimag:
            htmlcode += """   var cmpd_img = document.getElementById("cmpd_img");\n"""
            htmlcode += """   cmpd_img.src = %s;\n""" % (mouseimag)
        htmlcode += """}\n"""
        htmlcode += """</SCRIPT>\n"""

        htmlcode += """%s <BR>""" % (mousetitle)
        htmlcode += """Mouseover:<BR><a id="cid"></a><br>\n"""
        htmlcode += """<IMG SRC="%s" ismap usemap="#points" WIDTH="%d" HEIGHT="%d">\n""" % \
            (pltfile + '.png', img_width, img_height)
        htmlcode += """<IMG ID="cmpd_img" WIDTH="%d" HEIGHT="%d">\n""" % (200, 160)
        htmlcode += """<MAP name="points">\n"""

        # generating html image map code
        #   points sorted to avoid overlapping map areas that can overwhelm html for SSI
        #   y=0 on top for html and on bottom for mpl, so flip the numbers
        posnM.sort(key=lambda tup: tup[2])
        posnM.sort(key=lambda tup: tup[3])
        last = (0, 0)
        for dbse, rxn, val, x, y in posnM:
            if val is None:
                continue

            now = (int(x), int(y))
            if now == last:
                htmlcode += """<!-- map overlap! %s-%s %+.2f skipped -->\n""" % (dbse, rxn, val)
            else:
                htmlcode += """<AREA shape="rect" coords="%d,%d,%d,%d" onmouseover="javascript:mouseshow('%s', '%s', '%+.2f');">\n""" % \
                    (x - 2, img_height - y - 20,
                    x + 2, img_height - y + 20,
                    dbse, rxn, val)
                last = now

        htmlcode += """</MAP>\n"""

        return files_saved, htmlcode


#def thread_mouseover_web(pltfile, dbid, dbname, xmin, xmax, mcdats, labels, titles):
#    """Saves a plot with name *pltfile* with a slat representation of
#    the modelchems errors in *mcdat*. Mouseover shows geometry and error
#    from *labels* based on recipe of Andrew Dalke from
#    http://www.dalkescientific.com/writings/diary/archive/2005/04/24/interactive_html.html
#
#    """
#    from matplotlib.backends.backend_agg import FigureCanvasAgg
#    import matplotlib
#    import sapt_colors
#
#    cmpd_width = 200
#    cmpd_height = 160
#
#    nplots = len(mcdats)
#    fht = nplots * 0.8
#    fht = nplots * 0.8 * 1.4
#    fig = matplotlib.figure.Figure(figsize=(12.0, fht))
#    fig.subplots_adjust(left=0.01, right=0.99, hspace=0.3, top=0.8, bottom=0.2)
#    img_width = fig.get_figwidth() * 80
#    img_height = fig.get_figheight() * 80
#
#    htmlcode = """
#<SCRIPT>
#function mouseandshow(name, id, db, dbname) {
#  var cid = document.getElementById("cid");
#  cid.innerHTML = name;
#  cid.href = "fragmentviewer.py?name=" + id + "&dataset=" + db;
#  var cmpd_img = document.getElementById("cmpd_img");
#  cmpd_img.src = dbname + "/dimers/" + id + ".png";
#}
#</SCRIPT>
#
#Distribution of Fragment Errors in Interaction Energy (kcal/mol)<BR>
#Mouseover:<BR><a id="cid"></a><br>
#<IMG SRC="scratch/%s" ismap usemap="#points" WIDTH="%d" HEIGHT="%d">
#<IMG ID="cmpd_img" WIDTH="%d" HEIGHT="%d">
#<MAP name="points">
#""" % (pltfile, img_width, img_height, cmpd_width, cmpd_height)
#
#    for item in range(nplots):
#        print '<br><br><br><br><br><br>'
#        mcdat = mcdats[item]
#        label = labels[item]
#        tttle = titles[item]
#
#        erdat = np.array(mcdat)
#        # No masked_array because interferes with html map
#        #erdat = np.ma.masked_array(mcdat, mask=mask)
#        yvals = np.ones(len(mcdat))
#        y = np.array([sapt_colors.sapt_colors[dbname][i] for i in label])
#
#        ax = fig.add_subplot(nplots, 1, item + 1)
#        sc = ax.scatter(erdat, yvals, c=y, s=3000, marker="|", cmap=matplotlib.cm.jet, vmin=0, vmax=1)
#        ax.set_title(tttle, fontsize=8)
#        ax.set_yticks([])
#        lp = ax.plot([0, 0], [0.9, 1.1], color='#cccc00', lw=2)
#        ax.set_ylim([0.95, 1.05])
#        ax.text(xmin + 0.3, 1.0, stats(erdat), fontsize=7, family='monospace', verticalalignment='center')
#        if item + 1 == nplots:
#            ax.set_xticks([-12.0, -8.0, -4.0, -2.0, -1.0, 0.0, 1.0, 2.0, 4.0, 8.0, 12.0])
#            for tick in ax.xaxis.get_major_ticks():
#                tick.tick1line.set_markersize(0)
#                tick.tick2line.set_markersize(0)
#        else:
#            ax.set_xticks([])
#        ax.set_frame_on(False)
#        ax.set_xlim([xmin, xmax])
#
#        # Convert the data set points into screen space coordinates
#        #xyscreencoords = ax.transData.transform(zip(erdat, yvals))
#        xyscreencoords = ax.transData.transform(zip(erdat, yvals))
#        xcoords, ycoords = zip(*xyscreencoords)
#
#        # HTML image coordinates have y=0 on the top.  Matplotlib
#        # has y=0 on the bottom.  We'll need to flip the numbers
#        for cid, x, y, er in zip(label, xcoords, ycoords, erdat):
#            htmlcode += """<AREA shape="rect" coords="%d,%d,%d,%d" onmouseover="javascript:mouseandshow('%s %+.2f', '%s', %s, '%s');">\n""" % \
#                (x - 2, img_height - y - 20, x + 2, img_height - y + 20, cid, er, cid, dbid, dbname)
#
#    htmlcode += "</MAP>\n"
#    canvas = FigureCanvasAgg(fig)
#    canvas.print_figure('scratch/' + title, dpi=80, transparent=True)
#
#    #plt.savefig('mplflat_' + title + '.pdf', bbox_inches='tight', transparent=True, format='PDF')
#    #plt.savefig(os.environ['HOME'] + os.sep + 'mplflat_' + title + '.pdf', bbox_inches='tight', transparent=T    rue, format='PDF')
#
#    return htmlcode


def composition_tile(db, aa1, aa2):
    """Takes dictionary *db* of label, error pairs and amino acids *aa1*
    and *aa2* and returns a square array of all errors for that amino
    acid pair, buffered by zeros.

    """
    import re
    import numpy as np
    bfdbpattern = re.compile("\d\d\d([A-Z][A-Z][A-Z])-\d\d\d([A-Z][A-Z][A-Z])-\d")

    tiles = []
    for key, val in db.items():
        bfdbname = bfdbpattern.match(key)
        if (bfdbname.group(1) == aa1 and bfdbname.group(2) == aa2) or \
           (bfdbname.group(2) == aa1 and bfdbname.group(1) == aa2):
            tiles.append(val)

    dim = int(np.ceil(np.sqrt(len(tiles))))
    pad = dim * dim - len(tiles)
    tiles += [0] * pad

    return np.reshape(np.array(tiles), (dim, dim))


def iowa(mcdat, mclbl, title='', xlimit=2.0):
    """Saves a plot with (extensionless) name *pltfile* with an Iowa
    representation of the modelchems errors in *mcdat* for BBI/SSI-style
    *labels*.

    """
    import numpy as np
    import matplotlib
    import matplotlib.pyplot as plt

    aa = ['ARG', 'HIE', 'LYS', 'ASP', 'GLU', 'SER', 'THR', 'ASN', 'GLN', 'CYS', 'MET', 'GLY', 'ALA', 'VAL', 'ILE', 'LEU', 'PRO', 'PHE', 'TYR', 'TRP']
    #aa = ['ILE', 'LEU', 'ASP', 'GLU', 'PHE']
    err = dict(zip(mclbl, mcdat))

    # handle for frame, overall axis
    fig, axt = plt.subplots(figsize=(6, 6))

    axt.set_xticks(np.arange(len(aa)) + 0.3, minor=False)
    axt.set_yticks(np.arange(len(aa)) + 0.3, minor=False)
    axt.invert_yaxis()
    axt.xaxis.tick_top()
    axt.set_xticklabels(aa, minor=False, rotation=60, size='small')
    axt.set_yticklabels(aa, minor=False, size='small')
    axt.xaxis.set_tick_params(width=0, length=0)
    axt.yaxis.set_tick_params(width=0, length=0)
    #axt.set_title('%s' % (title), fontsize=16, verticalalignment='bottom')
    axt.text(10.0, -1.5, title, horizontalalignment='center', fontsize=16)

    # nill spacing between 20x20 heatmaps
    plt.subplots_adjust(hspace=0.001, wspace=0.001)

    index = 1
    for aa1 in aa:
        for aa2 in aa:
            cb = composition_tile(err, aa1, aa2)

            ax = matplotlib.axes.Subplot(fig, len(aa), len(aa), index)
            fig.add_subplot(ax)
            heatmap = ax.pcolor(cb, vmin=-xlimit, vmax=xlimit, cmap=plt.cm.PRGn)
            ax.set_xticks([])
            ax.set_yticks([])
            index += 1

    title = '_'.join(title.split())
    plt.savefig('iowa_' + title + '.pdf', bbox_inches='tight', transparent=True, format='PDF')
    plt.show()
    #plt.savefig(os.environ['HOME'] + os.sep + 'iowa_' + title + '.pdf', bbox_inches='tight', transparent=True, format='PDF')


if __name__ == "__main__":

    merge_dats = [
    {'db':'HSG', 'sys':'1', 'data':[0.3508, 0.1234, 0.0364, 0.0731, 0.0388]},
    {'db':'HSG', 'sys':'3', 'data':[0.2036, -0.0736, -0.1650, -0.1380, -0.1806]},
    {'db':'S22', 'sys':'14', 'data':[None, -3.2144, None, None, None]},
    {'db':'S22', 'sys':'15', 'data':[-1.5090, -2.5263, -2.9452, -2.8633, -3.1059]},
    {'db':'S22', 'sys':'22', 'data':[0.3046, -0.2632, -0.5070, -0.4925, -0.6359]}]

    threads(merge_dats, labels=['d', 't', 'dt', 'q', 'tq'], color='sapt',
        title='MP2-CPa[]z', mae=[0.25, 0.5, 0.5, 0.3, 1.0], mape=[20.1, 25, 15, 5.5, 3.6])

    more_dats = [
    {'mc':'MP2-CP-adz', 'data':[1.0, 0.8, 1.4, 1.6]},
    {'mc':'MP2-CP-adtz', 'data':[0.6, 0.2, 0.4, 0.6]},
    None,
    {'mc':'MP2-CP-adzagain', 'data':[1.0, 0.8, 1.4, 1.6]}]

    bars(more_dats, title='asdf')

    single_dats = [
    {'dbse':'HSG', 'sys':'1', 'data':[0.3508]},
    {'dbse':'HSG', 'sys':'3', 'data':[0.2036]},
    {'dbse':'S22', 'sys':'14', 'data':[None]},
    {'dbse':'S22', 'sys':'15', 'data':[-1.5090]},
    {'dbse':'S22', 'sys':'22', 'data':[0.3046]}]

    #flat(single_dats, color='sapt', title='fg_MP2_adz', mae=0.25, mape=20.1)

    flat([{'sys': '1', 'color': 0.6933450559423702, 'data': [0.45730000000000004]}, {'sys': '2', 'color': 0.7627027688599753, 'data': [0.6231999999999998]}, {'sys': '3', 'color': 0.7579958735528617, 'data': [2.7624999999999993]}, {'sys': '4', 'color': 0.7560883254421639, 'data': [2.108600000000001]}, {'sys': '5', 'color': 0.7515161912065955, 'data': [2.2304999999999993]}, {'sys': '6', 'color': 0.7235223893438876, 'data': [1.3782000000000014]}, {'sys': '7', 'color': 0.7120099024225569, 'data': [1.9519000000000002]}, {'sys': '8', 'color': 0.13721565059144678, 'data': [0.13670000000000004]}, {'sys': '9', 'color': 0.3087395095814767, 'data': [0.2966]}, {'sys': '10', 'color': 0.25493207637105103, 'data': [-0.020199999999999996]}, {'sys': '11', 'color': 0.24093814608979347, 'data': [-1.5949999999999998]}, {'sys': '12', 'color': 0.3304746631959777, 'data': [-1.7422000000000004]}, {'sys': '13', 'color': 0.4156050644764822, 'data': [0.0011999999999989797]}, {'sys': '14', 'color': 0.2667207259626991, 'data': [-2.6083999999999996]}, {'sys': '15', 'color': 0.3767053567641695, 'data': [-1.5090000000000003]}, {'sys': '16', 'color': 0.5572641509433963, 'data': [0.10749999999999993]}, {'sys': '17', 'color': 0.4788598239641578, 'data': [0.29669999999999996]}, {'sys': '18', 'color': 0.3799031371351281, 'data': [0.10209999999999964]}, {'sys': '19', 'color': 0.5053227185999078, 'data': [0.16610000000000014]}, {'sys': '20', 'color': 0.2967660584483015, 'data': [-0.37739999999999974]}, {'sys': '21', 'color': 0.38836460733750316, 'data': [-0.4712000000000005]}, {'sys': '22', 'color': 0.5585849893078809, 'data': [0.30460000000000065]}, {'sys': 'BzBz_PD36-1.8', 'color': 0.1383351040559965, 'data': [-1.1921]}, {'sys': 'BzBz_PD34-2.0', 'color': 0.23086034843049832, 'data': [-1.367]}, {'sys': 'BzBz_T-5.2', 'color': 0.254318060864096, 'data': [-0.32230000000000025]}, {'sys': 'BzBz_T-5.1', 'color': 0.26598486566733337, 'data': [-0.3428]}, {'sys': 'BzBz_T-5.0', 'color': 0.28011258347610224, 'data': [-0.36060000000000025]}, {'sys': 'PyPy_S2-3.9', 'color': 0.14520332101084785, 'data': [-0.9853000000000001]}, {'sys': 'PyPy_S2-3.8', 'color': 0.1690757103699542, 'data': [-1.0932]}, {'sys': 'PyPy_S2-3.5', 'color': 0.25615734567417053, 'data': [-1.4617]}, {'sys': 'PyPy_S2-3.7', 'color': 0.19566550224566906, 'data': [-1.2103999999999995]}, {'sys': 'PyPy_S2-3.6', 'color': 0.22476748600170826, 'data': [-1.3333]}, {'sys': 'BzBz_PD32-2.0', 'color': 0.31605681987208084, 'data': [-1.6637]}, {'sys': 'BzBz_T-4.8', 'color': 0.31533827331543723, 'data': [-0.38759999999999994]}, {'sys': 'BzBz_T-4.9', 'color': 0.2966146678069063, 'data': [-0.3759999999999999]}, {'sys': 'BzH2S-3.6', 'color': 0.38284814928043304, 'data': [-0.1886000000000001]}, {'sys': 'BzBz_PD32-1.7', 'color': 0.3128835191478639, 'data': [-1.8703999999999998]}, {'sys': 'BzMe-3.8', 'color': 0.24117892478245323, 'data': [-0.034399999999999986]}, {'sys': 'BzMe-3.9', 'color': 0.22230903086047088, 'data': [-0.046499999999999986]}, {'sys': 'BzH2S-3.7', 'color': 0.36724255203373696, 'data': [-0.21039999999999992]}, {'sys': 'BzMe-3.6', 'color': 0.284901522674611, 'data': [0.007099999999999884]}, {'sys': 'BzMe-3.7', 'color': 0.2621086166558813, 'data': [-0.01770000000000005]}, {'sys': 'BzBz_PD32-1.9', 'color': 0.314711251903219, 'data': [-1.7353999999999998]}, {'sys': 'BzBz_PD32-1.8', 'color': 0.3136181753200793, 'data': [-1.8039999999999998]}, {'sys': 'BzH2S-3.8', 'color': 0.3542001591399945, 'data': [-0.22230000000000016]}, {'sys': 'BzBz_PD36-1.9', 'color': 0.14128552184232473, 'data': [-1.1517]}, {'sys': 'BzBz_S-3.7', 'color': 0.08862098445220466, 'data': [-1.3414]}, {'sys': 'BzH2S-4.0', 'color': 0.33637540012259076, 'data': [-0.2265999999999999]}, {'sys': 'BzBz_PD36-1.5', 'color': 0.13203548045236127, 'data': [-1.3035]}, {'sys': 'BzBz_S-3.8', 'color': 0.0335358832178858, 'data': [-1.2022]}, {'sys': 'BzBz_S-3.9', 'color': 0.021704594689389095, 'data': [-1.0747]}, {'sys': 'PyPy_T3-5.1', 'color': 0.3207725129126432, 'data': [-0.2958000000000003]}, {'sys': 'PyPy_T3-5.0', 'color': 0.3254925304351165, 'data': [-0.30710000000000015]}, {'sys': 'BzBz_PD36-1.7', 'color': 0.13577087141986593, 'data': [-1.2333000000000003]}, {'sys': 'PyPy_T3-4.8', 'color': 0.3443704059902452, 'data': [-0.32010000000000005]}, {'sys': 'PyPy_T3-4.9', 'color': 0.3333442013628509, 'data': [-0.3158999999999996]}, {'sys': 'PyPy_T3-4.7', 'color': 0.35854000505665756, 'data': [-0.31530000000000014]}, {'sys': 'BzBz_PD36-1.6', 'color': 0.13364651314909243, 'data': [-1.2705000000000002]}, {'sys': 'BzMe-4.0', 'color': 0.20560117919562013, 'data': [-0.05389999999999984]}, {'sys': 'MeMe-3.6', 'color': 0.16934865900383142, 'data': [0.18420000000000003]}, {'sys': 'MeMe-3.7', 'color': 0.1422332591197123, 'data': [0.14680000000000004]}, {'sys': 'MeMe-3.4', 'color': 0.23032794290360467, 'data': [0.29279999999999995]}, {'sys': 'MeMe-3.5', 'color': 0.19879551978386897, 'data': [0.23260000000000003]}, {'sys': 'MeMe-3.8', 'color': 0.11744404936205816, 'data': [0.11680000000000001]}, {'sys': 'BzBz_PD34-1.7', 'color': 0.22537382457222138, 'data': [-1.5286999999999997]}, {'sys': 'BzBz_PD34-1.6', 'color': 0.22434088042760192, 'data': [-1.5754000000000001]}, {'sys': 'BzBz_PD32-2.2', 'color': 0.3189891685300601, 'data': [-1.5093999999999999]}, {'sys': 'BzBz_S-4.1', 'color': 0.10884135031532088, 'data': [-0.8547000000000002]}, {'sys': 'BzBz_S-4.0', 'color': 0.06911476296747143, 'data': [-0.9590000000000001]}, {'sys': 'BzBz_PD34-1.8', 'color': 0.22685419834431494, 'data': [-1.476]}, {'sys': 'BzBz_PD34-1.9', 'color': 0.2287079261672095, 'data': [-1.4223999999999997]}, {'sys': 'BzH2S-3.9', 'color': 0.3439077006047999, 'data': [-0.22739999999999982]}, {'sys': 'FaNNFaNN-4.1', 'color': 0.7512716174974567, 'data': [1.7188999999999997]}, {'sys': 'FaNNFaNN-4.0', 'color': 0.7531388297328865, 'data': [1.9555000000000007]}, {'sys': 'FaNNFaNN-4.3', 'color': 0.7478064149182957, 'data': [1.2514000000000003]}, {'sys': 'FaNNFaNN-4.2', 'color': 0.7493794908838113, 'data': [1.4758000000000013]}, {'sys': 'FaOOFaON-4.0', 'color': 0.7589275618320565, 'data': [2.0586]}, {'sys': 'FaOOFaON-3.7', 'color': 0.7619465815742713, 'data': [3.3492999999999995]}, {'sys': 'FaOOFaON-3.9', 'color': 0.7593958895631474, 'data': [2.4471000000000007]}, {'sys': 'FaOOFaON-3.8', 'color': 0.7605108059280967, 'data': [2.8793999999999986]}, {'sys': 'FaONFaON-4.1', 'color': 0.7577459277014137, 'data': [1.8697999999999997]}, {'sys': 'FaOOFaON-3.6', 'color': 0.7633298028299997, 'data': [3.847599999999998]}, {'sys': 'FaNNFaNN-3.9', 'color': 0.7548200901251662, 'data': [2.2089]}, {'sys': 'FaONFaON-3.8', 'color': 0.7582294603551467, 'data': [2.967699999999999]}, {'sys': 'FaONFaON-3.9', 'color': 0.7575285282217349, 'data': [2.578900000000001]}, {'sys': 'FaONFaON-4.2', 'color': 0.7594549221042256, 'data': [1.5579999999999998]}, {'sys': 'FaOOFaNN-3.6', 'color': 0.7661655616885379, 'data': [3.701599999999999]}, {'sys': 'FaOOFaNN-3.7', 'color': 0.7671068376007428, 'data': [3.156500000000001]}, {'sys': 'FaOOFaNN-3.8', 'color': 0.766947626251711, 'data': [2.720700000000001]}, {'sys': 'FaONFaNN-3.9', 'color': 0.7569836601896789, 'data': [2.4281000000000006]}, {'sys': 'FaONFaNN-3.8', 'color': 0.758024548462959, 'data': [2.7561999999999998]}, {'sys': 'FaOOFaOO-3.6', 'color': 0.7623422640217077, 'data': [3.851800000000001]}, {'sys': 'FaOOFaOO-3.7', 'color': 0.7597430792159379, 'data': [3.2754999999999974]}, {'sys': 'FaOOFaOO-3.4', 'color': 0.7672554950739594, 'data': [5.193299999999999]}, {'sys': 'FaOOFaOO-3.5', 'color': 0.764908813123865, 'data': [4.491900000000001]}, {'sys': 'FaONFaNN-4.2', 'color': 0.7549212942233738, 'data': [1.534699999999999]}, {'sys': 'FaONFaNN-4.0', 'color': 0.7559404310956357, 'data': [2.1133000000000024]}, {'sys': 'FaONFaNN-4.1', 'color': 0.7551574698775625, 'data': [1.813900000000002]}, {'sys': 'FaONFaON-4.0', 'color': 0.7572064604483282, 'data': [2.2113999999999994]}, {'sys': 'FaOOFaOO-3.8', 'color': 0.7573810956831686, 'data': [2.7634000000000007]}, {'sys': '1', 'color': 0.2784121805328983, 'data': [0.3508]}, {'sys': '2', 'color': 0.22013842798900166, 'data': [-0.034600000000000186]}, {'sys': '3', 'color': 0.12832496088281312, 'data': [0.20360000000000023]}, {'sys': '4', 'color': 0.6993695033529733, 'data': [1.9092000000000002]}, {'sys': '5', 'color': 0.7371192790053749, 'data': [1.656600000000001]}, {'sys': '6', 'color': 0.5367033190796172, 'data': [0.27970000000000006]}, {'sys': '7', 'color': 0.3014220615964802, 'data': [0.32289999999999974]}, {'sys': '8', 'color': 0.01605867807629261, 'data': [0.12199999999999994]}, {'sys': '9', 'color': 0.6106300539083558, 'data': [0.3075999999999999]}, {'sys': '10', 'color': 0.6146680031333968, 'data': [0.6436000000000002]}, {'sys': '11', 'color': 0.6139747851721759, 'data': [0.4551999999999996]}, {'sys': '12', 'color': 0.32122739401126593, 'data': [0.44260000000000005]}, {'sys': '13', 'color': 0.24678148099136055, 'data': [-0.11789999999999967]}, {'sys': '14', 'color': 0.23700950710597016, 'data': [0.42689999999999995]}, {'sys': '15', 'color': 0.23103396678138563, 'data': [0.3266]}, {'sys': '16', 'color': 0.1922070769654413, 'data': [0.0696000000000001]}, {'sys': '17', 'color': 0.19082151944747366, 'data': [0.11159999999999992]}, {'sys': '18', 'color': 0.2886200282444196, 'data': [0.4114]}, {'sys': '19', 'color': 0.23560171133945224, 'data': [-0.1392]}, {'sys': '20', 'color': 0.3268270751294533, 'data': [0.5593]}, {'sys': '21', 'color': 0.7324460869158442, 'data': [0.6806000000000001]}],
    color='sapt', title='MP2-CP-adz', mae=1.21356003247, mape=24.6665886087, xlimit=4.0)

    lin_dats = [-0.5, -0.4, -0.3, 0, .5, .8, 5]
    lin_labs = ['008ILE-012LEU-1', '012LEU-085ASP-1', '004GLU-063LEU-2',
        '011ILE-014PHE-1', '027GLU-031LEU-1', '038PHE-041ILE-1', '199LEU-202GLU-1']
    iowa(lin_dats, lin_labs, title='ttl', xlimit=0.5)

    disthist(lin_dats)
psi4-data 1:0.3-5 / usr / share / psi / python / qcdb / mpl.py
