text: import matplotlib as mpl
  from matplotlib import colors
  from matplotlib import pyplot as plt
  
  import numpy as np
  
  from pskf.tools.run import pythonmodule as pm
  from pskf.tools.plot import plotarrays as pa
  
  from pskf.scripts.numcomp import arrays as na
  
  
  ###############################################################################
  #             Plot bars of probabilites for numerical comparisons             #
  ###############################################################################
  def plot(
          ax,
          which_methods,
          which_methods_left,
          which_methods_right,
          which_res='endres',
          model='wavebc',
          method='meanc',
          ensemble_size=50,
          n_syn=1,
          n_syn_bold=1,
          n_comparisons=10000,
          pic_format='pdf',
          bar_colors=['black', 'white', 'grey']
  ):
  
      """
      Reads probability arrays which method is better,
      worse, or if they are even. Then plots those
      three possibilities in bars comparing the methods
      given in which_methods_left and which_methods_right.
  
      Parameters
      ----------
      ax : Axes
          The axes to draw to.
  
      which_methods : array int
          Array of integers containing the method specifiers
          from module plotarrays.
  
      which_methods_left : array int
          Array of integers containing the method specifiers
          for the left side of the comparisons.
  
      which_methods_right : array int
          Array of integers containing the method specifiers
          for the right side of the comparisons.
  
      which_res : string
          'endres' - use residuals after EnKF run
          'begres' - use residuals before EnKF run
  
      model : string
          'wavebc' - Model wavebc
          'wave' - Model wave
  
      method : string
          Which method to use for statistical comparison
          of the subset. If n_syn == 1, the comparison
          always defaults to comparing the residuals.
          'ttest' - Use the T-Test, testing if the
                    two samples belong to the same
                    Gaussian distribution.
          'gauss' - Calculate Gaussian distribution
                    of the difference and calculate
                    its probability to be larger
                    than zero.
          'meanc' - Calculate the means and compare.
  
      ensemble_size : integer
          Ensemble size of the job. Possibilities: 50,
          70, 100, 250, 500, 1000, 2000
  
      n_syn : integer
          Number of synthetic studies in subset.
  
      n_syn_bold : integer
          Number of synthetic studies in subset used for
          the bold ticklabels.
  
      n_comparisons : integer
          Number of comparisons calculated.
  
      pic_format : string
          Format of the picture
          'pdf' - pdf-format
          'eps' - eps-format
          'png' - png-format
          'jpg' - jpg-format
          'svg' - svg-format
  
      bar_colors : array of strings
          Three colors for the three patches of one bar.
  
      Returns
      -------
      ax : array
          Axes containing plot.
  
      pic_name : string
          Containing proposed saving location for Figure.
      """
  
      # Check
      for imethod in which_methods_left:
          if imethod not in which_methods:
              raise RuntimeError(
                  'Wrong methods in wrong_methods_left'
                  )
      for imethod in which_methods_left:
          if imethod not in which_methods:
              raise RuntimeError(
                  'Wrong methods in wrong_methods_right'
                  )
  
      if ensemble_size in [50, 70, 100, 250]:
          if n_syn > 1000:
              raise RuntimeError('n_syn wrong')
      elif ensemble_size in [500, 1000, 2000]:
          if n_syn > 100:
              raise RuntimeError('n_syn wrong')
      else:
          raise RuntimeError('ensemble size wrong')
  
      # Both methods in one array
      show_methods = [which_methods_left,
                      which_methods_right]
      # Number of bars and patches
      num_bars = len(show_methods[0])
      num_patches = 3*num_bars
  
      # Load probs
      probs = np.load(
          pm.py_output_filename(
              na.tag,
              'probs_'+which_res, model
              + '_'+method+'_'+str(ensemble_size)+'_'+str(n_syn)+'_'
              + str(n_comparisons)+'_'
              + '_'.join([str(i) for i in which_methods]),
              'npy'
          )
      )
  
      # Load probs for bold labels
      probs_bold = np.load(
          pm.py_output_filename(
              na.tagn,
              'probs_'+which_res,
              model+'_'
              + method+'_'+str(ensemble_size)+'_'+str(n_syn_bold)
              + '_'+str(n_comparisons)+'_'
              + '_'.join([str(i) for i in which_methods]),
              'npy'
          )
      )
  
      ax.set_position([0.3, 0.05, 0.4, 0.75])
      ax.set_frame_on(False)
  
      # Patch arrays for ax.barh()
      in_bottom = np.zeros(num_patches)
      in_height = np.zeros(num_patches)
      in_width = np.zeros(num_patches)
      in_left = np.zeros(num_patches)
      in_color = ['' for i in range(num_patches)]
      for i in range(num_patches):
          in_bottom[i] = num_bars-i/3
          in_height[i] = 0.8
          in_width[i] = probs[show_methods[0][i/3],
                              show_methods[1][i/3]][np.mod(i, 3)]
          in_left[i] = np.sum(probs[show_methods[0][i/3],
                                    show_methods[1][i/3]][0:np.mod(i, 3)])
          in_color[i] = bar_colors[np.mod(i, 3)]
  
      # Plot patches in bars
      ax.barh(bottom=in_bottom,
              height=in_height,
              width=in_width,
              left=in_left,
              color=in_color,
              edgecolor='k')
  
      # H_0 labels inside bar
      if method == "ttest":
          for i in range(1, num_bars+1):
              if in_left[3*i-1]-in_left[3*i-2] > 0.15:
                  ax.text(in_left[3*i-2]+0.4*(in_left[3*i-1]-in_left[3*i-2]),
                          in_bottom[3*i-2]+0.3,
                          "$H_0$",
                          fontsize=20)
  
      # Axis 1
      ax.tick_params(direction='out', length=0,
                     width=1, labelsize=20,
                     top=False, bottom=False,
                     labelright=False,
                     pad=8)
      ax.set_xlim([-0.01, 1.01])
      ax.set_ylim([0.9, num_bars+0.8])
      ax.set_xticks([])
      ax.set_yticks([num_bars-i+0.4 for i in range(num_bars)])
      ax.set_yticklabels(
          [pa.longnames_methods[show_methods[0][i]] for i in range(num_bars)]
          )
  
      # Twin Axis 2
      ax2 = ax.twinx()
      ax2.set_position([0.3, 0.05, 0.4, 0.75])
      ax2.set_frame_on(False)
      ax2.tick_params(direction='out', length=0,
                      width=1, labelsize=20,
                      top=False, bottom=False,
                      labelleft=False, labelright=True,
                      labelcolor='black',
                      pad=8)
      ax2.set_xlim([-0.01, 1.01])
      ax2.set_ylim([0.9, num_bars+0.8])
      ax2.set_xticks([])
      ax2.set_yticks([num_bars-i+0.4 for i in range(num_bars)])
      ax2.set_yticklabels(
          [pa.longnames_methods[show_methods[1][i]] for i in range(num_bars)]
          )
  
      # Boldness of axislabels
      for i in range(num_bars):
          if(probs_bold[show_methods[0][i], show_methods[1][i]][0] == 1):
              ax.yaxis.get_majorticklabels()[i].set_weight('bold')
          elif(probs_bold[show_methods[0][i], show_methods[1][i]][2] == 1):
              ax2.yaxis.get_majorticklabels()[i].set_weight('bold')
          else:
              ax.yaxis.get_majorticklabels()[i].set_style('italic')
              ax2.yaxis.get_majorticklabels()[i].set_style('italic')
  
      # Saving location
      pic_name = pm.py_output_filename(
          na.tag,
          'bars_'+which_res,
          model+'_'+method+'_'+str(ensemble_size)
          + '_'+str(n_syn)+str(n_comparisons)+'_'+'_'
          + '_'.join([str(i) for i in which_methods]),
          pic_format
          )
  
      return ax, pic_name
  
  
  ###############################################################################
  #           Plot matrix of probabilties for a numerical comparisons           #
  ###############################################################################
  def matrix(
          ax,
          which_methods=[0, 1, 2, 3, 4, 5, 6],
          which_res='endres',
          method='meanc',
          model='wavebc',
          ensemble_size=50,
          n_syn=1,  # number of synthetic studies
          n_comparisons=10000,
          is_text=False,
          pic_format='pdf',  # 'png' or 'eps' or 'svg' or 'pdf'
          figpos=[0.14, 0.14, 0.8, 0.8],  # xbeg, ybeg, xrange, yrange
          # ylims=[0.28,0.82],
          is_longnames=True,
          ticklabelfontsize=10,
          xtick_y=0.0,
          is_color_trafo=False,
          # num_pack=4,                     # Number of methods in pack
          # formatsos=['o','v','s','p','o','v','s','p'],
          # coleros=[(0.0,0.0,0.0),(0.0,0.0,0.0),(0.0,0.0,0.0),(0.0,0.0,0.0),
          #                (1.0,1.0,1.0),(1.0,1.0,1.0),(1.0,1.0,1.0),(1.0,1.0,1.0)],
          # markersize=10,
          # markeredgesize=1.5,
          # fontleg=30,                              #18
          # fonttit=40,
          # fontlab=40,
          # fonttic=30,
  ):
      """
      A plotting function for statistics of residual distributions.
  
      Parameters
      ----------
      ax : Axes
          The axes to draw to.
  
      which_methods : array of ints
          The methods to be loaded in ascending order.
  
      which_res : string
          'endres' - use residuals after EnKF run
          'begres' - use residuals before EnKF run
  
      method : string
          Which method to use for statistical comparison
          of the subset. If n_syn == 1, the comparison
          always defaults to comparing the residuals.
          'ttest' - Use the T-Test, testing if the
                    two samples belong to the same
                    Gaussian distribution.
          'gauss' - Calculate Gaussian distribution
                    of the difference and calculate
                    its probability to be larger
                    than zero.
          'meanc' - Calculate the means and compare.
  
      model : string
          'wavebc' - Model wavebc
          'wave' - Model wave
  
      ensemble_size : integer
          Ensemble size of the job. Possibilities: 50,
          70, 100, 250, 500, 1000, 2000
  
      n_syn : integer
          Number of synthetic studies in subset.
  
      n_comparisons : integer
          Number of comparisons calculated.
  
      pic_format : string
          Format of the picture
          'pdf' - pdf-format
          'eps' - eps-format
          'png' - png-format
          'jpg' - jpg-format
          'svg' - svg-format
  
      Returns
      -------
      ax : Axes
          Axes containing plot.
  
      pic_name : string
          Containing proposed saving location for Figure.
      """
  
      # Check
      if ensemble_size in [50, 70, 100, 250]:
          if n_syn > 1000:
              raise RuntimeError('n_syn wrong')
      elif ensemble_size in [500, 1000, 2000]:
          if n_syn > 100:
              raise RuntimeError('n_syn wrong')
      else:
          raise RuntimeError('ensemble size wrong')
  
      # Number of compared methods
      num_methods = len(which_methods)
  
      # Load probs
      probs = np.load(pm.py_output_filename(
          na.tag,
          'probs_'+which_res,
          model+'_'+method+'_'+str(ensemble_size)
          + '_'+str(n_syn)+'_'+str(n_comparisons)+'_'
          + '_'.join([str(i) for i in which_methods]),
          'npy'
          ))
  
      ax.set_position(figpos)
  
      # Rectangles in upper right half: Fraction of Undecided
      undecided = probs[:, :, 1]
      for ipm in range(num_methods):
          for jpm in range(num_methods):
              if ipm > jpm:
                  # Lower left half white
                  undecided[ipm, jpm] = None
              if ipm == jpm:
                  # Diagonal black
                  undecided[ipm, jpm] = 1.0
              if ipm < jpm:
                  # Single comparisons white
                  if n_syn == 1:
                      undecided[ipm, jpm] = None
                  # One comparison white
                  if n_syn == 1000:
                      undecided[ipm, jpm] = None
                  # For mean comparison white
                  if method == 'meanc':
                      undecided[ipm, jpm] = None
  
      ax.imshow(undecided, interpolation='nearest', cmap='Greys',
                norm=colors.Normalize(vmin=0, vmax=1, clip=False))
  
      # Triangles: Grid
      X, Y = np.meshgrid(np.arange(num_methods+1), np.arange(num_methods+1))
      X = X.flatten()-0.5
      Y = Y.flatten()-0.5
  
      # Triangles: Indices
      triangles = np.zeros([2*np.sum(range(1, num_methods)), 3])
  
      ix = 0
      # Upper triangles
      for i in range(1, num_methods):
          for j in range(i, num_methods):
              triangles[ix, :] = [(i-1)+j*(num_methods+1),
                                  i+j*(num_methods+1),
                                  i+num_methods+j*(num_methods+1)]
              ix = ix+1
      # Lower triangles
      for i in range(1, num_methods):
          for j in range(i, num_methods):
              triangles[ix, :] = [i+j*(num_methods+1),
                                  i+num_methods+j*(num_methods+1),
                                  i+num_methods+1+j*(num_methods+1)]
              ix = ix+1
  
      # Triangles: Triangulation instance
      tria = mpl.tri.Triangulation(X, Y, triangles)
  
      # Triangles: Colors
      coleros = np.array([[probs[i, j, 0] for i in range(j+1, num_methods)]
                          for j in range(num_methods-1)]
                         + [[probs[i, j, 2] for i in range(j+1, num_methods)]
                            for j in range(num_methods-1)])
      coleros = np.hstack(coleros)
  
      # Color Transformation
      if is_color_trafo:
          coleros = 0.5*(coleros-0.5)+0.5
          coleros[coleros < 0.275] = 0.0
          coleros[coleros > 0.725] = 1.0
  
      # Triangles: Plot with facecolor
      plt.tripcolor(
          tria,
          facecolors=coleros,
          cmap=mpl.cm.Greys,
          norm=colors.Normalize(vmin=0, vmax=1, clip=False),
          edgecolor='k'
      )
  
      # Plot: Mostly ticks
      ticklabelinput = ([pa.longnames_methods[which_methods[i]]
                         for i in range(num_methods)]
                        if is_longnames else
                        [pa.names_methods[which_methods[i]]
                         for i in range(num_methods)])
      ax.set_xticks([i for i in range(num_methods)])
      ax.set_xticklabels(ticklabelinput, fontsize=ticklabelfontsize,
                         rotation=90, y=xtick_y)
      ax.set_yticks([i for i in range(num_methods)])
      ax.set_yticklabels(ticklabelinput, fontsize=ticklabelfontsize)
      ax.tick_params(length=0)
      ax.set_frame_on(False)
  
      # Text: Upper triangles
      for i in range(3):
          for itext in range(num_methods):
              for jtext in range(num_methods):
                  if itext < jtext:
                      ntext = (np.around(100*probs[jtext, itext, i],
                                         decimals=1)
                               if i != 1 else
                               np.around(100*probs[itext, jtext, i],
                                         decimals=1))
                      ttext = (str(ntext)[0:4]
                               if 0 < ntext < 100 else
                               str(ntext)[0:0])
                      px = (itext-0.45
                            if i == 0 else
                            (jtext-0.125 if i == 1 else itext-0.10))
                      py = (jtext-0.15
                            if i == 0 else
                            (itext+0.05 if i == 1 else jtext+0.3))
                      colero = 'white' if ntext > 50 else 'black'
  
                      if i != 1 or (n_syn != 1000
                                    and n_syn != 1
                                    and method != "meanc"):
                          ax.text(px, py, ttext, color=colero, fontsize=20)
  
      # Text: n_syn and ensemble_size
      if is_text:
          model_spec = '  Tracer  ' if model == 'wavereal' else '   Well   '
          ax.text(
              3.5, 1.5,
              model_spec+'\n'
              + r'  $n_{e}$: '+str(ensemble_size).rjust(4)+'\n'
              + r' $n_{syn}$: '+str(n_syn).rjust(4),
              linespacing=1.5,
              fontsize=30,
              bbox={'facecolor': 'grey', 'alpha': 0.5, 'pad': 10},
          )
  
      # Saving location
      pic_name = pm.py_output_filename(
          na.tag,
          'matrix_'+which_res,
          model+'_'+method+'_'+str(ensemble_size)
          + '_'+str(n_syn)+'_'+str(n_comparisons)+'_'
          + '_'.join([str(i) for i in which_methods]),
          pic_format
      )
  
      return ax, pic_name
  