dlcplots.py

# -*- coding: utf-8 -*-
"""
Created on Tue Sep 16 10:21:11 2014

@author: dave
"""

from __future__ import division
from __future__ import print_function
#print(*objects, sep=' ', end='\n', file=sys.stdout)

import os
import socket
import gc

import numpy as np

import matplotlib.pyplot as plt
import matplotlib as mpl
#from matplotlib.figure import Figure
#from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigCanvas
#from scipy import interpolate as interp
#from scipy.optimize import fmin_slsqp
#from scipy.optimize import minimize
#from scipy.interpolate import interp1d
#import scipy.integrate as integrate
#http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
import pandas as pd

#import openpyxl as px
#import numpy as np

#import windIO
import mplutils
import Simulations as sim
import dlcdefs

plt.rc('font', family='serif')
plt.rc('xtick', labelsize=10)
plt.rc('ytick', labelsize=10)
plt.rc('axes', labelsize=12)
# do not use tex on Gorm
if not socket.gethostname()[:2] == 'g-':
    plt.rc('text', usetex=True)
plt.rc('legend', fontsize=11)
plt.rc('legend', numpoints=1)
plt.rc('legend', borderaxespad=0)

# =============================================================================
### STAT PLOTS
# =============================================================================

def plot_stats(sim_ids, post_dirs, fig_dir_base=None):
    """
    For each wind speed, take the max of the max.

    Only one or two sim_ids are supported. When they are from different
    models/projects, specify for each sim_id a different post_dir

    Parameters
    ----------

    sim_ids : list
        list of sim_id's, 1 or 2

    post_dirs
        list of post_dir's, 1 or 2. If 2, should correspond to sim_ids

    fig_dir_base : str, default=None

    """

    # if sim_id is a list, combine the two dataframes into one
    df_stats = pd.DataFrame()
    if type(sim_ids).__name__ == 'list':
        for ii, sim_id in enumerate(sim_ids):
            if isinstance(post_dirs, list):
                post_dir = post_dirs[ii]
            else:
                post_dir = post_dirs
            cc = sim.Cases(post_dir, sim_id, rem_failed=True)
            if ii == 0:
                df_stats, _, _ = cc.load_stats()
            else:
                # because there is no unique index, we will ignore it
                df_stats, _, _ = pd.concat([df_stats, cc.load_stats()],
                                            ignore_index=True)
    else:
        sim_id = sim_ids
        sim_ids = False
        post_dir = post_dirs
        cc = sim.Cases(post_dir, sim_id, rem_failed=True)
        df_stats, _, _ = cc.load_stats()

#    if force_dir:
#        cc.change_results_dir(resdir=force_dir)
#        for case in cc.cases:
#            sim_id = cc.cases[case]['[post_dir]']
#            cc.cases[case]['[post_dir]'] = post_dir

#    # add DLC category
#    f = lambda x: x.split('_')[0]
#    df_stats['DLC'] = df_stats['[Case id.]'].map(f)

