#!/usr/bin/python
# -*- coding: ascii -*-
#
from __future__ import print_function, division
import os, sys
from os import path
import numpy as np
import CoolProp
import glob
from warnings import warn
from time import clock

import CoolProp.constants
from CoolProp.CoolProp import PropsSI,generate_update_pair,get_parameter_index,set_debug_level,get_phase_index
from CoolProp import AbstractState as State

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import copy
from itertools import cycle
from matplotlib import gridspec, ticker
from jopy.dataPlotters import roundList, range_brace

try:
    from jopy.dataPlotters import BasePlotter
    bp = BasePlotter()
except:
    bp = None


# The basic settings for he plots
xypoints = 1000
loops    = 1
repeat   = 1
runs     = 0
maxruns  = 100
plot     = True
calc     = False
check    = False
folder   = "dataDESKTOP"
figures  = "figures"
#np.random.seed(1984)


fluids  = ["CO2","Pentane","R134a","Water","Air","LiBr-0%"]
#glskeys = [r"\glsentryshort{co2}",r"\glsentryshort{pentane}",r"\glsentryshort{r134a}",r"\glsentryshort{water}",r"\glsentryshort{air}",r"\glsentryshort{libr} \SI{0}{\percent}"]
#glskeys = [r"\ce{CO2}",r"n-Ppentane",r"R134a",r"Water",r"Air",r"\glsentryshort{libr} \SI{0}{\percent}"]
repList = []
#for i in range(len(fluids)):
#    repList.append(fluids[i])
#    repList.append(glskeys[i])

backends = ["INCOMP","HEOS","REFPROP"]
#repList.append("HEOS")
#repList.append(r"\glsentryshort{cp}")
#repList.append("REFPROP")
#repList.append(r"\glsentryshort{rp}")


#CoolProp.CoolProp.set_debug_level(51)

pStr = path.dirname(path.abspath(__file__))
fStr = path.splitext(path.basename(__file__))[0]
def getFolderName():
    folderName = path.join(pStr,folder)
    if not path.isdir(folderName):
        print("Creating data directory "+folderName)
        os.makedirs(folderName)
    return folderName

def getFigureFolder():
    folderName = path.join(pStr,figures)
    if not path.isdir(folderName):
        print("Creating data directory "+folderName)
        os.makedirs(folderName)
    return folderName

repList.append("TimeComp-")
repList.append("chapters/FluidProperties/"+path.basename(getFigureFolder())+"/TimeComp-")

def getFileName(qualifiers=[]):
    fileName = path.join(getFolderName(),"-".join(qualifiers))
    return fileName

# Some file handling
def loadNpzData(backend,fluid):
    dicts = {}
    globber = getFileName([backend,fluid])+'_[0-9][0-9][0-9].npz'
    for fname in glob.glob(globber):
        dataDict = dict(np.load(fname))
        dicts[str(dataDict["name"])] = dataDict
    #if len(dicts)<1:
    #    #print("No readable file found for {0}".format(globber))
    #    dataDict = dict(name=str(0).zfill(3))
    #    dicts[str(dataDict["name"])] = dataDict
    return dicts

def saveNpzData(backend,fluid,dicts,start=0,stop=-1):
    keys = dicts.keys()
    keys.sort()
    for k in keys[start:stop]:
        data = dicts[k]
        fname = getFileName([backend,fluid])+'_{0}.npz'.format(str(data['name']).zfill(3))
        np.savez(fname, **data)
    return True

def splitFluid(propsfluid):
    fld = propsfluid.split("::")
    if len(fld)==2:
        backend =  fld[0]
        fld = fld[1]
    else:
        backend = None
        fld = fld[0]
    fld = fld.split("-")
    if len(fld)==2:
        conc  = float(fld[1].strip('%'))/100.0
        fld = fld[0]
    else:
        conc  = None
        fld = fld[0]
    return backend,fld,conc

def getInpList(backend):
    if backend=="HEOS": return ["DT","HP"]
    elif backend=="REFPROP": return ["DT","HP"]
    elif backend=="INCOMP": return ["PT","HP"]
    else: raise ValueError("Unknown backend.")

def getOutList(inp=None):
    if inp == "HP":
        return [["Tmax"],["D"],["S"],["T"],["D","S","T"]]
    elif inp == "DT":
        return [["Tmax"],["H"],["P"],["S"],["H","P","S"]]
    elif inp == "PT":
        return [["Tmax"],["H"],["D"],["S"],["H","D","S"]]
    else:
        raise ValueError("Unknown inputs.")

#def getPhaseString(iPhase):
#    for i in range(11):
#        if getPhaseString(i)==sPhase:
#            return i
#     if iPhase==1: return "liquid"
#     elif iPhase==2: return "supercritical"
#     elif iPhase==3: return "supercritical_gas"
#     elif iPhase==4: return "supercritical_liquid"
#     elif iPhase==5: return "critical_point"
#     elif iPhase==6: return "gas"
#     elif iPhase==7: return "twophase"
#     elif iPhase==8: return "unknown"
#     elif iPhase==9: return "not_imposed"
#     else: raise ValueError("Couldn't find phase.")

def getPhaseNum(sPhase):
    return get_phase_index(sPhase)
#     for i in range(11):
#         if getPhaseString(i)==sPhase:
#             return i

def getOutKey(out): return "".join(out)

def getOutLabel(out): return ",".join(out)

def getTimeKey(inp,out): return "_".join([inp,getOutKey(out)])

def getVectorKey(inp,out): return getTimeKey(inp, out)+"_V"

def getCriticalProps(propsfluid):
    backend,_,_ = splitFluid(propsfluid)
    if backend!="INCOMP":
        p_crit_m = PropsSI('pcrit',"T",0,"D",0,propsfluid)  *0.995
        T_crit_m = PropsSI('Tcrit',"T",0,"D",0,propsfluid)  *1.005
        d_crit_m = PropsSI('rhocrit',"T",0,"D",0,propsfluid)*0.995
        h_crit_m = PropsSI('H',"T",T_crit_m,"D",d_crit_m,propsfluid)
        s_crit_m = PropsSI('H',"T",T_crit_m,"D",d_crit_m,propsfluid)
    else:
        p_crit_m = None
        T_crit_m = None
        d_crit_m = None
        h_crit_m = None
        s_crit_m = None
    return dict(P=p_crit_m,T=T_crit_m,D=d_crit_m,H=h_crit_m,S=s_crit_m)


def getPTRanges(propsfluid):
    backend,_,_ = splitFluid(propsfluid)
    # Setting the limits for enthalpy and pressure
    if backend == "REFPROP":
        T_min = PropsSI('Ttriple',"T",0,"D",0,propsfluid)
        T_max = PropsSI(   'Tmax',"T",0,"D",0,propsfluid)
        p_min = PropsSI(  'P',"T",T_min,"Q",0,propsfluid)
        p_max = PropsSI(   'pmax',"T",0,"D",0,propsfluid)
    elif backend == "INCOMP":
        T_min = PropsSI('Tmin',"T",0,"D",0,propsfluid)
        T_max = PropsSI('Tmax',"T",0,"D",0,propsfluid)
        p_min =   1.5*1e5
        p_max = 200.0*1e5
    else:
        p_min = PropsSI('ptriple',"T",0,"D",0,propsfluid)
        p_max = PropsSI(   'pmax',"T",0,"D",0,propsfluid)
        T_min = PropsSI('Ttriple',"T",0,"D",0,propsfluid)
        T_max = PropsSI(   'Tmax',"T",0,"D",0,propsfluid)

    p_range  = np.logspace(np.log10(p_min),np.log10(p_max),xypoints)
    T_range  = np.linspace(T_min,T_max,xypoints)
    return p_range,T_range

    #p_max = min(PropsSI('pcrit',"T",0,"D",0,fluid)*20, p_max)
    #T_max = min(PropsSI('Tcrit',"T",0,"D",0,fluid)* 3, T_max)


def getLists(propsfluid):
    backend,_,_ = splitFluid(propsfluid)
    """Returns randomised lists of all properties within the ranges"""
    p,T = getPTRanges(propsfluid)
    p_min = np.min(p)
    p_max = np.max(p)
    T_min = np.min(T)
    T_max = np.max(T)

    if backend=="INCOMP":
        h_min = PropsSI('H','T',T_min,'P',p_min,propsfluid)
        h_max = PropsSI('H','T',T_max,'P',p_max,propsfluid)
    else:
        critProps = getCriticalProps(propsfluid)
        h_min = PropsSI('H','T',T_min,'Q',0,propsfluid)
        h_max = PropsSI('H','T',T_min,'Q',1,propsfluid)
        h_max = max(PropsSI('H','T',critProps["T"],'D',critProps["D"],propsfluid),h_max)
        h_max = (h_max-h_min)*2.0+h_min

    loop = True
    count = 0
    while loop:
        count += 1

        h_list = np.random.uniform(          h_min,          h_max,int(xypoints*2.0))
        p_list = np.random.uniform(np.log10(p_min),np.log10(p_max),int(xypoints*2.0))
        p_list = np.power(10,p_list)

        out = ["T","D","S"]
        res = PropsSI(out,"P",p_list,"H",h_list,propsfluid)
        T_list = res[:,0]
        d_list = res[:,1]
        s_list = res[:,2]

        mask = np.isfinite(T_list) & np.isfinite(d_list) & np.isfinite(s_list)
        if np.sum(mask)<xypoints:
            if False:
                print(h_list); print(p_list); print(T_list); print(d_list); print(s_list)
            print("There were not enough valid entries in your result vector: {0:d} > {1:d} - rerunning".format(xypoints,np.sum(mask)))
            loop = True
        else:
            loop = False
            p_list = p_list[mask][0:xypoints]
            h_list = h_list[mask][0:xypoints]
            T_list = T_list[mask][0:xypoints]
            d_list = d_list[mask][0:xypoints]
            s_list = s_list[mask][0:xypoints]
            return dict(P=p_list,T=T_list,D=d_list,H=h_list,S=s_list)

        maxTries = 4
        if count > maxTries:
            loop = False
            raise ValueError("{0}: Could not fill the lists in {0} runs, aborting.".format(propsfluid,maxTries))

def getInpValues(inp,dataDict):
    in1 = inp[0]
    in2 = dataDict[in1]
    in3 = inp[1]
    in4 = dataDict[in3]
    return in1,in2,in3,in4

def getStateObj(propsfluid):
    backend,fld,conc = splitFluid(propsfluid)
    # fluidstr holds the full information and fluid is only the name
    # Initialise the state object
    if backend is not None:
        state = State(backend,fld)
    else:
        state = State(fld)
    #if backend=="INCOMP":
    #    state.set_mass_fractions([0.0])
    if conc is not None:
        state.set_mass_fractions([conc])
    return state

def getSpeedMeas(out,in1,in2,in3,in4,propsfluid,vector=False):
    pair, out1, _ = generate_update_pair(get_parameter_index(in1),in2[0],get_parameter_index(in3),in4[0])
    if out1==in2[0]: swap=False
    else: swap=True
    if swap:
        input1 = in4
        input2 = in2
    else:
        input1 = in2
        input2 = in4
    state = getStateObj(propsfluid)
    outList = [get_parameter_index(i) for i in out]
    resLst = np.empty((repeat,))
    timLst = np.empty((repeat,))
    if vector:
        for j in range(repeat):
            timLst.fill(np.inf)
            lrange = range(len(input1))
            resTmp = np.inf
            if len(outList)==1 and outList[0]==CoolProp.constants.iT_max:
                t1=clock()
                for l in lrange:
                    for o in outList:
                        resTmp = state.keyed_output(o)
                t2=clock()
                timLst[j] = (t2-t1)*1e6/float(len(input1))
            else: # We have to update before doing other things
                t1=clock()
                for l in lrange:
                    state.update(pair,input1[l],input2[l])
                    for o in outList:
                        resTmp = state.keyed_output(o)
                t2=clock()
                timLst[j] = (t2-t1)*1e6/float(len(input1))
        res = None
        tim = np.min(timLst) # Best of (repeat)
        return res,tim
    else:
        res = np.empty_like(input1)
        res.fill(np.inf)
        tim = np.empty_like(input1)
        tim.fill(np.inf)
        for i in range(len(input1)):
            resLst.fill(np.inf)
            timLst.fill(np.inf)
            for j in range(repeat):
                lrange = range(loops)
                resTmp = np.inf
                if len(outList)==1 and outList[0]==CoolProp.constants.iT_max:
                    t1=clock()
                    for _ in lrange:
                        for o in outList:
                            resTmp = state.keyed_output(o)
                    t2=clock()
                    timLst[j] = (t2-t1)*1e6/float(loops)
                    resLst[j] = resTmp
                else: # We have to update before doing other things
                    inV1 = input1[i]
                    inV2 = input2[i]#*(1.0+(l/1000.0)*pow(-1,l)) for l in lrange ]
                    t1=clock()
                    for l in lrange:
                        state.update(pair,inV1,inV2)
                        for o in outList:
                            resTmp = state.keyed_output(o)
                    t2=clock()
                    timLst[j] = (t2-t1)*1e6/float(loops)
                    resLst[j] = resTmp
            if not np.all(resLst==resLst[0]):
                raise ValueError("Not all results were the same.")
            res[i] = resLst[0]
            tim[i] = np.min(timLst) # Best of three (repeat)
        return res,tim


def checkDataSet(propsfluid,dataDict,fill=True,quiet=False):
    if not check: return
    backend,_,_ = splitFluid(propsfluid)
    if not quiet: print("\n\n-- {0:16s} --".format(propsfluid),end="")
    # Test for required inputs
    newInputs = False
    inLists = getLists(propsfluid)
    for inp in getInpList(backend):
        if not quiet: print("\n{0:2s}: ".format(inp),end="")
        for inVal in inp:
            if inVal not in dataDict: # A problem
                if not fill:
                    raise ValueError("The input {0:1s} is missing or faulty, cannot continue.".format(inVal))
                #dataDict[inVal] = inLists[inVal]
                dataDict.update(inLists)
                newInputs = True
                if not quiet: print("{0:s}*({1:d}),".format(inVal,len(dataDict[inVal])),end="")
            else:
                if not quiet: print("{0:s} ({1:d}),".format(inVal,len(dataDict[inVal])),end="")
        # All inputs are there
        in1,in2,in3,in4 = getInpValues(inp,dataDict)
        #in2 = in2[:3]
        #in4 = in4[:3]
        if in2.shape != in4.shape:
            raise ValueError("The stored data for {0:s} and {1:s} do not have the same shape.".format(in1,in3))
        if in2.shape != inLists[in1].shape:
            raise ValueError("The stored data for {0:s} and its list do not have the same shape {1} vs {2}.".format(in1,in2.shape,inLists[in1].shape))
        # Check for time data
        for out in getOutList(inp):
            key = getTimeKey(inp, out)
            okey = getOutKey(out)
            if key not in dataDict or newInputs or not np.all(np.isfinite(dataDict[key])):
                if not fill:
                    raise ValueError("The time data for {0:s} is missing or faulty, cannot continue.".format(key))
                res,tim = getSpeedMeas(out, in1, in2, in3, in4, propsfluid)
                dataDict[key] = tim
                dataDict[okey] = res # We calculated in, why not use it here...
                if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="")
            else:
                if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="")
            if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])):
                raise ValueError("The stored time data for {0:s} does not have the same shape as the inputs.".format(key))
        # Check for vectors
        for out in getOutList(inp):
            key = getVectorKey(inp, out)
            if key not in dataDict or not np.all(np.isfinite(dataDict[key])):
                if not fill:
                    raise ValueError("The fluid data for {0:s} is missing or faulty, cannot continue.".format(key))
                res,tim = getSpeedMeas(out, in1, in2, in3, in4, propsfluid, vector=True)
                dataDict[key] = tim
                if not quiet: print("{0:s}*({1:d}),".format(key,dataDict[key].size),end="")
            else:
                if not quiet: print("{0:s} ({1:d}),".format(key,dataDict[key].size),end="")
            if dataDict[key].size!=1 or not np.all(np.isfinite(dataDict[key])):
                raise ValueError("The vector data for {0:s} does not have the correct size {1}..".format(key,dataDict[key].size))