#    fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(12,8), num=1)

    # define the number of positions you want to have the color for
    N = 22
    # select a color map
    cmap = mpl.cm.get_cmap('jet', N)
    # convert to array
    cmap_arr = cmap(np.arange(N))
    # color=cmap_arr[icol][0:3]

    # make a stastics plot for each channel
    gb_ch = df_stats.groupby(df_stats.channel)

    # channel selection
    plot_chans = {}
    plot_chans['DLL-2-inpvec-2'] = 'P_e'
    plot_chans['bearing-shaft_rot-angle_speed-rpm'] = 'RPM'

    plot_chans['tower-tower-node-001-momentvec-x'] = 'M_x T_B'
    plot_chans['tower-tower-node-001-momentvec-y'] = 'M_y T_B'
    plot_chans['tower-tower-node-001-momentvec-z'] = 'M_z T_B'

    plot_chans['tower-tower-node-008-momentvec-z'] = 'M_x T_T'
    plot_chans['tower-tower-node-008-momentvec-z'] = 'M_y T_T'
    plot_chans['tower-tower-node-008-momentvec-z'] = 'M_z T_T'

    plot_chans['shaft-shaft-node-004-momentvec-x'] = 'M_x Shaft_{MB}'
    plot_chans['shaft-shaft-node-004-momentvec-y'] = 'M_y Shaft_{MB}'
    plot_chans['shaft-shaft-node-004-momentvec-z'] = 'M_z Shaft_{MB}'

    plot_chans['blade1-blade1-node-003-momentvec-x'] = 'M_x B1_{root}'
    plot_chans['blade1-blade1-node-003-momentvec-y'] = 'M_y B1_{root}'
    plot_chans['blade1-blade1-node-003-momentvec-z'] = 'M_z B1_{root}'
    plot_chans['blade2-blade2-node-003-momentvec-x'] = 'M_x B2_{root}'
    plot_chans['blade2-blade2-node-003-momentvec-y'] = 'M_y B2_{root}'
    plot_chans['blade2-blade2-node-003-momentvec-z'] = 'M_z B2_{root}'
    plot_chans['blade3-blade3-node-003-momentvec-x'] = 'M_x B3_{root}'
    plot_chans['blade3-blade3-node-003-momentvec-y'] = 'M_y B3_{root}'
    plot_chans['blade3-blade3-node-003-momentvec-z'] = 'M_z B3_{root}'

    plot_chans['DLL-5-inpvec-1'] = 'Min tower clearance'

    plot_chans['bearing-pitch1-angle-deg'] = 'B1_{pitch}'
    plot_chans['bearing-pitch2-angle-deg'] = 'B2_{pitch}'
    plot_chans['bearing-pitch3-angle-deg'] = 'B3_{pitch}'

    plot_chans['setbeta-bladenr-1-flapnr-1'] = 'B1_{flap}'
    plot_chans['setbeta-bladenr-2-flapnr-1'] = 'B2_{flap}'
    plot_chans['setbeta-bladenr-3-flapnr-1'] = 'B3_{flap}'

    mfcs1 = ['k', 'w']
    mfcs2 = ['b', 'w']
    mfcs3 = ['r', 'w']
    stds = ['r', 'b']

    for nr, (ch_name, gr_ch) in enumerate(gb_ch):
        if ch_name not in plot_chans:
            continue
        for dlc_name, gr_ch_dlc in gr_ch.groupby(df_stats['[DLC]']):
            print('start plotting:  %s %s' % (str(dlc_name).ljust(7), ch_name))

            fig, axes = mplutils.make_fig(nrows=1, ncols=1, figsize=(7,5))
            ax = axes[0,0]
            # seperate figure for the standard deviations
            fig2, axes2 = mplutils.make_fig(nrows=1, ncols=1, figsize=(7,5))
            ax2 = axes2[0,0]

            if fig_dir_base is None and not sim_ids:
                res_dir = gr_ch_dlc['[res_dir]'][:1].values[0]
                run_dir = gr_ch_dlc['[run_dir]'][:1].values[0]
                fig_dir = os.path.join(fig_dir_base, res_dir)
            elif fig_dir_base is None and isinstance(sim_ids, list):
                fig_dir = os.path.join(fig_dir_base, '-'.join(sim_ids))
            elif fig_dir_base and not sim_ids:
                res_dir = gr_ch_dlc['[res_dir]'][:1].values[0]
                fig_dir = os.path.join(fig_dir_base, res_dir)
            elif sim_ids and fig_dir_base is not None:
                # create the compare directory if not defined
                fig_dir = fig_dir_base

            # if we have a list of different cases, we also need to group those
            # because the sim_id wasn't saved before in the data frame,
            # we need to derive that from the run dir
            # if there is only one run dir nothing changes
            ii = 0
            sid_names = []
            for run_dir, gr_ch_dlc_sid in gr_ch_dlc.groupby(df_stats['[run_dir]']):
                sid_name = run_dir.split('/')[-2]
                sid_names.append(sid_name)
                print(sid_name)
                wind = gr_ch_dlc_sid['[Windspeed]'].values
                dmin = gr_ch_dlc_sid['min'].values
                dmean = gr_ch_dlc_sid['mean'].values
                dmax = gr_ch_dlc_sid['max'].values
                dstd = gr_ch_dlc_sid['std'].values
                if not sim_ids:
                    lab1 = 'mean'
                    lab2 = 'min'
                    lab3 = 'max'
                    lab4 = 'std'
                else:
                    lab1 = 'mean %s' % sid_name
                    lab2 = 'min %s' % sid_name
                    lab3 = 'max %s' % sid_name
                    lab4 = 'std %s' % sid_name
                mfc1 = mfcs1[ii]
                mfc2 = mfcs2[ii]
                mfc3 = mfcs3[ii]
                ax.plot(wind, dmean, mec='k', marker='o', mfc=mfc1, ls='',
                        label=lab1, alpha=0.7)
                ax.plot(wind, dmin, mec='b', marker='^', mfc=mfc2, ls='',
                        label=lab2, alpha=0.7)
                ax.plot(wind, dmax, mec='r', marker='v', mfc=mfc3, ls='',
                        label=lab3, alpha=0.7)

                ax2.plot(wind, dstd, mec=stds[ii], marker='s', mfc=stds[ii], ls='',
                        label=lab4, alpha=0.7)

                ii += 1

#            for wind, gr_wind in  gr_ch_dlc.groupby(df_stats['[Windspeed]']):
#                wind = gr_wind['[Windspeed]'].values
#                dmin = gr_wind['min'].values#.mean()
#                dmean = gr_wind['mean'].values#.mean()
#                dmax = gr_wind['max'].values#.mean()
##                dstd = gr_wind['std'].mean()
#                ax.plot(wind, dmean, 'ko', label='mean', alpha=0.7)
#                ax.plot(wind, dmin, 'b^', label='min', alpha=0.7)
#                ax.plot(wind, dmax, 'rv', label='max', alpha=0.7)
##                ax.errorbar(wind, dmean, c='k', ls='', marker='s', mfc='w',
##                        label='mean and std', yerr=dstd)
            ax.grid()
            ax.set_xlim([3, 27])
            leg = ax.legend(loc='best', ncol=2)
            leg.get_frame().set_alpha(0.7)
            ax.set_title(r'{DLC%s} $%s$' % (dlc_name, plot_chans[ch_name]))
            ax.set_xlabel('Wind speed [m/s]')
            fig.tight_layout()
            fig.subplots_adjust(top=0.92)
            if not sim_ids:
                fig_path = os.path.join(fig_dir,
                                        ch_name.replace(' ', '_') + '.png')
            else:
                sids = '_'.join(sid_names)
#                fig_dir = run_dir.split('/')[:-1] + 'figures/'
                fname = '%s_%s.png' % (ch_name.replace(' ', '_'), sids)
                fig_path = os.path.join(fig_dir, 'dlc%s/' % dlc_name)
                if not os.path.exists(fig_path):
                    os.makedirs(fig_path)
                fig_path = fig_path + fname
            fig.savefig(fig_path)#.encode('latin-1')
#            canvas.close()
            fig.clear()
            print('saved: %s' % fig_path)


            ax2.grid()
            ax2.set_xlim([3, 27])
            leg = ax2.legend(loc='best', ncol=2)
            leg.get_frame().set_alpha(0.7)
            ax2.set_title(r'{DLC%s} $%s$' % (dlc_name, plot_chans[ch_name]))
            ax2.set_xlabel('Wind speed [m/s]')
            fig2.tight_layout()
            fig2.subplots_adjust(top=0.92)
            if not sim_ids:
                fig_path = os.path.join(fig_dir,
                                        ch_name.replace(' ', '_') + '_std.png')
            else:
                sids = '_'.join(sid_names)
                fname = '%s_std_%s.png' % (ch_name.replace(' ', '_'), sids)
                fig_path = os.path.join(fig_dir, 'dlc%s/' % dlc_name)
                if not os.path.exists(fig_path):
                    os.makedirs(fig_path)
                fig_path = fig_path + fname
            fig2.savefig(fig_path)#.encode('latin-1')
#            canvas.close()
            fig2.clear()
            print('saved: %s' % fig_path)


def plot_stats2(sim_ids, post_dirs, fig_dir_base=None, labels=None,
                post_dir_save=False, dlc_ignore=['00']):
    """
    Map which channels have to be compared
    """

    plot_chans = {}

    plot_chans['B1_{flap}'] = ['setbeta-bladenr-1-flapnr-1']
    plot_chans['B2_{flap}'] = ['setbeta-bladenr-2-flapnr-1']
    plot_chans['B3_{flap}'] = ['setbeta-bladenr-3-flapnr-1']
    plot_chans['M_x B1_{root}'] = ['blade1-blade1-node-003-momentvec-x',
                                   'blade1-blade1-node-004-momentvec-x']
    plot_chans['M_y B1_{root}'] = ['blade1-blade1-node-003-momentvec-y',
                                   'blade1-blade1-node-004-momentvec-y']
    plot_chans['M_z B1_{root}'] = ['blade1-blade1-node-003-momentvec-z',
                                   'blade1-blade1-node-004-momentvec-z']
    plot_chans['B3_{pitch}'] = ['bearing-pitch3-angle-deg']
    plot_chans['RPM'] = ['bearing-shaft_rot-angle_speed-rpm']
    plot_chans['P_e'] = ['DLL-2-inpvec-2']
    plot_chans['P_{mech}'] = ['stats-shaft-power']
    plot_chans['M_x B3_{root}'] = ['blade3-blade3-node-003-momentvec-x',
                                   'blade3-blade3-node-004-momentvec-x']
    plot_chans['M_y B3_{root}'] = ['blade3-blade3-node-003-momentvec-y',
                                   'blade3-blade3-node-004-momentvec-y']
    plot_chans['M_z B3_{root}'] = ['blade3-blade3-node-003-momentvec-z',
                                   'blade3-blade3-node-004-momentvec-z']
    plot_chans['B2_{pitch}'] = ['bearing-pitch2-angle-deg']

    plot_chans['B3 U_y'] = ['global-blade3-elem-018-zrel-1.00-State pos-y']
    plot_chans['M_z B2_{root}'] = ['blade2-blade2-node-003-momentvec-z',
                                   'blade2-blade2-node-004-momentvec-z']
    plot_chans['M_x B2_{root}'] = ['blade2-blade2-node-003-momentvec-x',
                                   'blade2-blade2-node-004-momentvec-x']
    plot_chans['M_y B2_{root}'] = ['blade2-blade2-node-003-momentvec-y',
                                   'blade2-blade2-node-004-momentvec-y']
    plot_chans['B1_{pitch}'] = ['bearing-pitch1-angle-deg']
    plot_chans['M_x T_B'] = ['tower-tower-node-001-momentvec-x']
    plot_chans['M_y T_B'] = ['tower-tower-node-001-momentvec-y']
    plot_chans['M_z T_B'] = ['tower-tower-node-001-momentvec-z']
    plot_chans['tower clearance'] = ['DLL-5-inpvec-1']
    plot_chans['M_z T_T'] = ['tower-tower-node-008-momentvec-z']
    plot_chans['M_y Shaft_{MB}'] = ['shaft-shaft-node-004-momentvec-y']
    plot_chans['M_x Shaft_{MB}'] = ['shaft-shaft-node-004-momentvec-x']
    plot_chans['M_z Shaft_{MB}'] = ['shaft-shaft-node-004-momentvec-z']

    # reduce required memory, only use following columns
    cols = ['[run_dir]', '[DLC]', 'channel', '[res_dir]', '[Windspeed]',
            'mean', 'max', 'min', 'std', '[wdir]']

    # if sim_id is a list, combine the two dataframes into one
    df_stats = pd.DataFrame()
    if type(sim_ids).__name__ == 'list':
        for ii, sim_id in enumerate(sim_ids):
            if isinstance(post_dirs, list):
                post_dir = post_dirs[ii]
            else:
                post_dir = post_dirs
            cc = sim.Cases(post_dir, sim_id, rem_failed=True)
            df_stats, _, _ = cc.load_stats(columns=cols, leq=False)
            print('%s Cases loaded.' % sim_id)

            # if specified, save the merged sims elsewhere
            if post_dir_save:
                fpath = os.path.join(post_dir_save, '-'.join(sim_ids) + '.h5')
                try:
                    os.makedirs(post_dir_save)
                except OSError:
                    pass
            else:
                fpath = os.path.join(post_dir, '-'.join(sim_ids) + '.h5')
            if ii == 0:
                # and save somewhere so we can add the second data frame on
                # disc
                df_stats.to_hdf(fpath, 'table', mode='w', format='table',
                                complevel=9, complib='blosc')
                print('%s merged stats written to: %s' % (sim_id, fpath))
            else:
                # instead of doing a concat in memory, add to the hdf store
                df_stats.to_hdf(fpath, 'table', mode='r+', format='table',
                                complevel=9, complib='blosc', append=True)
                print('%s merging stats into:      %s' % (sim_id, fpath))