#     inp = getInpList(backend)[0] # Explicit calls
#      # Check for properties
#      for out in getOutList(inp)[:-1]:
#          key = getOutKey(out)
#          if key not in dataDict or not np.all(np.isfinite(dataDict[key])):
#              if not fill:
#                  raise ValueError("The fluid data for {0:s} is missing or faulty, cannot continue.".format(key))
#              res = PropsSI(out,in1,in2,in3,in4,propsfluid)
#              dataDict[key] = res
#              if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="")
#          else:
#              if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="")
#          if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])):
#              raise ValueError("The stored data for {0:s} does not have the same shape as the inputs {1} vs {2}..".format(key,dataDict[key].shape,in2.shape))
#     # Check for phase
#     for out in ["Phase"]:
#         if backend!="HEOS":
#             dataDict[key] = np.zeros_like(a, dtype, order, subok)
#         key = getOutKey(out)
#         if key not in dataDict or newInputs or not np.all(np.isfinite(dataDict[key])):
#             if not fill:
#                 raise ValueError("The phase data for {0:s} is missing or faulty, cannot continue.".format(key))
#             res = np.empty_like(in2)
#             res.fill(np.inf)
#             for i in range(len(in2)):
#                 res[i] = PropsSI(out,in1,in2[i],in3,in4[i],propsfluid)
#             dataDict[key] = res
#             if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="")
#         else:
#             if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="")
#         if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])):
#             raise ValueError("The stored data for {0:s} does not have the same shape as the inputs {1} vs {2}..".format(key,dataDict[key].shape,in2.shape))
#
#             # Now we use the vector data
#             key = getVectorKey(inp, out)
#             if key not in dataDict or not np.all(np.isfinite(dataDict[key])):
#                 if not fill:
#                     raise ValueError("The vector data for {0:s} is missing or faulty, cannot continue.".format(key))
#                 dataDict[key] = np.empty_like(in2)
#                 dataDict[key].fill(np.inf)
#                 res = []
#                 for _ in range(repeat):
#                     t1=clock()
#                     PropsSI(out,in1,in2,in3,in4,propsfluid)
#                     t2=clock()
#                     res.append((t2-t1)/float(len(in2)))
#                 dataDict[key] = np.min(res)
#                 if not quiet: print("{0:s}*({1}),".format(key,dataDict[key]),end="")
#             else:
#                 if not quiet: print("{0:s} ({1}),".format(key,dataDict[key]),end="")
#             try:
#                 float(dataDict[key])
#             except:
#                 raise ValueError("The stored vector data for {0:s} cannot be casted to float.".format(key))
#        if not quiet: print("")

# All data is loaded and checked, we can calculate more now
def getEntryCount(dicts,backend,fld):
    return len(fluidData[fld][backend].keys())

def getUKey(fld,bck,inp,out):
    return "-".join([fld,bck,inp,"".join(out)])

def getData(fld,backend,inp,out,fluidData):
    inputs1 = []
    inputs2 = []
    values  = []
    times   = []
    i1key = inp[0]
    i2key = inp[1]
    vkey  = getOutKey(out)
    tkey  = getTimeKey(inp,out)
    for dkey in fluidData[fld][backend]:
        cData = fluidData[fld][backend][dkey]
        inputs1.append(cData[i1key])
        inputs2.append(cData[i2key])
        values.append( cData[vkey] )
        times.append(  cData[tkey] )
    ret = {}
    if len(inputs1)>0:
        ret[i1key] = np.concatenate(inputs1)
        ret[i2key] = np.concatenate(inputs2)
        ret[vkey]  = np.concatenate(values )
        ret[tkey]  = np.concatenate(times  )
    return ret

def getSingleData(fld,backend,key,fluidData):
    #print("Getting: "+fld+", "+backend+", "+key)
    values  = []
    for dkey in fluidData[fld][backend]:
        if key in fluidData[fld][backend][dkey]:
            if "P" in fluidData[fld][backend][dkey]:
                # TODO: Fix this, do we need the mask?
                #mask = fluidData[fld][backend][dkey]["P"]>0.3e5
                mask = fluidData[fld][backend][dkey]["P"]>0.3e5
                try:
                    values.append( fluidData[fld][backend][dkey][key][mask] )
                except Exception as e:
                    values.append( fluidData[fld][backend][dkey][key] )
                    print(e)
                    pass
            else:
                values.append( fluidData[fld][backend][dkey][key] )
    if len(values)>0:
        if np.size(values[0])>1:
            return np.concatenate(values)
        else:
            return np.array(values)
    return None

def fillDict(fld,backend,fluidData,curDict,curKeys):
    if curDict is None: curDict = {}
    for key in curKeys:
        vals = getSingleData(fld, backend, key, fluidData)
        if vals is not None: curDict[key] = vals
    return curDict

################################################

fluidData = {}

for fluidstr in fluids:
    _,fld,_ = splitFluid(fluidstr)
    if fld not in fluidData: fluidData[fld] = {}
    for backend in backends:
        if backend not in fluidData[fld]: # Try to add it
            propsfluid = "::".join([backend,fluidstr])
            try:
                PropsSI('Tmax',"T",0,"D",0,propsfluid)
                fluidData[fld][backend] = loadNpzData(backend,fld)
            except:
                pass
            if backend in fluidData[fld]:
                for k in fluidData[fld][backend]:
                    checkDataSet(propsfluid, fluidData[fld][backend][k], fill=False, quiet=True)
        else: # Already added backend for fluid
            pass

lastSave = 0
while runs<maxruns and calc:
    check = True # force checking for new records
    runs +=1
    # Now we have the data structure with the precalculated data
    for fluidstr in fluids:
        _,fld,_ = splitFluid(fluidstr)
        for backend in fluidData[fld]:
            propsfluid = "::".join([backend,fluidstr])
            dicts      = fluidData[fld][backend]
            keys = list(dicts.keys())
            keys.sort()
            while len(keys)<runs:
                if len(keys)<1: newKey = 0
                else: newKey = int(keys[-1])+1
                if newKey>999:
                    raise ValueError("Too many dicts: {0}>999".format(newKey))
                k = str(newKey).zfill(3)
                dicts[k] = {}
                dicts[k]['name']=k
                try:
                    checkDataSet(propsfluid, dicts[k], fill=True, quiet=False)
                except Exception as e:
                    print("There was an error, dataset {0} from {1}.might be faulty:\n{2}".format(k,propsfluid,str(e)))
                    pass
                todel = []
                #for k in dicts:
                try:
                    checkDataSet(propsfluid, dicts[k], fill=False, quiet=True)
                except Exception as e:
                    print("There was an error, removing dataset {0} from {1}.\n{2}".format(k,propsfluid,str(e)))
                    todel.append(k)
                for t in todel: del dicts[t]
                keys = list(dicts.keys())
                keys.sort()

    # Updated all dicts for this backend, saving data
    if runs>=maxruns or (lastSave+4)<runs:
        print("\n\nDone processing fluids, saving data: ")
        for fld in fluidData:
            for backend in fluidData[fld]:
                saveNpzData(backend, fld, fluidData[fld][backend], start=lastSave, stop=runs)
                print("{0} ({1})".format(backend+"::"+fld,len(fluidData[fld][backend].keys()[lastSave:runs])),end=", ")
        print("")
        lastSave = runs