#                df_stats = pd.concat([df_stats, df_stats2], ignore_index=True)
#                df_stats2 = None
            # we might run into memory issues
            del df_stats, _, cc
            gc.collect()
        # and load the reduced combined set
        print('loading merged stats:            %s' % fpath)
        df_stats = pd.read_hdf(fpath, 'table')
    else:
        sim_id = sim_ids
        sim_ids = [sim_id]
        post_dir = post_dirs
        cc = sim.Cases(post_dir, sim_id, rem_failed=True)
        df_stats, _, _ = cc.load_stats(leq=False)

    mfcs1 = ['k', 'w']
    mfcs2 = ['b', 'w']
    mfcs3 = ['r', 'w']
    stds = ['r', 'b']

    # first, take each DLC appart
    for dlc_name, gr_dlc in df_stats.groupby(df_stats['[DLC]']):
        # do not plot the stats for dlc00
        if dlc_name in dlc_ignore:
            continue
        # cycle through all the target plot channels
        for ch_dscr, ch_names in plot_chans.iteritems():
            # second, group per channel. Note that when the channel names are not
            # identical, we need to manually pick them.
            # figure file name will be the first channel
            if isinstance(ch_names, list):
                df_chan = gr_dlc[gr_dlc.channel == ch_names[0]]
                fname_base = ch_names[0].replace(' ', '_')
                try:
                    df2 = gr_dlc[gr_dlc.channel == ch_names[1]]
                    df_chan = pd.concat([df_chan, df2], ignore_index=True)
                except IndexError:
                    pass
            else:
                ch_name = ch_names
                ch_names = [ch_name]
                df_chan = gr_dlc[gr_dlc.channel == ch_names]
                fname_base = ch_names.replace(' ', '_')

            # if not, than we are missing a channel description, or the channel
            # is simply not available in the given result set
#            if not len(df_chan.channel.unique()) == len(ch_names):
#                continue
            lens = []
            for key, gr_ch_dlc_sid in df_chan.groupby(df_chan['[run_dir]']):
                lens.append(len(gr_ch_dlc_sid))
            # when the channel is simply not present
            if len(lens) == 0:
                continue
            # when only one of the channels was present, but the set is still
            # complete.
            # FIXME: what if both channels are present?
            if len(ch_names) > 1 and (lens[0] < 1 or lens[1] < 1):
                continue

            print('start plotting:  %s %s' % (str(dlc_name).ljust(7), ch_dscr))

            fig, axes = mplutils.make_fig(nrows=1, ncols=1,
                                           figsize=(11,7.15), dpi=120)
            ax = axes[0,0]
            # seperate figure for the standard deviations
            fig2, axes2 = mplutils.make_fig(nrows=1, ncols=1,
                                             figsize=(11,7.15), dpi=120)
            ax2 = axes2[0,0]

            if fig_dir_base is None and len(sim_ids) < 2:
                res_dir = df_chan['[res_dir]'][:1].values[0]
                fig_dir = os.path.join(fig_dir_base, res_dir)
            elif fig_dir_base is None and isinstance(sim_ids, list):
                fig_dir = os.path.join(fig_dir_base, '-'.join(sim_ids))