if not plot: sys.exit(0)

# exclLst = [["Tmax"],["H"],["D"],["S"],["H","D","S"],["D"],["S"],["T"],["D","S","T"]]
#
### Start with a temporary dictionary that holds all the data we need
# for fld in fluidData:
#     cstData = {} # Data from calls to constants (overhead)
#     expData = {} # Data from explicit EOS calls
#     impData = {} # Data from EOS calls that require iterations
#     for backend in fluidData[fld]:
#
#         curKeys = []
#         for inp in getInpList(backend):
#             for out in getOutList(inp)[:1]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         cstData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[:1]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         expData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[1:]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         impData[backend] = curDict


#         curDict[backend] = {}
#         for key in :
#             vals = getSingleData(fld, backend, key, fluidData)
#             if vals is not None: curDict[backend][key] = vals
#             if curDict
#
#         cstData[backend] = {}
#         for out in getOutList(inp[0])[0]:
#             res = getData(fld,backend,inp[0],out,fluidData)
#             cstData[backend].update(res)
#         for out in getOutList(inp[1])[0]:
#             res = getData(fld,backend,inp[1],out,fluidData)
#             cstData[backend].update
#
#         expData[backend] = {}
#         for out in getOutList(inp[0])[1:]:
#             res = getData(fld,backend,inp[0],out,fluidData)
#             expData[backend].update(res)
#
#         impData[backend] = {}
#         for out in getOutList(inp[1])[1:]:
#             res = getData(fld,backend,inp[1],out,fluidData)
#             impData[backend].update(res)

    #############################################################
    # All data is available in the dicts now.
    #############################################################
    # The first thuing to do is to print some statistical
    # measures to give you an idea about the data.
#     try:
#         #dataOHCP = [cstData["HEOS"]["DT_Tmax"]   , cstData["HEOS"]["HP_Tmax"]   ]
#         #dataOHRP = [cstData["REFPROP"]["DT_Tmax"], cstData["REFPROP"]["HP_Tmax"]]
#         print("\n{0} - {1} points ".format(fld,np.size(cstData["HEOS"]["DT_Tmax"])))
#         print("Overhead CoolProp: {0:5.3f} us".format(np.mean(cstData["HEOS"]["DT_Tmax"])))#,   np.mean(cstData["HEOS"]["HP_Tmax"]))
#         print("Overhead REFPROP : {0:5.3f} us".format(np.mean(cstData["REFPROP"]["DT_Tmax"])))#,np.mean(cstData["REFPROP"]["HP_Tmax"]))
#         print("Mean EOS in  CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["HEOS"]["DT_HPS"]),np.mean(impData["HEOS"]["HP_DST"])))
#         print("Std. dev. in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["HEOS"]["DT_HPS"]) ,np.std(impData["HEOS"]["HP_DST"])))
#         print("Minimum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["HEOS"]["DT_HPS"]) ,np.min(impData["HEOS"]["HP_DST"])))
#         print("Maximum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["HEOS"]["DT_HPS"]) ,np.max(impData["HEOS"]["HP_DST"])))
#         print("Mean EOS in   REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["REFPROP"]["DT_HPS"]),np.mean(impData["REFPROP"]["HP_DST"])))
#         print("Std. dev. in  REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["REFPROP"]["DT_HPS"]) ,np.std(impData["REFPROP"]["HP_DST"])))
#         print("Minimum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["REFPROP"]["DT_HPS"]) ,np.min(impData["REFPROP"]["HP_DST"])))
#         print("Maximum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["REFPROP"]["DT_HPS"]) ,np.max(impData["REFPROP"]["HP_DST"])))
#         print("")
#
#         print("\n{0} - {1} points ".format(fld,np.size(cstData["HEOS"]["DT_Tmax_V"])))
#         print("Overhead CoolProp: {0:5.3f} us".format(np.mean(cstData["HEOS"]["DT_Tmax_V"])))#,   np.mean(cstData["HEOS"]["HP_Tmax"]))
#         print("Overhead REFPROP : {0:5.3f} us".format(np.mean(cstData["REFPROP"]["DT_Tmax_V"])))#,np.mean(cstData["REFPROP"]["HP_Tmax"]))
#         print("Mean EOS in  CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["HEOS"]["DT_HPS_V"]),np.mean(impData["HEOS"]["HP_DST_V"])))
#         print("Std. dev. in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["HEOS"]["DT_HPS_V"]) ,np.std(impData["HEOS"]["HP_DST_V"])))
#         print("Minimum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["HEOS"]["DT_HPS_V"]) ,np.min(impData["HEOS"]["HP_DST_V"])))
#         print("Maximum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["HEOS"]["DT_HPS_V"]) ,np.max(impData["HEOS"]["HP_DST_V"])))
#         print("Mean EOS in   REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["REFPROP"]["DT_HPS_V"]),np.mean(impData["REFPROP"]["HP_DST_V"])))
#         print("Std. dev. in  REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["REFPROP"]["DT_HPS_V"]) ,np.std(impData["REFPROP"]["HP_DST_V"])))
#         print("Minimum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["REFPROP"]["DT_HPS_V"]) ,np.min(impData["REFPROP"]["HP_DST_V"])))
#         print("Maximum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["REFPROP"]["DT_HPS_V"]) ,np.max(impData["REFPROP"]["HP_DST_V"])))
#         print("")
#
#     except Exception as e:
#         print(str(e.message))
#         pass
#
# try:
#     fld = 'Water'
#     cstData = {} # Data from calls to constants (overhead)
#     expData = {} # Data from explicit EOS calls
#     impData = {} # Data from EOS calls that require iterations
#     for backend in fluidData[fld]:
#         curKeys = []
#         for inp in getInpList(backend):
#             for out in getOutList(inp)[:1]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         cstData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[:1]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         expData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[1:]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         impData[backend] = curDict
#
#
#     print("Done")

def autolabel(ax,rects,means,stds,lens=0,fac=1):
    return
    # attach some text labels
    yerr = (stds*100.0)/means
    ypos = np.max(means) #+ np.max(stds)
    for i in range(len(rects)):
        xpos = rects[i].get_x()+rects[i].get_width()/2.
        #ax.text(xpos, 1.05*ypos, '{0:s}{1:4.2f}{2:s}{3:3.1f}{4:s}'.format(r'',means[i]*fac,r'us +- ',yerr[i],r'%'), rotation=90, ha='center', va='bottom', fontsize='smaller')
        ax.text(xpos, 1.05*ypos, '{0:s}{1:4.2f}{2:s}{3:3.1f}{4:s}'.format(r'\SI{',means[i]*fac,r'}{\us} (\SI{\pm',yerr[i],r'}{\percent})'), rotation=90, ha='center', va='bottom', fontsize='smaller')
        #ax.text(xpos, 0.25*ypos, str(lens[i]), rotation=90, ha='center', va='bottom', fontsize='smaller')