#            elif fig_dir_base and len(sim_ids) < 2:
#                res_dir = df_chan['[res_dir]'][:1].values[0]
#                fig_dir = os.path.join(fig_dir_base, res_dir)
            elif sim_ids and fig_dir_base is not None:
                # create the compare directory if not defined
                fig_dir = fig_dir_base

            # if we have a list of different cases, we also need to group those
            # because the sim_id wasn't saved before in the data frame,
            # we need to derive that from the run dir
            # if there is only one run dir nothing changes
            ii = 0
            sid_names = []
            for run_dir, gr_ch_dlc_sid in df_chan.groupby(df_chan['[run_dir]']):
                if labels is None:
                    sid_name = run_dir.split(os.path.sep)[-2]
                else:
                    sid_name = labels[ii]
                sid_names.append(sid_name)
                print('   sim_id/label:', sid_name)
                # FIXME: will this go wrong in PY3?
                if str(dlc_name) in ['61', '62']:
                    xdata = gr_ch_dlc_sid['[wdir]'].values
                    xlabel = 'wind direction [deg]'
                    xlims = [0, 360]
                else:
                    xdata = gr_ch_dlc_sid['[Windspeed]'].values
                    xlabel = 'Wind speed [m/s]'
                    xlims = [3, 27]
                dmin = gr_ch_dlc_sid['min'].values
                dmean = gr_ch_dlc_sid['mean'].values
                dmax = gr_ch_dlc_sid['max'].values
                dstd = gr_ch_dlc_sid['std'].values
                if not sim_ids:
                    lab1 = 'mean'
                    lab2 = 'min'
                    lab3 = 'max'
                    lab4 = 'std'
                else:
                    lab1 = 'mean %s' % sid_name
                    lab2 = 'min %s' % sid_name
                    lab3 = 'max %s' % sid_name
                    lab4 = 'std %s' % sid_name
                mfc1 = mfcs1[ii]
                mfc2 = mfcs2[ii]
                mfc3 = mfcs3[ii]
                ax.errorbar(xdata, dmean, mec='k', marker='o', mfc=mfc1, ls='',
                            label=lab1, alpha=0.7, yerr=dstd)
                ax.plot(xdata, dmin, mec='b', marker='^', mfc=mfc2, ls='',
                        label=lab2, alpha=0.7)
                ax.plot(xdata, dmax, mec='r', marker='v', mfc=mfc3, ls='',
                        label=lab3, alpha=0.7)

                ax2.plot(xdata, dstd, mec=stds[ii], marker='s', mfc=stds[ii],
                        ls='', label=lab4, alpha=0.7)

                ii += 1

#            for wind, gr_wind in  gr_ch_dlc.groupby(df_stats['[Windspeed]']):
#                wind = gr_wind['[Windspeed]'].values
#                dmin = gr_wind['min'].values#.mean()
#                dmean = gr_wind['mean'].values#.mean()
#                dmax = gr_wind['max'].values#.mean()
##                dstd = gr_wind['std'].mean()
#                ax.plot(wind, dmean, 'ko', label='mean', alpha=0.7)
#                ax.plot(wind, dmin, 'b^', label='min', alpha=0.7)
#                ax.plot(wind, dmax, 'rv', label='max', alpha=0.7)
##                ax.errorbar(wind, dmean, c='k', ls='', marker='s', mfc='w',
##                        label='mean and std', yerr=dstd)
            ax.grid()
            ax.set_xlim(xlims)
            leg = ax.legend(loc='best', ncol=3)
            leg.get_frame().set_alpha(0.7)
            ax.set_title(r'{DLC%s} $%s$' % (dlc_name, ch_dscr))
            ax.set_xlabel(xlabel)
            fig.tight_layout()
            fig.subplots_adjust(top=0.92)
            if not sim_ids:
                fig_path = os.path.join(fig_dir, fname_base + '.png')
            else:
                sids = '_'.join(sid_names)
                fname = '%s_%s.png' % (fname_base, sids)
                fig_path = os.path.join(fig_dir, 'dlc%s/' % dlc_name)
                if not os.path.exists(fig_path):
                    os.makedirs(fig_path)
                fig_path = fig_path + fname
            fig.savefig(fig_path)#.encode('latin-1')
            fig.clear()
            print('saved: %s' % fig_path)

            ax2.grid()
            ax2.set_xlim(xlims)
            leg = ax2.legend(loc='best', ncol=3)
            leg.get_frame().set_alpha(0.7)
            ax2.set_title(r'{DLC%s} $%s$' % (dlc_name, ch_dscr))
            ax2.set_xlabel('Wind speed [m/s]')
            fig2.tight_layout()
            fig2.subplots_adjust(top=0.92)
            if not sim_ids:
                fig_path = os.path.join(fig_dir, fname_base + '_std.png')
            else:
                sids = '_'.join(sid_names)
                fname = '%s_std_%s.png' % (fname_base, sids)
                fig_path = os.path.join(fig_dir, 'dlc%s/' % dlc_name)
                if not os.path.exists(fig_path):
                    os.makedirs(fig_path)
                fig_path = fig_path + fname
            fig2.savefig(fig_path)#.encode('latin-1')
            fig2.clear()
            print('saved: %s' % fig_path)