def axislabel(txt,ax,xPos,yPos=-1):
    ax.text(xPos, yPos, txt, rotation=45, ha='center', va='top', fontsize='xx-small')


#############################################################
# All data is available in the dicts now.
#############################################################
# The first plot contains the time data, this is averaged and
# plotted as a bar graph with the standard deviation.
hatchLst    = ["", "///"  ,"\\\\\\" ]
backendsLst = ["INCOMP","HEOS","REFPROP"]
for fluidstr in fluids[:-1]:
    _,fld,_ = splitFluid(fluidstr)
    outCst = []
    labCst = []
    outExp = []
    labExp = []
    outImp = []
    labImp = []
    DEBUG=False
    for backend in backendsLst:
        # Backend exists in fluid data?
        try:
            for i,inp in enumerate(getInpList(backend)):
                for j,out in enumerate(getOutList(inp)):
                    if j==0: # First output is Tmax
                        if backend not in fluidData[fld]:
                            outCst.append([0])
                            labCst.append("Dummy")
                            if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld,backend,inp,out))
                        else:
                            outCst.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
                            labCst.append(getTimeKey(inp, out))
                        continue
                    elif i==0: # First input is explicit
                        if backend not in fluidData[fld]:
                            outExp.append([0])
                            labExp.append("Dummy")
                            if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld,backend,inp,out))
                        else:
                            outExp.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
                            labExp.append(getTimeKey(inp, out))
                        continue
                    elif i==1:
                        if backend not in fluidData[fld]:
                            outImp.append([0])
                            labImp.append("Dummy")
                            if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld,backend,inp,out))
                        else:
                            outImp.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
                            labImp.append(getTimeKey(inp, out))
                        continue
                    else:
                        raise ValueError("Wrong number.")
        except Exception as e:
            print(e)
            sys.exit(1)

    # Do the plotting
    if bp is not None:
        bp.figure = None
        fg = bp.getFigure()
        ccycle = bp.getColorCycle(length=3)
    else:
        fg = plt.figure()
        ccycle = cycle(["b","g","r"])

    #fg = plt.figure()
    ax1 = fg.add_subplot(111)
    ax2 = ax1.twinx()
    rects1  = []
    labels1 = []
    rects2  = []
    labels2 = []
    rects3  = []
    labels3 = []
    col1 = ccycle.next()
    col2 = ccycle.next()
    col3 = ccycle.next()


    numBackends = len(backendsLst)

    step   = 1
    width  = step/(numBackends+1)  # the width of the bars
    offset = -0.5*numBackends*width

    entries = int((len(outCst)+len(outExp)+len(outImp))/numBackends)

#     for o in range(entries):
#         ids = np.empty((numBackends,))
#         for b in range(numBackends):
#             i = o*numBackends+b
#             ids[b] = offset + o*step + b*width
#             j = i - 0
#             if j < len(outCst):
#                 rects1.extend(ax1.bar(ids[b], np.mean(outCst[j]), width, color=col1, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#
#
#             j = i - 0
#             if j < len(outCst):
#                 rects1.extend(ax1.bar(ids[b], np.mean(outCst[j]), width, color=col1, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#             else:
#                 j = i-len(outCst)
#                 if j < len(outExp):
#                     rects2.extend(ax1.bar(ids[b], np.mean(outExp[j]), width, color=col2, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#                 else:
#                     j = i-len(outCst)-len(outExp)
#                     if j < len(outImp):
#                         rects3.extend(ax2.bar(ids[b], np.mean(outImp[j]), width, color=col3, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#                     else:
#                         raise ValueError("Do not go here!")
    DEBUG=False

    entries = 2
    for o in range(entries):
        ids = np.empty((numBackends,))
        for b in range(numBackends):
            i = b*entries+o
            try:
                ids[b] = offset + o*step + b*width
                rects1.extend(ax1.bar(ids[b], np.mean(outCst[i]), width, color=col1, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
                if DEBUG:
                    print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld,np.mean(outCst[i]),"cst.",b))
                    #print(ids[b],labCst[i])
                    #axislabel(labCst[i], ax1, ids[b]+0.5*width)
            except:
                pass

    offset += entries*step
    #entries = int(len(outExp)/numBackends)
    entries = 4
    for o in range(entries):
        ids = np.empty((numBackends,))
        for b in range(numBackends):
            i = b*entries+o
            try:
                ids[b] = offset + o*step + b*width
                rects2.extend(ax1.bar(ids[b], np.mean(outExp[i]), width, color=col2, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
                if DEBUG:
                    print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld,np.mean(outExp[i]),"exp.",b))
                    #print(ids[b],labExp[i])
                    #axislabel(labExp[i], ax1, ids[b]+0.5*width)
            except:
                pass

    x_newaxis = np.max(ids)+1.5*width
    plt.axvline(x_newaxis, color=bp._black, linestyle='dashed')

    offset += entries*step
    entries = 4
    for o in range(entries):
        ids = np.empty((numBackends,))
        for b in range(numBackends):
            i = b*entries+o
            try:
                ids[b] = offset + o*step + b*width
                rects3.extend(ax2.bar(ids[b], np.mean(outImp[i]), width, color=col3, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
                if DEBUG:
                    print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld,np.mean(outImp[i]),"imp.",b))
                    #print(ids[b],labImp[i])
                    #axislabel(labImp[i], ax1, ids[b]+0.5*width)
            except:
                pass


    #ax1.set_xlim([ids.min()-2.5*width,ids.max()+2.5*width])
    ax1.spines['top'].set_visible(False)
    ax2.spines['top'].set_visible(False)
    ax1.xaxis.set_ticks_position('bottom')
    ax2.xaxis.set_ticks_position('bottom')

    labels=[r"ex.",r"im.",r"$h$",r"$\rho|p$",r"$s$",r"all",r"$\rho$",r"$s$",r"$T$",r"all"]
    ax1.set_xticks(range(len(labels)))
    ax1.set_xticklabels(labels)
    #ax1.yaxis.get_label().set_verticalalignment("baseline")

    x_min = rects1[0].get_x()
    dx = rects1[0].get_width()
    x_max = rects3[-1].get_x()

    x_min = x_min - 1*dx
    x_max = x_max + 2*dx

    ax1.set_xlim([x_min,x_max])

    y_min = 0

    y_max_c = np.nanmax([a.get_height() for a in rects1])
    y_max_e = np.nanmax([a.get_height() for a in rects2])
    y_max_i = np.nanmax([a.get_height() for a in rects3])

    y_max = np.max([y_max_c,y_max_e,y_max_i/10.0])

    y_max = np.ceil(1.3*y_max/10.0)*10.0

    ax1.set_ylim([y_min,y_max])
    ax2.set_ylim([y_min,y_max*10.0])

    ratio = 10.0/4.0 * y_max / 250.0 # height of 10 for 4 points if y_max==250

    x_min = rects1[ 0].get_x()
    x_max = rects1[-1].get_x() + dx
    x,y = range_brace(x_min, x_max)
    dy  = np.ceil(y_max_c/10.0)*10.0
    y   = dy+y*ratio*(x[-1]-x[0])
    ax1.plot(x, y, ls='-',color=bp._black)
    ax1.text(np.mean(x), np.max(y), "const.", rotation=0, ha='center', va='bottom', fontsize='medium')

    x_min = rects2[ 0].get_x()
    x_max = rects2[-1].get_x() + dx
    x,y = range_brace(x_min, x_max)
    dy  = np.ceil(y_max_e/10.0)*10.0
    y   = dy+y*ratio*(x[-1]-x[0])
    ax1.plot(x, y, ls='-',color=bp._black)
    ax1.text(np.mean(x), np.max(y), "explicit", rotation=0, ha='center', va='bottom', fontsize='medium')

    x_min = rects3[ 0].get_x()
    x_max = rects3[-1].get_x() + dx
    x,y = range_brace(x_min, x_max)
    dy  = np.ceil(y_max_i/100.0)*10
    y   = dy+y*ratio*(x[-1]-x[0])
    ax1.plot(x, y, ls='-',color=bp._black)
    ax1.text(np.mean(x), np.max(y), "implicit", rotation=0, ha='center', va='bottom', fontsize='medium')

    #ax1.text(x_newaxis*0.9, y_max*0.9, "<- left axis", rotation=0, ha='right', va='bottom', fontsize='medium')
    #ax1.text(x_newaxis*1.1, y_max*0.9, "right axis ->", rotation=0, ha='left', va='bottom', fontsize='medium')

    handles = []
    for h in (rects1[0], rects2[1], rects3[2]):
        handles.append(copy.copy(h))
        handles[-1].set_facecolor('white')
        handles.append(copy.copy(h))
        handles[-1].set_hatch('')
    labels = (r'$p,T$-fit', r'constant', r'CoolProp', r'explicit, $f(p|\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$')

    if bp is not None:
        bp.drawLegend(ax=ax1,
          loc='lower center',
          bbox_to_anchor=(0.5, 1.05),
          ncol=3,
          handles=handles,
          labels=labels)
    else:
        ax1.legend(handles,labels,
          loc='lower center',
          bbox_to_anchor=(0.5, 1.05),
          ncol=3)


    ax1.set_ylabel(r'Time per explicit call (us)')
    ax2.set_ylabel(r'Time per implicit call (us)')
    fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+".pdf"))
    if bp is not None:
        ax1.set_ylabel(r'Time per explicit call (\si{\us})')
        ax2.set_ylabel(r'Time per implicit call (\si{\us})')
        mpl.rcParams['text.usetex'] = True
        fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+"-tex.pdf"))
        mpl.rcParams['text.usetex'] = False
        # Fix the wrong baseline
        for tick in ax1.get_xaxis().get_major_ticks():
            tick.set_pad(2*tick.get_pad())
            tick.label1 = tick._get_text1()
        for lab in ax1.xaxis.get_ticklabels():
            lab.set_verticalalignment("baseline")
            #lab.set_pad(1.5*lab.get_pad())


        #ax1.set_xticklabels(labels)
        #
        #for tick in ax1.xaxis.get_major_ticks():
        #    tick.label1.set_horizontalalignment('center')
        bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+".pgf"),fg,repList)

    plt.close()