class PlotStats(object):

    def __init__(self):
        pass

    def load_stats(self, sim_ids, post_dirs, post_dir_save=False):

        self.sim_ids = sim_ids
        self.post_dirs = post_dirs

        # reduce required memory, only use following columns
        cols = ['[run_dir]', '[DLC]', 'channel', '[res_dir]', '[Windspeed]',
                'mean', 'max', 'min', 'std', '[wdir]']

        # if sim_id is a list, combine the two dataframes into one
        df_stats = pd.DataFrame()
        if type(sim_ids).__name__ == 'list':
            for ii, sim_id in enumerate(sim_ids):
                if isinstance(post_dirs, list):
                    post_dir = post_dirs[ii]
                else:
                    post_dir = post_dirs
                cc = sim.Cases(post_dir, sim_id, rem_failed=True)
                df_stats, _, _ = cc.load_stats(columns=cols, leq=False)
                print('%s Cases loaded.' % sim_id)

                # if specified, save the merged sims elsewhere
                if post_dir_save:
                    fpath = os.path.join(post_dir_save, '-'.join(sim_ids) + '.h5')
                    try:
                        os.makedirs(post_dir_save)
                    except OSError:
                        pass
                else:
                    fpath = os.path.join(post_dir, '-'.join(sim_ids) + '.h5')
                if ii == 0:
                    # and save somewhere so we can add the second data frame on
                    # disc
                    df_stats.to_hdf(fpath, 'table', mode='w', format='table',
                                    complevel=9, complib='blosc')
                    print('%s merged stats written to: %s' % (sim_id, fpath))
                else:
                    # instead of doing a concat in memory, add to the hdf store
                    df_stats.to_hdf(fpath, 'table', mode='r+', format='table',
                                    complevel=9, complib='blosc', append=True)
                    print('%s merging stats into:      %s' % (sim_id, fpath))
                # we might run into memory issues
                del df_stats, _, cc
                gc.collect()
            # and load the reduced combined set
            print('loading merged stats:            %s' % fpath)
            df_stats = pd.read_hdf(fpath, 'table')
        else:
            sim_id = sim_ids
            sim_ids = [sim_id]
            post_dir = post_dirs
            cc = sim.Cases(post_dir, sim_id, rem_failed=True)
            df_stats, _, _ = cc.load_stats(leq=False)

        return df_stats

    def select_extremes_blade_radial(self, df):
        """
        For each radial position of the blade, find the extremes
        """

        def selector(x):
            """
            select following channels:
            'local-blade%1i-node-%03i-momentvec-x'
            'blade2-blade2-node-003-momentvec-y'
            """
            if x[:11] == 'local-blade' and x[22:31] == 'momentvec':
                return True
            else:
                return False

        # find all blade local load channels
        criteria = df.channel.map(lambda x: selector(x))
        df = df[criteria]
        # split channel columns so we can select channels properly
        df = df.join(df.channel.apply(lambda x: pd.Series(x.split('-'))))

        df_ext = {'dlc':[], 'case':[], 'node':[], 'max':[], 'min':[], 'comp':[]}

        def fillvalues(x, ii, maxmin):
            x['node'].append(m_group[3].ix[ii])
            x['dlc'].append(m_group['[DLC]'].ix[ii])
            x['case'].append(m_group['[case_id]'].ix[ii])
            x['comp'].append(m_group[5].ix[ii])
            if maxmin == 'max':
                x['max'].append(m_group['max'].ix[ii])
                x['min'].append(np.nan)
            else:
                x['max'].append(np.nan)
                x['min'].append(m_group['min'].ix[ii])
            return x

        # we take blade1, blade2, and blade3
        df_b2 = df[df[0]=='local']
#        df_b2 = df_b2[df_b2[1]=='blade2']
        df_b2 = df_b2[df_b2[4]=='momentvec']
#        df_b2 = df_b2[df_b2[5]=='x']
        # make sure we only have blade1, 2 and 3
        assert set(df_b2[1].unique()) == set(['blade3', 'blade2', 'blade1'])
#        # number of nodes
#        nrnodes = len(df_b2[3].unique())
        # group by node number, and take the max
        for nodenr, group in df_b2.groupby(df_b2[3]):
            print(nodenr, end='   ')
            for comp, m_group in df_b2.groupby(group[5]):
                print(comp)
                imax = m_group['max'].argmax()
                imin = m_group['min'].argmin()
                df_ext = fillvalues(df_ext, imax, 'max')
                df_ext = fillvalues(df_ext, imin, 'min')

        df_ext = pd.DataFrame(df_ext)
        df_ext.sort(columns='node', inplace=True)

        return df_ext

    def plot_extremes_blade_radial(self, df_ext, fpath):
        nrows = 2
        ncols = 2
        figsize = (11,7.15)
        fig, axes = mplutils.make_fig(nrows=nrows, ncols=ncols, figsize=figsize)

#        self.col = ['r', 'k']
#        self.alf = [1.0, 0.7]
#        self.i = 0

        mx_max = df_ext['max'][df_ext.comp == 'x'].dropna()
        mx_min = df_ext['min'][df_ext.comp == 'x'].dropna()
        my_max = df_ext['max'][df_ext.comp == 'y'].dropna()
        my_min = df_ext['min'][df_ext.comp == 'y'].dropna()
#        nodes = df_ext.node.ix[mx_max.index]
        axes[0,0].plot(mx_max, 'r^', label='$M_{x_{max}}$')
        axes[0,1].plot(mx_min, 'bv', label='$M_{x_{min}}$')
        axes[1,0].plot(my_max, 'r^', label='$M_{y_{max}}$')
        axes[1,1].plot(my_min, 'bv', label='$M_{y_{min}}$')

        axs = axes.ravel()
        for ax in axs:
            ax.grid()
            ax.legend(loc='best')

#        axs[0].set_xticklabels([])
#        axs[1].set_xticklabels([])
#        axs[2].set_xticklabels([])
#        axs[-1].set_xlabel('time [s]')

        fig.tight_layout()
        fig.subplots_adjust(hspace=0.06)
        fig.subplots_adjust(top=0.98)

        fdir = os.path.dirname(fpath)
#        fname = os.path.basename(fpath)

        if not os.path.exists(fdir):
            os.makedirs(fdir)
        print('saving: %s ...' % fpath, end='')
        fig.savefig(fpath)#.encode('latin-1')
        print('done')
        fig.clear()

        # save as tables
        df_ext.ix[mx_max.index].to_excel(fpath.replace('.png', '_mx_max.xls'))
        df_ext.ix[mx_min.index].to_excel(fpath.replace('.png', '_mx_min.xls'))
        df_ext.ix[my_max.index].to_excel(fpath.replace('.png', '_my_max.xls'))
        df_ext.ix[my_min.index].to_excel(fpath.replace('.png', '_my_min.xls'))

    def extract_leq_blade_radial(self, df_leq, fpath):

        def selector(x):
            """
            select following channels:
            'local-blade%1i-node-%03i-momentvec-x'
            'blade2-blade2-node-003-momentvec-y'
            """
            if x[:11] == 'local-blade' and x[22:31] == 'momentvec':
                return True
            else:
                return False

        # find all blade local load channels
        criteria = df_leq.channel.map(lambda x: selector(x))
        df = df_leq[criteria]
        # split channel columns so we can select channels properly
        df = df.join(df.channel.apply(lambda x: pd.Series(x.split('-'))))
        df.sort(columns=3, inplace=True)
        assert set(df[1].unique()) == set(['blade3', 'blade2', 'blade1'])

        leqs = df.keys()[1:10]
        df_ext = {leq:[] for leq in leqs}
        df_ext['node'] = []
        df_ext['comp'] = []

        for nodenr, group in df.groupby(df[3]):
            print(nodenr, end='   ')
            for comp, m_group in df.groupby(group[5]):
                print(comp)
                for leq in leqs:
                    df_ext[leq].append(m_group[leq].max())
                df_ext['node'].append(nodenr)
                df_ext['comp'].append(comp)

        df_ext = pd.DataFrame(df_ext)
        df_ext.sort(columns='node', inplace=True)

        df_ext[df_ext.comp=='x'].to_excel(fpath.replace('.xls', '_x.xls'))
        df_ext[df_ext.comp=='y'].to_excel(fpath.replace('.xls', '_y.xls'))
        df_ext[df_ext.comp=='z'].to_excel(fpath.replace('.xls', '_z.xls'))

        return df_ext


class PlotPerf(object):

    def __init__(self, nrows=4, ncols=1, figsize=(14,11)):

        self.fig, self.axes = mplutils.make_fig(nrows=nrows, ncols=ncols,
                                                 figsize=figsize)
#        self.axs = self.axes.ravel()
        self.col = ['r', 'k']
        self.alf = [1.0, 0.7]
        self.i = 0

    def plot(self, res, label_id):
        """
        """
        i = self.i

        sim_id = label_id
        time = res.sig[:,0]
        self.t0, self.t1 = time[0], time[-1]

        # find the wind speed
        for channame, chan in res.ch_dict.iteritems():
            if channame.startswith('windspeed-global-Vy-0.00-0.00'):
                break
        wind = res.sig[:,chan['chi']]

        chi = res.ch_dict['bearing-shaft_rot-angle_speed-rpm']['chi']
        rpm = res.sig[:,chi]

        chi = res.ch_dict['bearing-pitch1-angle-deg']['chi']
        pitch = res.sig[:,chi]

        chi = res.ch_dict['tower-tower-node-001-momentvec-x']['chi']
        tx = res.sig[:,chi]

        chi = res.ch_dict['tower-tower-node-001-momentvec-y']['chi']
        ty = res.sig[:,chi]

        chi = res.ch_dict['DLL-2-inpvec-2']['chi']
        power = res.sig[:,chi]

        try:
            chi = res.ch_dict['Tors_e-1-100.11']['chi']
        except KeyError:
            chi = res.ch_dict['Tors_e-1-86.50']['chi']
        tors_1 = res.sig[:,chi]

#        try:
#            chi = res.ch_dict['Tors_e-1-96.22']['chi']
#        except:
#            chi = res.ch_dict['Tors_e-1-83.13']['chi']
#        tors_2 = res.sig[:,chi]

        try:
            chi = res.ch_dict['Tors_e-1-84.53']['chi']
        except:
            chi = res.ch_dict['Tors_e-1-73.21']['chi']
        tors_3 = res.sig[:,chi]

        ax = self.axes.ravel()
        ax[0].plot(time, wind, self.col[i]+'--', label='%s wind speed' % sim_id,
                   alpha=self.alf[i])
        ax[0].plot(time, pitch, self.col[i]+'-.', label='%s pitch' % sim_id,
                   alpha=self.alf[i])
        ax[0].plot(time, rpm, self.col[i]+'-', label='%s RPM' % sim_id,
                   alpha=self.alf[i])

        ax[1].plot(time, tx, self.col[i]+'--', label='%s Tower FA' % sim_id,
                   alpha=self.alf[i])
        ax[1].plot(time, ty, self.col[i]+'-', label='%s Tower SS' % sim_id,
                   alpha=self.alf[i])

        ax[2].plot(time, power/1e6, self.col[i]+'-', alpha=self.alf[i],
                   label='%s El Power' % sim_id)

        ax[3].plot(time, tors_1, self.col[i]+'--', label='%s torsion tip' % sim_id,
                   alpha=self.alf[i])