# for fld in fluids:
#     try:
#         if bp is not None:
#             bp.figure = None
#             fg = bp.getFigure()
#             ccycle = bp.getColorCycle(length=3)
#         else:
#             fg = plt.figure()
#             ccycle = cycle(["b","g","r"])
#
#         #fg = plt.figure()
#         ax1 = fg.add_subplot(111)
#         ax2 = ax1.twinx()
#
#         if "INCOMP" in fluidData[fld]:
#             el = 3
#             hatch = ["","//","x"]
#         else:
#             el = 2
#             hatch = ["//","x"]
#
#         #one standard deviation above and below the mean of the data
#         width  = 0.25        # the width of the bars
#         step   = 1
#         offset = -step-0.5*el*width
#
#         lab   = []
#         rects1 = []
#         rects2 = []
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_Tmax", fluidData)*1e2)
#         curList.append(getSingleData(fld, "HEOS"   , "DT_Tmax", fluidData)*1e2)
#         curList.append(getSingleData(fld, "REFPROP", "DT_Tmax", fluidData)*1e2)
#
#         lab.extend(["Tmax"])
#
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         curCol = ccycle.next()
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[0:],np.mean(curList,axis=1),np.std(curList,axis=1),fac=1e-2)
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_H", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_H", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_H", fluidData))
#         lab.extend(["H"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         curCol = ccycle.next()
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[el:],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_D", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_P", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_P", fluidData))
#         if el==3: lab.extend(["D/P"])
#         else: lab.extend(["P"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "g"
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[int(2*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_S", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_S", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_S", fluidData))
#         lab.extend(["S"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "g"
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[int(3*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_HDS", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_HPS", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_HPS", fluidData))
#         lab.extend(["all"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "g"
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[int(4*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_D", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_D", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_D", fluidData))
#         lab.extend(["D"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         curCol = ccycle.next()
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[0:],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_T", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_T", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_T", fluidData))
#         lab.extend(["T"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "r"
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[int(1*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_S", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_S", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_S", fluidData))
#         lab.extend(["S"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "r"
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[int(2*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_DST", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_DST", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_DST", fluidData))
#         lab.extend(["all"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "r"
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[int(3*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#
#         ids = np.arange(len(lab))
#         # add some text for labels, title and axes ticks
#         #if backend=="INCOMP": ax1.set_ylabel(r'Time per $f(p,T)$ call (\si{\us})')
#         if el==3: ax1.set_ylabel(r'Time per \texttt{Tmax} call (\SI{0.01}{\us}) and'+"\n"+r'per $f(p,T)$ and $f(\rho,T)$ call (\si{\us})')
#         else: ax1.set_ylabel(r'Time per \texttt{Tmax} call (\SI{0.01}{\us})'+"\n"+r'and per $f(\rho,T)$ call (\si{\us})')
#         ax2.set_ylabel(r'Time per $f(h,p)$ call (\si{\us})')
#
#         ax1.set_xticks(ids)
#         ax1.set_xticklabels([r"\texttt{"+i+r"}" for i in lab], rotation=0)
#         ax1.set_xlim([ids.min()-2.5*width,ids.max()+2.5*width])
#
#         ax1.spines['top'].set_visible(False)
#         ax2.spines['top'].set_visible(False)
#         ax1.xaxis.set_ticks_position('bottom')
#         ax2.xaxis.set_ticks_position('bottom')
#
#         handles = []
#         if el==3:
#             for h in (rects1[0], rects1[4], rects2[2]):
#                 handles.append(copy.copy(h))
#                 handles[-1].set_facecolor('white')
#                 handles.append(copy.copy(h))
#                 handles[-1].set_hatch('')
#             labels = (r'$p,T$-fit', r'\texttt{Tmax}', r'CoolProp', r'explicit, $f(p|\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$')
#         else:
#             for h in (rects1[0], rects1[2], rects2[1]):
#                 handles.append(copy.copy(h))
#                 handles[-1].set_facecolor('white')
#                 handles.append(copy.copy(h))
#                 handles[-1].set_hatch('')
#             labels = (r'', r'\texttt{Tmax}', r'CoolProp', r'explicit, $f(\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$')
#             handles[0] = mpatches.Patch(visible=False)
#
#         if bp is not None:
#             bp.drawLegend(ax=ax1,
#               loc='upper center',
#               bbox_to_anchor=(0.5, 1.4),
#               ncol=3,
#               handles=handles,
#               labels=labels)
#         else:
#             ax1.legend(handles,labels,
#               loc='upper center',
#               bbox_to_anchor=(0.5, 1.4),
#               ncol=3)
#
#         fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld+".pdf"))
#         if bp is not None: bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+fld+".pgf"),fg,repList)
#         plt.close()
#
#     except Exception as e:
#         print(e)
#         pass

#############################################################
# The second figure compares the backend for the full calculation
#############################################################
backendExp = []
backendImp = []
fluidLabel = []
hatchLst    = ["///"  ,"\\\\\\"]
backendsLst = ["HEOS","REFPROP"]
for fluidstr in fluids[:-1]:
    _,fld,_ = splitFluid(fluidstr)
    #if fld=="CO2": fluidLabel.append("\ce{CO2}")
    #else:
    fluidLabel.append(fld)
    for backend in backendsLst:
        if backend not in fluidData[fld]:
            backendExp.append([0])
            backendImp.append([0])
            continue
        # Backend exists in fluid data
        try:
            inp = getInpList(backend)
            outExp = getTimeKey(inp[0], getOutList(inp[0])[-1])
            outImp = getTimeKey(inp[1], getOutList(inp[1])[-1])
            backendExp.append(getSingleData(fld, backend, outExp, fluidData))
            backendImp.append(getSingleData(fld, backend, outImp, fluidData))
        except Exception as e:
            backendExp.append([0])
            backendImp.append([0])
            print(e)
            pass

# Data is prepared, we can plot now.
if bp is not None:
    bp.figure = None
    fg1 = bp.getFigure()
    bp2 = BasePlotter()
    fg2 = bp2.getFigure()
    ccycle = bp.getColorCycle(length=3)
else:
    fg1 = plt.figure()
    fg2 = plt.figure()
    ccycle = cycle(["b","g","r"])

fg1.set_size_inches( (fg1.get_size_inches()[0]*1, fg1.get_size_inches()[1]*0.75) )
fg2.set_size_inches( (fg2.get_size_inches()[0]*1, fg2.get_size_inches()[1]*0.75) )
ccycle.next() # No incomp
#
#ax1 = fg.add_subplot(111)
#ax2 = ax1.twinx()
ax1 = fg1.add_subplot(111)
ax2 = fg2.add_subplot(111)

#entries = int(len(backendExp)/len(fluidLabel))
#one standard deviation above and below the mean of the data
rects1 = []
rects2 = []
col1 = ccycle.next()
col2 = ccycle.next()

numFluids   = len(fluidLabel)
numBackends = len(backendsLst)

step   = 1
width  = step/(numBackends+1)  # the width of the bars
offset = -0.5*numBackends*width

for f in range(numFluids):
    ids = np.empty((numBackends,))
    for b in range(numBackends):
        i = f*numBackends+b
        ids[b] = offset + f*step + b*width
        rects1.extend(ax1.bar(ids[b], np.mean(backendExp[i]), width, color=col1, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
        rects2.extend(ax2.bar(ids[b], np.mean(backendImp[i]), width, color=col2, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))

y_max = np.max(np.concatenate((np.ravel(ax1.get_ylim()),np.ravel(ax2.get_ylim())/10.0)))

ax1.set_ylim([0,y_max])
ax2.set_ylim([0,y_max*10.0])

for ax in [ax1,ax2]:
    ax.set_xticks(range(numFluids))
    ax.set_xticklabels(fluidLabel, rotation=25)
    ax.set_xlim([0.0-0.5*step,numFluids-1+0.5*step])
    ax.spines['top'].set_visible(False)
    ax.xaxis.set_ticks_position('bottom')

    if   ax==ax1: rects = rects1
    elif ax==ax2: rects = rects2
    handles = (rects[0], rects[1])
    labels  = (r'CoolProp', r'REFPROP')
    anchor  = (0.5, 1.2)

    if bp is not None:
        bp.drawLegend(ax=ax,
          loc='upper center',
          bbox_to_anchor=anchor,
          ncol=numBackends,
          handles=handles,
          labels=labels)
    else:
        ax.legend(handles,labels,
          loc='upper center',
          bbox_to_anchor=anchor,
          ncol=numBackends)

ax1.set_ylabel(r'Time per $f(\rho,T)$ call (us)')
ax2.set_ylabel(r'Time per $f(h,p)$ call (us)')
fg1.savefig(path.join(getFigureFolder(),"TimeComp-backends-exp.pdf"))
fg2.savefig(path.join(getFigureFolder(),"TimeComp-backends-imp.pdf"))
if bp is not None:
    ax1.set_ylabel(r'Time per $f(\rho,T)$ call (\si{\us})')
    ax2.set_ylabel(r'Time per $f(h,p)$ call (\si{\us})')
    mpl.rcParams['text.usetex'] = True
    fg1.savefig(path.join(getFigureFolder(),"TimeComp-backends-exp-tex.pdf"))
    fg2.savefig(path.join(getFigureFolder(),"TimeComp-backends-imp-tex.pdf"))
    mpl.rcParams['text.usetex'] = False
    bp.savepgf(path.join(getFigureFolder(),"TimeComp-backends-exp.pgf"),fg1,repList)
    bp.savepgf(path.join(getFigureFolder(),"TimeComp-backends-imp.pgf"),fg2,repList)
plt.close('all')

#############################################################
# The third figure is a heat map of the execution times in
# log p h diagram
#############################################################
for fluidstr in fluids:
    try:
        _,fld,_ = splitFluid(fluidstr)
        for backend in fluidData[fld]:
            propsfluid = "::".join([backend,fluidstr])

            if backend!="INCOMP":
                TP = {}
                points = int(xypoints/2)
                T_range_TP = np.linspace(PropsSI('Ttriple',"T",0,"D",0,propsfluid),PropsSI(  'Tcrit',"T",0,"D",0,propsfluid),points)
                T_TP = np.append(T_range_TP,T_range_TP[::-1])
                Q_TP = np.zeros_like(T_TP)
                Q_TP[points:] = 1
                points *= 2

                out = ["D","H","P","S"]
                res = PropsSI(out,"T",T_TP,"Q",Q_TP,propsfluid)
                D_TP = res[:,0]
                H_TP = res[:,1]
                P_TP = res[:,2]
                S_TP = res[:,3]

                mask = np.isfinite(D_TP)
                if np.sum(mask)<points:
                    warn("There were not enough valid entries in your result vector. Reducing the number of points from {0:d} to {1:d}.".format(points,np.sum(mask)))
                    points = np.sum(mask)

                TP["T"] = T_TP[mask]
                TP["D"] = D_TP[mask]
                TP["H"] = H_TP[mask]
                TP["P"] = P_TP[mask]
                TP["S"] = S_TP[mask]
                #saveNpzData(TP)
            else:
                TP = None

            state = getStateObj(propsfluid)
            if backend=="HEOS" and state.has_melting_line():
                p_melt = np.logspace(np.log10(state.melting_line(CoolProp.constants.iP_min, CoolProp.constants.iT, 0)),np.log10(state.melting_line(CoolProp.constants.iP_max, CoolProp.constants.iT, 0)),xypoints)
                #p_melt = p_range
                ML = dict(T=[],D=[],H=[],S=[],P=p_melt)
                for p in p_melt:
                    try:
                        ML["T"].append(state.melting_line(CoolProp.constants.iT, CoolProp.constants.iP, p))
                    except Exception as ve:
                        ML["T"].append(np.inf)
                res = PropsSI(["D","H","P","S","T"],"T",ML["T"],"P",ML["P"],propsfluid)
                ML["D"] = res[:,0]
                ML["H"] = res[:,1]
                ML["P"] = res[:,2]
                ML["S"] = res[:,3]
                ML["T"] = res[:,4]
                mask = np.isfinite(ML["T"])
                ML["P"] = ML["P"][mask]
                ML["T"] = ML["T"][mask]
                ML["D"] = ML["D"][mask]
                ML["H"] = ML["H"][mask]
                ML["S"] = ML["S"][mask]
            else:
                ML = None
            #ML = {}

            IP = {}
            p_range,T_range = getPTRanges(propsfluid)
            critProps = getCriticalProps(propsfluid)
            try:
                IP["T"] = T_range
                IP["P"] = np.zeros_like(T_range) + critProps["P"]
                res = PropsSI(["D","H"],"T",IP["T"],"P",IP["P"],propsfluid)
                IP["D"] = res[:,0]
                IP["H"] = res[:,1]
            except Exception as ve:
                IP = None

            IT = {}
            try:
                IT["P"] = p_range
                IT["T"] = np.zeros_like(p_range) + critProps["T"]
                res = PropsSI(["D","H"],"T",IT["T"],"P",IT["P"],propsfluid)
                IT["D"] = res[:,0]
                IT["H"] = res[:,1]
            except Exception as ve:
                IT = None

            ID = {}
            try:
                ID["T"] = T_range
                ID["D"] = np.zeros_like(p_range) + critProps["D"]
                res = PropsSI(["P","H"],"T",ID["T"],"D",ID["D"],propsfluid)
                ID["P"] = res[:,0]
                ID["H"] = res[:,1]
            except Exception as ve:
                ID = None

            IH = {}
            try:
                IH["P"] = p_range
                IH["H"] = np.zeros_like(p_range) + critProps["H"]
                res = PropsSI(["D","T"],"P",IH["P"],"H",IH["H"],propsfluid)
                IH["D"] = res[:,0]
                IH["T"] = res[:,1]
            except Exception as ve:
                IH = None

            IS = {}
            try:
                IS["P"] = p_range
                IS["S"] = np.zeros_like(p_range) + critProps["S"]
                res = PropsSI(["D","H","T"],"P",IS["P"],"S",IS["S"],propsfluid)
                IS["D"] = res[:,0]
                IS["H"] = res[:,1]
                IS["T"] = res[:,2]
            except Exception as ve:
                IS = None

            for I in [IP,IT,ID,IH,IS]:
                if I is not None:
                    mask = np.isfinite(I["D"]) & np.isfinite(I["H"])
                    if np.sum(mask)<20: I = None
                    else:
                        for k in I:
                            I[k] = I[k][mask]






            for inp in getInpList(backend):
                if bp is not None:
                    bp.figure = None
                    fg = bp.getFigure()
                else:
                    fg = plt.figure()

                kind = getTimeKey(inp, getOutList(inp)[-1])
                t_data = getSingleData(fld, backend, kind, fluidData)
                x_data = getSingleData(fld, backend, "H", fluidData)
                y_data = getSingleData(fld, backend, "P", fluidData)

                gs  = gridspec.GridSpec(1, 2, wspace=None, hspace=None, width_ratios=[10,1])
                ax1 = fg.add_subplot(gs[0,0])
                ax1.set_yscale('log')
                #ax2 = ax1.twinx()

                minHP = np.min(t_data)
                maxHP = np.max(t_data)
                minIT = np.percentile(t_data,10)
                maxIT = np.percentile(t_data,90)
                difIT = np.log10(maxIT/minIT)*0.25

                print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
                if kind == "DT":
                    if bp is not None:
                        cx1=bp.getColourMap(reverse=True)
                    else:
                        cx1   = mpl.cm.get_cmap('cubehelix_r')
                    minHP = minIT
                    maxHP = np.power(10,np.log10(maxIT)+difIT)
                    #minHP = np.power(10,np.log10(np.percentile(t_data,10)*1e6))
                    #maxHP = np.power(10,np.log10(np.percentile(t_data,90)*1e6)*1.10)
                    #maxHP = np.power(10,1.10*np.log10(maxHP))
                    #minHP = np.percentile(t_data,10)*1e6
                    #maxHP = np.percentile(t_data,99)*1e6
                    #print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
                    #minHP = 100
                    #maxHP = 20000
                else:
                    if bp is not None:
                        cx1 = bp.getColourMap()
                    else:
                        cx1   = mpl.cm.get_cmap('cubehelix')
                    minHP = np.power(10,np.log10(minIT)-difIT)
                    maxHP = maxIT
                    #minHP = np.power(10,np.log10(np.percentile(t_data,10)*1e6)*0.90)
                    #maxHP = np.power(10,np.log10(np.percentile(t_data,90)*1e6))
                    #minHP = np.percentile(t_data,01)*1e6
                    #maxHP = np.percentile(t_data,90)*1e6
                    #print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
                    #minHP = 100
                    #maxHP = 20000




                #cx1_r = reverse_colourmap(cx1)
                cNorm  = mpl.colors.LogNorm(vmin=minHP, vmax=maxHP)
                #cNorm  = mpl.colors.LogNorm(vmin=ceil(minHP/1e1)*1e1, vmax=floor(maxHP/1e2)*1e2)
                #cNorm  = mpl.colors.Normalize(vmin=round(minHP,-2), vmax=round(maxHP,-2))
                colourSettings  = dict(c=t_data, edgecolors='none', cmap=cx1, norm=cNorm)
                pointSettings   = dict(s=4)
                scatterSettings = dict(rasterized=True, alpha=0.5)
                #scatterSettings = dict(rasterized=False, alpha=0.5)
                scatterSettings.update(colourSettings)
                scatterSettings.update(pointSettings)
                SC = ax1.scatter(x_data/1e6,y_data/1e5, **scatterSettings)

                for I in [TP,ML]:
                    if I is not None:
                        ax1.plot(I["H"]/1e6,I["P"]/1e5,lw=1.5,c=bp._black)

                for I in [IP,IT,ID,IS,IH]:
                    if I is not None:
                        ax1.plot(I["H"]/1e6,I["P"]/1e5,lw=1.0,c=bp._black,alpha=1)


                #ax1.set_xlim([0e+0,6e1])
                #ax1.set_ylim([5e-1,2e4])
                ax1.set_xlim([np.percentile(x_data/1e6,0.1),np.percentile(x_data/1e6,99.9)])
                ax1.set_ylim([np.percentile(y_data/1e5,0.1),np.percentile(y_data/1e5,99.9)])

                formatter = ticker.LogFormatter(base=10.0, labelOnlyBase=False)
                #formatter = ticker.ScalarFormatter()
                #ticks     = roundList(np.logspace(np.log10(ax1.get_ylim()[0]), np.log10(ax1.get_ylim()[1]), 5))
                #locator   = ticker.FixedLocator(ticks)
                #ax1.yaxis.set_major_locator(locator)
                ax1.yaxis.set_major_formatter(formatter)

                cax = fg.add_subplot(gs[0,1])
                formatter = ticker.ScalarFormatter()
                CB  = fg.colorbar(SC, cax=cax,format=formatter)
                CB.set_alpha(1)
                #CB.locator = ticker.MaxNLocator(nbins=7)
                ticks = roundList(np.logspace(np.log10(minHP), np.log10(maxHP), 5))
                CB.locator = ticker.FixedLocator(ticks)
                CB.update_ticks()
                CB.draw_all()


                CB.set_label(r'Execution time per call (us)')
                ax1.set_xlabel(r'Specific enthalpy (MJ/kg)')
                ax1.set_ylabel(r'Pressure (bar)')
                fg.tight_layout()
                fg.savefig(path.join(getFigureFolder(),"TimeComp-"+inp+"-"+backend.lower()+"-"+fld.lower()+".pdf"))

                CB.set_label(r'Execution time per call (\si{\us})')
                ax1.set_xlabel(r'Specific enthalpy (\si{\mega\J\per\kg})')
                ax1.set_ylabel(r'Pressure (\si{\bar})')
                fg.tight_layout()
                bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+inp+"-"+backend.lower()+"-"+fld.lower()+".pgf"),fg,repList)
                plt.close()

    except Exception as e:
        print(e)
        pass
    plt.close('all')




