#        ax[3].plot(time, tors_2, self.col[i]+'-.', label='%s torsion 96 pc' % sim_id,
#                   alpha=self.alf[i])
#        ax[3].plot(time, tors_3, self.col[i]+'-', label='%s torsion 84 pc' % sim_id,
#                   alpha=self.alf[i])

        self.i += 1

    def final(self, fig_path, fig_name):

        axs = self.axes.ravel()

        for ax in axs:
            ax.set_xlim([self.t0, self.t1])
            ax.grid()
            ax.legend(loc='best')

        axs[0].set_xticklabels([])
        axs[1].set_xticklabels([])
        axs[2].set_xticklabels([])
        axs[-1].set_xlabel('time [s]')

        self.fig.tight_layout()
        self.fig.subplots_adjust(hspace=0.06)
        self.fig.subplots_adjust(top=0.98)

        if not os.path.exists(fig_path):
            os.makedirs(fig_path)
        fname = os.path.join(fig_path, fig_name)
        print('saving: %s ...' % fname, end='')
        self.fig.savefig(fname)#.encode('latin-1')
        print('done')
        self.fig.clear()

def plot_dlc01_powercurve(sim_ids, post_dirs, run_dirs, fig_dir_base):
    """
    Create power curve based on steady DLC01 results
    Use the same format as for HS2 for easy comparison!
    """



def plot_dlc00(sim_ids, post_dirs, run_dirs, fig_dir_base=None, labels=None,
               cnames=['dlc00_stair_wsp04_25_noturb.htc',
                       'dlc00_ramp_wsp04_25_04_noturb.htc'], figsize=(14,11)):
    """
    This version is an update over plot_staircase.
    """

    stairs = []
    # if sim_id is a list, combine the two dataframes into one
    if type(sim_ids).__name__ == 'list':
        for ii, sim_id in enumerate(sim_ids):
            if isinstance(post_dirs, list):
                post_dir = post_dirs[ii]
            else:
                post_dir = post_dirs
            stairs.append(sim.Cases(post_dir, sim_id, rem_failed=True))
    else:
        post_dir = post_dirs
        stairs.append(sim.Cases(post_dir, sim_id, rem_failed=True))

    for cname in cnames:
        fp = PlotPerf(figsize=figsize)
        for i, cc in enumerate(stairs):
            if isinstance(cname, list):
                _cname = cname[i]
            else:
                _cname = cname
            if _cname in cc.cases_fail:
                print('no result for %s' % cc.sim_id)
                continue
            cc.change_results_dir(run_dirs[i])
            try:
                res = cc.load_result_file(cc.cases[_cname])
            except KeyError:
                for k in sorted(cc.cases.keys()):
                    print(k)
                print('-'*79)
                print(cc.sim_id, _cname)
                print('-'*79)
                raise KeyError
            if labels is not None:
                label = labels[i]
            else:
                label = cc.sim_id
            fp.plot(res, label)
        dlcf = 'dlc' + cc.cases[_cname]['[DLC]']
        fig_path = os.path.join(fig_dir_base, dlcf)
        fp.final(fig_path, _cname.replace('.htc', '.png'))

def plot_staircase(sim_ids, post_dirs, run_dirs, fig_dir_base=None,
                   cname='dlc00_stair_wsp04_25_noturb.htc'):
    """
    Default stair and ramp names:

    dlc00_stair_wsp04_25_noturb
    dlc00_ramp_wsp04_25_04_noturb
    """

    stairs = []

    col = ['r', 'k']
    alf = [1.0, 0.7]

    # if sim_id is a list, combine the two dataframes into one
    if type(sim_ids).__name__ == 'list':
        for ii, sim_id in enumerate(sim_ids):
            if isinstance(post_dirs, list):
                post_dir = post_dirs[ii]
            else:
                post_dir = post_dirs
            stairs.append(sim.Cases(post_dir, sim_id, rem_failed=True))
    else:
        sim_id = sim_ids
        sim_ids = [sim_id]
        post_dir = post_dirs
        stairs.append(sim.Cases(post_dir, sim_id, rem_failed=True))

    fig, axes = mplutils.make_fig(nrows=3, ncols=1, figsize=(14,10))
    ax = axes.ravel()

    for i, cc in enumerate(stairs):
        if cname in cc.cases_fail:
            print('no result for %s' % cc.sim_id)
            continue
        cc.change_results_dir(run_dirs[i])
        res = cc.load_result_file(cc.cases[cname])
        respath = cc.cases[cname]['[run_dir]']
        fname = os.path.join(respath, cname)
        df_respost = pd.read_hdf(fname + '_postres.h5', 'table')
        sim_id = cc.sim_id
        time = res.sig[:,0]
        t0, t1 = time[0], time[-1]

        # find the wind speed
        for channame, chan in res.ch_dict.iteritems():
            if channame.startswith('windspeed-global-Vy-0.00-0.00'):
                break
        wind = res.sig[:,chan['chi']]
        chi = res.ch_dict['bearing-pitch1-angle-deg']['chi']
        pitch = res.sig[:,chi]

        chi = res.ch_dict['bearing-shaft_rot-angle_speed-rpm']['chi']
        rpm = res.sig[:,chi]

        chi = res.ch_dict['bearing-pitch1-angle-deg']['chi']
        pitch = res.sig[:,chi]

        chi = res.ch_dict['tower-tower-node-001-momentvec-x']['chi']
        tx = res.sig[:,chi]

        chi = res.ch_dict['tower-tower-node-001-momentvec-y']['chi']
        ty = res.sig[:,chi]