wetb/hawc2/tests/test_turbulence_file.py

+'''
+Created on 17/07/2014
+
+@author: MMPE
+'''
+from wetb.hawc2.turbulence_file import TurbulenceFile
+import os
+import unittest
+from wetb.hawc2.pc_file import PCFile
+
+
+import numpy as np
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
+
+
+class TestPCFile(unittest.TestCase):
+
+    def setUp(self):
+        unittest.TestCase.setUp(self)
+
+    # def test_TurbulenceFile(self):
+    #    TurbulenceFile(tfp + "turb/mann_l29.4_ae1.00_g3.9_h1_64x8x8_1.000x2.00x2.00_s1001u.turb")
+
+    def test_load_from_htc(self):
+        u, v, w = TurbulenceFile.load_from_htc(tfp + "htcfiles/test_turb.htc")
+        self.assertEqual(u.data.shape, (512, 64))
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()

wetb/hawc2/turbulence_file.py

+'''
+Created on 22. jun. 2017
+
+@author: mmpe
+'''
+import os
+
+import numpy as np
+from wetb.hawc2.htc_file import HTCFile
+from wetb.wind.turbulence import mann_turbulence
+
+
+class TurbulenceFile(object):
+    def __init__(self, filename, Nxyz, dxyz, transport_speed=10, mean_wsp=0, center_position=(0, 0, -20)):
+        self.filename = filename
+        self.Nxyz = Nxyz
+        self.dxyz = dxyz
+        self.transport_speed = transport_speed
+        self.mean_wsp = mean_wsp
+        self.center_position = center_position
+        self.data = mann_turbulence.load(filename, Nxyz)
+
+    @property
+    def data3d(self):
+        return self.data.reshape(self.Nxyz)
+
+    @staticmethod
+    def load_from_htc(htcfilename, modelpath=None, type='mann'):
+        htc = HTCFile(htcfilename, modelpath)
+
+        Nxyz = np.array([htc.wind[type]['box_dim_%s' % uvw][0] for uvw in 'uvw'])
+        dxyz = np.array([htc.wind[type]['box_dim_%s' % uvw][1] for uvw in 'uvw'])
+        center_position = htc.wind.center_pos0.values
+        wsp = htc.wind.wsp
+        return [TurbulenceFile(os.path.join(htc.modelpath, htc.wind[type]['filename_%s' % uvw][0]), Nxyz, dxyz, wsp, (0, wsp)[uvw == 'u'], center_position) for uvw in 'uvw']
+
+
+if __name__ == '__main__':
+    pass

wetb/hawc2flow/Evaluate.py

+import pandas as pd
+import numpy as np
+import re, os
+import click
+from scipy import signal
+from functools import wraps 
+from wetb.hawc2.Hawc2io import ReadHawc2
+from wetb.fatigue_tools.fatigue import eq_load
+
+    
+def isFloat(string):
+    try:
+        float(string)
+        return True
+    except ValueError:
+        return False
+
+
+class HAWC2DataFrame(pd.DataFrame):
+    """
+    A subclass of a pandas.DataFrame. The HAWC2DataFrame is able to link to a
+    directory of HAWC2 result files. The HAWC2DataFrame contains an Accessor,
+    named 'wetb', which contains functions to  populate the HAWC2DataFrame with
+    post-processed statistics, and to aggregate the statistics.
+    """
+
+    _metadata = ['_fields', 'channels', '_directory', '_filenames']
+    
+    @property
+    def _constructor(self):
+        return HAWC2DataFrame
+        
+        
+    def __init__(self, *args, dir=None, pattern=None, channels=None, **kwargs):
+        
+        if all(x is not None for x in [dir, channels]):            
+            if pattern is None:
+                res = self.link_directory(dir, channels)
+            else:
+                res = self.link_results(dir, pattern, channels)
+            super(HAWC2DataFrame, self).__init__(res)
+
+        else:
+            super(HAWC2DataFrame, self).__init__(*args, **kwargs)
+
+
+    
+    def __finalize__(self, other, method=None, **kwargs):
+        """propagate metadata from other to self """
+        # merge operation: using metadata of the left object
+        if method == 'merge':
+            for name in self._metadata:
+                object.__setattr__(self, name, getattr(other.left, name, None))
+        # concat operation: using metadata of the first object
+        elif method == 'concat':
+            for name in self._metadata:
+                object.__setattr__(self, name, getattr(other.objs[0], name, None))
+        else:
+            for name in self._metadata:
+                object.__setattr__(self, name, getattr(other, name, None))
+        return self
+
+    
+    def link_results(self, directory, pattern_string, channels):
+        self._directory = directory
+        pattern, self._fields = self._compile_pattern(pattern_string)
+
+        # get all filenames that fit the pattern
+        self._filenames = [x[:-4] for x in os.listdir(directory) if x.endswith('.sel')]
+        self._filenames = [x for x in self._filenames if pattern.match(x)]
+        
+        self.channels = channels
+        # Extract input attributes and put in dataframe
+        dat = []
+        for fn in self._filenames:
+            dat.append([float(x) if isFloat(x) else x for x in pattern.findall(fn)[0]])
+        dat = pd.DataFrame(dat)
+        
+        # set column multi index
+        if not dat.empty:
+            column_tuples = list(zip(*[self._fields, ['']*len(self._fields)]))
+            dat.columns = pd.MultiIndex.from_tuples(column_tuples, names=['channel', 'stat'])
+        return dat
+    
+    
+    def link_directory(self, directory, channels):
+        '''
+        Links all result files in a directory. Used if no pattern is given.
+        '''
+        self._directory = directory
+        self.channels = channels
+
+        self._fields = ['filename']
+        # get all filenames of result files
+        self._filenames = [x[:-4] for x in os.listdir(directory) if x.endswith('.sel')]
+        dat = pd.DataFrame(self._filenames)
+
+        # set column multi index
+        if not dat.empty:
+            column_tuples = list(zip(*[self._fields, ['']*len(self._fields)]))
+            dat.columns = pd.MultiIndex.from_tuples(column_tuples, names=['channel', 'stat'])
+        return dat
+        
+            
+    @staticmethod
+    def _compile_pattern(pattern_string):
+        brackets = re.compile('{(.*?)}')
+        fields = brackets.findall(pattern_string)
+        for field in fields:
+            pattern_string = pattern_string.replace('{'+ field +'}', '(.*)')
+        return re.compile(pattern_string), fields   
+        
+                 
+    def __call__(self, **kwargs):
+        return self[self._mask( **kwargs)]
+
+
+    def _mask(self, **kwargs):
+        '''
+        Returns a mask for refering to a dataframe, or self.Data, or self.Data_f, etc.
+        example. dlc.mask(wsp=[12, 14], controller='noIPC')
+        '''
+        N = len(self)
+        mask = [True] * N
+        for key, value in kwargs.items():
+            if isinstance(value, (list, tuple, np.ndarray)):
+                mask_temp = [False] * N
+                for v in value:
+                    mask_temp = mask_temp | (self[key] == v)
+                mask = mask & mask_temp
+            else: #scalar, or single value
+                mask = mask & (self[key] == value)
+        return mask
+
+
+
+@pd.api.extensions.register_dataframe_accessor("wetb")
+class wetbAccessor(object):
+    def __init__(self, pandas_obj):
+        self._validate(pandas_obj)
+        self._obj = pandas_obj
+    
+        
+    @staticmethod
+    def _validate(obj):
+        # verify that a HAWC2DataFrame is using the wetb accessor.
+        if not isinstance(obj, HAWC2DataFrame):
+            raise TypeError(f"wetb namespace is not accessible by type {type(obj)}. Expected: { type(HAWC2DataFrame())}")
+            
+            
+    @classmethod    
+    def populate_method(cls, method_name, label=None):
+        label = label or method_name
+        def decorator(func):
+            @wraps(func) 
+            def wrapper(self, channels, *args, **kwargs): 
+                return self._add_stat(func, label, channels, *args, **kwargs)
+            setattr(cls, method_name, wrapper)
+            # Note we are not binding func, but wrapper which accepts self but does exactly the same as func
+            return func # returning func means func can still be used normally
+        return decorator   
+        
+                 
+    def fetch_result(self, idx, channels=None):
+        '''
+        Returns a time series dataframe of a single result file given an index number.
+        '''
+        fn = self._obj._filenames[idx]
+        channels = channels or self.channels
+
+        raw = ReadHawc2(os.path.join(self._obj._directory, fn)).ReadAll(ChVec=[i-1 for i in channels.values()])
+        #raw = readHawc2Res(os.path.join(self._directory, fn), channels)
+        return pd.DataFrame(raw, columns=channels.keys())
+        # except:
+        #     print(f'File {fn} could not be loaded.')
+        
+            
+    def iter_sim(self, **kwargs):
+        '''
+        Iterates over simulation result files given the input filter defined in kwargs.
+        '''
+        sims_to_iterate = self._obj(**kwargs)
+        for idx, row in sims_to_iterate.iterrows():
+            raw = self.fetch_result(idx)
+            yield row[self._obj._fields], raw
+
+    
+    def _add_stat(self, func, stat_name, channels=None, *args, **kwargs):
+        '''
+        Adds a column of statistics for the given channels using the given function.
+        The function should take a pandas series (1d array) and return a float.
+        '''
+        values = []
+        if channels is None:
+            channels = self._obj.channels
+        else:
+            channels = {k:v for k,v in self._obj.channels.items() if k in channels}
+        
+        channel_string = ', '.join(channels)
+        print(f'Calculating {stat_name} for {channel_string}...')
+        
+        with click.progressbar(self._obj.index) as bar:
+            for idx in bar:
+                raw = self.fetch_result(idx, channels)
+                values.append(func(raw, *args, **kwargs))
+
+        df = pd.DataFrame(values)
+        # add multi index columns
+        col_ch     = list(channels.keys())
+        col_stat   = [stat_name]*len(col_ch)
+        col_tuples = list(zip(*[col_ch, col_stat]))
+        df.columns = pd.MultiIndex.from_tuples(col_tuples, names=['channel', 'stat'])
+
+        return self._obj.join(df)
+
+
+    def aggregate_rows(self, key):
+        # used for taking the mean over all seeds
+        print(f'Calculating mean over key={key}...')
+        in_fields = [x for x in self._obj._fields if x != key]
+        out_fields = [x for x in list(self._obj) if x[0] not in self._obj._fields]
+        in_atts = self._obj[in_fields].drop_duplicates()
+        new_dat = []
+        with click.progressbar(in_atts.iterrows()) as bar:
+            for _, x in bar:
+                filt = {k[0]: v for k, v in dict(x).items()}
+                new_dat.append(list(x) + list(self._obj(**filt)[out_fields].mean().values))
+            
+        # Make new HAWC2DataFrame and copy metadata (is there a better way?)
+        new_df = HAWC2DataFrame(new_dat)
+        for attr in self._obj._metadata:
+            setattr(new_df, attr, getattr(self._obj, attr))
+        new_df._fields = in_fields
+
+        # add multi index columns
+        col_ch     = in_fields + [x[0] for x in out_fields]
+        col_stat   = ['']*len(in_fields) + [x[1] for x in out_fields]
+        col_tuples = list(zip(*[col_ch, col_stat]))
+        new_df.columns = pd.MultiIndex.from_tuples(col_tuples,
+                    names=['channel', 'stat'])
+        return new_df
+
+
+    def aggregate_columns(self, source, dest, drop=False):
+        # used for taking the mean over all blades.
+        channel_string = ', '.join(source)
+        print(f'Calculating mean over {channel_string}...')
+        new_df = self._obj.copy()
+        # get all unique stat indices
+        col_stat = list(set(x[1] for x in list(self._obj)))
+        for stat in col_stat:
+            keys = [x for x in list(self._obj) if x[0] in source and x[1] == stat]
+            if not keys:
+                continue
+            new_df[(dest, stat)] = self._obj[keys].mean(axis=1)
+            
+        if drop:
+            # delete source columns
+            new_df = new_df.drop(columns=source)
+        return new_df
+
+    
+    @staticmethod
+    def read_csv(filename):
+        '''
+        Reads a postproc csv file which was generated by the HAWC2Res class. Returns
+        a DataFrame
+        '''
+        df           = HAWC2DataFrame(pd.read_csv(filename, header =[0, 1]))
+        multicol     = list(df)
+        multicol_new = []
+        for col in multicol:
+            if 'Unnamed' in col[1]:
+                multicol_new.append((col[0], ''))
+            else:
+                multicol_new.append(col)
+        df.columns = pd.MultiIndex.from_tuples(multicol_new)
+        return df
+
+
+@wetbAccessor.populate_method('DEL')
+def DEL(x, m=4):
+    '''
+    Calculates the damage equivalent load of a set of time series loads.
+    args:
+        x (pd.DataFrame): load time series, where each column is a different load channel
+        m (float): wohler constant (default: 4)
+    returns:
+        DEL (list of floats): Damage equivalent load for each load channel.
+    '''
+    DEL = []
+    for k in list(x):
+        DEL.append(eq_load(x[k].values, m=m)[0][0])
+    return DEL
+    
+    
+@wetbAccessor.populate_method('mean', label='Mean')
+def _mean(x):
+    '''
+    Calculates the mean of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): Mean of each time series.
+    '''
+    return x.mean().values
+    
+    
+@wetbAccessor.populate_method('var', label='Var')
+def _var(x):
+    '''
+    Calculates the variance of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): variance of each time series.
+    '''
+    return x.var().values
+
+
+@wetbAccessor.populate_method('std', label='Std')
+def _std(x):
+    '''
+    Calculates the standard deviation of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): standard deviation of each time series.
+    '''
+    return x.std().values
+
+
+@wetbAccessor.populate_method('final')
+def _final(x):
+    '''
+    Returns the final value of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+    returns:
+        (list of floats): final value of each time series.
+    '''
+    return x.iloc[-1].values
+
+
+@wetbAccessor.populate_method('PSD')
+def _psd(x, fs=100):
+    '''
+    Calculates power spectral density (PSD) of a set of time series.
+    args:
+        x (pd.DataFrame): time series, where each column contains a different time series.
+        fs (float): sampling frequency in Hz (default: 100Hz).
+    returns:
+        (list of arrays): PSD of each time series.
+    '''
+    PSD = []
+    for k in list(x):
+        f, Py = signal.welch(x[k].values, fs=fs, nperseg=1024*8)
+        PSD.append(Py)
+    return PSD
+            

wetb/hawc2flow/Generate.py

+from wetb.hawc2.hawc2_input_writer import HAWC2InputWriter, JinjaWriter

wetb/prepost/DataChecks.py

@@ -4,15 +4,6 @@ Created on Mon Mar  5 16:00:02 2012
 
 @author: dave
 """
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from future import standard_library
-standard_library.install_aliases()
-
-
-
 
 # time and data should be 1 dimensional arrays
 def array_1d(array):

wetb/prepost/GenerateDLCs.py

@@ -4,17 +4,12 @@ Created on Fri Nov 20 10:11:06 2015
 
 @author: tlbl
 """
-from __future__ import print_function
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import absolute_import
-
 # arctan and pi are required because they are in the formulas that are
 # evaluated
 from numpy import (floor, arctan, pi, log, log10, sin, cos, tan, e, arcsin,
                    arccos)
 import pandas as pd
-import xlrd
+from openpyxl import load_workbook
 from argparse import ArgumentParser
 import os
 
@@ -99,6 +94,8 @@ class GeneralDLC(object):
                     counter = floor(irow/len_wsp) + 1
                     value = '%4.4i' % (1000*counter + row[i_wsp] + 1)
                 elif variables_order[icol] == '[wave_seed]':
+                    # FIXME: shouldn't we also have unique wave seeds??
+                    # that is not the case with this implementation
                     value = '%4.4i' % (100*(row[i_wsp]+1) + row[i_wave_seed] + 1)
 #                    value = '%4.4i' % (1000*counter + row[i_wsp] + 101)
 #                    value = '%4.4i' % (irow+1)
@@ -120,7 +117,6 @@ class GeneralDLC(object):
 
     def sort_formulas(self, formulas):
         # sort formulas based on their dependency
-
         keys_list = sorted(formulas)
         for i in range(len(keys_list)):
             for ikey, key in enumerate(keys_list):
@@ -218,72 +214,71 @@ class GenerateDLCCases(GeneralDLC):
 
     def execute(self, filename='DLCs.xlsx', folder='', isheets=None):
 
-        book = xlrd.open_workbook(filename)
+        book = load_workbook(filename, data_only=True)
+
+        # Read all the initialization constants and functions of the main sheet
+        # The main sheet is assumed to be the first one (i=0)
+        general_constants = {}
+        general_functions = {}
+        sheet = book.worksheets[0]
+        # note that sheet.cell(i,j) is using 1-based indexing (row nr 1)
+        # first column is just for readability
+        for colnr in range(2, sheet.max_column+1):
+            if sheet.cell(10, colnr).value != None:
+                general_constants[str(sheet.cell(10, colnr).value)] = \
+                    sheet.cell(11, colnr).value
+            if sheet.cell(14, colnr).value != None:
+                general_functions[str(sheet.cell(14, colnr).value)] = \
+                    sheet.cell(15, colnr).value
 
         if isheets is None:
-            isheets = list(range(1, book.nsheets))
+            # refer to sheet number, so 1-based indexing
+            isheets = list(range(1,len(book.sheetnames)))
 
         # Loop through all the sheets. Each sheet correspond to a DLC.
         for isheet in isheets:
 
-            # Read all the initialization constants and functions in the
-            # first sheet
-            general_constants = {}
-            general_functions = {}
-            sheet = book.sheets()[0]
-            for i in range(1, sheet.ncols):
-                if sheet.cell_value(9, i) != '':
-                    general_constants[str(sheet.cell_value(9, i))] = \
-                        sheet.cell_value(10, i)
-                if sheet.cell_value(13, i) != '':
-                    general_functions[str(sheet.cell_value(13, i))] = \
-                        sheet.cell_value(14, i)
-
-            sheet = book.sheets()[isheet]
-
-            print('Sheet #%i' % isheet, sheet.name)
+            sheet = book.worksheets[isheet]
 
+            print('Sheet #%i' % isheet, book.sheetnames[isheet])
             # Read the actual sheet.
             constants = {}
             variables = {}
             formulas = {}
             variables_order = []
             # Loop through the columns
-            for i in range(sheet.ncols):
-                if sheet.cell_value(1, i) is not None:
-                    tag = str(sheet.cell_value(1, i))
-                    if tag is not '':
+            for i in range(1,sheet.max_column+1):
+                if sheet.cell(1, i).value is not None:
+                    tag = str(sheet.cell(1+1, i).value)
+                    if len(tag) > 0:
                         # FIXME: only works if [wsp] is defined as variable
                         # and [seed] tags are present
-                        if sheet.cell_value(0, i) == 'C':
-                            constants[tag] = sheet.cell_value(2, i)
-                        if sheet.cell_value(0, i) == 'V':
+                        if sheet.cell(1, i).value == 'C':
+                            constants[tag] = sheet.cell(3, i).value
+                        if sheet.cell(1, i).value == 'V':
                             variables_order.append(tag)
                             variables[tag] = \
-                                [sheet.cell_value(j, i) for j in range(2, sheet.nrows)]
-                        if sheet.cell_value(0, i) == 'F':
-                            formulas[tag] = str(sheet.cell_value(2, i))
-
+                                [sheet.cell(j, i).value for j in range(3, sheet.max_row+1) if sheet.cell(j, i).value != None]
+                        if sheet.cell(1, i).value == 'F':
+                            formulas[tag] = str(sheet.cell(3, i).value)
             dlc = {}
 
-            general_constants = self.remove_from_dict(variables,
-                                                      general_constants)
-            general_constants = self.remove_from_dict(constants,
-                                                      general_constants)
-            general_functions = self.remove_from_dict(formulas,
-                                                      general_functions)
+            # make copies of the general constants and functions since remove
+            # will otherwise remove them for the following sheets as well
+            sheet_gen_con = self.remove_from_dict(variables, general_constants.copy())
+            sheet_gen_con = self.remove_from_dict(constants, sheet_gen_con)
+            sheet_gen_fun = self.remove_from_dict(formulas, general_functions.copy())
 
             self.add_variables_tag(dlc, variables, variables_order)
-            self.add_constants_tag(dlc, general_constants)
+            self.add_constants_tag(dlc, sheet_gen_con)
             self.add_constants_tag(dlc, constants)
             self.add_formulas(dlc, formulas)
-            self.add_formulas(dlc, general_functions)
+            self.add_formulas(dlc, sheet_gen_fun)
             # TODO: before eval, check if all tags in formula's are present
             self.eval_formulas(dlc)
             df = pd.DataFrame(dlc)
-            if not os.path.exists(folder):
-                os.makedirs(folder)
-            df.to_excel(os.path.join(folder, sheet.name+'.xlsx'), index=False)
+            os.makedirs(folder, exist_ok=True)
+            df.to_excel(os.path.join(folder, book.sheetnames[isheet]+'.xlsx'), index=False)
 
 
 if __name__ == '__main__':

wetb/prepost/GenerateHydro.py

@@ -10,13 +10,35 @@ Description: This script is used for writing the hydro input files for HAWC2
 
 import os
 
+template = """
+begin wkin_input ;
+    wavetype {wavetype} ;
+    wdepth  {wdepth} ;
+    begin reg_airy ;
+        stretching {stretching} ; 0=none, 1=wheeler
+        wave {Hs} {Tp} {wave_phase_shift} ; (Hs, Tp, phase shift)
+;        ignore_water_surface ;
+    end reg_airy ;
+    begin ireg_airy ;
+        stretching {stretching} ; 0=none, 1=wheeler
+        spectrum   {spectrum_id} ; (1=jonswap, 2=pm)
+        mccamyfuchs {mccamyfuchs} {mccamyfuchs_radius};
+        jonswap {Hs} {Tp} {gamma} ; (Hs, Tp, gamma)
+        pm {Hs} {Tp} ; (Hs, Tp)
+;        wn {wn_variance} {wn_f0} {wn_f1};
+        coef {coef_nrs} {coef_seed} {coef_phase_shift} ; (number of coefficients, seed, phase shift)
+        spreading {spreading} {spreading_param};
+    end;
+end;
+"""
+
 class hydro_input(object):
 
     """ Basic class aiming for write the hydrodynamics input file"""
 
     def __init__(self, wavetype, wdepth, spectrum, Hs, Tp, seed, gamma=3.3,
                  stretching=1, wn=None, coef=200, spreading=None,
-                 embed_sf=None, embed_sf_t0=None):
+                 embed_sf=None, embed_sf_t0=None, mccamyfuchs=1):
 
         self.wdepth = wdepth
         if self.wdepth < 0:

wetb/prepost/Simulations.py

@@ -6,21 +6,6 @@ Created on Tue Nov  1 15:16:34 2011
 __author__ = "David Verelst <dave@dtu.dk>"
 __license__ = "GPL-2+"
 """
-
-from __future__ import print_function
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import dict
-from io import open
-from builtins import zip
-from builtins import range
-from builtins import str
-from builtins import int
-from future import standard_library
-standard_library.install_aliases()
-from builtins import object
-
 # standard python library
 import os
 import subprocess as sproc
@@ -41,22 +26,27 @@ from time import time
 #import threading
 #from multiprocessing import Pool
 
-# numpy and scipy only used in HtcMaster._all_in_one_blade_tag
 import numpy as np
 import scipy
-import scipy.interpolate as interpolate
 #import matplotlib.pyplot as plt
 import pandas as pd
 import tables as tbl
 
 # custom libraries
+from wetb.fatigue_tools.bearing_damage import bearing_damage
 from wetb.prepost import misc
 from wetb.prepost import windIO
 from wetb.prepost import prepost
 from wetb.dlc import high_level as dlc
-from wetb.prepost.GenerateHydro import hydro_input
+#from wetb.prepost.GenerateHydro import hydro_input
 from wetb.utils.envelope import compute_envelope
 
+#def join_path(*args):
+#    return os_path_join(*args).replace("\\","/")
+#os.path.join = join_path
+
+
+
 def load_pickled_file(source):
     FILE = open(source, 'rb')
     result = pickle.load(FILE)
@@ -398,7 +388,6 @@ def run_local_ram(cases, check_log=True):
 
     return cases
 
-
 def run_local(cases, silent=False, check_log=True):
     """
     Run all HAWC2 simulations locally from cases
@@ -464,7 +453,7 @@ def run_local(cases, silent=False, check_log=True):
         # create the required directories
         dirkeys = ['[data_dir]', '[htc_dir]', '[res_dir]', '[log_dir]',
                    '[eigenfreq_dir]', '[animation_dir]', '[turb_dir]',
-                   '[wake_dir]', '[meander_dir]', '[opt_dir]', '[control_dir]',
+                   '[micro_dir]', '[meander_dir]', '[opt_dir]', '[control_dir]',
                    '[mooring_dir]', '[hydro_dir]', '[externalforce]']
         for dirkey in dirkeys:
             if tags[dirkey]:
@@ -552,17 +541,17 @@ def run_local(cases, silent=False, check_log=True):
 
     return cases
 
-
 def prepare_launch(iter_dict, opt_tags, master, variable_tag_func,
                 write_htc=True, runmethod='none', verbose=False,
                 copyback_turb=True, msg='', silent=False, check_log=True,
                 update_cases=False, ignore_non_unique=False,
-                run_only_new=False, windows_nr_cpus=2, wine_64bit=False,
+                run_only_new=False, windows_nr_cpus=2,
                 pbs_fname_appendix=True, short_job_names=True, qsub='',
-                update_model_data=True, maxcpu=1, pyenv='wetb_py3',
+                update_model_data=False, maxcpu=1, pyenv='py36-wetb',
                 m=[3,4,6,8,9,10,12], postpro_node_zipchunks=True,
                 postpro_node=False, exesingle=None, exechunks=None,
-                wine_arch='win32', wine_prefix='~/.wine32'):
+                wine_arch='win32', wine_prefix='~/.wine32', prelude='',
+                pyenv_cmd='source /home/python/miniconda3/bin/activate'):
     """
     Create the htc files, pbs scripts and replace the tags in master file
     =====================================================================
@@ -634,10 +623,9 @@ def prepare_launch(iter_dict, opt_tags, master, variable_tag_func,
         line is added, and when launching the job this dependency needs to
         specified.
 
-    update_model_data : default=True
+    update_model_data : default=False
         If set to False, the zip file will not be created, and the data files
-        are not copied to the run_dir. Use this when only updating the htc
-        files.
+        are not copied to the run_dir.
 
     Returns
     -------
@@ -722,7 +710,8 @@ def prepare_launch(iter_dict, opt_tags, master, variable_tag_func,
             # returns a dictionary with all the tags used for this
             # specific case
             htc = master.createcase(write_htc=write_htc)
-            master.create_run_dir()
+            if update_model_data:
+                master.create_run_dir()
             #htc=master.createcase_check(cases_repo,write_htc=write_htc)
 
             # make sure the current cases is unique!
@@ -802,213 +791,20 @@ def prepare_launch(iter_dict, opt_tags, master, variable_tag_func,
            copyback_turb=copyback_turb, qsub=qsub,
            windows_nr_cpus=windows_nr_cpus, short_job_names=short_job_names,
            pbs_fname_appendix=pbs_fname_appendix, silent=silent, maxcpu=maxcpu,
-           pyenv=pyenv, wine_64bit=wine_64bit, m=[3,4,6,8,9,10,12],
-           postpro_node_zipchunks=postpro_node_zipchunks,
+           pyenv=pyenv, m=[3,4,6,8,9,10,12],
+           postpro_node_zipchunks=postpro_node_zipchunks, prelude=prelude,
            postpro_node=postpro_node, exesingle=exesingle, exechunks=exechunks,
-           wine_arch=wine_arch, wine_prefix=wine_prefix)
+           wine_arch=wine_arch, wine_prefix=wine_prefix, pyenv_cmd=pyenv_cmd)
 
     return cases
 
-def prepare_relaunch(cases, runmethod='gorm', verbose=False, write_htc=True,
-                     copyback_turb=True, silent=False, check_log=True):
-    """
-    Instead of redoing everything, we know recreate the HTC file for those
-    in the given cases dict. Nothing else changes. The data and zip files
-    are not updated, the convience tagfile is not recreated. However, the
-    saved (pickled) cases dict corresponding to the sim_id is updated!
-
-    This method is usefull to correct mistakes made for some cases.
-
-    It is adviced to not change the case_id, sim_id, from the cases.
-    """
-
-    # initiate the HtcMaster object, load the master file
-    master = HtcMaster()
-    # for invariant tags, load random case. Necessary before we can load
-    # the master file, otherwise we don't know which master to load
-    master.tags = cases[list(cases.keys())[0]]
-    master.loadmaster()
-
-    # load the original cases dict
-    post_dir = master.tags['[post_dir]']
-    FILE = open(post_dir + master.tags['[sim_id]'] + '.pkl', 'rb')
-    cases_orig = pickle.load(FILE)
-    FILE.close()
-
-    sim_nr = 0
-    sim_total = len(cases)
-    for case, casedict in cases.items():
-        sim_nr += 1
-
-        # set all the tags in the HtcMaster file
-        master.tags = casedict
-        # returns a dictionary with all the tags used for this
-        # specific case
-        htc = master.createcase(write_htc=write_htc)
-        #htc=master.createcase_check(cases_repo,write_htc=write_htc)
-
-        if not silent:
-            print('htc progress: ' + format(sim_nr, '3.0f') + '/' + \
-                   format(sim_total, '3.0f'))
-
-        if verbose:
-            print('===master.tags===\n', master.tags)
-
-        # make sure the current cases already exists, otherwise we are not
-        # relaunching!
-        if case not in cases_orig:
-            msg = 'relaunch only works for existing cases: %s' % case
-            raise KeyError(msg)
-
-        # save in the big cases. Note that values() gives a copy!
-        # remark, what about the copying done at the end of master.createcase?
-        # is that redundant then?
-        cases[list(htc.keys())[0]] = list(htc.values())[0]
-
-        if verbose:
-            print('created cases for: %s.htc\n' % master.tags['[case_id]'])
-
-    launch(cases, runmethod=runmethod, verbose=verbose, check_log=check_log,
-           copyback_turb=copyback_turb, silent=silent)
-
-    # update the original file: overwrite the newly set cases
-    FILE = open(post_dir + master.tags['[sim_id]'] + '.pkl', 'wb')
-    cases_orig.update(cases)
-    pickle.dump(cases_orig, FILE, protocol=2)
-    FILE.close()
-
-def prepare_launch_cases(cases, runmethod='gorm', verbose=False,write_htc=True,
-                         copyback_turb=True, silent=False, check_log=True,
-                         variable_tag_func=None, sim_id_new=None):
-    """
-    Same as prepare_launch, but now the input is just a cases object (cao).
-    If relaunching some earlier defined simulations, make sure to at least
-    rename the sim_id, otherwise it could become messy: things end up in the
-    same folder, sim_id post file get overwritten, ...
-
-    In case you do not use a variable_tag_fuc, make sure all your tags are
-    defined in cases. First and foremost, this means that the case_id does not
-    get updated to have a new sim_id, the path's are not updated, etc
-
-    When given a variable_tag_func, make sure it is properly
-    defined: do not base a variable tag's value on itself to avoid value chains
-
-    The master htc file will be loaded and alls tags defined in the cases dict
-    will be applied to it as is.
-    """
-
-    # initiate the HtcMaster object, load the master file
-    master = HtcMaster()
-    # for invariant tags, load random case. Necessary before we can load
-    # the master file, otherwise we don't know which master to load
-    master.tags = cases[list(cases.keys())[0]]
-    # load the master htc file as a string under the master.tags
-    master.loadmaster()
-    # create the execution folder structure and copy all data to it
-    # but reset to the correct launch dirs first
-    sim_id = master.tags['[sim_id]']
-    if runmethod in ['local', 'local-script', 'none']:
-        path = '/home/dave/PhD_data/HAWC2_results/ojf_post/%s/' % sim_id
-        master.tags['[run_dir]'] = path
-    elif runmethod == 'jess':
-        master.tags['[run_dir]'] = '/mnt/jess/HAWC2/ojf_post/%s/' % sim_id
-    elif runmethod == 'gorm':
-        master.tags['[run_dir]'] = '/mnt/gorm/HAWC2/ojf_post/%s/' % sim_id
-    else:
-        msg='unsupported runmethod, options: none, local, thyra, gorm, opt'
-        raise ValueError(msg)
-
-    master.create_run_dir()
-    master.copy_model_data()
-    # create the zip file
-    master.create_model_zip()
-
-    sim_nr = 0
-    sim_total = len(cases)
-
-    # for safety, create a new cases dict. At the end of the ride both cases
-    # and cases_new should be identical!
-    cases_new = {}
-
-    # cycle thourgh all the combinations
-    for case, casedict in cases.items():
-        sim_nr += 1
-
-        sim_id = casedict['[sim_id]']
-        # reset the launch dirs
-        if runmethod in ['local', 'local-script', 'none']:
-            path = '/home/dave/PhD_data/HAWC2_results/ojf_post/%s/' % sim_id
-            casedict['[run_dir]'] = path
-        elif runmethod == 'thyra':
-            casedict['[run_dir]'] = '/mnt/thyra/HAWC2/ojf_post/%s/' % sim_id
-        elif runmethod == 'gorm':
-            casedict['[run_dir]'] = '/mnt/gorm/HAWC2/ojf_post/%s/' % sim_id
-        else:
-            msg='unsupported runmethod, options: none, local, thyra, gorm, opt'
-            raise ValueError(msg)
-
-        # -----------------------------------------------------------
-        # set all the tags in the HtcMaster file
-        master.tags = casedict
-        # apply the variable tags if applicable
-        if variable_tag_func:
-            master = variable_tag_func(master)
-        elif sim_id_new:
-            # TODO: finish this
-            # replace all the sim_id occurences with the updated one
-            # this means also the case_id tag changes!
-            pass
-        # -----------------------------------------------------------
-
-        # returns a dictionary with all the tags used for this specific case
-        htc = master.createcase(write_htc=write_htc)
-
-        if not silent:
-            print('htc progress: ' + format(sim_nr, '3.0f') + '/' + \
-                   format(sim_total, '3.0f'))
-
-        if verbose:
-            print('===master.tags===\n', master.tags)
-
-        # make sure the current cases is unique!
-        if list(htc.keys())[0] in cases_new:
-            msg = 'non unique case in cases: %s' % list(htc.keys())[0]
-            raise KeyError(msg)
-        # save in the big cases. Note that values() gives a copy!
-        # remark, what about the copying done at the end of master.createcase?
-        # is that redundant then?
-        cases_new[list(htc.keys())[0]] = list(htc.values())[0]
-
-        if verbose:
-            print('created cases for: %s.htc\n' % master.tags['[case_id]'])
-
-    post_dir = master.tags['[post_dir]']
-
-    # create directory if post_dir does not exists
-    try:
-        os.makedirs(post_dir)
-    except OSError:
-        pass
-    FILE = open(post_dir + master.tags['[sim_id]'] + '.pkl', 'wb')
-    pickle.dump(cases_new, FILE, protocol=2)
-    FILE.close()
-
-    if not silent:
-        print('\ncases saved at:')
-        print(post_dir + master.tags['[sim_id]'] + '.pkl')
-
-    launch(cases_new, runmethod=runmethod, verbose=verbose,
-           copyback_turb=copyback_turb, check_log=check_log)
-
-    return cases_new
-
-
 def launch(cases, runmethod='none', verbose=False, copyback_turb=True,
            silent=False, check_log=True, windows_nr_cpus=2, qsub='time',
            pbs_fname_appendix=True, short_job_names=True,
-           maxcpu=1, pyenv='wetb_py3', wine_64bit=False, m=[3,4,6,8,9,10,12],
+           maxcpu=1, pyenv='py36-wetb', m=[3,4,6,8,9,10,12], prelude='',
            postpro_node_zipchunks=True, postpro_node=False, exesingle=None,
-           exechunks=None, wine_arch='win32', wine_prefix='~/.wine32'):
+           exechunks=None, wine_arch='win32', wine_prefix='~/.wine32',
+           pyenv_cmd='source /home/python/miniconda3/bin/activate'):
     """
     The actual launching of all cases in the Cases dictionary. Note that here
     only the PBS files are written and not the actuall htc files.
@@ -1049,11 +845,12 @@ def launch(cases, runmethod='none', verbose=False, copyback_turb=True,
         # create the pbs object
         pbs = PBS(cases, short_job_names=short_job_names, pyenv=pyenv,
                   pbs_fname_appendix=pbs_fname_appendix, qsub=qsub,
-                  verbose=verbose, silent=silent, wine_64bit=wine_64bit,
+                  verbose=verbose, silent=silent, prelude=prelude,
                   m=m, postpro_node_zipchunks=postpro_node_zipchunks,
                   postpro_node=postpro_node, exesingle=exesingle,
                   exechunks=exechunks, wine_arch=wine_arch,
                   wine_prefix=wine_prefix)
+        pbs.pyenv_cmd = pyenv_cmd
         pbs.copyback_turb = copyback_turb
         pbs.pbs_out_dir = pbs_out_dir
         pbs.maxcpu = maxcpu
@@ -1069,7 +866,6 @@ def launch(cases, runmethod='none', verbose=False, copyback_turb=True,
               'windows-script, local-ram, none, pbs'
         raise ValueError(msg)
 
-
 def post_launch(cases, save_iter=False, silent=False, suffix=None,
                 path_errorlog=None):
     """
@@ -1206,7 +1002,6 @@ def post_launch(cases, save_iter=False, silent=False, suffix=None,
 
     return cases_fail
 
-
 def copy_pbs_in_failedcases(cases_fail, path='pbs_in_fail', silent=True):
     """
     Copy all the pbs_in files from failed cases to a new directory so it
@@ -1230,7 +1025,6 @@ def copy_pbs_in_failedcases(cases_fail, path='pbs_in_fail', silent=True):
             os.makedirs(os.path.dirname(dst))
         shutil.copy2(src, dst)
 
-
 def logcheck_case(errorlogs, cases, case, silent=False):
     """
     Check logfile of a single case
@@ -1320,6 +1114,7 @@ class Log(object):
         for k in self.log:
             print(k)
 
+
 class HtcMaster(object):
     """
     """
@@ -1369,24 +1164,22 @@ class HtcMaster(object):
         self.tags['[iter_dir]'] = 'iter/'
         self.tags['[log_dir]'] = 'logfiles/'
         self.tags['[turb_dir]'] = 'turb/'
-        self.tags['[wake_dir]'] = None
-        self.tags['[meand_dir]'] = None
+        self.tags['[micro_dir]'] = None
+        self.tags['[meander_dir]'] = None
         self.tags['[turb_db_dir]'] = None
-        self.tags['[wake_db_dir]'] = None
-        self.tags['[meand_db_dir]'] = None
+        self.tags['[micro_db_dir]'] = None
+        self.tags['[meander_db_dir]'] = None
         self.tags['[control_dir]'] = 'control/'
         self.tags['[externalforce]'] = 'externalforce/'
         self.tags['[animation_dir]'] = 'animation/'
         self.tags['[eigenfreq_dir]'] = 'eigenfreq/'
-        self.tags['[wake_dir]'] = 'wake/'
-        self.tags['[meander_dir]'] = 'meander/'
         self.tags['[htc_dir]'] = 'htc/'
         self.tags['[mooring_dir]'] = 'mooring/'
         self.tags['[hydro_dir]'] = 'htc_hydro/'
         self.tags['[pbs_out_dir]'] = 'pbs_out/'
         self.tags['[turb_base_name]'] = None
-        self.tags['[wake_base_name]'] = None
-        self.tags['[meand_base_name]'] = None
+        self.tags['[micro_base_name]'] = None
+        self.tags['[meander_base_name]'] = None
         self.tags['[zip_root_files]'] = []
 
         self.tags['[fname_source]'] = []
@@ -1394,7 +1187,6 @@ class HtcMaster(object):
 
         self.tags['[eigen_analysis]'] = False
 
-        self.tags['[pbs_queue_command]'] = '#PBS -q workq'
         # the express que has 2 thyra nodes with max walltime of 1h
 #        self.tags['[pbs_queue_command]'] = '#PBS -q xpresq'
         # walltime should have following format: hh:mm:ss
@@ -1403,8 +1195,8 @@ class HtcMaster(object):
 #        self.queue = Queue.Queue()
 
         self.output_dirs = ['[res_dir]', '[log_dir]', '[turb_dir]',
-                            '[case_id]', '[wake_dir]', '[animation_dir]',
-                            '[meand_dir]', '[eigenfreq_dir]']
+                            '[case_id]', '[micro_dir]', '[animation_dir]',
+                            '[meander_dir]', '[eigenfreq_dir]']
 
     def create_run_dir(self):
         """
@@ -1413,13 +1205,15 @@ class HtcMaster(object):
 
         dirkeys = ['[data_dir]', '[htc_dir]', '[res_dir]', '[log_dir]',
                    '[eigenfreq_dir]', '[animation_dir]', '[turb_dir]',
-                   '[wake_dir]', '[meander_dir]', '[opt_dir]', '[control_dir]',
+                   '[micro_dir]', '[meander_dir]', '[opt_dir]', '[control_dir]',
                    '[mooring_dir]', '[hydro_dir]', '[externalforce]']
 
         # create all the necessary directories
         for dirkey in dirkeys:
             if isinstance(self.tags[dirkey], str):
                 path = os.path.join(self.tags['[run_dir]'], self.tags[dirkey])
+                if self.tags[dirkey].lower() == 'none':
+                    continue
                 if not os.path.exists(path):
                     os.makedirs(path)
 
@@ -1544,7 +1338,7 @@ class HtcMaster(object):
         htcmaster_dir = self.tags['[master_htc_dir]']
         data_dir = self.tags['[data_dir]']
         turb_dir = self.tags['[turb_dir]']
-        wake_dir = self.tags['[wake_dir]']
+        wake_dir = self.tags['[micro_dir]']
         meander_dir = self.tags['[meander_dir]']
         mooring_dir = self.tags['[mooring_dir]']
         hydro_dir = self.tags['[hydro_dir]']
@@ -1571,7 +1365,7 @@ class HtcMaster(object):
 
         # the master file
         src = os.path.join(htcmaster_dir, htcmaster)
-        dst = os.path.join('htc', '_master', htcmaster)
+        dst = os.path.join('htc', '_master', os.path.basename(htcmaster))
         zf.write(src, dst, zipfile.ZIP_DEFLATED)
 
         # manually add all that resides in control, mooring and hydro
@@ -1629,150 +1423,6 @@ class HtcMaster(object):
         #dst += self.tags['[model_zip]']
         #shutil.copy2(src, dst)
 
-    def _sweep_tags(self):
-        """
-        The original way with all tags in the htc file for each blade node
-        """
-        # set the correct sweep cruve, these values are used
-        a = self.tags['[sweep_amp]']
-        b = self.tags['[sweep_exp]']
-        z0 = self.tags['[sweep_curve_z0]']
-        ze = self.tags['[sweep_curve_ze]']
-        nr = self.tags['[nr_nodes_blade]']
-        # format for the x values in the htc file
-        ff = ' 1.03f'
-        for zz in range(nr):
-            it_nosweep = '[x'+str(zz+1)+'-nosweep]'
-            item = '[x'+str(zz+1)+']'
-            z = self.tags['[z'+str(zz+1)+']']
-            if z >= z0:
-                curve = eval(self.tags['[sweep_curve_def]'])
-                # new swept position = original + sweep curve
-                self.tags[item]=format(self.tags[it_nosweep]+curve,ff)
-            else:
-                self.tags[item]=format(self.tags[it_nosweep], ff)
-
-    def _staircase_windramp(self, nr_steps, wind_step, ramptime, septime):
-        """Create a stair case wind ramp
-
-
-        """
-
-        pass
-
-    def _all_in_one_blade_tag(self, radius_new=None):
-        """
-        Create htc input based on a HAWTOPT blade result file
-
-        Automatically get the number of nodes correct in master.tags based
-        on the number of blade nodes
-
-        WARNING: initial x position of the half chord point is assumed to be
-        zero
-
-        zaxis_fact : int, default=1.0 --> is member of default tags
-            Factor for the htc z-axis coordinates. The htc z axis is mapped to
-            the HAWTOPT radius. If the blade radius develops in negative z
-            direction, set to -1
-
-        Parameters
-        ----------
-
-        radius_new : ndarray(n), default=False
-            z coordinates of the nodes. If False, a linear distribution is
-            used and the tag [nr--of-nodes-per-blade] sets the number of nodes
-
-
-        """
-        # TODO: implement support for x position to be other than zero
-
-        # TODO: This is not a good place, should live somewhere else. Or
-        # reconsider inputs etc so there is more freedom in changing the
-        # location of the nodes, set initial x position of the blade etc
-
-        # and save under tag [blade_htc_node_input] in htc input format
-
-        nr_nodes = self.tags['[nr_nodes_blade]']
-
-        blade = self.tags['[blade_hawtopt]']
-        # in the htc file, blade root =0 and not blade hub radius
-        blade[:,0] = blade[:,0] - blade[0,0]
-
-        if type(radius_new).__name__ == 'NoneType':
-            # interpolate to the specified number of nodes
-            radius_new = np.linspace(blade[0,0], blade[-1,0], nr_nodes)
-
-        # Data checks on radius_new
-        elif not type(radius_new).__name__ == 'ndarray':
-            raise ValueError('radius_new has to be either NoneType or ndarray')
-        else:
-            if not len(radius_new.shape) == 1:
-                raise ValueError('radius_new has to be 1D')
-            elif not len(radius_new) == nr_nodes:
-                msg = 'radius_new has to have ' + str(nr_nodes) + ' elements'
-                raise ValueError(msg)
-
-        # save the nodal positions in the tag cloud
-        self.tags['[blade_nodes_z_positions]'] = radius_new
-
-        # make sure that radius_hr is just slightly smaller than radius low res
-        radius_new[-1] = blade[-1,0]-0.00000001
-        twist_new = interpolate.griddata(blade[:,0], blade[:,2], radius_new)
-        # blade_new is the htc node input part:
-        # sec 1   x     y     z   twist;
-        blade_new = scipy.zeros((len(radius_new),4))
-        blade_new[:,2] = radius_new*self.tags['[zaxis_fact]']
-        # twist angle remains the same in either case (standard/ojf rotation)
-        blade_new[:,3] = twist_new*-1.
-
-        # set the correct sweep cruve, these values are used
-        a = self.tags['[sweep_amp]']
-        b = self.tags['[sweep_exp]']
-        z0 = self.tags['[sweep_curve_z0]']
-        ze = self.tags['[sweep_curve_ze]']
-        tmp = 'nsec ' + str(nr_nodes) + ';'
-        for k in range(nr_nodes):
-            tmp += '\n'
-            i = k+1
-            z = blade_new[k,2]
-            y = blade_new[k,1]
-            twist = blade_new[k,3]
-            # x position, sweeping?
-            if z >= z0:
-                x = eval(self.tags['[sweep_curve_def]'])
-            else:
-                x = 0.0
-
-            # the node number
-            tmp += '        sec ' + format(i, '2.0f')
-            tmp += format(x, ' 11.03f')
-            tmp += format(y, ' 11.03f')
-            tmp += format(z, ' 11.03f')
-            tmp += format(twist, ' 11.03f')
-            tmp += ' ;'
-
-        self.tags['[blade_htc_node_input]'] = tmp
-
-        # and create the ae file
-        #5	Blade Radius [m] 	Chord[m]  T/C[%]  Set no. of pc file
-        #1 25 some comments
-        #0.000     0.100    21.000   1
-        nr_points = blade.shape[0]
-        tmp2 = '1  Blade Radius [m] Chord [m] T/C [%] pc file set nr\n'
-        tmp2 += '1  %i auto generated by _all_in_one_blade_tag()' % nr_points
-
-        for k in range(nr_points):
-            tmp2 += '\n'
-            tmp2 += '%9.3f %9.3f %9.3f' % (blade[k,0], blade[k,1], blade[k,3])
-            tmp2 += ' %4i' % (k+1)
-        # end with newline
-        tmp2 += '\n'
-
-        # TODO: finish writing file, implement proper handling of hawtopt path
-        # and save the file
-        #if self.tags['aefile']
-        #write_file(file_path, tmp2, 'w')
-
     def loadmaster(self):
         """
         Load the master file, path to master file is defined in
@@ -1795,8 +1445,7 @@ class HtcMaster(object):
         if not self.silent:
             print('loading master: ' + fpath)
 
-        with open(fpath, 'r') as f:
-            lines = f.readlines()
+        lines = misc.readlines_try_encodings(fpath)
 
         # regex for finding all tags in a line
         regex = re.compile('(\\[.*?\\])')
@@ -1867,7 +1516,14 @@ class HtcMaster(object):
             else:
                 value = self.tags[k]
             # if string is not found, it will do nothing
-            htc = htc.replace(str(k), str(value))
+            # case sensitive tag names
+#            htc = htc.replace(str(k), str(value))
+            # tag names are case insensitive
+            try:
+                htc = re.sub(re.escape(str(k)), str(value), htc, flags=re.IGNORECASE)
+            except re.error:
+                print('ignore replace tag, value: ', k, value)
+            # htc = re.sub(str(k), str(value), htc, flags=re.IGNORECASE)
 
         # and save the the case htc file:
         cname = self.tags['[case_id]'] + '.htc'
@@ -1940,11 +1596,11 @@ class PBS(object):
     such as the turbulence file and folder, htc folder and others
     """
 
-    def __init__(self, cases, qsub='time', silent=False, pyenv='wetb_py3',
+    def __init__(self, cases, qsub='time', silent=False, pyenv='py36-wetb',
                  pbs_fname_appendix=True, short_job_names=True, verbose=False,
-                 wine_64bit=False, m=[3,4,6,8,9,10,12], exesingle=None,
-                 postpro_node_zipchunks=True, postpro_node=False,
-                 exechunks=None, wine_arch='win32', wine_prefix='~/.wine32'):
+                 m=[3,4,6,8,9,10,12], exesingle=None, prelude='',
+                 postpro_node_zipchunks=True, postpro_node=False, queue='workq',
+                 exechunks=None, wine_arch='win32', wine_prefix='.wine32'):
         """
         Define the settings here. This should be done outside, but how?
         In a text file, paramters list or first create the object and than set
@@ -1976,43 +1632,49 @@ class PBS(object):
         self.verbose = verbose
         self.silent = silent
         self.pyenv = pyenv
-        self.pyenv_cmd = 'source /home/python/miniconda3/bin/activate'
+        self.pyenv_cmd = "source /home/python/miniconda3/bin/activate"
         self.postpro_node_zipchunks = postpro_node_zipchunks
         self.postpro_node = postpro_node
+        # Prelude is executed just before HAWC2 is
+        self.prelude = prelude
 
         self.m = m
 
+        wine_prefix = misc.sanitize_wine_prefix(wine_prefix)
+
         # run in 32-bit or 64-bit mode. Note this uses the same assumptions
         # on how to configure wine in toolbox/pbsutils/config-wine-hawc2.sh
         wineparam = (wine_arch, wine_prefix)
-        if wine_64bit:
-            wineparam = ('win64', '~/.wine')
-        self.winebase = 'time WINEARCH=%s WINEPREFIX=%s ' % wineparam
-
+        if wine_arch==None or wine_prefix==None:
+            self.winebase = "time "
+        else:
+            self.winebase = "time WINEARCH='%s' WINEPREFIX=\"%s\" " % wineparam
         self.wine = self.winebase + 'wine'
-        self.winenumactl = self.winebase + 'numactl --physcpubind=$CPU_NR wine'
+        self.winenumactl = self.winebase + "numactl --physcpubind=$CPU_NR wine"
 
         # in case you want to redirect stdout to /dev/nul, append as follows:
         # '> /dev/null 2>&1'
         self.exesingle = exesingle
         if exesingle is None:
-            self.exesingle = "{wine:} {hawc2_exe:} {fname_htc:}"
+            self.exesingle = "{wine:} '{hawc2_exe:}' '{fname_htc:}'"
         # in zipchunks mode we will output std out and err to a separate
         # pbs_out file, and that in addition to the pbs_out_zipchunks file
         self.exechunks = exechunks
         if exechunks is None:
-            self.exechunks = "({winenumactl:} {hawc2_exe:} {fname_htc:}) "
-            self.exechunks += "2>&1 | tee {fname_pbs_out:}"
+            self.exechunks = "({winenumactl:} '{hawc2_exe:}' '{fname_htc:}') "
+            self.exechunks += "2>&1 | tee '{fname_pbs_out:}'"
 
         # TODO: based on a certain host/architecture you can change these
         self.maxcpu = 1
         self.secperiter = 0.012
 
         # determine at runtime if winefix has to be ran
-        self.winefix = '  _HOSTNAME_=`hostname`\n'
-        self.winefix += '  if [[ ${_HOSTNAME_:0:1} == "j" ]] ; then\n'
-        self.winefix += '    WINEARCH=%s WINEPREFIX=%s winefix\n' % wineparam
-        self.winefix += '  fi\n'
+        self.winefix = ''
+        if wine_arch!=None or wine_prefix!=None:
+            self.winefix = "  _HOSTNAME_=`hostname`\n"
+            self.winefix += "  if [[ ${_HOSTNAME_:0:1} == \"j\" ]] ; then\n"
+            self.winefix += "    WINEARCH='%s' WINEPREFIX=\"%s\" winefix\n" % wineparam
+            self.winefix += "  fi\n"
 
         # the output channels comes with a price tag. Each time step
         # will have a penelty depending on the number of output channels
@@ -2034,6 +1696,12 @@ class PBS(object):
 #        self.node_run_root = '/dev/shm'
         self.node_run_root = '/scratch'
 
+        # only allow valid queues for jess
+        queues = ('windq', 'workq', 'fatq', 'windfatq', 'xpresq', 'xpresfatq')
+        if not queue in queues:
+            raise ValueError('Invalid queue name: %s' % queue)
+        self.pbs_queue_command = '#PBS -q %s' % queue
+
         self.cases = cases
 
         # location of the output messages .err and .out created by the node
@@ -2083,21 +1751,6 @@ class PBS(object):
         self.t0 = []
         # '[time_stop]' '[dt_sim]'
 
-        # REMARK: this i not realy consistent with how the result and log file
-        # dirs are allowed to change for each individual case...
-        # first check if the pbs_out_dir exists, this dir is considered to be
-        # the same for all cases present in cases
-        # self.tags['[run_dir]']
-        case0 = list(self.cases.keys())[0]
-        path = self.cases[case0]['[run_dir]'] + self.pbs_out_dir
-        if not os.path.exists(path):
-            os.makedirs(path)
-
-        # create pbs_in base dir
-        path = self.cases[case0]['[run_dir]'] + self.pbs_in_dir
-        if not os.path.exists(path):
-            os.makedirs(path)
-
         # number the pbs jobs:
         count2 = self.pbs_start_number
         # initial cpu count is zero
@@ -2111,6 +1764,9 @@ class PBS(object):
             # get a shorter version for the current cases tag_dict:
             tag_dict = self.cases[case]
 
+            # sanitize paths from a list of predefined tags in tag_dict
+            self.sanitize_paths(case, tag_dict)
+
             # group all values loaded from the tag_dict here, to keep overview
             # the directories to SAVE the results/logs/turb files
             # load all relevant dir settings: the result/logfile/turbulence/zip
@@ -2124,9 +1780,9 @@ class PBS(object):
             self.animation_dir = tag_dict['[animation_dir]']
             self.TurbDirName = tag_dict['[turb_dir]']
             self.TurbDb = tag_dict['[turb_db_dir]']
-            self.wakeDb = tag_dict['[wake_db_dir]']
-            self.meandDb = tag_dict['[meand_db_dir]']
-            self.WakeDirName = tag_dict['[wake_dir]']
+            self.wakeDb = tag_dict['[micro_db_dir]']
+            self.meandDb = tag_dict['[meander_db_dir]']
+            self.WakeDirName = tag_dict['[micro_dir]']
             self.MeanderDirName = tag_dict['[meander_dir]']
             self.ModelZipFile = tag_dict['[model_zip]']
             self.htc_dir = tag_dict['[htc_dir]']
@@ -2134,16 +1790,15 @@ class PBS(object):
             self.mooring_dir = tag_dict['[mooring_dir]']
             self.model_path = tag_dict['[run_dir]']
             self.turb_base_name = tag_dict['[turb_base_name]']
-            self.wake_base_name = tag_dict['[wake_base_name]']
-            self.meand_base_name = tag_dict['[meand_base_name]']
-            self.pbs_queue_command = tag_dict['[pbs_queue_command]']
+            self.wake_base_name = tag_dict['[micro_base_name]']
+            self.meand_base_name = tag_dict['[meander_base_name]']
             self.walltime = tag_dict['[walltime]']
             self.dyn_walltime = tag_dict['[auto_walltime]']
             self.case_duration = tag_dict['[duration]']
 
             # create the pbs_out_dir if necesary
             try:
-                path = tag_dict['[run_dir]'] + tag_dict['[pbs_out_dir]']
+                path = os.path.join(tag_dict['[run_dir]'], tag_dict['[pbs_out_dir]'])
                 if not os.path.exists(path):
                     os.makedirs(path)
                 self.pbs_out_dir = tag_dict['[pbs_out_dir]']
@@ -2152,7 +1807,7 @@ class PBS(object):
 
             # create pbs_in subdirectories if necessary
             try:
-                path = tag_dict['[run_dir]'] + tag_dict['[pbs_in_dir]']
+                path = os.path.join(tag_dict['[run_dir]'], tag_dict['[pbs_in_dir]'])
                 if not os.path.exists(path):
                     os.makedirs(path)
                 self.pbs_in_dir = tag_dict['[pbs_in_dir]']
@@ -2168,6 +1823,10 @@ class PBS(object):
             try:
                 self.copyto_generic = tag_dict['[copyto_generic]']
                 self.copyto_files = tag_dict['[copyto_files]']
+                if not isinstance(self.copyto_generic, list):
+                    raise ValueError('[copyto_generic] should be a list')
+                if not isinstance(self.copyto_files, list):
+                    raise ValueError('[copyto_files] should be a list')
             except KeyError:
                 pass
 
@@ -2181,6 +1840,10 @@ class PBS(object):
             try:
                 self.copyto_generic = [tag_dict['[copyto_generic_f1]']]
                 self.copyto_files = [tag_dict['[copyto_f1]']]
+                if not isinstance(self.copyto_generic, list):
+                    raise ValueError('[copyto_generic] should be a list')
+                if not isinstance(self.copyto_files, list):
+                    raise ValueError('[copyto_files] should be a list')
             except KeyError:
                 pass
 
@@ -2214,10 +1877,11 @@ class PBS(object):
                     jobid = tag_dict['[case_id]']
                 if self.pbs_fname_appendix and self.short_job_names:
                     # define the path for the new pbs script
-                    pbs_in_fname = '%s_%s.p' % (tag_dict['[case_id]'], jobid)
+                    pbs_in_fname = "%s_%s.p" % (tag_dict['[case_id]'], jobid)
                 else:
-                    pbs_in_fname = '%s.p' % (tag_dict['[case_id]'])
-                pbs_path = self.model_path + self.pbs_in_dir + pbs_in_fname
+                    pbs_in_fname = "%s.p" % (tag_dict['[case_id]'])
+                pbs_path = os.path.join(self.model_path, self.pbs_in_dir,
+                                        pbs_in_fname)
                 # Start a new pbs script, we only need the tag_dict here
                 self.starting(tag_dict, jobid)
                 ended = False
@@ -2228,16 +1892,16 @@ class PBS(object):
             # browse to the current scratch directory
             self.pbs += "\n\n"
             # mark start of single PBS mode
-            self.pbs += '# ' + '='*78 + '\n'
+            self.pbs += "# " + "="*78 + "\n"
             # evaluates to true if LAUNCH_PBS_MODE is NOT set
-            self.pbs += '# single PBS mode: one case per PBS job\n'
-            self.pbs += '# evaluates to true if LAUNCH_PBS_MODE is NOT set\n'
+            self.pbs += "# single PBS mode: one case per PBS job\n"
+            self.pbs += "# evaluates to true if LAUNCH_PBS_MODE is NOT set\n"
             self.pbs += "if [ -z ${LAUNCH_PBS_MODE+x} ] ; then\n"
             self.pbs += "  echo\n"
             self.pbs += "  echo 'Execute commands on scratch nodes'\n"
-            self.pbs += "  cd %s/$USER/$PBS_JOBID\n" % self.node_run_root
+            self.pbs += "  cd \"%s/$USER/$PBS_JOBID\"\n" % self.node_run_root
             self.pbs += "  # create unique dir for each CPU\n"
-            self.pbs += '  mkdir "%i"; cd "%i"\n' % (count1, count1)
+            self.pbs += "  mkdir '%i'; cd '%i'\n" % (count1, count1)
 
             # output the current scratch directory
             self.pbs += "  pwd\n"
@@ -2246,24 +1910,30 @@ class PBS(object):
             self.pbs += "  /usr/bin/unzip ../" + self.ModelZipFile + '\n'
             # create all directories, especially relevant if there are case
             # dependent sub directories that are not present in the ZIP file
-            self.pbs += "  mkdir -p " + self.htc_dir + '\n'
-            self.pbs += "  mkdir -p " + self.results_dir + '\n'
-            self.pbs += "  mkdir -p " + self.logs_dir + '\n'
-            if self.TurbDirName is not None or self.TurbDirName != 'None':
-                self.pbs += "  mkdir -p " + self.TurbDirName + '\n'
+            self.pbs += "  mkdir -p '" + self.htc_dir + "'\n"
+            self.pbs += "  mkdir -p '" + self.results_dir + "'\n"
+            self.pbs += "  mkdir -p '" + self.logs_dir + "'\n"
+            if self.TurbDirName is not None or self.TurbDirName.lower()!='none':
+                self.pbs += "  mkdir -p '" + self.TurbDirName + "'\n"
+
             if self.WakeDirName and self.WakeDirName != self.TurbDirName:
-                self.pbs += "  mkdir -p " + self.WakeDirName + '\n'
+                if str(self.WakeDirName).lower() != 'none':
+                    self.pbs += "  mkdir -p '" + self.WakeDirName + "'\n"
+
             if self.MeanderDirName and self.MeanderDirName != self.TurbDirName:
-                self.pbs += "  mkdir -p " + self.MeanderDirName + '\n'
+                if str(self.MeanderDirName).lower() != 'none':
+                    self.pbs += "  mkdir -p '" + self.MeanderDirName + "'\n"
+
             if self.hydro_dir:
-                self.pbs += "  mkdir -p " + self.hydro_dir + '\n'
+                self.pbs += "  mkdir -p '" + self.hydro_dir + "'\n"
             # create the eigen analysis dir just in case that is necessary
             if self.eigenfreq_dir:
-                self.pbs += '  mkdir -p %s\n' % self.eigenfreq_dir
+                self.pbs += "  mkdir -p '%s'\n" % self.eigenfreq_dir
 
             # and copy the htc file to the node
-            self.pbs += "  cp -R $PBS_O_WORKDIR/" + self.htc_dir \
-                + case +" ./" + self.htc_dir + '\n'
+            # use double quotes here because we need to expand $PBS_O_WORKDIR
+            case_source = os.path.join('$PBS_O_WORKDIR', self.htc_dir, case)
+            self.pbs += "  cp -R \"%s\" '%s'\n" % (case_source, self.htc_dir)
 
             # if there is a turbulence file data base dir, copy from there
             if self.TurbDb:
@@ -2275,16 +1945,17 @@ class PBS(object):
             # names: turb_base_name_xxx_u.bin, turb_base_name_xxx_v.bin
             if self.turb_base_name is not None:
                 turb_src = os.path.join(turb_dir_src, self.turb_base_name)
-                self.pbs += "  cp -R %s*.bin %s\n" % (turb_src, self.TurbDirName)
+                # double quotes for expanding variables at run time
+                self.pbs += "  cp -R \"%s\"*.bin \"%s\"\n" % (turb_src, self.TurbDirName)
             # more generally, literally define the names of the boxes for u,v,w
             # components
             elif '[turb_fname_u]' in tag_dict:
                 turb_u = os.path.join(turb_dir_src, tag_dict['[turb_fname_u]'])
                 turb_v = os.path.join(turb_dir_src, tag_dict['[turb_fname_v]'])
                 turb_w = os.path.join(turb_dir_src, tag_dict['[turb_fname_w]'])
-                self.pbs += "  cp %s %s\n" % (turb_u, self.TurbDirName)
-                self.pbs += "  cp %s %s\n" % (turb_v, self.TurbDirName)
-                self.pbs += "  cp %s %s\n" % (turb_w, self.TurbDirName)
+                self.pbs += "  cp \"%s\" \"%s\"\n" % (turb_u, self.TurbDirName)
+                self.pbs += "  cp \"%s\" \"%s\"\n" % (turb_v, self.TurbDirName)
+                self.pbs += "  cp \"%s\" \"%s\"\n" % (turb_w, self.TurbDirName)
 
             # if there is a turbulence file data base dir, copy from there
             if self.wakeDb and self.WakeDirName:
@@ -2293,7 +1964,7 @@ class PBS(object):
                 wake_dir_src = os.path.join('$PBS_O_WORKDIR', self.WakeDirName)
             if self.wake_base_name is not None:
                 wake_src = os.path.join(wake_dir_src, self.wake_base_name)
-                self.pbs += "  cp -R %s*.bin %s\n" % (wake_src, self.WakeDirName)
+                self.pbs += "  cp -R \"%s\"*.bin \"%s\"\n" % (wake_src, self.WakeDirName)
 
             # if there is a turbulence file data base dir, copy from there
             if self.meandDb and self.MeanderDirName:
@@ -2302,55 +1973,53 @@ class PBS(object):
                 meand_dir_src = os.path.join('$PBS_O_WORKDIR', self.MeanderDirName)
             if self.meand_base_name is not None:
                 meand_src = os.path.join(meand_dir_src, self.meand_base_name)
-                self.pbs += "  cp -R %s*.bin %s\n" % (meand_src, self.MeanderDirName)
+                self.pbs += "  cp -R \"%s\"*.bin \"%s\"\n" % (meand_src, self.MeanderDirName)
 
             # copy and rename input files with given versioned name to the
             # required non unique generic version
             for fname, fgen in zip(self.copyto_files, self.copyto_generic):
-                self.pbs += "  cp -R $PBS_O_WORKDIR/%s ./%s\n" % (fname, fgen)
+                fname_full = os.path.join('$PBS_O_WORKDIR', fname)
+                self.pbs += "  cp -R \"%s\" \"%s\"\n" % (fname_full, fgen)
 
             # only apply the wine fix in PBS mode
             self.pbs += self.winefix
             # TODO: activate python env, calculate post-processing
 #            self.pbs += 'echo `python -c "import wetb; print(wetb.__version__)"`\n'
             # mark end of single PBS mode
-            self.pbs += '# ' + '='*78 + '\n\n'
+            self.pbs += "# " + "="*78 + "\n\n"
 
             # end of the file copying in PBS mode
             # mark start of find+xargs mode
-            self.pbs += '# ' + '-'*78 + '\n'
-            self.pbs += '# find+xargs mode: 1 PBS job, multiple cases\n'
+            self.pbs += "# " + "-"*78 + "\n"
+            self.pbs += "# find+xargs mode: 1 PBS job, multiple cases\n"
             self.pbs += "else\n"
             # when in find+xargs mode, browse to the relevant CPU
-            self.pbs += '  # with find+xargs we first browse to CPU folder\n'
-            self.pbs += '  cd "$CPU_NR"\n'
+            self.pbs += "  # with find+xargs we first browse to CPU folder\n"
+            self.pbs += "  cd \"$CPU_NR\"\n"
             self.pbs += "fi\n"
             # mark end of find+xargs mode
-            self.pbs += '# ' + '-'*78 + '\n\n'
+            self.pbs += "# " + "-"*78 + "\n\n"
 
-            self.pbs += 'echo ""\n'
+            self.pbs += "echo ''\n"
             # mark start of single PBS mode
-            self.pbs += '# ' + '='*78 + '\n'
-            self.pbs += '# single PBS mode: one case per PBS job\n'
-            self.pbs += '# evaluates to true if LAUNCH_PBS_MODE is NOT set\n'
+            self.pbs += "# " + "="*78 + "\n"
+            self.pbs += "# single PBS mode: one case per PBS job\n"
+            self.pbs += "# evaluates to true if LAUNCH_PBS_MODE is NOT set\n"
             self.pbs += "if [ -z ${LAUNCH_PBS_MODE+x} ] ; then\n"
 
-            self.pbs += '  echo "execute HAWC2, fork to background"\n'
+            self.pbs += "  echo 'execute HAWC2, fork to background'\n"
             # the hawc2 execution commands via wine, in PBS mode fork and wait
             # METHOD MORE GENERAL
             # case contains the htc file name extension, self.case doesn't
-            fname_htc = "./" + os.path.join(self.htc_dir, case)
+            # we include a leading ./ to help HAWC2
+            fname_htc = os.path.join(self.htc_dir, case)
             fname_log = os.path.join(self.logs_dir, self.case)
             ext = '.err.out'
             fname_pbs_out = os.path.join(self.pbs_out_dir, self.case + ext)
             execstr = self.exesingle.format(wine=self.wine, case=case,
                                             fname_htc=fname_htc,
-                                            hawc2_exe=hawc2_exe,
-                                            pbs_out_dir=self.pbs_out_dir,
-                                            logs_dir=self.logs_dir,
-                                            fname_log=fname_log,
-                                            fname_pbs_out=fname_pbs_out,
-                                            winenumactl=self.winenumactl)
+                                            hawc2_exe=hawc2_exe)
+            self.pbs += self.prelude
             self.pbs += "  %s &\n" % execstr
             # # OLD METHOD
             # param = (self.wine, hawc2_exe, self.htc_dir+case)
@@ -2360,24 +2029,24 @@ class PBS(object):
             # per PBS file, otherwise you have to wait for each job to finish
             # first and then run the post-processing for all those cases
             if self.maxcpu == 1:
-                self.pbs += '  wait\n'
+                self.pbs += "  wait\n"
                 if self.pyenv is not None and self.postpro_node:
-                    self.pbs += '  echo "POST-PROCESSING"\n'
-                    self.pbs += '  %s %s\n' % (self.pyenv_cmd, self.pyenv)
+                    self.pbs += "  echo 'POST-PROCESSING'\n"
+                    self.pbs += "  %s %s\n" % (self.pyenv_cmd, self.pyenv)
                     self.pbs += "  "
                     self.checklogs()
                     self.pbs += "  "
                     self.postprocessing()
-                    self.pbs += '  source deactivate\n'
+                    self.pbs += "  source deactivate\n"
             # mark end of single PBS mode
-            self.pbs += '# ' + '='*78 + '\n\n'
+            self.pbs += "# " + "="*78 + "\n\n"
             # mark start of find+xargs mode
-            self.pbs += '# ' + '-'*78 + '\n'
-            self.pbs += '# find+xargs mode: 1 PBS job, multiple cases\n'
+            self.pbs += "# " + "-"*78 + "\n"
+            self.pbs += "# find+xargs mode: 1 PBS job, multiple cases\n"
             self.pbs += "else\n"
             # numactl --physcpubind=$CPU_NR
 
-            fname_htc = "./" + os.path.join(self.htc_dir, case)
+            fname_htc = os.path.join(self.htc_dir, case)
             fname_log = os.path.join(self.logs_dir, self.case)
             ext = '.err.out'
             fname_pbs_out = os.path.join(self.pbs_out_dir, self.case + ext)
@@ -2389,22 +2058,22 @@ class PBS(object):
                                             fname_log=fname_log,
                                             fname_pbs_out=fname_pbs_out,
                                             winenumactl=self.winenumactl)
-            self.pbs += '  echo "execute HAWC2, do not fork and wait"\n'
-            self.pbs += "  %s \n" % execstr
+            self.pbs += "  echo 'execute HAWC2, do not fork and wait'\n"
+            self.pbs += "  %s\n" % execstr
 
             # param = (self.winenumactl, hawc2_exe, self.htc_dir+case,
             #          self.wine_appendix)
             # self.pbs += '  echo "execute HAWC2, do not fork and wait"\n'
             # self.pbs += "  " + ("%s %s ./%s %s" % param).strip() + "\n"
             if self.pyenv is not None and self.postpro_node_zipchunks:
-                self.pbs += '  echo "POST-PROCESSING"\n'
+                self.pbs += "  echo 'POST-PROCESSING'\n"
                 self.pbs += "  "
                 self.checklogs()
                 self.pbs += "  "
                 self.postprocessing()
             self.pbs += "fi\n"
             # mark end of find+xargs mode
-            self.pbs += '# ' + '-'*78 + '\n'
+            self.pbs += "# " + "-"*78 + "\n"
 
             #self.pbs += "wine get_mac_adresses" + '\n'
             # self.pbs += "cp -R ./*.mac  $PBS_O_WORKDIR/." + '\n'
@@ -2444,28 +2113,10 @@ class PBS(object):
         """
 
         # a new clean pbs script!
-        self.pbs = ''
-        self.pbs += "### Standard Output" + '\n'
 
         case_id = tag_dict['[case_id]']
-
-        # PBS job name
-        self.pbs += "#PBS -N %s\n" % (jobid)
-        self.pbs += "#PBS -o ./" + self.pbs_out_dir + case_id + ".out" + '\n'
-        # self.pbs += "#PBS -o ./pbs_out/" + jobid + ".out" + '\n'
-        self.pbs += "### Standard Error" + '\n'
-        self.pbs += "#PBS -e ./" + self.pbs_out_dir + case_id + ".err" + '\n'
-        # self.pbs += "#PBS -e ./pbs_out/" + jobid + ".err" + '\n'
-        self.pbs += '#PBS -W umask=0003\n'
-        self.pbs += "### Maximum wallclock time format HOURS:MINUTES:SECONDS\n"
-#        self.pbs += "#PBS -l walltime=" + self.walltime + '\n'
-        self.pbs += "#PBS -l walltime=[walltime]\n"
-        if self.qsub == 'time':
-            self.pbs += "#PBS -a [start_time]" + '\n'
-        elif self.qsub == 'depend':
-            # set job dependencies, job_id refers to PBS job_id, which is only
-            # assigned to a job at the moment it gets qsubbed into the que
-            self.pbs += "[nodeps]PBS -W depend=afterany:[job_id]\n"
+        std_out = os.path.join(self.pbs_out_dir, case_id + ".out")
+        std_err = os.path.join(self.pbs_out_dir, case_id + ".err")
 
 #        if self.que_jobdeps:
 #            self.pbs += "#PBS -W depend=afterany:%s\n" % jobid_dep
@@ -2482,18 +2133,49 @@ class PBS(object):
             # Number of nodes and cpus per node (ppn)
             lnodes = int(math.ceil(len(self.cases)/float(self.maxcpu)))
             lnodes = 1
-            self.pbs += "#PBS -l nodes=%i:ppn=%i\n" % (lnodes, self.maxcpu)
+            ppn = self.maxcpu
         else:
-            self.pbs += "#PBS -l nodes=1:ppn=1\n"
+            ppn = 1
+            lnodes = 1
             # Number of nodes and cpus per node (ppn)
 
-        self.pbs += "### Queue name" + '\n'
-        # queue names for Thyra are as follows:
-        # short walltime queue (shorter than an hour): '#PBS -q xpresq'
-        # or otherwise for longer jobs: '#PBS -q workq'
-        self.pbs += self.pbs_queue_command + '\n'
+        header_jess = """### Standard Output
+#PBS -N {jobid}
+#PBS -o '{std_out}'
+### Standard Error
+#PBS -e '{std_err}'
+#PBS -W umask=0003
+### Maximum wallclock time format HOURS:MINUTES:SECONDS
+#PBS -l walltime=[walltime]
+#PBS -l nodes={lnodes}:ppn={ppn}
+### Queue name
+{pbs_queue_command}
+"""
+
+        header_sophia = """#!/bin/bash
+#==============================================================================
+# ONLY RELEVANT/USED FOR/BY SLURM ON SOPHIA
+#==============================================================================
+#SBATCH --job-name={jobid}
+#SBATCH --output={std_out}
+#SBATCH --error={std_err}
+#SBATCH --nodes={lnodes}
+#SBATCH --ntasks=4
+#SBATCH --time=[walltime]
+#==============================================================================
+"""
+
+        self.pbs = header_jess.format(jobid=jobid, std_out=std_out,
+                                      std_err=std_err, ppn=ppn,
+                                      pbs_queue_command=self.pbs_queue_command,
+                                      lnodes=lnodes)
+        # self.pbs += header_sophia.format(jobid=jobid, std_out=std_out,
+        #                               std_err=std_err, ppn=ppn,
+        #                               pbs_queue_command=self.pbs_queue_command,
+        #                               lnodes=lnodes)
+
         # mark start of single PBS mode
-        self.pbs += '\n' + '# ' + '='*78 + '\n'
+        self.pbs += "\n" + "# " + "="*78 + "\n"
 
         # ignore all the file copying when running in xargs mode:
         # when varibale LAUNCH_PBS_MODE is set, file copying is ignored
@@ -2501,29 +2183,26 @@ class PBS(object):
         # we do this so the same launch script can be used either with the node
         # scheduler and the PBS system (for example when re-running cases)
         # evaluates to true if LAUNCH_PBS_MODE is NOT set
-        self.pbs += '# single PBS mode: one case per PBS job\n'
-        self.pbs += '# evaluates to true if LAUNCH_PBS_MODE is NOT set\n'
+        self.pbs += "# single PBS mode: one case per PBS job\n"
+        self.pbs += "# evaluates to true if LAUNCH_PBS_MODE is NOT set\n"
         self.pbs += "if [ -z ${LAUNCH_PBS_MODE+x} ] ; then\n"
 
         self.pbs += "  ### Create scratch directory and copy data to it\n"
         # output the current directory
-        self.pbs += "  cd $PBS_O_WORKDIR" + '\n'
-        self.pbs += '  echo "current working dir (pwd):"\n'
+        self.pbs += "  cd \"$PBS_O_WORKDIR\"\n"
+        self.pbs += "  echo 'current working dir (pwd):'\n"
         self.pbs += "  pwd\n"
-        # The batch system on Gorm allows more than one job per node.
-        # Because of this the scratch directory name includes both the
-        # user name and the job ID, that is /scratch/$USER/$PBS_JOBID
-        # if not scratch, make the dir
-        if self.node_run_root != '/scratch':
-            self.pbs += '  mkdir -p %s/$USER\n' % self.node_run_root
-            self.pbs += '  mkdir -p %s/$USER/$PBS_JOBID\n' % self.node_run_root
+        # # in very special cases you might want the scratch directory be made
+        # if self.node_run_root != '/scratch':
+        #     self.pbs += '  mkdir -p %s/$USER\n' % self.node_run_root
+        #     self.pbs += '  mkdir -p %s/$USER/$PBS_JOBID\n' % self.node_run_root
 
         # copy the zip files to the scratch dir on the node
-        self.pbs += "  cp -R ./" + self.ModelZipFile + \
-            ' %s/$USER/$PBS_JOBID\n' % (self.node_run_root)
+        rpl = (self.ModelZipFile, self.node_run_root)
+        self.pbs += "  cp -R '%s' \"%s/$USER/$PBS_JOBID\"\n" % rpl
         self.pbs += "fi\n"
         # mark end of single PBS mode
-        self.pbs += '# ' + '='*78 + '\n'
+        self.pbs += "# " + "="*78 + "\n"
 
     def ending(self, pbs_path):
         """
@@ -2533,15 +2212,15 @@ class PBS(object):
         self.pbs += "\n\n"
         self.pbs += "### Epilogue\n"
         # mark start of single PBS mode
-        self.pbs += '# ' + "="*78 + "\n"
+        self.pbs += "# " + "="*78 + "\n"
         # evaluates to true if LAUNCH_PBS_MODE is NOT set
-        self.pbs += '# single PBS mode: one case per PBS job\n'
-        self.pbs += '# evaluates to true if LAUNCH_PBS_MODE is NOT set\n'
+        self.pbs += "# single PBS mode: one case per PBS job\n"
+        self.pbs += "# evaluates to true if LAUNCH_PBS_MODE is NOT set\n"
         self.pbs += "if [ -z ${LAUNCH_PBS_MODE+x} ] ; then\n"
         self.pbs += "  ### wait for jobs to finish\n"
         self.pbs += "  wait\n"
-        self.pbs += '  echo ""\n'
-        self.pbs += '  echo "Copy back from scratch directory"\n'
+        self.pbs += "  echo ''\n"
+        self.pbs += "  echo 'Copy back from scratch directory'\n"
         for i in range(1, self.maxcpu+1, 1):
 
             # navigate to the cpu dir on the node
@@ -2553,23 +2232,11 @@ class PBS(object):
             # for this mode is handled elsewhere
             if self.maxcpu == 1:
                 # mark start of find+xargs mode
-                self.pbs += '# ' + "-"*78 + "\n"
-                self.pbs += '# find+xargs mode: 1 PBS job, multiple cases\n'
-                self.pbs += 'else\n'
+                self.pbs += "# " + "-"*78 + "\n"
+                self.pbs += "# find+xargs mode: 1 PBS job, multiple cases\n"
+                self.pbs += "else\n"
                 self.copyback_all_files("find+xargs", None)
-                # mark end of find+xargs mode
-#                self.pbs += '# ' + "-"*78 + "\n"
-
-#            # and delete it all (but that is not allowed)
-#            self.pbs += 'cd ..\n'
-#            self.pbs += 'ls -lah\n'
-#            self.pbs += 'echo $PBS_JOBID\n'
-#            self.pbs += 'rm -r $PBS_JOBID\n'
-
-            # Delete the batch file at the end. However, is this possible since
-            # the batch file is still open at this point????
-            # self.pbs += "rm "
-        self.pbs += 'fi\n'
+        self.pbs += "fi\n"
 
         # base walltime on the longest simulation in the batch
         nr_time_steps = max(self.nr_time_steps)
@@ -2599,7 +2266,10 @@ class PBS(object):
         self.pbs += "exit\n"
 
         if self.verbose:
-            print('writing pbs script to path: ' + pbs_path)
+            print("writing pbs script to path: " + pbs_path)
+
+        if self.pbs.find('rm -r') > -1 or self.pbs.find('rm -f') > -1:
+            raise UserWarning('Anything that looks like rm -rf is prohibited.')
 
         # and write the script to a file:
         write_file(pbs_path, self.pbs, 'w')
@@ -2611,20 +2281,22 @@ class PBS(object):
         or from CPU sub-directory back to main directory in find+xargs mode.
         """
         if mode=="find+xargs":
-            foper = "rsync -a --remove-source-files" # move files instead of copy
+            # we use --remove-source-files otherwise all the results will be
+            # here twice: in the main directory and in the cpu-sub directories
+            foper = "rsync -a --remove-source-files"
             dst = os.path.join('..', self.sim_id, '')
             # copy back to DB dir, and not the scratch dir root
-            dst_db = "$PBS_O_WORKDIR/"
-            cd2model = "  cd %s\n" % os.path.join(self.node_run_root, '$USER',
-                                                  '$PBS_JOBID', '$CPU_NR', '')
+            dst_db = "$PBS_O_WORKDIR"
+            cd2model = "  cd \"%s\"\n" % os.path.join(self.node_run_root, '$USER',
+                                                      '$PBS_JOBID', '$CPU_NR', '')
             pbs_mode = False
         else:
             foper = "cp -R"
-            dst = "$PBS_O_WORKDIR/"
+            dst = "$PBS_O_WORKDIR"
             dst_db = dst
             pbs_mode = True
-            cd2model = "  cd %s\n" % os.path.join(self.node_run_root, '$USER',
-                                                  '$PBS_JOBID', '%i' % cpu_nr, '')
+            cd2model = "  cd \"%s\"\n" % os.path.join(self.node_run_root, '$USER',
+                                                      '$PBS_JOBID', '%i' % cpu_nr, '')
 
         # navigate to the cpu dir on the node
         # The batch system on Gorm/Jess allows more than one job per node.
@@ -2635,121 +2307,129 @@ class PBS(object):
         # create the log, res etc dirs in case they do not exist. Only relevant
         # for pbs_mode, they are created in advance in find+xargs
         if pbs_mode:
-            mk = '  mkdir -p'
-            self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.results_dir))
-            self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.logs_dir))
+            mk = "  mkdir -p"
+            self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.results_dir))
+            self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.logs_dir))
             if self.animation_dir:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.animation_dir))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.animation_dir))
             if self.copyback_turb and self.TurbDb:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.TurbDb))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.TurbDb))
             elif self.copyback_turb:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.TurbDirName))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.TurbDirName))
             if self.copyback_turb and self.wakeDb:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.wakeDb))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.wakeDb))
             elif self.WakeDirName and self.WakeDirName != self.TurbDirName:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.WakeDirName))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.WakeDirName))
             if self.copyback_turb and self.meandDb:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.meandDb))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.meandDb))
             elif self.MeanderDirName and self.MeanderDirName != self.TurbDirName:
-                self.pbs += "%s %s\n" % (mk, os.path.join(dst, self.MeanderDirName))
+                self.pbs += "%s \"%s\"\n" % (mk, os.path.join(dst, self.MeanderDirName))
 
         # and copy the results and log files frome the scratch to dst
         res_dst = os.path.join(dst, self.results_dir, ".")
-        self.pbs += "  %s %s. %s\n" % (foper, self.results_dir, res_dst)
+        # strict on the first path (single quote), allow expanding on the
+        # second path (contains variables)
+        self.pbs += "  %s '%s.' \"%s\"\n" % (foper, self.results_dir, res_dst)
         log_dst = os.path.join(dst, self.logs_dir, ".")
-        self.pbs += "  %s %s. %s\n" % (foper, self.logs_dir, log_dst)
+        self.pbs += "  %s '%s.' \"%s\"\n" % (foper, self.logs_dir, log_dst)
         # in zipchunks mode by default we also copy the time+std out/err to
         # an additional file that is in pbs_out for consistancy with the
         # pbs_mode approach
         if not pbs_mode:
             pbs_dst = os.path.join(dst, self.pbs_out_dir, ".")
-            self.pbs += "  %s %s. %s\n" % (foper, self.pbs_out_dir, pbs_dst)
+            self.pbs += "  %s '%s.' \"%s\"\n" % (foper, self.pbs_out_dir, pbs_dst)
         if self.animation_dir:
             ani_dst = os.path.join(dst, self.animation_dir, ".")
-            self.pbs += "  %s %s. %s\n" % (foper, self.animation_dir, ani_dst)
+            self.pbs += "  %s '%s.' \"%s\"\n" % (foper, self.animation_dir, ani_dst)
 
         if self.eigenfreq_dir:
             # just in case the eig dir has subdirs for the results, only
             # select the base path and cp -r will take care of the rest
             p1 = self.eigenfreq_dir.split('/')[0]
             p2 = os.path.join(dst, p1, ".")
-            self.pbs += "  cp -R %s/. %s\n" % (p1, p2)
+            self.pbs += "  cp -R '%s/.' \"%s\"\n" % (p1, p2)
             # for eigen analysis with floater, modes are in root
             eig_dir_sys = os.path.join(dst, self.eigenfreq_dir, 'system/', '.')
-            self.pbs += '  mkdir -p %s\n' % eig_dir_sys
-            self.pbs += "  cp -R mode* %s\n" % eig_dir_sys
-            self.pbs += "  %s mode* %s\n" % (foper, eig_dir_sys)
+            self.pbs += "  mkdir -p \"%s\"\n" % eig_dir_sys
+            self.pbs += "  cp -R mode* \"%s\"\n" % eig_dir_sys
+            self.pbs += "  %s mode* \"%s\"\n" % (foper, eig_dir_sys)
 
         # only copy the turbulence files back if they do not exist
         # for all *.bin files on the node
-        cmd = '  for i in `ls *.bin`; do  if [ -e %s$i ]; '
+        cmd = '  for i in `ls *.bin`; do  if [ -e "%s" ]; '
         cmd += 'then echo "$i exists no copyback"; else echo "$i copyback"; '
-        cmd += 'cp $i %s; fi; done\n'
+        cmd += 'cp "$i" "%s"; fi; done\n'
         # copy back turbulence file?
         # browse to the node turb dir
-        self.pbs += '\n  echo ""\n'
-        self.pbs += '  echo "COPY BACK TURB IF APPLICABLE"\n'
+        self.pbs += "\n  echo ''\n"
+        self.pbs += "  echo 'COPY BACK TURB IF APPLICABLE'\n"
         if self.copyback_turb and self.TurbDb:
-            self.pbs += '  cd %s\n' % self.TurbDirName
-            tmp = (os.path.join(dst_db, self.TurbDb, ''),)*2
+            self.pbs += "  cd '%s'\n" % self.TurbDirName
+            tmp = (os.path.join(dst_db, self.TurbDb, '$i'),
+                   os.path.join(dst_db, self.TurbDb, ''))
             self.pbs += cmd % tmp
             # and back to normal model root
             self.pbs += cd2model
         elif self.copyback_turb:
-            self.pbs += '  cd %s\n' % self.TurbDirName
-            tmp = (os.path.join(dst, self.TurbDirName, ''),)*2
+            self.pbs += "  cd '%s'\n" % self.TurbDirName
+            tmp = (os.path.join(dst, self.TurbDirName, '$i'),
+                   os.path.join(dst, self.TurbDirName, ''))
             self.pbs += cmd % tmp
             # and back to normal model root
             self.pbs += cd2model
 
         if self.copyback_turb and self.wakeDb:
-            self.pbs += '  cd %s\n' % self.WakeDirName
-            tmp = (os.path.join(dst_db, self.wakeDb, ''),)*2
+            self.pbs += "  cd '%s'\n" % self.WakeDirName
+            tmp = (os.path.join(dst_db, self.wakeDb, '$i'),
+                   os.path.join(dst_db, self.wakeDb, ''))
             self.pbs += cmd % tmp
             # and back to normal model root
             self.pbs += cd2model
         elif self.copyback_turb and self.WakeDirName:
-            self.pbs += '  cd %s\n' % self.WakeDirName
-            tmp = (os.path.join(dst, self.WakeDirName, ''),)*2
+            self.pbs += "  cd '%s'\n" % self.WakeDirName
+            tmp = (os.path.join(dst, self.WakeDirName, '$i'),
+                   os.path.join(dst, self.WakeDirName, ''))
             self.pbs += cmd % tmp
             # and back to normal model root
             self.pbs += cd2model
 
         if self.copyback_turb and self.meandDb:
-            self.pbs += '  cd %s\n' % self.MeanderDirName
-            tmp = (os.path.join(dst_db, self.meandDb, ''),)*2
+            self.pbs += "  cd '%s'\n" % self.MeanderDirName
+            tmp = (os.path.join(dst_db, self.meandDb, '$i'),
+                   os.path.join(dst_db, self.meandDb, ''))
             self.pbs += cmd % tmp
             # and back to normal model root
             self.pbs += cd2model
         elif self.copyback_turb and self.MeanderDirName:
-            self.pbs += '  cd %s\n' % self.MeanderDirName
-            tmp = (os.path.join(dst, self.MeanderDirName, ''),)*2
+            self.pbs += "  cd '%s'\n" % self.MeanderDirName
+            tmp = (os.path.join(dst, self.MeanderDirName, '$i'),
+                   os.path.join(dst, self.MeanderDirName, ''))
             self.pbs += cmd % tmp
             # and back to normal model root
             self.pbs += cd2model
 
-        self.pbs += '  echo "END COPY BACK TURB"\n'
-        self.pbs += '  echo ""\n\n'
+        self.pbs += "  echo 'END COPY BACK TURB'\n"
+        self.pbs += "  echo ''\n\n"
 
         # copy back any other kind of files, as specified in copyback_files
-        self.pbs += '  echo "COPYBACK [copyback_files]/[copyback_frename]"\n'
+        self.pbs += "  echo 'COPYBACK [copyback_files]/[copyback_frename]'\n"
         if len(self.copyback_frename) == 0:
             self.copyback_frename = self.copyback_files
         for fname, fnew in zip(self.copyback_files, self.copyback_frename):
             dst_fnew = os.path.join(dst, fnew)
-            self.pbs += "  %s %s %s\n" % (foper, fname, dst_fnew)
-        self.pbs += '  echo "END COPYBACK"\n'
-        self.pbs += '  echo ""\n'
+            self.pbs += "  %s '%s' \"%s\"\n" % (foper, fname, dst_fnew)
+        self.pbs += "  echo 'END COPYBACK'\n"
+        self.pbs += "  echo ''\n"
 
         if pbs_mode:
             # check what is left
-            self.pbs += '  echo ""\n'
-            self.pbs += '  echo "following files are on '
-            self.pbs += 'node/cpu %i (find .):"\n' % cpu_nr
-            self.pbs += '  find .\n'
-            self.pbs += '# ' + '='*78 + '\n'
+            self.pbs += "  echo ''\n"
+            self.pbs += "  echo 'following files are on "
+            self.pbs += "node/cpu %i (find .):'\n" % cpu_nr
+            self.pbs += "  find .\n"
+            self.pbs += "# " + "="*78 + "\n"
         else:
-            self.pbs += '# ' + '-'*78 + '\n'
+            self.pbs += "# " + "-"*78 + "\n"
 
     def checklogs(self):
         """
@@ -2797,7 +2477,11 @@ class PBS(object):
             except OSError:
                 size_sel = 0
                 size_dat = 0
-            if size_sel < 5 or size_dat < 5:
+            try:
+                size_hdf = os.stat(f_res + '.hdf5').st_size
+            except OSError:
+                size_hdf = 0
+            if (size_sel < 5 or size_dat < 5) and size_hdf < 5:
                 cases_fail[cname] = copy.copy(cases[cname])
 
         if not self.silent:
@@ -2806,6 +2490,53 @@ class PBS(object):
         # length will be zero if there are no failures
         return cases_fail
 
+    def sanitize_paths(self, case, tag_dict):
+        """Do some checks on the user defined paths
+        """
+        dirs = ['[hawc2_exe]', '[sim_id]', '[res_dir]', '[eigenfreq_dir]',
+                '[log_dir]', '[animation_dir]', '[turb_dir]', '[micro_dir]',
+                '[meander_dir]', '[model_zip]', '[htc_dir]', '[hydro_dir]',
+                '[mooring_dir]', '[turb_base_name]', '[case_id]'
+                '[micro_base_name]', '[meander_base_name]', '[pbs_out_dir]',
+                '[pbs_in_dir]']
+
+        # are allowed to have up to 2 ..
+        dirs_db = ['[turb_db_dir]', '[micro_db_dir]', '[meander_db_dir]']
+
+        dirlists = ['[copyto_generic]', '[copyto_files]', '[copyback_f1]',
+                    '[copyback_f1_rename]', '[copyto_generic_f1]',
+                    '[copyto_f1]', '[copyback_files]', '[copyback_frename]']
+
+        # is the only allowed absolute path
+        misc.path_sanitize(tag_dict['[run_dir]'], allowabs=True)
+
+        for pathtag in dirs:
+            if pathtag not in tag_dict:
+                continue
+            path = tag_dict[pathtag]
+            # Booleans are allowed
+            if isinstance(path, bool) or isinstance(path, type(None)):
+                path = str(path)
+            misc.path_sanitize(path)
+
+        for pathtag in dirs_db:
+            if pathtag not in tag_dict:
+                continue
+            path = tag_dict[pathtag]
+            # Booleans are allowed
+            if isinstance(path, bool) or isinstance(path, type(None)):
+                path = str(path)
+            misc.path_sanitize(path, allowdd=True)
+
+        for pathtag in dirlists:
+            if pathtag not in tag_dict:
+                continue
+            for path in tag_dict[pathtag]:
+                # Booleans are allowed
+                if isinstance(path, bool) or isinstance(path, type(None)):
+                    path = str(path)
+                misc.path_sanitize(path)
+
 
 # TODO: rewrite the error log analysis to something better. Take different
 # approach: start from the case and see if the results are present. Than we
@@ -3864,7 +3595,8 @@ class Cases(object):
         """Given the log file analysis and the Cases tag list, generate a list
         of failed cases. This is usefull when either some cases have been
         re-run or when the post-processing is done at the same time as the
-        simulations (e.g. zipchunks approach).
+        simulations (e.g. zipchunks approach). Cases for which the elapsted_time
+        column of the error logs is 0 or smaller are also considered as failed
 
         Parameters
         ----------
@@ -3900,6 +3632,9 @@ class Cases(object):
         # convert case_id to log file names
         # logids = pd.DataFrame(columns=[''])
         df_cases['logid'] = df_cases['[log_dir]'] + df_cases['[case_id]'] + '.log'
+        # remove those cases for which the logfile has not ended with the last
+        # statement "Elapsed time", results in value 0
+        df_err = df_err[df_err['elapsted_time'] > 0]
         # we only need to merge with errorlogs using a portion of the data
         # join error logs and df_cases on the logid
         df = pd.merge(df_cases[['logid', '[case_id]']], df_err[['file_name']],
@@ -3919,23 +3654,28 @@ class Cases(object):
             save_pickle(os.path.join(self.post_dir, self.sim_id + '_fail.pkl'),
                         self.cases_fail)
 
-    def remove_failed(self):
+    def remove_failed(self, verbose=False):
 
         # don't do anything if there is nothing defined
         if self.cases_fail == None:
             print('no failed cases to remove')
             return
 
+        nr_cases = len(self.cases)
         # ditch all the failed cases out of the htc_dict
         # otherwise we will have fails when reading the results data files
         for k in self.cases_fail:
             try:
                 self.cases_fail[k] = copy.copy(self.cases[k])
                 del self.cases[k]
-                print('removed from htc_dict due to error: ' + k)
+                if verbose:
+                    print('removed from htc_dict due to error: ' + k)
             except KeyError:
-                print('WARNING: failed case does not occur in cases')
-                print('   ', k)
+                if verbose:
+                    print('WARNING: failed case does not occur in cases')
+                    print('   ', k)
+        rpl = (len(self.cases_fail), nr_cases)
+        print('removed %i failed cases (out of %i)' % rpl)
 
     def load_failed(self, sim_id):
 
@@ -4007,7 +3747,8 @@ class Cases(object):
                    ch_fatigue={}, update=False, add_sensor=None,
                    chs_resultant=[], i0=0, i1=None, saveinterval=1000,
                    csv=True, suffix=None, A=None, add_sigs={},
-                   ch_wind=None, save_new_sigs=False, xlsx=False):
+                   ch_wind=None, save_new_sigs=False, xlsx=False,
+                   bearing_damage_lst=()):
         """
         Calculate statistics and save them in a pandas dataframe. Save also
         every 500 cases the statistics file.
@@ -4031,7 +3772,11 @@ class Cases(object):
 
         add_sigs : dict, default={}
             channel name, expression key/value paires. For example,
-            '[p1-p1-node-002-forcevec-z]*3 + [p1-p1-node-002-forcevec-y]'
+            'p1-p1-node-002-forcevec-z*3 + p1-p1-node-002-forcevec-y'
+
+        bearing_damage_lst : iterable, default=()
+            Input for wetb.fatigue_tools.bearing_damage: angle and moment
+            channels of the bearing of interest.
 
         chs_resultant
 
@@ -4344,6 +4089,15 @@ class Cases(object):
             # calculate the statistics values
             stats = self.res.calc_stats(self.sig, i0=i0, i1=i1)
 
+            # calculate any bearing damage
+            for name, angle_moment_lst in bearing_damage_lst:
+                angle_moment_timeseries_lst = []
+                for aa, mm in angle_moment_lst:
+                    angle = self.sig[:,self.res.ch_dict[aa]['chi']]
+                    moment = self.sig[:,self.res.ch_dict[mm]['chi']]
+                    angle_moment_timeseries_lst.append((angle, moment))
+                stats[name] = bearing_damage(angle_moment_timeseries_lst)
+
             # Because each channel is a new row, it doesn't matter how many
             # data channels each case has, and this approach does not brake
             # when different cases have a different number of output channels
@@ -4582,7 +4336,7 @@ class Cases(object):
             # if we have a list, convert to string
             if type(col[0]).__name__ == 'list':
                 for ii, item in enumerate(col):
-                    col[ii] = '**'.join(item)
+                    col[ii] = '*;*'.join(item)
             # if we already have an array (statistics) or a list of numbers
             # do not try to cast into another data type, because downcasting
             # in that case will not raise any exception
@@ -4602,11 +4356,11 @@ class Cases(object):
                 try:
                     df_dict2[str(colkey)] = np.array(col, dtype=np.float64)
                 except ValueError:
-                    df_dict2[str(colkey)] = np.array(col, dtype=np.str)
+                    df_dict2[str(colkey)] = np.array(col, dtype=str)
             except TypeError:
                 # in all other cases, make sure we have converted them to
                 # strings and NOT unicode
-                df_dict2[str(colkey)] = np.array(col, dtype=np.str)
+                df_dict2[str(colkey)] = np.array(col, dtype=str)
             except Exception as e:
                 print('failed to convert column %s to single data type' % colkey)
                 raise(e)
@@ -4823,7 +4577,7 @@ class Cases(object):
 
     def AEP(self, dfs, fh_lst=None, ch_powe='DLL-2-inpvec-2', extra_cols=[],
             res_dir='res/', dlc_folder="dlc%s_iec61400-1ed3/", csv=False,
-            new_sim_id=False, save=False, years=20.0, update=False, xlsx=False):
+            new_sim_id=False, save=False, update=False, xlsx=False):
 
         """
         Calculate the Annual Energy Production (AEP) for DLC1.2 cases.
@@ -4923,6 +4677,12 @@ class Cases(object):
                                 save=save, check_datatypes=True, xlsx=xlsx,
                                 complib=self.complib)
 
+        # check if the power channel actually exists first!
+        if ch_powe not in dfs[chan_col_name].unique():
+            msg = 'The defined channel for the electrical power does not '
+            msg =+ 'exist: %s' % ch_powe
+            raise UserWarning(msg)
+
         # and select only the power channels
         dfs_powe = dfs[dfs[chan_col_name]==ch_powe]
 
@@ -5208,13 +4968,26 @@ class Cases(object):
 
         return envelope
 
-    def envelopes(self, silent=False, ch_list=[], append=''):
+    def envelopes(self, silent=False, ch_list=[], append='', int_env=False,
+                  Nx=300):
         """
         Calculate envelopes and save them in a table.
 
         Parameters
         ----------
 
+        silent
+
+        ch_list
+
+        append
+
+        int_env : boolean, default=False
+            If the logic parameter is True, the function will interpolate the
+            envelope on a given number of points
+
+        Nx : int, default=300
+            Number of points for the envelope interpolation
 
         Returns
         -------
@@ -5254,7 +5027,7 @@ class Cases(object):
 
             self.load_result_file(case)
 
-            envelope = self.compute_envelopes(ch_list, int_env=False, Nx=300)
+            envelope = self.compute_envelopes(ch_list, int_env=int_env, Nx=Nx)
 
             for ch_id in ch_list:
                 title = str(ch_id[0].replace('-', '_'))
@@ -5311,83 +5084,6 @@ class EnvelopeClass(object):
         Fz = tbl.Float32Col()
 
 
-# TODO: implement this
-class Results(object):
-    """
-    Move all Hawc2io to here? NO: this should be the wrapper, to interface
-    the htc_dict with the io functions
-
-    There should be a bare metal module/class for those who only want basic
-    python support for HAWC2 result files and/or launching simulations.
-
-    How to properly design this module? Change each class into a module? Or
-    leave like this?
-    """
-
-    # OK, for now use this to do operations on HAWC2 results files
-
-    def __init___(self):
-        """
-        """
-        pass
-
-    def m_equiv(self, st_arr, load, pos):
-        r"""Centrifugal corrected equivalent moment
-
-        Convert beam loading into a single equivalent bending moment. Note that
-        this is dependent on the location in the cross section. Due to the
-        way we measure the strain on the blade and how we did the calibration
-        of those sensors.
-
-        .. math::
-
-            \epsilon = \frac{M_{x_{equiv}}y}{EI_{xx}} = \frac{M_x y}{EI_{xx}}
-            + \frac{M_y x}{EI_{yy}} + \frac{F_z}{EA}
-
-            M_{x_{equiv}} = M_x + \frac{I_{xx}}{I_{yy}} M_y \frac{x}{y}
-            + \frac{I_{xx}}{Ay} F_z
-
-        Parameters
-        ----------
-
-        st_arr : np.ndarray(19)
-            Only one line of the st_arr is allowed and it should correspond
-            to the correct radial position of the strain gauge.
-
-        load : list(6)
-            list containing the load time series of following components
-            .. math:: load = F_x, F_y, F_z, M_x, M_y, M_z
-            and where each component is an ndarray(m)
-
-        pos : np.ndarray(2)
-            x,y position wrt neutral axis in the cross section for which the
-            equivalent load should be calculated
-
-        Returns
-        -------
-
-        m_eq : ndarray(m)
-            Equivalent load, see main title
-
-        """
-
-        F_z = load[2]
-        M_x = load[3]
-        M_y = load[4]
-
-        x, y = pos[0], pos[1]
-
-        A = st_arr[ModelData.st_headers.A]
-        I_xx = st_arr[ModelData.st_headers.Ixx]
-        I_yy = st_arr[ModelData.st_headers.Iyy]
-
-        M_x_equiv = M_x + ( (I_xx/I_yy)*M_y*(x/y) ) + ( F_z*I_xx/(A*y) )
-        # or ignore edgewise moment
-        #M_x_equiv = M_x + ( F_z*I_xx/(A*y) )
-
-        return M_x_equiv
-
-
 class MannTurb64(prepost.PBSScript):
     """
     alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
@@ -5422,6 +5118,64 @@ class MannTurb64(prepost.PBSScript):
         self.silent = silent
         self.pbs_in_dir = 'pbs_in_turb/'
 
+    def create_turb(self, base_name, out_base, turb_dir, turb_db_dir, param):
+        """
+
+        Parameters
+        ----------
+
+        base_name
+
+        out_base
+
+        turb_dir
+
+        turb_db_dir
+
+        param : dictionary
+            Should contain the following keys: [MannAlfaEpsilon], [MannL]
+            [MannGamma], [seed], [turb_nr_u], [turb_nr_u], [turb_nr_w],
+            [turb_dx], [turb_dy], [turb_dz], [high_freq_comp]
+        """
+
+        self.path_pbs_e = os.path.join(out_base, turb_dir, base_name + '.err')
+        self.path_pbs_o = os.path.join(out_base, turb_dir, base_name + '.out')
+        self.path_pbs_i = os.path.join(self.pbs_in_dir, turb_dir, base_name + '.p')
+
+        # apply winefix
+        self.prelude = self.winefix
+        # browse to scratch dir
+        self.prelude += 'cd {}\n'.format(self.scratchdir)
+
+        self.coda = '# COPY BACK FROM SCRATCH AND RENAME, remove _ at end\n'
+        # copy back to turb dir at the end
+        if turb_db_dir is not None:
+            dst = os.path.join('$PBS_O_WORKDIR', turb_db_dir, base_name)
+        else:
+            dst = os.path.join('$PBS_O_WORKDIR', turb_dir, base_name)
+        # FIXME: Mann64 turb exe creator adds an underscore to output
+        for comp in list('uvw'):
+            src = '{}_{}.bin'.format(base_name, comp)
+            dst2 = '{}{}.bin'.format(dst, comp)
+            # dst contains $PBS_O_WORKDIR, hence double quote
+            self.coda += "cp '{}' \"{}\"\n".format(src, dst2)
+
+        # alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
+        rpl = (float(param['[MannAlfaEpsilon]']),
+               float(param['[MannL]']),
+               float(param['[MannGamma]']),
+               int(float(param['[seed]'])),
+               int(float(param['[turb_nr_u]'])),
+               int(float(param['[turb_nr_v]'])),
+               int(float(param['[turb_nr_w]'])),
+               float(param['[turb_dx]']),
+               float(param['[turb_dy]']),
+               float(param['[turb_dz]']),
+               int(float(param['[high_freq_comp]'])))
+        params = '%1.6f %1.6f %1.6f %i %i %i %i %1.4f %1.4f %1.4f %i' % rpl
+        self.execution = "%s '%s' %s" % (self.exe, base_name, params)
+        self.create(check_dirs=True)
+
     def gen_pbs(self, cases):
         """
         Parameters
@@ -5436,57 +5190,51 @@ class MannTurb64(prepost.PBSScript):
         self.pbsworkdir = os.path.join(case0['[run_dir]'], '')
         if not self.silent:
             print('\nStart creating PBS files for turbulence with Mann64...')
-        for cname, case in cases.items():
 
-            # only relevant for cases with turbulence
-            if '[tu_model]' in case and int(case['[tu_model]']) == 0:
-                continue
-            if '[turb_base_name]' not in case:
-                continue
+        mann_turb = ['[MannAlfaEpsilon]', '[MannL]', '[MannGamma]',
+                     '[turb_nr_u]', '[turb_nr_v]', '[turb_nr_w]', '[seed]',
+                     '[turb_dx]', '[turb_dy]', '[turb_dz]',
+                     '[high_freq_comp]']
+        mann_micro = {k:k.replace(']', '_micro]') for k in mann_turb}
+        mann_meander = {k:k.replace(']', '_meander]') for k in mann_turb}
 
-            base_name = case['[turb_base_name]']
-            # pbs_in/out dir can contain subdirs, only take the inner directory
+        for cname, case in cases.items():
+
+            # pbs_in/out dir can contain subdirs, only take the base directory
             out_base = misc.path_split_dirs(case['[pbs_out_dir]'])[0]
-            turb = case['[turb_dir]']
-
-            self.path_pbs_e = os.path.join(out_base, turb, base_name + '.err')
-            self.path_pbs_o = os.path.join(out_base, turb, base_name + '.out')
-            self.path_pbs_i = os.path.join(self.pbs_in_dir, base_name + '.p')
-
-            # apply winefix
-            self.prelude = self.winefix
-            # browse to scratch dir
-            self.prelude += 'cd {}\n'.format(self.scratchdir)
-
-            self.coda = '# COPY BACK FROM SCRATCH AND RENAME, remove _ at end\n'
-            # copy back to turb dir at the end
-            if case['[turb_db_dir]'] is not None:
-                dst = os.path.join('$PBS_O_WORKDIR', case['[turb_db_dir]'],
-                                   base_name)
-            else:
-                dst = os.path.join('$PBS_O_WORKDIR', case['[turb_dir]'],
-                                   base_name)
-            # FIXME: Mann64 turb exe creator adds an underscore to output
-            for comp in list('uvw'):
-                src = '{}_{}.bin'.format(base_name, comp)
-                dst2 = '{}{}.bin'.format(dst, comp)
-                self.coda += 'cp {} {}\n'.format(src, dst2)
-
-            # alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
-            rpl = (float(case['[MannAlfaEpsilon]']),
-                   float(case['[MannL]']),
-                   float(case['[MannGamma]']),
-                   int(case['[seed]']),
-                   int(case['[turb_nr_u]']),
-                   int(case['[turb_nr_v]']),
-                   int(case['[turb_nr_w]']),
-                   float(case['[turb_dx]']),
-                   float(case['[turb_dy]']),
-                   float(case['[turb_dz]']),
-                   int(case['[high_freq_comp]']))
-            params = '%1.6f %1.6f %1.6f %i %i %i %i %1.4f %1.4f %1.4f %i' % rpl
-            self.execution = '%s %s %s' % (self.exe, base_name, params)
-            self.create(check_dirs=True)
+
+            # NORMAL ATMOSPHERIC TURBULENCE
+            # only relevant for cases with turbulence
+            req = '[tu_model]' in case and '[turb_base_name]' in case
+            if req and int(case['[tu_model]'])==1:
+                base = case['[turb_base_name]']
+                # pbs_in/out dir can contain subdirs, only take the base directory
+                turb_dir = case['[turb_dir]']
+                turb_db_dir = case['[turb_db_dir]']
+                # more fail safe switches in case user did not change defaults
+                if turb_dir and base and base.lower()!='none' and base!='':
+                    self.create_turb(base, out_base, turb_dir, turb_db_dir,
+                                     {key:case[key] for key in mann_turb})
+
+            # MICRO TURBULENCE
+            if ('[micro_dir]' in case) and ('[micro_base_name]' in case):
+                base = case['[micro_base_name]']
+                turb_dir = case['[micro_dir]']
+                turb_db_dir = case['[micro_db_dir]']
+                # more fail safe switches in case user did not change defaults
+                if turb_dir and base and base.lower()!='none' and base!='':
+                    p = {key:case[mann_micro[key]] for key in mann_turb}
+                    self.create_turb(base, out_base, turb_dir, turb_db_dir, p)
+
+            # MEANDER TURBULENCE
+            if ('[meander_dir]' in case) and ('[meander_base_name]' in case):
+                base = case['[meander_base_name]']
+                turb_dir = case['[meander_dir]']
+                turb_db_dir = case['[meander_db_dir]']
+                # more fail safe switches in case user did not change defaults
+                if turb_dir and base and base.lower()!='none' and base!='':
+                    p = {key:case[mann_meander[key]] for key in mann_turb}
+                    self.create_turb(base, out_base, turb_dir, turb_db_dir, p)
 
 
 def eigenbody(cases, debug=False):
@@ -5535,6 +5283,7 @@ def eigenbody(cases, debug=False):
 
     return cases
 
+
 def eigenstructure(cases, debug=False):
     """
     Read HAWC2 structure eigenalysis result file

wetb/prepost/aep_del.py

+# -*- coding: utf-8 -*-
+"""
+Created on Fri Oct 30 13:40:44 2020
+
+@author: dave
+"""
+
+import os
+from argparse import ArgumentParser
+import math
+
+import numpy as np
+import pandas as pd
+
+from wetb.dlc.high_level import Weibull_IEC
+
+
+def aep_del(df_stats, resdir):
+    """Works for regular wetb statistics df.
+    """
+
+    # note about the unique channel names: they are saved under
+    # prepost-data/{sim_id}_unique-channel-names.csv
+    # after running parpost.py
+
+    # TODO: make option to have user define list of channels
+    # channels = df_stats['channel'].unique()
+
+    # DLC12: assumes that case_id (the file name) has the following format:
+    # dlc10_wsp04_wdir000_s1001.htc
+
+    # -----------------------------------------------------
+    # USER INPUT
+    neq_life = 1e7
+    vref = 50 # used for the IEC Weibull distribution
+    fact_1year = 365 * 24 * 0.975
+    years_del = 20
+    # specify which channel to use for the wind speed
+    chwind = 'windspeed-global-Vy-0.00-0.00--119.00'
+    # specify which channel to use for the electrical power output
+    chpowe = 'DLL-generator_servo-inpvec-2'
+    yawdist = {'wdir000':0.5, 'wdir010':0.25, 'wdir350':0.25}
+    dlc10_dir = 'dlc10_powercurve'
+    dlc12_dir = 'dlc12_iec61400-1ed3'
+    # -----------------------------------------------------
+
+    df_stats['hour_weight'] = 0
+
+    # -----------------
+    # DLC10
+    # -----------------
+    df_pc_dlc10 = pd.DataFrame()
+    resdir_dlc10 = os.path.join(resdir, dlc10_dir)
+    df_dlc10 = df_stats[df_stats['res_dir']==resdir_dlc10]
+    if len(df_dlc10) > 0:
+        df_pivot = df_dlc10.pivot(index='case_id', columns='channel', values='mean')
+        df_pivot.sort_values(chwind, inplace=True)
+        # wind speed and power
+        wsp = df_pivot[chwind].values
+        powe = df_pivot[chpowe].values
+        # wind speed probability distribution for 1 year
+        wsp_dist = Weibull_IEC(vref, wsp) * fact_1year
+        # save data into the DataFrame table
+        df_pc_dlc10['wsp'] = wsp
+        df_pc_dlc10['powe'] = powe
+        df_pc_dlc10['hours'] = wsp_dist
+        df_pc_dlc10['aep'] = (powe * wsp_dist).sum()
+
+    # -----------------
+    # DLC12
+    # -----------------
+    resdir_dlc12 = os.path.join(resdir, dlc12_dir)
+    df_dlc12 = df_stats[df_stats['res_dir']==resdir_dlc12].copy()
+    df_pivot = df_dlc12.pivot(index='case_id', columns='channel', values='mean')
+    # now add the wsp, yaw, seed columns again
+    tmp = pd.DataFrame(df_pivot.index.values)
+    df_pivot['dlc'] = ''
+    df_pivot['wsp'] = ''
+    df_pivot['yaw'] = ''
+    df_pivot['seed'] = ''
+    # assumes that case_id (the file name) has the following format:
+    # dlc12_wsp04_wdir000_s1001.htc
+    df_pivot[['dlc', 'wsp', 'yaw', 'seed']] = tmp[0].str.split('_', expand=True).values
+    df_pivot['wsp'] = (df_pivot['wsp'].str[3:]).astype(int)
+
+    # sorted on wsp, yaw, seed
+    df_pivot.sort_index(inplace=True)
+
+    # figure out how many seeds there are per yaw and wsp
+    seeds = {'wsp_yaw':[], 'nr_seeds':[]}
+    for ii, grii in df_pivot.groupby('wsp'):
+        for jj, grjj in grii.groupby('yaw'):
+            seeds['wsp_yaw'].append('%02i_%s' % (ii, jj))
+            seeds['nr_seeds'].append(len(grjj))
+
+    nr_seeds = np.array(seeds['nr_seeds']).max()
+    seeds_min = np.array(seeds['nr_seeds']).min()
+    assert nr_seeds == seeds_min
+
+    # yaw error probabilities
+    for yaw, dist in yawdist.items():
+        df_pivot.loc[df_pivot['yaw']==yaw, 'hour_weight'] = dist/nr_seeds
+    df_pivot['annual_hours_wsp'] = 0
+
+    # Weibull hour distribution as function of wind speeds
+    wsps = df_pivot['wsp'].unique()
+    wsps.sort()
+    wsp_dist = Weibull_IEC(vref, wsps) * fact_1year
+    for i, wsp in enumerate(wsps):
+        sel = df_pivot['wsp']==wsp
+        df_pivot.loc[sel, 'annual_hours_wsp'] = wsp_dist[i]
+
+    df_pivot['annual_hour_dist'] = df_pivot['annual_hours_wsp'] * df_pivot['hour_weight']
+
+    # check we still have the same amount of hours in a year
+    hours1 = df_pivot['annual_hour_dist'].sum()
+    hours2 = wsp_dist.sum()
+    assert np.allclose(hours1, hours2)
+
+    aep = (df_pivot['annual_hour_dist'] * df_pivot[chpowe]).sum()
+
+    df_pc_dlc12 = df_pivot[[chwind, chpowe, 'dlc', 'wsp', 'yaw', 'seed',
+                            'hour_weight', 'annual_hours_wsp', 'annual_hour_dist']].copy()
+    df_pc_dlc12['aep'] = aep
+
+    # -----------------
+    # DLC12 DEL
+    ms = [k for k in df_dlc12.columns if k.find('m=') > -1]
+    dict_Leq = {m:[] for m in ms}
+    dict_Leq['channel'] = []
+
+    # TODO: make option to have user define list of channels
+    # for channel in channels:
+        #statsel = df_dlc12[df_dlc12['channel']==channel].copy()
+
+    for channel, grch in df_dlc12.groupby('channel'):
+        # sort the case_id values in the same way as the pivot table is to
+        # align the hours with the corresponding case
+        statsel = grch.copy()
+        statsel.set_index('case_id', inplace=True)
+        statsel = statsel.loc[df_pivot.index,:]
+
+        dict_Leq['channel'].append(channel)
+
+        for m in ms:
+            m_ = float(m.split('=')[1])
+            eq1hz_mod = np.power(statsel[m].values, m_)
+            # R_eq_mod will have to be scaled from its simulation length
+            # to 1 hour (hour distribution is in hours...). Since the
+            # simulation time has not been multiplied out of R_eq_mod yet,
+            # we can just multiply with 3600 (instead of doing 3600/neq)
+            tmp = (eq1hz_mod * df_pivot['annual_hour_dist'] * years_del * 3600).sum()
+            # the effective Leq for each of the material constants
+            leq = math.pow(tmp/neq_life, 1.0/m_)
+            dict_Leq[m].append(leq)
+
+    df_leq = pd.DataFrame(dict_Leq)
+    df_leq.set_index('channel', inplace=True)
+
+    return df_pc_dlc10, df_pc_dlc12, df_leq
+
+
+if __name__ == '__main__':
+
+    parser = ArgumentParser(description = "Calculate AEP and LEQ")
+    parser.add_argument('--res', type=str, default='res', action='store',
+                        dest='resdir', help='Directory containing result files')
+    opt = parser.parse_args()
+
+    sim_id = os.path.basename(os.getcwd())
+    fname = f'prepost-data/{sim_id}_statistics.h5'
+    df_stats = pd.read_hdf(fname, key='table')
+
+    df_pc_dlc10, df_pc_dlc12, df_leq= aep_del(df_stats, opt.resdir)
+
+    # save fot Excel files
+    df_leq.to_excel(f'prepost-data/{sim_id}_del.xlsx')
+    df_pc_dlc10.to_excel(f'prepost-data/{sim_id}_aep_10.xlsx')
+    df_pc_dlc12.to_excel(f'prepost-data/{sim_id}_aep_12.xlsx')

wetb/prepost/dlcdefs.py

@@ -4,14 +4,6 @@ Created on Wed Nov  5 14:01:25 2014
 
 @author: dave
 """
-from __future__ import print_function
-from __future__ import unicode_literals
-from __future__ import division
-from __future__ import absolute_import
-from builtins import str
-from future.utils import viewitems
-from future import standard_library
-standard_library.install_aliases()
 
 import os
 import unittest
@@ -132,34 +124,6 @@ def variable_tag_func(master, case_id_short=False):
     return master
 
 
-def vartag_dlcs(master):
-
-    mt = master.tags
-
-    dlc_case = mt['[Case folder]']
-    mt['[data_dir]'] = 'data/'
-    mt['[res_dir]'] = 'res/%s/' % dlc_case
-    mt['[log_dir]'] = 'logfiles/%s/' % dlc_case
-    mt['[htc_dir]'] = 'htc/%s/' % dlc_case
-    mt['[case_id]'] = mt['[Case id.]']
-    mt['[time_stop]'] = mt['[time stop]']
-    mt['[turb_base_name]'] = mt['[Turb base name]']
-    mt['[DLC]'] = mt['[Case id.]'].split('_')[0][3:]
-    mt['[pbs_out_dir]'] = 'pbs_out/%s/' % dlc_case
-    mt['[pbs_in_dir]'] = 'pbs_in/%s/' % dlc_case
-    mt['[iter_dir]'] = 'iter/%s/' % dlc_case
-    if mt['[eigen_analysis]']:
-        rpl = (dlc_case, mt['[Case id.]'])
-        mt['[eigenfreq_dir]'] = 'res_eigen/%s/%s/' % rpl
-    mt['[duration]'] = str(float(mt['[time_stop]']) - float(mt['[t0]']))
-    # replace nan with empty
-    for ii, jj in mt.items():
-        if jj == 'nan':
-            mt[ii] = ''
-
-    return master
-
-
 def vartag_excel_stabcon(master):
     """Variable tag function type that generates a hydro input file for the
     wave kinematics dll if [hydro input name] is defined properly.
@@ -172,9 +136,25 @@ def vartag_excel_stabcon(master):
     print('creating hydro input file for: %s.inp\n' % mt['[hydro input name]'])
 
     mt['[wdepth]'] = float(mt['[wdepth]'])
-    mt['[Hs]'] = float(mt['[Hs]'])
-    mt['[Tp]'] = float(mt['[Tp]'])
 
+    if '[Hs]' in mt:
+        hs_key = '[Hs]'
+    elif '[hs]' in mt:
+        hs_key = '[hs]'
+    else:
+        msg = 'Significant wave tag is expected to be called: [Hs]'
+        raise KeyError(msg)
+
+    if '[Tp]' in mt:
+        tp_key = '[Tp]'
+    elif '[tp]' in mt:
+        tp_key = '[tp]'
+    else:
+        msg = 'Wave period tag is expected to be called: [Tp]'
+        raise KeyError(msg)
+
+    mt[hs_key] = float(mt[hs_key])
+    mt[tp_key] = float(mt[tp_key])
 
     if '[wave_gamma]' not in mt or not mt['[wave_gamma]']:
         mt['[wave_gamma]'] = 3.3
@@ -203,7 +183,7 @@ def vartag_excel_stabcon(master):
         embed_sf = None
         embed_sf_t0 = None
 
-    hio = hydro_input(wavetype=mt['[wave_type]'], Hs=mt['[Hs]'], Tp=mt['[Tp]'],
+    hio = hydro_input(wavetype=mt['[wave_type]'], Hs=mt[hs_key], Tp=mt[tp_key],
                       gamma=mt['[wave_gamma]'], wdepth=mt['[wdepth]'],
                       spectrum=mt['[wave_spectrum]'], seed=mt['[wave_seed]'],
                       stretching=mt['[stretching]'], coef=mt['[wave_coef]'],
@@ -296,7 +276,7 @@ def tags_defaults(master):
     master.tags['[iter_dir]']      = 'iter/'
     master.tags['[turb_dir]']      = 'turb/'
     master.tags['[turb_db_dir]']   = '../turb/'
-    master.tags['[wake_dir]']      = False
+    master.tags['[micro_dir]']      = False
     master.tags['[hydro_dir]']     = False
     master.tags['[mooring_dir]']   = False
     master.tags['[externalforce]'] = False
@@ -391,13 +371,13 @@ def excel_stabcon(proot, fext='xlsx', pignore=None, pinclude=None, sheet=0,
 
     if not silent:
         k = 0
-        for df in dict_dfs:
+        for dlc, df in dict_dfs.items():
             k += len(df)
         print('in which a total of %i cases are defined.' % k)
 
     opt_tags = []
 
-    for (dlc, df) in sorted(viewitems(dict_dfs)):
+    for (dlc, df) in sorted(dict_dfs.items()):
         # replace ';' with False, and Nan(='') with True
         # this is more easy when testing for the presence of stuff compared
         # to checking if a value is either True/False or ''/';'
@@ -411,12 +391,12 @@ def excel_stabcon(proot, fext='xlsx', pignore=None, pinclude=None, sheet=0,
         for count, row in df2.iterrows():
             tags_dict = {}
             # construct to dict, convert unicode keys/values to strings
-            for key, value in row.iteritems():
+            for key, value in row.items():
                 if isinstance(value, str):
                     tags_dict[str(key)] = str(value)
                 else:
                     tags_dict[str(key)] = value
-                # convert ; and empty to False/True
+                # convert ; and empty to False/True, "none" to None
                 if isinstance(tags_dict[str(key)], str):
                     if tags_dict[str(key)] == ';':
                         tags_dict[str(key)] = False
@@ -424,6 +404,8 @@ def excel_stabcon(proot, fext='xlsx', pignore=None, pinclude=None, sheet=0,
                         tags_dict[str(key)] = True
                     elif tags_dict[str(key)].lower() == 'nan':
                         tags_dict[str(key)] = True
+                    elif tags_dict[str(key)].lower() == 'none':
+                        tags_dict[str(key)] = None
 
             # FIXME: this horrible mess requires a nice and clearly defined
             # tag spec/naming convention, and with special tag prefix
@@ -443,6 +425,14 @@ def excel_stabcon(proot, fext='xlsx', pignore=None, pinclude=None, sheet=0,
             else:
                 raise KeyError('[seed] should be used as tag for turb. seed')
 
+            # for backwards compatibility
+            if '[wake_dir]' in tags_dict and not '[micro_dir]':
+                tags_dict['[micro_dir]'] = tags_dict['[wake_dir]']
+            if '[wake_db_dir]' in tags_dict and not '[micro_db_dir]':
+                tags_dict['[micro_db_dir]'] = tags_dict['[wake_db_dir]']
+            if '[wake_base_name]' in tags_dict and not '[micro_base_name]':
+                tags_dict['[micro_base_name]'] = tags_dict['[wake_base_name]']
+
             tags_dict['[Case folder]'] = tags_dict['[Case folder]'].lower()
             tags_dict['[Case id.]'] = tags_dict['[Case id.]'].lower()
             dlc_case = tags_dict['[Case folder]']
@@ -460,11 +450,22 @@ def excel_stabcon(proot, fext='xlsx', pignore=None, pinclude=None, sheet=0,
                 tags_dict['[time_stop]'] = tags_dict['[time stop]']
             else:
                 tags_dict['[time stop]'] = tags_dict['[time_stop]']
-            try:
-                tags_dict['[turb_base_name]'] = tags_dict['[Turb base name]']
-            except KeyError:
+
+            # [tu_model] is in the docs, but some of dave's examples have
+            # [turb_format]...so that is great...
+            if '[turb_format]' in tags_dict:
+                tags_dict['[tu_model]'] = tags_dict['[turb_format]']
+
+            if '[Turb base name]' in tags_dict:
+                if not '[turb_base_name]' in tags_dict:
+                    tags_dict['[turb_base_name]'] = tags_dict['[Turb base name]']
+            elif not '[turb_base_name]' in tags_dict:
                 tags_dict['[turb_base_name]'] = None
                 tags_dict['[Turb base name]'] = None
+            # for backwards compatibility: user has only defined [Turb base name]
+            elif not '[Turb base name]' in tags_dict:
+                tags_dict['[Turb base name]'] = tags_dict['[turb_base_name]']
+
             tags_dict['[DLC]'] = tags_dict['[Case id.]'].split('_')[0][3:]
             if '[pbs_out_dir]' not in tags_dict:
                 tags_dict['[pbs_out_dir]'] = 'pbs_out/%s/' % dlc_case

wetb/prepost/dlcplots.py

@@ -4,16 +4,6 @@ Created on Tue Sep 16 10:21:11 2014
 
 @author: dave
 """
-from __future__ import print_function
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import str
-from future import standard_library
-standard_library.install_aliases()
-
-#print(*objects, sep=' ', end='\n', file=sys.stdout)
-
 import os
 # import socket
 import gc
@@ -52,14 +42,17 @@ plt.rc('legend', numpoints=1)
 plt.rc('legend', borderaxespad=0)
 
 
-def merge_sim_ids(sim_ids, post_dirs, post_dir_save=False):
+def merge_sim_ids(sim_ids, post_dirs, post_dir_save=False, columns=None):
     """
     """
 
     cols_extra = ['[run_dir]', '[res_dir]', '[wdir]', '[DLC]', '[Case folder]']
+    min_itemsize={'channel':100, '[run_dir]':100, '[res_dir]':100, '[DLC]':10,
+                  '[Case folder]':100}
 
     # map the run_dir to the same order as the post_dirs, labels
     run_dirs = []
+
     # avoid saving merged cases if there is only one!
     if type(sim_ids).__name__ == 'list' and len(sim_ids) == 1:
         sim_ids = sim_ids[0]
@@ -73,7 +66,9 @@ def merge_sim_ids(sim_ids, post_dirs, post_dir_save=False):
             else:
                 post_dir = post_dirs
             cc = sim.Cases(post_dir, sim_id, rem_failed=True)
-            df_stats, _, _ = cc.load_stats(columns=None, leq=False)
+            df_stats, _, _ = cc.load_stats(leq=False)
+            if columns is not None:
+                df_stats = df_stats[columns]
 
             # stats has only a few columns identifying the different cases
             # add some more for selecting them
@@ -114,19 +109,27 @@ def merge_sim_ids(sim_ids, post_dirs, post_dir_save=False):
             if ii == 0:
                 # and save somewhere so we can add the second data frame on
                 # disc
-                df_stats.to_hdf(fmerged, 'table', mode='w', format='table',
-                                complevel=9, complib='blosc')
+                store = pd.HDFStore(fmerged, mode='w', complevel=9, complib='zlib')
+                store.append('table', df_stats, min_itemsize=min_itemsize)
+                print(store.get_storer('table').table.description)
+                # df_stats.to_hdf(fmerged, '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(fmerged, 'table', mode='r+', format='table',
-                                complevel=9, complib='blosc', append=True)
+                store.append('table', df_stats)
+                # will fail if there are longer string columns compared to ii=0
+                # df_stats.to_hdf(fmerged, '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()
+
+        store.close()
+
         # and load the reduced combined set
         print('loading merged stats:            %s' % fmerged)
         df_stats = pd.read_hdf(fmerged, 'table')
@@ -137,8 +140,14 @@ def merge_sim_ids(sim_ids, post_dirs, post_dir_save=False):
         if isinstance(post_dirs, list):
             post_dir = post_dirs[0]
         cc = sim.Cases(post_dir, sim_id, rem_failed=True)
-        df_stats, _, _ = cc.load_stats(leq=False)
-        run_dirs = [df_stats['[run_dir]'].unique()[0]]
+        df_stats, _, _ = cc.load_stats(columns=columns, leq=False)
+        if columns is not None:
+            df_stats = df_stats[columns]
+
+        try:
+            run_dirs = [df_stats['[run_dir]'].unique()[0]]
+        except KeyError:
+            run_dirs = []
 
         # stats has only a few columns identifying the different cases
         # add some more for selecting them
@@ -153,7 +162,7 @@ def merge_sim_ids(sim_ids, post_dirs, post_dir_save=False):
         add_cols = list(cols_cc - set(df_stats.columns))
         add_cols.append('[case_id]')
         dfc = dfc[add_cols]
-        df_stats = pd.merge(df_stats, dfc, on='[case_id]')
+        df_stats = pd.merge(df_stats, dfc, right_on='[case_id]', left_on='case_id')
         if '[Windspeed]' in df_stats.columns and '[wsp]' in df_stats.columns:
             df_stats.drop('[wsp]', axis=1, inplace=True)
         if wsp != '[Windspeed]':
@@ -317,13 +326,13 @@ def plot_dlc_stats(df_stats, plot_chans, fig_dir_base, labels=None,
             # figure file name will be the first channel
             if isinstance(ch_names, list):
                 ch_name = ch_names[0]
-                fname_base = ch_names[0]#.replace(' ', '_')
+                fname_base = ch_names[0].replace('/', '_')
                 df_chan = gr_dlc[gr_dlc.channel.isin(ch_names)]
             else:
                 ch_name = ch_names
                 ch_names = [ch_names]
                 df_chan = gr_dlc[gr_dlc.channel == ch_names]
-                fname_base = ch_names#.replace(' ', '_')
+                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
@@ -351,11 +360,11 @@ def plot_dlc_stats(df_stats, plot_chans, fig_dir_base, labels=None,
 
             fig, axes = mplutils.make_fig(nrows=1, ncols=1,
                                           figsize=figsize, dpi=120)
-            ax = axes[0,0]
+            ax = axes
             # seperate figure for the mean of the 1Hz equivalent loads
             fig2, axes2 = mplutils.make_fig(nrows=1, ncols=1,
                                             figsize=figsize, dpi=120)
-            ax2 = axes2[0,0]
+            ax2 = axes2
 
             if fig_dir_base is None and len(sim_ids) < 2:
                 res_dir = df_chan['[res_dir]'][:1].values[0]

wetb/prepost/dlctemplate.py

@@ -4,20 +4,8 @@ Created on Thu Sep 18 13:00:25 2014
 
 @author: dave
 """
-from __future__ import print_function
-from __future__ import unicode_literals
-from __future__ import division
-from __future__ import absolute_import
-from builtins import dict
-from builtins import str
-from builtins import range
-from future import standard_library
-standard_library.install_aliases()
-
 import os
-import socket
 from argparse import ArgumentParser
-from sys import platform
 
 import numpy as np
 import pandas as pd
@@ -33,18 +21,18 @@ plt.rc('font', family='serif')
 plt.rc('xtick', labelsize=10)
 plt.rc('ytick', labelsize=10)
 plt.rc('axes', labelsize=12)
-# on Gorm tex printing doesn't work
-if socket.gethostname()[:2] in ['g-', 'je', 'j-']:
-    RUNMETHOD = 'pbs'
-else:
-    plt.rc('text', usetex=True)
-    # set runmethod based on the platform host
-    if platform in ["linux", "linux2", "darwin"]:
-        RUNMETHOD = 'linux-script'
-    elif platform == "win32":
-        RUNMETHOD = 'windows-script'
-    else:
-        RUNMETHOD = 'none'
+# plt.rc('text', usetex=True)
+
+# import socket
+# from sys import platform
+# set runmethod based on the platform host
+# if platform in ["linux", "linux2", "darwin"]:
+#     RUNMETHOD = 'linux-script'
+# elif platform == "win32":
+#     RUNMETHOD = 'windows-script'
+# else:
+#     RUNMETHOD = 'none'
+RUNMETHOD = 'none'
 plt.rc('legend', fontsize=11)
 plt.rc('legend', numpoints=1)
 plt.rc('legend', borderaxespad=0)
@@ -266,9 +254,10 @@ def variable_tag_func_mod1(master, case_id_short=False):
 def launch_dlcs_excel(sim_id, silent=False, verbose=False, pbs_turb=False,
                       runmethod=None, write_htc=True, zipchunks=False,
                       walltime='04:00:00', postpro_node=False, compress=False,
-                      dlcs_dir='htc/DLCs', wine_64bit=False,
-                      m=[3,4,6,8,9,10,12], postpro_node_zipchunks=True,
-                      wine_arch='win32', wine_prefix='~/.wine32'):
+                      dlcs_dir='htc/DLCs', postpro_node_zipchunks=True,
+                      wine_arch='win32', wine_prefix='.wine32', ppn=17,
+                      m=[3,4,6,8,9,10,12], prelude='', linux=False,
+                      update_model_data=False):
     """
     Launch load cases defined in Excel files
     """
@@ -277,10 +266,16 @@ def launch_dlcs_excel(sim_id, silent=False, verbose=False, pbs_turb=False,
     iter_dict['[empty]'] = [False]
 
     if postpro_node or postpro_node_zipchunks:
-        pyenv = 'wetb_py3'
+        pyenv = 'py36-wetb'
     else:
         pyenv = None
 
+    # FIXME: THIS IS VERY MESSY, we have wine_prefix/arch and exesingle/chunks
+    if linux:
+        wine_arch = None
+        wine_prefix = None
+        prelude = 'module load mpi/openmpi_1.6.5_intelv14.0.0\n'
+
     # see if a htc/DLCs dir exists
     # Load all DLC definitions and make some assumptions on tags that are not
     # defined
@@ -331,11 +326,11 @@ def launch_dlcs_excel(sim_id, silent=False, verbose=False, pbs_turb=False,
                                ignore_non_unique=False, run_only_new=False,
                                pbs_fname_appendix=False, short_job_names=False,
                                silent=silent, verbose=verbose, pyenv=pyenv,
-                               wine_64bit=wine_64bit, m=[3,4,6,8,9,10,12],
+                               m=[3,4,6,8,9,10,12], postpro_node=postpro_node,
+                               exechunks=None, exesingle=None, prelude=prelude,
                                postpro_node_zipchunks=postpro_node_zipchunks,
-                               postpro_node=postpro_node, exesingle=None,
-                               exechunks=None, wine_arch=wine_arch,
-                               wine_prefix=wine_prefix)
+                               wine_arch=wine_arch, wine_prefix=wine_prefix,
+                               update_model_data=update_model_data)
 
     if pbs_turb:
         # to avoid confusing HAWC2 simulations and Mann64 generator PBS files,
@@ -351,25 +346,20 @@ def launch_dlcs_excel(sim_id, silent=False, verbose=False, pbs_turb=False,
         # note that walltime here is for running all cases assigned to the
         # respective nodes. It is not walltime per case.
         sorts_on = ['[DLC]', '[Windspeed]']
-        create_chunks_htc_pbs(cases, sort_by_values=sorts_on, ppn=20,
-                              nr_procs_series=3, walltime='20:00:00',
+        create_chunks_htc_pbs(cases, sort_by_values=sorts_on, queue='workq',
+                              ppn=ppn, nr_procs_series=3, walltime='09:00:00',
                               chunks_dir='zip-chunks-jess', compress=compress,
-                              queue='workq', wine_64bit=wine_64bit,
-                              wine_arch=wine_arch, wine_prefix=wine_prefix)
-        create_chunks_htc_pbs(cases, sort_by_values=sorts_on, ppn=12,
-                              nr_procs_series=3, walltime='20:00:00',
-                              chunks_dir='zip-chunks-gorm', compress=compress,
-                              queue='workq', wine_64bit=wine_64bit,
-                              wine_arch=wine_arch, wine_prefix=wine_prefix)
+                              wine_arch=wine_arch, wine_prefix=wine_prefix,
+                              prelude=prelude, ppn_pbs=20)
 
     df = sim.Cases(cases).cases2df()
-    df.to_excel(os.path.join(POST_DIR, sim_id + '.xls'))
+    df.to_excel(os.path.join(POST_DIR, sim_id + '.xlsx'))
 
 
 def post_launch(sim_id, statistics=True, rem_failed=True, check_logs=True,
                 force_dir=False, update=False, saveinterval=2000, csv=False,
-                m=[3, 4, 6, 8, 9, 10, 12], neq=1e7, no_bins=46,
-                years=20.0, fatigue=True, A=None, AEP=False,
+                m=[3, 4, 6, 8, 9, 10, 12], neq=1e7, no_bins=46, int_env=False,
+                years=20.0, fatigue=True, A=None, AEP=False, nx=300,
                 save_new_sigs=False, envelopeturbine=False, envelopeblade=False,
                 save_iter=False, pbs_failed_path=False):
 
@@ -412,6 +402,7 @@ def post_launch(sim_id, statistics=True, rem_failed=True, check_logs=True,
 #        add_sigs = {name:expr}
 
         # in addition, sim_id and case_id are always added by default
+        # FIXME: HAS TO BE THE SAME AS required IN postpro_node_merge
         tags = ['[Case folder]', '[run_dir]', '[res_dir]', '[DLC]',
                 '[wsp]', '[Windspeed]', '[wdir]']
         add = None
@@ -434,8 +425,18 @@ def post_launch(sim_id, statistics=True, rem_failed=True, check_logs=True,
         # load the statistics in case they are missing
         if not statistics:
             df_stats, Leq_df, AEP_df = cc.load_stats()
+
+        # CAUTION: depending on the type output, electrical power can be two
+        # different things with the DTU Wind Energy Controller.
+        # Either manually set ch_powe to the correct value or use simple
+        # mechanism to figure out which one of two expected values it is.
+        if 'DLL-2-inpvec-2' in df_stats['channel'].unique():
+            ch_powe = 'DLL-2-inpvec-2'
+        elif 'DLL-generator_servo-inpvec-2' in df_stats['channel'].unique():
+            ch_powe = 'DLL-generator_servo-inpvec-2'
+
         df_AEP = cc.AEP(df_stats, csv=csv, update=update, save=True,
-                        ch_powe='DLL-2-inpvec-2')
+                        ch_powe=ch_powe)
 
     if envelopeblade:
         ch_list = []
@@ -448,7 +449,7 @@ def post_launch(sim_id, statistics=True, rem_failed=True, check_logs=True,
                                 'blade%i-blade%i-node-%3.3i-forcevec-x' % rpl,
                                 'blade%i-blade%i-node-%3.3i-forcevec-y' % rpl,
                                 'blade%i-blade%i-node-%3.3i-forcevec-z' % rpl])
-        cc.envelopes(ch_list=ch_list, append='_blade')
+        cc.envelopes(ch_list=ch_list, append='_blade', int_env=int_env, Nx=nx)
 
     if envelopeturbine:
         ch_list = [['tower-tower-node-001-momentvec-x',
@@ -463,7 +464,7 @@ def post_launch(sim_id, statistics=True, rem_failed=True, check_logs=True,
                    ['hub1-hub1-node-001-momentvec-x',
                    'hub1-hub1-node-001-momentvec-y',
                    'hub1-hub1-node-001-momentvec-z']]
-        cc.envelopes(ch_list=ch_list, append='_turbine')
+        cc.envelopes(ch_list=ch_list, append='_turbine', int_env=int_env, Nx=nx)
 
     if fatigue:
         # load the statistics in case they are missing
@@ -476,7 +477,7 @@ def post_launch(sim_id, statistics=True, rem_failed=True, check_logs=True,
     return df_stats, df_AEP, df_Leq
 
 
-def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
+def postpro_node_merge(tqdm=False, zipchunks=False):
     """With postpro_node each individual case has a .csv file for the log file
     analysis and a .csv file for the statistics tables. Merge all these single
     files into one table/DataFrame.
@@ -508,10 +509,22 @@ def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
     mdf = AppendDataFrames(tqdm=tqdm)
     # individual log file analysis does not have header, make sure to include
     # a line for the header
-    mdf.txt2txt(fcsv, path_pattern, tarmode='r:xz', header=None,
-                header_fjoined=lf._header(), recursive=True)
+    cols = mdf.txt2txt(fcsv, path_pattern, tarmode='r:xz', header=None,
+                       header_fjoined=lf._header(), recursive=True)
+
+    # FIXME: this is due to bug in log file analysis. What is going on here??
+    # fix that some cases do not have enough columns
+    with open(fcsv.replace('.csv', '2.csv'), 'w') as f1:
+        with open(fcsv) as f2:
+            for line in f2.readlines():
+                # older versions had 95, columns, newer 103
+                if len(line.split(';'))==96 or len(line.split(';'))==104:
+                    line = line.replace(';0.00000000000;nan;-0.0000;',
+                                        '0.00000000000;nan;-0.0000;')
+                f1.write(line)
+
     # convert from CSV to DataFrame
-    df = lf.csv2df(fcsv)
+    df = lf.csv2df(fcsv.replace('.csv', '2.csv'))
     df.to_hdf(fcsv.replace('.csv', '.h5'), 'table')
     # -------------------------------------------------------------------------
     path_pattern = os.path.join(P_RUN, 'res', '*', '*.csv')
@@ -522,15 +535,16 @@ def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
     mdf = AppendDataFrames(tqdm=tqdm)
     # individual log file analysis does not have header, make sure to include
     # a line for the header
-    mdf.txt2txt(fcsv, path_pattern, tarmode='r:xz', header=0, sep=',',
-                header_fjoined=None, recursive=True, fname_col='[case_id]')
+    cols = mdf.txt2txt(fcsv, path_pattern, tarmode='r:xz', header=0, sep=',',
+                       header_fjoined=None, recursive=True, fname_col='[case_id]')
     # and convert to df: takes 2 minutes
     fdf = fcsv.replace('.csv', '.h5')
     store = pd.HDFStore(fdf, mode='w', format='table', complevel=9,
                         complib='zlib')
-    colnames = ['channel', 'max', 'min', 'mean', 'std', 'range',
-                'absmax', 'rms', 'int', 'intabs', '[case_id]']
-    colnames.extend(['m=%1.0f' % k for k in m])
+    colnames = cols.split(',')
+    # the first column is the channel name but the header is emtpy
+    assert colnames[0] == ''
+    colnames[0] = 'channel'
     dtypes = {col:np.float64 for col in colnames}
     dtypes['channel'] = str
     dtypes['[case_id]'] = str
@@ -541,6 +555,7 @@ def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
     store.close()
     # -------------------------------------------------------------------------
     # merge missing cols onto stats
+    # FIXME: HAS TO BE THE SAME AS tags IN post_launch
     required = ['[DLC]', '[run_dir]', '[wdir]', '[Windspeed]', '[res_dir]',
                 '[case_id]', '[Case folder]']
     df = pd.read_hdf(fdf, 'table')
@@ -556,8 +571,8 @@ def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
     cc = sim.Cases(POST_DIR, sim_id)
     df_tags = cc.cases2df()
     df_stats = pd.merge(df, df_tags[required], on=['[case_id]'])
-    # if the merge didn't work due to other misaligned case_id tags, do not
-    # overwrite our otherwise ok tables!
+    # find out if we have some misalignment between generated cases and results
+    # this could happen when we added new cases and removed others
     if len(df_stats) != len(df):
         print('failed to merge required tags, something is wrong!')
         # find out which cases we lost and why
@@ -568,8 +583,12 @@ def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
         msg = 'nr of case_ids lost:'
         print(msg, (len(df)-len(df_stats))/len(df['channel'].unique()))
         print('following case_ids have mysteriously disappeared:')
-        print(s_df-s_stats)
-        return
+        missing = list(s_df-s_stats)
+        print(missing)
+        # save misalligned cases
+        fname = os.path.join(POST_DIR, '%s_misallgined_cases.tsv' % sim_id)
+        pd.DataFrame(missing).to_csv(fname, sep='\t')
+
     df_stats.to_hdf(fdf, 'table', mode='w')
     df_stats.to_csv(fdf.replace('.h5', '.csv'))
 
@@ -580,11 +599,27 @@ def postpro_node_merge(tqdm=False, zipchunks=False, m=[3,4,6,8,9,10,12]):
     fname = os.path.join(POST_DIR, '%s_unique-channel-names.csv' % sim_id)
     pd.DataFrame(chans).to_csv(fname)
 
+
+def prepare_failed(compress=False, wine_arch='win32', wine_prefix='.wine32',
+                   prelude='', zipchunks=False):
+
+    cc = sim.Cases(POST_DIR, sim_id, rem_failed=False)
+    df_tags = cc.cases2df()
+
     # -------------------------------------------------------------------------
     # find failed cases and create pbs_in_failed dir
     cc.find_failed(df_cases=df_tags)
     sim.copy_pbs_in_failedcases(cc.cases_fail, path=opt.pbs_failed_path)
 
+    if zipchunks and len(cc.cases_fail) > 0:
+        # and for chunks as well
+        sorts_on = ['[DLC]', '[Windspeed]']
+        create_chunks_htc_pbs(cc.cases_fail, sort_by_values=sorts_on,
+                              ppn=17, nr_procs_series=3, walltime='09:00:00',
+                              chunks_dir='zip-chunks-jess-fail', compress=compress,
+                              wine_arch=wine_arch, wine_prefix=wine_prefix,
+                              prelude=prelude, queue='windq', i0=1000)
+
 
 if __name__ == '__main__':
 
@@ -593,6 +628,10 @@ if __name__ == '__main__':
                         dest='prep', help='create htc, pbs, files')
     parser.add_argument('--check_logs', action='store_true', default=False,
                         dest='check_logs', help='check the log files')
+    parser.add_argument('--failed', action='store_true', default=False,
+                        dest='failed', help='Create new pbs_in files for all '
+                        'failed cases. Combine with --zipchunks to also create '
+                        'new zipchunks for the failed cases.')
     parser.add_argument('--pbs_failed_path', default='pbs_in_fail', type=str,
                         action='store', dest='pbs_failed_path',
                         help='Copy pbs launch files of the failed cases to a '
@@ -625,6 +664,12 @@ if __name__ == '__main__':
                         dest='envelopeblade', help='Compute envelopeblade')
     parser.add_argument('--envelopeturbine', default=False, action='store_true',
                         dest='envelopeturbine', help='Compute envelopeturbine')
+    parser.add_argument('--int_env', default=False, action='store_true',
+                        dest='int_env', help='Interpolate load envelopes to '
+                        'a fixed number of nx angles.')
+    parser.add_argument('--nx', type=int, default=300, action='store',
+                        dest='nx', help='Number of angles to interpolate the '
+                        'load envelopes to. Default=300.')
     parser.add_argument('--zipchunks', default=False, action='store_true',
                         dest='zipchunks', help='Create PBS launch files for'
                         'running in zip-chunk find+xargs mode.')
@@ -674,19 +719,20 @@ if __name__ == '__main__':
                         'generated DLC exchange files, default: htc/DLCs/')
     parser.add_argument('--wine_64bit', default=False, action='store_true',
                         dest='wine_64bit', help='Run wine in 64-bit mode. '
-                        'Only works on Jess.')
+                        'Only works on Jess. Sets --wine_arch and '
+                        '--wine_prefix to win64 and .wine respectively.')
     parser.add_argument('--wine_arch', action='store', default='win32', type=str,
                         dest='wine_arch', help='Set to win32 for 32-bit, and '
                         'win64 for 64-bit. 64-bit only works on Jess. '
                         'Defaults to win32.')
-    parser.add_argument('--wine_prefix', action='store', default='~/.wine32',
+    parser.add_argument('--wine_prefix', action='store', default='.wine32',
                         type=str, dest='wine_prefix', help='WINEPREFIX: '
-                        'Directory used by wineserver. Default ~/.wine32')
+                        'Directory used by wineserver. Default .wine32')
+    parser.add_argument('--linux', action='store_true', default=False,
+                        dest='linux', help='Do not use wine. Implies that '
+                        'wine_prefix and wine_arch is set to None.')
     opt = parser.parse_args()
 
-    # Wholer coefficients to be considered for the fatigue analysis
-    m = [3, 4, 6, 8, 9, 10, 12]
-
     # -------------------------------------------------------------------------
 #    # manually configure paths, HAWC2 model root path is then constructed as
 #    # p_root_remote/PROJECT/sim_id, and p_root_local/PROJECT/sim_id
@@ -716,6 +762,9 @@ if __name__ == '__main__':
     P_RUN, P_SOURCE, PROJECT, sim_id, P_MASTERFILE, MASTERFILE, POST_DIR \
         = dlcdefs.configure_dirs(verbose=True)
 
+    if opt.wine_64bit:
+        opt.wine_arch, opt.wine_prefix = ('win64', '.wine')
+
     if opt.gendlcs:
         DLB = GenerateDLCCases()
         DLB.execute(filename=os.path.join(P_SOURCE, opt.dlcmaster),
@@ -727,10 +776,10 @@ if __name__ == '__main__':
         launch_dlcs_excel(sim_id, silent=False, zipchunks=opt.zipchunks,
                           pbs_turb=opt.pbs_turb, walltime=opt.walltime,
                           postpro_node=opt.postpro_node, runmethod=RUNMETHOD,
-                          dlcs_dir=os.path.join(P_SOURCE, 'htc', 'DLCs'),
-                          compress=opt.compress, wine_64bit=opt.wine_64bit,
+                          dlcs_dir=os.path.join(P_SOURCE, opt.dlcfolder),
                           postpro_node_zipchunks=opt.no_postpro_node_zipchunks,
-                          wine_arch=opt.wine_arch, wine_prefix=opt.wine_prefix)
+                          wine_arch=opt.wine_arch, wine_prefix=opt.wine_prefix,
+                          compress=opt.compress, linux=opt.linux)
     # post processing: check log files, calculate statistics
     if opt.check_logs or opt.stats or opt.fatigue or opt.envelopeblade \
         or opt.envelopeturbine or opt.AEP:
@@ -740,11 +789,22 @@ if __name__ == '__main__':
                     fatigue=opt.fatigue, A=opt.rotarea, AEP=opt.AEP,
                     no_bins=opt.no_bins, pbs_failed_path=opt.pbs_failed_path,
                     save_new_sigs=opt.save_new_sigs, save_iter=False,
-                    envelopeturbine=opt.envelopeturbine,
-                    envelopeblade=opt.envelopeblade)
+                    envelopeturbine=opt.envelopeturbine, int_env=opt.int_env,
+                    envelopeblade=opt.envelopeblade, nx=opt.nx)
     if opt.postpro_node_merge:
-        postpro_node_merge(zipchunks=opt.zipchunks, m=m)
+        postpro_node_merge(zipchunks=opt.zipchunks)
+    if opt.failed:
+        prepare_failed(zipchunks=opt.zipchunks, compress=opt.compress,
+                       wine_arch=opt.wine_arch, wine_prefix=opt.wine_prefix)
     if opt.dlcplot:
+
+        # simply plot all channels
+        # fname = os.path.join(POST_DIR, '%s_unique-channel-names.csv' % sim_id)
+        # df = pd.read_csv(fname, header=0, names=['c1', 'c2'])
+        # plot_chans = {}
+        # for i, k in enumerate(df['c2'].values.tolist()):
+        #     plot_chans['i=%i' % i] = [k]
+
         plot_chans = {}
         plot_chans['$B1_{flap}$'] = ['setbeta-bladenr-1-flapnr-1']
         plot_chans['$B2_{flap}$'] = ['setbeta-bladenr-2-flapnr-1']
@@ -767,14 +827,12 @@ if __name__ == '__main__':
         plot_chans['$B123_{pitch}$'] = chans
 
         plot_chans['RPM'] = ['bearing-shaft_rot-angle_speed-rpm']
-        plot_chans['$P_e$'] = ['DLL-2-inpvec-2']
+        plot_chans['$P_e$'] = ['DLL-generator_servo-inpvec-2']
         plot_chans['$P_{mech}$'] = ['stats-shaft-power']
         plot_chans['$B3 U_y$'] = ['global-blade3-elem-018-zrel-1.00-State pos-y']
         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['TC blade to tower'] = ['DLL-5-inpvec-2']
-        plot_chans['TC tower to blade'] = ['DLL-5-inpvec-3']
         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']

wetb/prepost/filters.py

@@ -4,17 +4,6 @@ Created on Sun Jan 20 18:14:02 2013
 
 @author: dave
 """
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import absolute_import
-from future import standard_library
-standard_library.install_aliases()
-from builtins import object
-
-
-
-
 import numpy as np
 import scipy as sp
 
@@ -464,4 +453,4 @@ class Filters(object):
         ax3.set_xlabel('t')
         ax3.grid(True)
 
-        plot.save_fig()
\ No newline at end of file
+        plot.save_fig()

wetb/prepost/h2_vs_hs2.py

@@ -4,20 +4,8 @@ Created on Mon Nov  2 15:23:15 2015
 
 @author: dave
 """
-from __future__ import division
-from __future__ import print_function
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import range
-from builtins import zip
-from builtins import dict
-from builtins import str
-from builtins import int
-from future import standard_library
-standard_library.install_aliases()
-from builtins import object
-
 import os
+from copy import copy
 
 import numpy as np
 #import scipy.interpolate as interpolate
@@ -100,7 +88,7 @@ class ConfigBase(object):
         """
         opt_tags = [self.opt_h2.copy()]
         opt_tags[0].update(self.eigenan.copy())
-        opt_tags[0]['[Case id.]'] = '%s_hawc2_eigenanalysis' % basename
+        opt_tags[0]['[case_id]'] = '%s_hawc2_eigenanalysis' % basename
         opt_tags[0]['[blade_damp_x]'] = 0.0
         opt_tags[0]['[blade_damp_y]'] = 0.0
         opt_tags[0]['[blade_damp_z]'] = 0.0
@@ -111,7 +99,7 @@ class ConfigBase(object):
         opt_tags[0]['[eigen_analysis]'] = True
         opt_tags[0]['[output]'] = False
         opt_tags[0]['[t0]'] = 0.0
-        opt_tags[0]['[time stop]'] = 0.0
+        opt_tags[0]['[time_stop]'] = 0.0
 
         return opt_tags
 
@@ -120,7 +108,7 @@ class ConfigBase(object):
         analysis, at 0 RPM.
         """
         opt_tags = [self.opt_hs2.copy()]
-        opt_tags[0]['[Case id.]'] = '%s_hs2_eigenanalysis' % basename
+        opt_tags[0]['[case_id]'] = '%s_hs2_eigenanalysis' % basename
         opt_tags[0]['[blade_damp_x]'] = 0.0
         opt_tags[0]['[blade_damp_y]'] = 0.0
         opt_tags[0]['[blade_damp_z]'] = 0.0
@@ -136,7 +124,7 @@ class ConfigBase(object):
     def opt_tags_hs2(self, basename):
 
         opt_tags = [self.opt_hs2.copy()]
-        opt_tags[0]['[Case id.]'] = '%s_hawcstab2' % basename
+        opt_tags[0]['[case_id]'] = '%s_hawcstab2' % basename
         return opt_tags
 
     def set_hs2opdata(self, master, basename):
@@ -150,10 +138,10 @@ class ConfigBase(object):
         fpath = os.path.join(master.tags['[data_dir]'],
                              master.tags['[operational_data]'])
         hs2_res = hs2.results()
-        hs2_res.load_operation(fpath)
-        omegas = hs2_res.operation.rotorspeed_rpm.values*np.pi/30.0
-        winds = hs2_res.operation.windspeed.values
-        pitchs = -1.0*hs2_res.operation.pitch_deg.values
+        operation = hs2_res.load_operation(fpath)
+        omegas = operation.rotorspeed_rpm.values*np.pi/30.0
+        winds = operation.windspeed.values
+        pitchs = -1.0*operation.pitch_deg.values
 
         return self.set_opdata(winds, pitchs, omegas, basename=basename)
 
@@ -174,9 +162,9 @@ class ConfigBase(object):
             rotor speed at given operating point [rad/s]
 
         basename : str, default=None
-            If not None, the [Case id.] tag is composed out of the basename,
+            If not None, the [case_id] tag is composed out of the basename,
             wind speed, pitch angle and rotor speed. If set to None, the
-            [Case id.] tag is not set.
+            [case_id] tag is not set.
 
         Returns
         -------
@@ -192,13 +180,13 @@ class ConfigBase(object):
             rpl = (basename, wind, pitch, omega)
             if basename is not None:
                 tmp = '%s_%02.0fms_%04.01fdeg_%04.02frads_hawc2' % rpl
-                opt_dict['[Case id.]'] = tmp
+                opt_dict['[case_id]'] = tmp
             opt_dict['[Windspeed]'] = wind
             opt_dict['[blade_pitch_deg]'] = pitch
             opt_dict['[fixspeed_rotor_rads]'] = omega
             opt_dict['[initspeed_rotor_rads]'] = omega
 #            opt_dict['[t0]'] = int(2000.0/opt_dict['[Windspeed]']) # or 2000?
-#            opt_dict['[time stop]'] = opt_dict['[t0]']+100
+#            opt_dict['[time_stop]'] = opt_dict['[t0]']+100
 #            opt_dict['[time_stop]'] = opt_dict['[t0]']+100
             opt_tags.append(opt_dict.copy())
         return opt_tags
@@ -286,28 +274,29 @@ class Sims(object):
 
         """
 
-        mt = master.tags
+        master_copy = copy(master)
+        mt = master_copy.tags
 
         dlc_case = mt['[Case folder]']
         mt['[data_dir]'] = 'data/'
         mt['[res_dir]'] = 'res/%s/' % dlc_case
         mt['[log_dir]'] = 'logfiles/%s/' % dlc_case
         mt['[htc_dir]'] = 'htc/%s/' % dlc_case
-        mt['[case_id]'] = mt['[Case id.]']
+        mt['[case_id]'] = mt['[case_id]']
         mt['[DLC]'] = dlc_case
         mt['[pbs_out_dir]'] = 'pbs_out/%s/' % dlc_case
         mt['[pbs_in_dir]'] = 'pbs_in/%s/' % dlc_case
         mt['[iter_dir]'] = 'iter/%s/' % dlc_case
 
         if mt['[eigen_analysis]']:
-            rpl = os.path.join(dlc_case, mt['[Case id.]'])
+            rpl = os.path.join(dlc_case, mt['[case_id]'])
             mt['[eigenfreq_dir]'] = 'res_eigen/%s/' % rpl
 
         # for HAWCStab2 certain things have to be done differently
         if mt['[hs2]']:
             mt['[htc_dir]'] = ''
             mt['[t0]'] = 0
-            mt['[time stop]'] = 1
+            mt['[time_stop]'] = 1
             mt['[hawc2]'] = False
             mt['[output]'] = False
             mt['[copyback_files]'] = ['./*.ind', './*.pwr', './*.log',
@@ -336,10 +325,10 @@ class Sims(object):
                 mt['[hs2_gradients]'] = 'nogradients'
 
         mt['[windspeed]'] = mt['[Windspeed]']
-        mt['[time_stop]'] = mt['[time stop]']
+        mt['[time_stop]'] = mt['[time_stop]']
         mt['[duration]'] = str(float(mt['[time_stop]']) - float(mt['[t0]']))
 
-        return master
+        return master_copy
 
     def _set_path_auto_config(self, verbose=True):
         """
@@ -414,7 +403,7 @@ class Sims(object):
             opt['[zip_root_files]'] = f_ziproot
 
         self.master.output_dirs.extend('[Case folder]')
-        self.master.output_dirs.extend('[Case id.]')
+        self.master.output_dirs.extend('[case_id]')
 
         return iter_dict, opt_tags
 
@@ -425,7 +414,8 @@ class Sims(object):
                            write_htc=True, runmethod=runmethod, verbose=False,
                            copyback_turb=False, msg='', update_cases=False,
                            ignore_non_unique=False, run_only_new=False,
-                           pbs_fname_appendix=False, short_job_names=False)
+                           pbs_fname_appendix=False, short_job_names=False,
+                           windows_nr_cpus=2, update_model_data=True)
 
     def get_control_tuning(self, fpath):
         """
@@ -540,9 +530,9 @@ class MappingsH2HS2(object):
 
     def _powercurve_hs2(self, fname):
 
-        mappings = {'P [kW]'  :'P_aero',
-                    'T [kN]'  :'T_aero',
-                    'V [m/s]' :'windspeed'}
+        mappings = {'P'  :'P_aero',
+                    'T'  :'T_aero',
+                    'V' :'windspeed'}
 
         df_pwr, units = self.hs2_res.load_pwr_df(fname)
 
@@ -553,21 +543,28 @@ class MappingsH2HS2(object):
     def blade_distribution(self, fname_h2, fname_hs2, h2_df_stats=None,
                            fname_h2_tors=None):
 
-        self.hs2_res.df_ind = self.hs2_res.load_ind(fname_hs2)
+        self.df_ind = self.hs2_res.load_ind(fname_hs2)
         self.h2_res = sim.windIO.ReadOutputAtTime(fname_h2)
         self._distribution_hs2()
         self._distribution_h2()
         if h2_df_stats is not None:
             self.h2_df_stats = h2_df_stats
             if fname_h2_tors is not None:
-                self.distribution_stats_h2(fname_h2_tors, 'Tors_e', 'torsion')
+                # self.distribution_stats_h2(fname_h2_tors, 'Tors_e', 'torsion',
+                #                            xaxis='radius')
+                self.distribution_stats_h2(fname_h2_tors, 'statevec_new', 'torsion',
+                                           component='Rz', xaxis='s')
+                self.distribution_stats_h2(fname_h2_tors, 'statevec_new', 'def_x_svn',
+                                            component='Dx', xaxis='s')
+                self.distribution_stats_h2(fname_h2_tors, 'statevec_new', 'def_y_svn',
+                                            component='Dy', xaxis='s')
 
     def _distribution_hs2(self):
         """Read a HAWCStab2 *.ind file (blade distribution loading)
 
         rot_angle and rot_vec_123 in HS2 should be in rotor polar coordinates
         """
-
+        # mapping_hs2[hawcstab2_key] = general_key
         mapping_hs2 =  {'s [m]'       :'curved_s',
                         'CL0 [-]'     :'Cl',
                         'CD0 [-]'     :'Cd',
@@ -601,7 +598,7 @@ class MappingsH2HS2(object):
             hs2_cols = list(mapping_hs2)
             # select only the HS channels that will be used for the mapping
             std_cols = list(mapping_hs2.values())
-            self.hs_aero = self.hs2_res.df_ind[hs2_cols].copy()
+            self.hs_aero = self.df_ind[hs2_cols].copy()
         except KeyError:
             # some results have been created with older HAWCStab2 that did not
             # include CT and CP columns
@@ -610,7 +607,7 @@ class MappingsH2HS2(object):
             hs2_cols = list(mapping_hs2)
             std_cols = list(mapping_hs2.values())
             # select only the HS channels that will be used for the mapping
-            self.hs_aero = self.hs2_res.df_ind[hs2_cols].copy()
+            self.hs_aero = self.df_ind[hs2_cols].copy()
 
         # change column names to the standard form that is shared with H2
         self.hs_aero.columns = std_cols
@@ -623,6 +620,7 @@ class MappingsH2HS2(object):
         self.hs_aero['twist'] *= (180.0/np.pi)
 
     def _distribution_h2(self):
+        # mapping_h2[hawc2_key] = general_key
         mapping_h2 =  { 'Radius_s'  :'curved_s',
                         'Cl'        :'Cl',
                         'Cd'        :'Cd',
@@ -633,9 +631,9 @@ class MappingsH2HS2(object):
                         'Vrel'      :'vrel',
                         'Inflow_ang':'inflow_angle',
                         'alfa'      :'AoA',
-                        'pos_RP_x'  :'pos_x',
-                        'pos_RP_y'  :'pos_y',
-                        'pos_RP_z'  :'pos_z',
+                        'Pos_RP_x'  :'pos_x',
+                        'Pos_RP_y'  :'pos_y',
+                        'Pos_RP_z'  :'pos_z',
                         'Chord'     :'chord',
                         'Secfrc_RPx':'F_x',
                         'Secfrc_RPy':'F_y',
@@ -647,9 +645,13 @@ class MappingsH2HS2(object):
         h2_aero = self.h2_res[h2_cols].copy()
         # change column names to the standard form that is shared with HS
         h2_aero.columns = std_cols
+
+        # using the distributed aero outputs
         h2_aero['def_x'] = self.h2_res['Pos_B_x'] - self.h2_res['Inipos_x_x']
         h2_aero['def_y'] = self.h2_res['Pos_B_y'] - self.h2_res['Inipos_y_y']
         h2_aero['def_z'] = self.h2_res['Pos_B_z'] - self.h2_res['Inipos_z_z']
+        # note that distribution_stats_h2() will set def_x from statevec_new
+
         h2_aero['ax_ind_vel'] *= (-1.0)
 #        h2_aero['pos_x'] += (self.h2_res['Chord'] / 2.0)
         h2_aero['F_x'] *= (1e3)
@@ -661,7 +663,8 @@ class MappingsH2HS2(object):
 #        h2_aero = h2_aero[1:-1]
         self.h2_aero = h2_aero
 
-    def distribution_stats_h2(self, fname_h2, sensortype, newname):
+    def distribution_stats_h2(self, fname_h2, sensortype, newname, xaxis='s',
+                              component=None):
         """Determine blade distribution sensor from the HAWC2 statistics.
         This requires that for each aerodynamic calculation point there should
         be an output sensor defined manually in the output section.
@@ -675,6 +678,14 @@ class MappingsH2HS2(object):
 
         newname
 
+        xaxis : string
+            column name in LoadResults.ch_df identifying the blade radial
+            coordinate. Valid values are 's', 'radius', 'radius_actual',
+            'srel'
+
+        component : string
+            For statevec_new we also need to specify which component, Dx, Rx, etc.
+
         """
         if not hasattr(self, 'h2_aero'):
             raise UserWarning('first run blade_distribution')
@@ -682,36 +693,41 @@ class MappingsH2HS2(object):
         # load the HAWC2 .sel file for the channels
         fpath = os.path.dirname(fname_h2)
         fname = os.path.basename(fname_h2)
-        res = sim.windIO.LoadResults(fpath, fname, readdata=False)
+        res = sim.windIO.LoadResults(fpath, fname, readdata=True)
         sel = res.ch_df[res.ch_df.sensortype == sensortype].copy()
+        if component is not None:
+            sel = sel[sel['component']=='Rz']
         if len(sel) == 0:
             msg = 'HAWC2 sensor type "%s" is missing, are they defined?'
             raise ValueError(msg % sensortype)
-        sel.sort_values(['radius'], inplace=True)
-        tors_e_channels = sel.unique_ch_name.tolist()
+        # for backward compatibilitiy with tors_e
+        sel.sort_values([xaxis], inplace=True)
+        chan_sel = sel.unique_ch_name.tolist()
 
         # find the current case in the statistics DataFrame
         case = fname.replace('.htc', '')
         df_case = self.h2_df_stats[self.h2_df_stats['[case_id]']==case].copy()
         # and select all the torsion channels
-        df_tors_e = df_case[df_case.channel.isin(tors_e_channels)].copy()
+        df_chan_sel = df_case[df_case.channel.isin(chan_sel)].copy()
         # join the stats with the channel descriptions DataFrames, have the
         # same name on the joining column
-        df_tors_e.set_index('channel', inplace=True)
+        df_chan_sel.set_index('channel', inplace=True)
         sel.set_index('unique_ch_name', inplace=True)
 
         # joining happens on the index, and for which the same channel has been
         # used: the unique HAWC2 channel naming scheme
-        df_tors_e = pd.concat([df_tors_e, sel], axis=1)
-        df_tors_e.radius = df_tors_e.radius.astype(np.float64)
-        # sorting on radius, combine with ch_df
-        df_tors_e.sort_values(['radius'], inplace=True)
+        df_chan_sel = pd.concat([df_chan_sel, sel], axis=1)
+        df_chan_sel[xaxis] = df_chan_sel[xaxis].astype(np.float64)
+        # sorting on radius/s/etc, combine with ch_df
+        df_chan_sel.sort_values([xaxis], inplace=True)
+
+        # TODO: check if the outputs are at the same positions as the aero outputs
 
         # FIXME: what if number of torsion outputs is less than aero
         # calculation points?
-        self.h2_aero['%s' % newname] = df_tors_e['mean'].values.copy()
-        self.h2_aero['%s_std' % newname] = df_tors_e['std'].values.copy()
-        self.h2_aero['%s_radius_s' % newname] = df_tors_e['radius'].values.copy()
+        self.h2_aero['%s' % newname] = df_chan_sel['mean'].values.copy()
+        self.h2_aero['%s_std' % newname] = df_chan_sel['std'].values.copy()
+        self.h2_aero['%s_%s_s' % (newname, xaxis)] = df_chan_sel[xaxis].values.copy()
 
     def body_structure_modes(self, fname_h2, fname_hs):
         self._body_structure_modes_h2(fname_h2)
@@ -801,7 +817,14 @@ class Plots(object):
         if h2_df_stats is not None:
             res.h2_df_stats = h2_df_stats
             if fname_h2_tors is not None:
-                res.distribution_stats_h2(fname_h2_tors, 'Tors_e', 'torsion')
+                # res.distribution_stats_h2(fname_h2_tors, 'Tors_e', 'torsion',
+                #                           xaxis='radius')
+                res.distribution_stats_h2(fname_h2_tors, 'statevec_new', 'torsion',
+                                          component='Rz', xaxis='s')
+                res.distribution_stats_h2(fname_h2_tors, 'statevec_new', 'def_x_svn',
+                                          component='Dx', xaxis='s')
+                res.distribution_stats_h2(fname_h2_tors, 'statevec_new', 'def_y_svn',
+                                          component='Dy', xaxis='s')
 
         return res
 

wetb/prepost/hawcstab2.py

-# -*- coding: utf-8 -*-
 """
 Created on Tue Jan 14 14:12:58 2014
 
 @author: dave
 """
-from __future__ import division
-from __future__ import print_function
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import range
-from io import open
-from builtins import int
-from future import standard_library
-standard_library.install_aliases()
-from builtins import object
-
 import os
 import re
 
@@ -33,22 +21,34 @@ regex_units = re.compile('(\\[.*?\\])')
 
 def ReadFileHAWCStab2Header(fname):
     """
-    Read a file with a weird HAWCStab2 header that starts with a #, and
-    includes the column number and units between square brackets.
+    Read a file with a weird HAWCStab2 header that starts with a # (or not),
+    and includes the column number and units between square brackets. Can also
+    strip away the first integer number on the header line (e.g. .opt files).
     """
 
-    def get_lines(fname):
+    def get_lines():
         # get the line that contains the header/column names and the first
         # line that holds the data
-        with open(fname) as f:
-            line_header = f.readline()
+        line_header = f.readline()
+        seek_data = f.tell()
+        line_data = f.readline()
+        # the opt file has a number of data points as the first element
+        nrl = line_header.split()[0]
+        try:
+            int(nrl)
+            # replace with empty text to no destroy the nice line-up of
+            # the header with the data columns, assumed later
+            line_header = line_header.replace(nrl, ' '*len(nrl))
+        except ValueError:
+            pass
+        # sometimes there are more header lines. The header is always on the
+        # last of the lines marked with #
+        while line_data[:2].strip() == '#':
+            line_header = line_data
+            seek_data = f.tell()
             line_data = f.readline()
-            # sometimes there are more header lines. The header is always on the
-            # last of the lines marked with #
-            while line_data[:2].strip() == '#':
-                line_header = line_data
-                line_data = f.readline()
-        return line_header, line_data
+
+        return line_header, line_data, seek_data
 
     def get_col_widths(line):
         # it is very annoying that various files can have various column widths
@@ -68,28 +68,41 @@ def ReadFileHAWCStab2Header(fname):
         ci[1:] = ci[1:] - 1
         columns = []
         for i in range(len(ci)-1):
-            # also lose the index in the header
-            colname = line[ci[i]:ci[i+1]][:-2].replace('#', '').strip()
+            # also lose the index in the header, if there is one
+            try:
+                int(line.split()[-1])
+                colname = line[ci[i]:ci[i+1]][:-2].replace('#', '').strip()
+            except ValueError:
+                colname = line[ci[i]:ci[i+1]].replace('#', '').strip()
             columns.append(colname)
         return columns
 
-    line_header, line_data = get_lines(fname)
-    colwidths = get_col_widths(line_data)
-    columns = get_col_names(line_header, colwidths)
-    # gradients have duplicate columns: set for with wake updated
-    # and another with frozen wake assumption, append _fw to the columns
-    # used to indicate frozen wake gradients
-    if 'dQ/dt [kNm/deg]' in columns:
-        i1 = columns.index('dQ/dt [kNm/deg]')
-        if i1 > -1:
-            i2 = columns.index('dQ/dt [kNm/deg]', i1+1)
-        if i2 > i1:
-            for i in range(i2, len(columns)):
-                columns[i] = columns[i].replace(' [', '_fw [')
-
-    df = pd.read_fwf(fname, widths=colwidths, comment='#', header=None,
-                     names=columns)
-    units = regex_units.findall(''.join(columns))
+    with open(fname) as f:
+        line_header, line_data, seek_data = get_lines()
+        colwidths = get_col_widths(line_data)
+        columns = get_col_names(line_header, colwidths)
+        units = regex_units.findall(''.join(columns))
+        # strip units from the column name, units is dictionary
+        columns = [k.split('[')[0].strip() for k in columns]
+        units = {col:unit for col, unit in zip(columns, units)}
+
+        # gradients have duplicate columns: set for with wake updated
+        # and another with frozen wake assumption, append _fw to the columns
+        # used to indicate frozen wake gradients
+        # the doubles occur at the end and are the following:
+        # 'dQ/dt', 'dQ/dV', 'dQ/dO', 'dT/dt', 'dT/dV', 'dT/dO'
+        if 'dQ/dt' in columns:
+            i1 = columns.index('dQ/dt')
+            if i1 > -1:
+                i2 = columns.index('dQ/dt', i1+1)
+            if i2 > i1:
+                for i in range(i2, len(columns)):
+                    columns[i] = columns[i] + '_fw'
+
+        # we are at the second data line, go back to the start of the data
+        f.seek(seek_data)
+        df = pd.read_fwf(f, widths=colwidths, comment='#', header=None,
+                         names=columns)
 
     return df, units
 
@@ -150,7 +163,7 @@ class results(object):
         pass
 
     def load_pwr(self, fname):
-        pwr = np.loadtxt(fname)
+        pwr = np.atleast_2d(np.loadtxt(fname))
 
         res = dummy()
 
@@ -233,12 +246,30 @@ class results(object):
 
         with open(fname) as f:
             line = f.readline()
+            f.readline()
+            line3 = f.readline()
+            line4 = f.readline() # first data entry
+
+        # assuming first we have: "# Mode number:"
+        mode_nrs = line.split()[3:]
+
+        elements_line1 = line4.split()
+#        nrcols = len(elements_line1)
+        element1 = elements_line1[0]
+        width_col1 = line4.find(element1) + len(element1)
+
+#        width = 14
+#        nrcols = int((len(line)-1)/width)
+#        # first columns has one extra character
+#        # col nr1: rotor speed, col nr2: radius
+#        widths = [width+1] + [width]*(nrcols-1)
+
+        # the rotor/wind speed column does not have a unit
+        ilocs = [0, width_col1-1] + [m.start() for m in re.finditer('\]', line3)]
+        widths = np.diff(np.array(ilocs))
+        # the length of the first element should be +1 due to zero based index
+        widths[0] += 1
 
-        width = 14
-        nrcols = int((len(line)-1)/width)
-        # first columns has one extra character
-        # col nr1: rotor speed, col nr2: radius
-        widths = [width+1] + [width]*(nrcols-1)
         # last line is empty
         df = pd.read_fwf(fname, header=2, widths=widths, skipfooter=1)
         units = regex_units.findall(''.join(df.columns))
@@ -246,11 +277,15 @@ class results(object):
         # since U_x, u_y, phase and theta will be repeated as many times as
         # there are modes, add the mode number in the column name
         columns = [k.replace('#', '').strip() for k in df.columns]
-        nrmodes = int((len(columns) - 2 )/6)
-        for k in range(nrmodes):
-            for i in range(6):
-                j = 2+k*6+i
-                columns[j] = columns[j].split('.')[0] + ' nr%i' % (k+1)
+#        nrmodes = int((len(columns) - 2 )/6)
+#        for k in range(nrmodes):
+#            for i in range(6):
+#                j = 2+k*6+i
+#                columns[j] = columns[j].split('.')[0] + ' nr%i' % (k+1)
+
+        for i, k in enumerate(mode_nrs):
+            columns[i+2] = columns[i+2].split('.')[0] + ' nr%s' % k
+
         df.columns = columns
 
         return df, units
@@ -261,19 +296,19 @@ class results(object):
         # when the array is empty, set operation to an empty DataFrame
         if len(operation) == 0:
             cols = ['windspeed', 'pitch_deg', 'rotorspeed_rpm']
-            self.operation = pd.DataFrame(columns=cols)
-            return
+            return pd.DataFrame(columns=cols)
         # when there is only one data point, the array is 1D, we allways need
         # a 2D array otherwise the columns become rows in the DataFrame
         elif len(operation.shape) == 1:
             operation = operation.reshape((1, operation.shape[0]))
         try:
             cols = ['windspeed', 'pitch_deg', 'rotorspeed_rpm']
-            self.operation = pd.DataFrame(operation, columns=cols)
+            operation = pd.DataFrame(operation, columns=cols)
         except ValueError:
             cols = ['windspeed', 'pitch_deg', 'rotorspeed_rpm', 'P_aero',
                     'T_aero']
-            self.operation = pd.DataFrame(operation, columns=cols)
+            operation = pd.DataFrame(operation, columns=cols)
+        return operation
 
     def load_matrices(self, fpath, basename, operating_point=1,
                       control_mat=False, local_wind_mat=False):
@@ -335,84 +370,40 @@ class results(object):
 
         return matrices
 
+    def load_cntrl_tuning(self, fname):
+        """Load contruller tuning file.
+        """
+
+        tune = ReadControlTuning()
+        tune.read_parameters(fname)
+        tuning = tune.parameters2tags()
+        return tuning
+
     def write_ae_sections_h2(self):
         """
         Get the aerosection positions from the HS2 ind result file and
         write them as outputs for HAWC2
+
+        NOT IMPLEMENTED YET
         """
         self.ind
 
-    def plot_pwr(self, figname, fnames, labels=[], figsize=(11,7.15), dpi=120):
-
-        results = []
-        if isinstance(fnames, list):
-            if len(fnames) > 4:
-                raise ValueError('compare up to maximum 4 HawcStab2 cases')
-            for fname in fnames:
-                results.append(self.load_pwr(fname))
-                # if the labels are not defined, take the file name
-                if len(labels) < len(fnames):
-                    labels.append(os.path.basename(fname))
-        else:
-            results.append(self.load_pwr(fname))
-
-        colors = list('krbg')
-        symbols = list('o<+x')
-        alphas = [1.0, 0.9, 0.8, 0.75]
-
-        fig, axes = mplutils.subplots(nrows=2, ncols=2, figsize=figsize,
-                                       dpi=dpi, num=0)
-        for i, res in enumerate(results):
-            ax = axes[0,0]
-            ax.plot(res.wind, res.power, color=colors[i],
-                    label='Power %s ' % labels[i],
-                    marker=symbols[i], ls='-', alpha=alphas[i])
-            ax.set_title('Aerodynamic Power [kW]')#, RPM')
-
-            ax = axes[0,1]
-            ax.plot(res.wind, res.pitch_deg, color=colors[i],
-                    label='Pitch %s' % labels[i],
-                    marker=symbols[i], ls='-', alpha=alphas[i])
-            ax.plot(res.wind, res.rpm, color=colors[i],
-                    label='RPM %s ' % labels[i],
-                    marker=symbols[i], ls='--', alpha=alphas[i])
-            ax.set_title('Pitch [deg], RPM')
-
-            ax = axes[1,0]
-            ax.plot(res.wind, res.thrust, color=colors[i], label=labels[i],
-                    marker=symbols[i], ls='-', alpha=alphas[i])
-            ax.set_title('Thrust [kN]')
-
-            ax = axes[1,1]
-            ax.plot(res.wind, res.cp, label='$C_p$ %s ' % labels[i], ls='-',
-                    color=colors[i], marker=symbols[i], alpha=alphas[i])
-            ax.plot(res.wind, res.ct, label='$C_t$ %s ' % labels[i], ls='--',
-                    color=colors[i], marker=symbols[i], alpha=alphas[i])
-            ax.set_title('Power and Thrust coefficients [-]')
-
-        for ax in axes.ravel():
-            ax.legend(loc='best')
-            ax.grid(True)
-            ax.set_xlim([res.wind[0], res.wind[-1]])
-        fig.tight_layout()
-
-        print('saving figure: %s ... ' % figname, end='')
-        figpath = os.path.dirname(figname)
-        if not os.path.exists(figpath):
-            os.makedirs(figpath)
-        fig.savefig(figname)
-        fig.clear()
-        print('done!')
-
 
 class ReadControlTuning(object):
 
     def __init__(self):
         """
         """
-        pass
+        self._aerogains = False
 
     def parse_line(self, line, controller):
+        """Parses the output lines with the controller tuning parameters.
+        Does not parse the aerodynamic gain lines.
+        """
+
+        if line.startswith('Aerodynamic gains'):
+            self._aerogains = True
+            return
 
         split1 = line.split('=')
         var1 = split1[0].strip()
@@ -449,8 +440,23 @@ class ReadControlTuning(object):
                 elif i == 10:
                     controller = 'aero_damp'
                     setattr(self, controller, dummy())
-                else:
+                elif not self._aerogains:
                     self.parse_line(line, controller)
+                # do not break since we want to see if aero gains are included
+                # elif self._aerogains:
+                #     break
+
+        self.aero_gains_units = ['[deg]', '[kNm/deg]', '[kNm/deg]',
+                                 '[kNm/(rad/s)]', '[kNm/(rad/s)]']
+        self.aero_gains = pd.DataFrame()
+        # in case the gains are missing from the file, don't try to read it
+        if i > 17:
+            arr = np.loadtxt(fpath, skiprows=17)
+            columns = ['theta', 'dq/dtheta', 'dq/dtheta_fit', 'dq/domega',
+                       'dq/domega_fit']
+            # sometimes not all columns are present here, only take those that
+            # are present
+            self.aero_gains = pd.DataFrame(arr, columns=columns[:arr.shape[1]])
 
         # set some parameters to zero for the linear case, or when aerodynamic
         # gain scheduling is not used
@@ -494,6 +500,369 @@ class ReadControlTuning(object):
         return tune_tags
 
 
+def read_modid(fname):
+    """Separate text file describing the modes.
+    """
+
+    df_modes = pd.read_csv(fname, comment='#', delimiter=';',
+                           header=None, names=['mode_nr', 'description'],
+                           converters={0:lambda x:x.strip(),
+                                       1:lambda x:x.strip()})
+    return df_modes['description'].values.tolist()
+
+
+def read_cmb_all(f_cmb, f_pwr=None, f_modid=None, f_save=None, ps=[1,3,6]):
+    """Convenience method to load Campbell and Performance data in one go
+
+    Parameters
+    ----------
+
+    f_cmb : str
+
+    f_pwr : str, default=None
+
+    f_modid : str, default=None
+
+    f_save : str, default=None
+
+    ps : list, default=[1,2,3]
+
+
+    Returns
+    -------
+
+    df_perf, df_freq, df_damp
+
+    """
+
+    hs2 = results()
+
+    # create DataFrames with the freq/damp on the column name
+    wind, freqs, damps, real_eig = hs2.load_cmb(f_cmb)
+    nr_oper, nr_modes = freqs.shape
+
+    # performance indicators, optional
+    df_perf = None
+    if f_pwr is not None:
+        df_perf, units = hs2.load_pwr_df(f_pwr)
+        nr_perf = df_perf.shape[0]
+        if nr_oper != nr_perf or not np.allclose(wind, df_perf['V'].values):
+            raise UserWarning('pwr and cmb files must have same operating points')
+
+    # strip characters if there is a comment after the description
+    if f_modid is None or not os.path.isfile(f_modid):
+        modes_descr = ['{:02d}'.format(k+1) for k in range(nr_modes)]
+    else:
+        modes_descr = read_modid(f_modid)[:nr_modes]
+
+    df_freq = pd.DataFrame(freqs, columns=modes_descr, index=wind)
+    df_damp = pd.DataFrame(damps, columns=modes_descr, index=wind)
+
+    if f_save is not None:
+        tmp = pd.DataFrame(freqs, columns=modes_descr, index=wind)
+        tmp.index.name = 'windspeed'
+        # add p in frequencies
+        if f_pwr is not None:
+            for p in ps:
+                tmp[f'{p}P'] = p*df_perf['Speed'].values / 60
+        # sort columns on mean frequeny over wind speeds
+        icolsort = tmp.values.mean(axis=0).argsort()
+        tmp = tmp[tmp.columns[icolsort]]
+        tmp.to_excel(f_save + '_freqs.xlsx')
+
+        tmp = pd.DataFrame(damps, columns=modes_descr, index=wind)
+        tmp.index.name = 'windspeed'
+        tmp.to_excel(f_save + '_damps.xlsx')
+
+    return df_perf, df_freq, df_damp
+
+
+def plot_pwr(figname, fnames, labels=[], figsize=(11,7.15), dpi=120):
+
+    hs2res = results()
+    reslist = []
+    if isinstance(fnames, list):
+        if len(fnames) > 4:
+            raise ValueError('compare up to maximum 4 HawcStab2 cases')
+        for fname in fnames:
+            reslist.append(hs2res.load_pwr(fname))
+            # if the labels are not defined, take the file name
+            if len(labels) < len(fnames):
+                labels.append(os.path.basename(fname))
+    else:
+        reslist.append(hs2res.load_pwr(fname))
+
+    colors = list('krbg')
+    symbols = list('o<+x')
+    alphas = [1.0, 0.9, 0.8, 0.75]
+
+    fig, axes = mplutils.subplots(nrows=2, ncols=2, figsize=figsize,
+                                   dpi=dpi, num=0)
+    for i, res in enumerate(reslist):
+        ax = axes[0,0]
+        ax.plot(res.wind, res.power, color=colors[i],
+                label='Power %s ' % labels[i],
+                marker=symbols[i], ls='-', alpha=alphas[i])
+        ax.set_title('Aerodynamic Power [kW]')#, RPM')
+
+        ax = axes[0,1]
+        ax.plot(res.wind, res.pitch_deg, color=colors[i],
+                label='Pitch %s' % labels[i],
+                marker=symbols[i], ls='-', alpha=alphas[i])
+        ax.plot(res.wind, res.rpm, color=colors[i],
+                label='RPM %s ' % labels[i],
+                marker=symbols[i], ls='--', alpha=alphas[i])
+        ax.set_title('Pitch [deg], RPM')
+
+        ax = axes[1,0]
+        ax.plot(res.wind, res.thrust, color=colors[i], label=labels[i],
+                marker=symbols[i], ls='-', alpha=alphas[i])
+        ax.set_title('Thrust [kN]')
+
+        ax = axes[1,1]
+        ax.plot(res.wind, res.cp, label='$C_p$ %s ' % labels[i], ls='-',
+                color=colors[i], marker=symbols[i], alpha=alphas[i])
+        ax.plot(res.wind, res.ct, label='$C_t$ %s ' % labels[i], ls='--',
+                color=colors[i], marker=symbols[i], alpha=alphas[i])
+        ax.set_title('Power and Thrust coefficients [-]')
+
+    for ax in axes.ravel():
+        ax.legend(loc='best')
+        ax.grid(True)
+        ax.set_xlim([res.wind[0], res.wind[-1]])
+    fig.tight_layout()
+
+    print('saving figure: %s ... ' % figname, end='')
+    figpath = os.path.dirname(figname)
+    if len(figpath)>0 and not os.path.exists(figpath):
+        os.makedirs(figpath)
+    fig.savefig(figname)
+    fig.clear()
+    print('done!')
+
+
+class PlotCampbell(object):
+    """
+    Generic base class for a Campbell diagram plot. This class assumes a
+    pandas DataFrame df_freq and df_damp is set by a subclass that deals
+    with the tool dependen file IO. The following DataFrames are requried:
+        * df_freq[operating points, modes] holds the respective mode frequencies
+        * df_damp[operating points, modes] holds the respective mode damping
+        * df_perf['wind'] holding the relevant wind speeds
+
+    The column names are used as labels in the Campbell diagram.
+
+    """
+
+    alpha_box = 0.5
+
+    def __init__(self, wind, df_freq, df_damp):
+        self.wind = wind
+        self.df_freq = df_freq
+        self.df_damp = df_damp
+
+    def _inplot_label_pos(self, nr_xpos, nr_series, xpos):
+        """
+        Generate sensible label positions
+        """
+
+        if xpos == 'random':
+            pos = np.random.randint(1, nr_xpos-5, nr_series)
+        elif xpos == 'centre':
+            pos = np.zeros((nr_series,), dtype=int)
+            pos[0:len(pos):2] = np.ceil(nr_xpos/4.0)
+            pos[1:len(pos):2] = np.ceil(2.0*nr_xpos/4.0)
+        elif xpos == 'borders':
+            pos = np.zeros((nr_series,), dtype=int)
+            pos[0:len(pos):2] = 2
+            pos[2:len(pos):4] += 1
+            pos[1:len(pos):2] = np.floor(3.3*nr_xpos/4.0)
+            # and +1 alternating on the right
+            pos[1:len(pos):4] += 1
+            pos[3:len(pos):4] -= 1
+        elif xpos == 'right':
+            pos = np.zeros((nr_series,), dtype=int)
+            pos[0:len(pos):2] = 2
+            pos[1:len(pos):2] = np.ceil(1.0*nr_xpos/4.0)
+        elif xpos == 'left':
+            pos = np.zeros((nr_series,), dtype=int)
+            pos[0:len(pos):2] = np.ceil(2.0*nr_xpos/4.0)
+            pos[1:len(pos):2] = np.ceil(3.0*nr_xpos/4.0)
+        elif xpos == 'spread':
+            # TODO: we should determine the spacing based on wind speeds since
+            # it might not be equally spaced
+            pos = np.zeros((nr_series,), dtype=int)
+            pos[0:len(pos):4] = 0
+            pos[1:len(pos):4] = np.ceil(1.0*nr_xpos/4.0)
+            pos[2:len(pos):4] = np.ceil(2.0*nr_xpos/4.0)
+            pos[3:len(pos):4] = nr_xpos - 2
+
+        return pos
+
+    def plot_freq(self, ax, xpos='random', col='k', mark='^', ls='-',
+                  modes='all'):
+        """
+
+
+        Parameters
+        ----------
+        ax : TYPE
+            An Axes matplotlib instance that will be used for plotting.
+        xpos : TYPE, str
+            Specify strategy how to place the labels indifying the different modes.
+            Valid string values are: 'random', 'centre', 'borders', 'right', 'left', 'spread'.
+            The default is 'random'. If a list is passed, it should contain at
+            least nr_modes elements. Items in the list refer to the index of the
+            operating point (wind speed).
+        col : string, optional
+            Color used for the line plot. The default is 'k'. Can also be a list
+            of colors in case modes should have different colors.
+        mark : string, optional
+            Marker to use for the line plot. The default is '^'. Can also be a list
+            of mark symbols in case modes should have different symbols.
+        ls : TYPE, optional
+            Line style to use. The default is '-'. Can also be a list
+            of line styles in case modes should have different line styles.
+        modes : string, optional
+            Identify which modes to plot. The default is 'all'. Optionally
+            specify an integer to select the first x modes, or a list of indices
+            to cherry pick modes of interest.
+
+        Returns
+        -------
+        ax : matplotlib.axes
+            An updated Axes matplotlib instance of the plot.
+
+        """
+
+        if isinstance(modes, str) and modes == 'all':
+            df_freq = self.df_freq
+        elif isinstance(modes, int):
+            df_freq = self.df_freq[:modes]
+        else:
+            df_freq = self.df_freq[modes]
+
+        nr_winds = df_freq.shape[0]
+        nr_modes = df_freq.shape[1]
+
+        pos = xpos
+        if isinstance(xpos, str):
+            pos = self._inplot_label_pos(nr_winds, nr_modes, xpos)
+        elif isinstance(xpos, list) and len(xpos) < nr_modes:
+            raise ValueError(f'xpos has only {len(xpos)}, while it needs to be >= {nr_modes}')
+
+        if isinstance(col, str):
+            col = [col]
+        if isinstance(mark, str):
+            mark = [mark]
+        if isinstance(ls, str):
+            ls = [ls]
+        # just to make sure we always have enough colors/marks,lss
+        col = col*nr_modes
+        mark = mark*nr_modes
+        ls = ls*nr_modes
+
+        for i, (name, row) in enumerate(df_freq.items()):
+            colmark = '%s%s%s' % (col[i], ls[i], mark[i])
+            ax.plot(self.wind, row.values, colmark)#, mfc='w')
+            x, y = self.wind[pos[i]], row.values[pos[i]]
+            bbox = dict(boxstyle="round", alpha=self.alpha_box,
+                        edgecolor=col[i], facecolor=col[i])
+            ax.annotate(name, xy=(x, y), xycoords='data',
+                        xytext=(-6, +20), textcoords='offset points',
+                        fontsize=12, bbox=bbox, arrowprops=dict(arrowstyle="->",
+                        connectionstyle="arc3,rad=.05"), color='w')
+
+        self.nr_modes = nr_modes
+
+        return ax
+
+    def plot_damp(self, ax, xpos='random', col='r', mark='o' ,ls='--',
+                  modes=14):
+        """see plot_freq
+
+        """
+
+        # reduce the number of modes we are going to plot
+        if isinstance(modes, str) and modes == 'all':
+            df_damp = self.df_damp
+        elif isinstance(modes, int):
+            nr_modes = modes
+            # sort the columns according to damping: lowest damped modes first
+            # sort according to damping at lowest wind speed
+            isort = self.df_damp.iloc[0,:].values.argsort()
+            modes_sort_reduced = self.df_damp.columns[isort][:nr_modes]
+            df_damp = self.df_damp[modes_sort_reduced]
+        else:
+            df_damp = self.df_damp[modes]
+
+        nr_winds = df_damp.shape[0]
+        nr_modes = df_damp.shape[1]
+
+        if isinstance(col, str):
+            col = [col]
+        if isinstance(mark, str):
+            mark = [mark]
+        if isinstance(ls, str):
+            ls = [ls]
+        # just to make sure we always have enough colors/marks,lss
+        col = col*nr_modes
+        mark = mark*nr_modes
+        ls = ls*nr_modes
+
+        # put the labels in sensible places
+        pos = xpos
+        if isinstance(xpos, str):
+            pos = self._inplot_label_pos(nr_winds, nr_modes, xpos)
+        elif isinstance(xpos, list) and len(xpos) < nr_modes:
+            raise ValueError(f'xpos has only {len(xpos)}, while it needs to be >= {nr_modes}')
+
+        for i, (name, row) in enumerate(df_damp.items()):
+            colmark = '%s%s%s' % (col[i], ls[i], mark[i])
+            ax.plot(self.wind, row.values, colmark, alpha=0.8)#, mfc='w')
+            x, y = self.wind[pos[i]], row.values[pos[i]]
+            bbox = dict(boxstyle="round", alpha=self.alpha_box,
+                        edgecolor=col[i], facecolor=col[i])
+            ax.annotate(name, xy=(x, y), xycoords='data',
+                        xytext=(-6, 20), textcoords='offset points',
+                        fontsize=12, bbox=bbox, arrowprops=dict(arrowstyle="->",
+                        connectionstyle="arc3,rad=.05"), color='w')
+
+        self.nr_modes = nr_modes
+
+        return ax
+
+    def add_legend(self, ax, labels, on='freq'):
+        if on == 'freq':
+            i0s = self.istart_freqs
+        else:
+            i0s = self.istart_damps
+
+        return ax.legend([ax.lines[k] for k in i0s[:-1]], labels, loc='best')
+
+
+def plot_add_ps(ax, wind, rpm, col='g', fmax=10, ps=None):
+    # plot 1P, 2P, ..., 9P
+    bbox = dict(boxstyle="round", alpha=0.4, edgecolor=col, facecolor=col,)
+    if ps is None:
+        pmax = int(60*fmax/rpm.mean())
+        ps = list(range(1, pmax))
+    for p in ps:
+        if p%3 == 0:
+            alpha=0.6
+            ax.plot(wind, rpm*p/60, '%s--' % col, )
+        else:
+            alpha = 0.4
+            ax.plot(wind, rpm*p/60, '%s-.' % col)#, alpha=alpha)
+        x, y = wind[10], rpm[10]*p/60
+        p_str = '%iP' % p
+        bbox['alpha'] = alpha
+        ax.text(x, y, p_str, fontsize=9, verticalalignment='bottom',
+                horizontalalignment='center', bbox=bbox, color='w')
+
+    return ax
+
+
 if __name__ == '__main__':
 
     pass

wetb/prepost/misc.py

@@ -6,22 +6,9 @@ Library for general stuff
 
 @author: dave
 """
-from __future__ import print_function
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import range
-from builtins import dict
-from builtins import int
-from io import open
-from builtins import str
-from future import standard_library
-standard_library.install_aliases()
-from builtins import object
 
 import os
-import sys
-import shutil
+import io
 import unittest
 import pickle
 import re
@@ -77,13 +64,69 @@ def path_split_dirs(path):
     return dirs
 
 
-def print_both(f, text, end='\n'):
+def path_sanitize(path, allowdd=False, allowabs=False):
+    """Raises a ValueError if not considered safe. A leading ../ is allowed
+    when allowdd=True.
     """
-    Print both to a file and the console
+
+    if not isinstance(path, str):
+        raise ValueError('path_sanitize requires a string as input.')
+
+    keepcharacters = set(['/', '.', '_', '-'])
+
+    # in special cases we allow one leading ../
+    if allowdd and path[:3] == '../':
+        path = path[3:]
+
+    # no absolute paths if not allowabs
+    if path[0] == '/' and not allowabs:
+        raise ValueError('Absolute paths not allowed: "%s"' % path)
+
+    # no empty paths either
+    if path == '':
+        raise ValueError('Empty paths are not allowed')
+
+    # only alphanummerical characters (includes unicode)
+    if not all([c.isalnum() or c in keepcharacters for c in path]):
+        raise ValueError('Invalid or unsafe path: "%s"' % path)
+
+    # additional checks on sub-directories
+    items = path.split('/')
+    # leading/trailing '/' leads to first/last element being emtpy
+    if path[0] == "/":
+        items = items[1:]
+    if path[-1] == '/':
+        items = items[:-1]
+
+    no_lt = set(['.', '-'])  # characters not allowed leading/trailing a sub-dir
+    for item in items:
+        # require a sub-dir to be at least 2 characters or more
+        if len(item) < 2:
+            msg = 'Directories and filenames need to be longer than 1 character.'
+            raise ValueError('Invalid or unsafe path: "%s". %s' % (path, msg))
+        if item[0] in no_lt or item[-1] in no_lt:
+            msg = 'No leading/trailing . or - allowed.'
+            raise ValueError('Invalid or unsafe path: "%s". %s' % (path, msg))
+
+
+def sanitize_wine_prefix(wine_prefix):
+    """special case to sanitize
     """
-    print(text)
-    if isinstance(f, file):
-        f.write(text + end)
+
+    # sanitize wine_prefix
+    # for backwards compatibility
+    if wine_prefix[:2] == '~/':
+        wine_prefix = wine_prefix[2:]
+    # first character can also be a dot
+    if not wine_prefix[0].isalnum():
+        if wine_prefix[0]=='.':
+            pass
+        else:
+            raise ValueError('Invalid wine_prefix value: %s' % wine_prefix)
+    if not wine_prefix[1:].isalnum():
+        raise ValueError('Invalid wine_prefix value: %s' % wine_prefix)
+    return os.path.join('$HOME', wine_prefix).replace("\\", "/")
+
 
 def unique(s):
     """
@@ -157,6 +200,7 @@ def unique(s):
             u.append(x)
     return u
 
+
 def CoeffDeter(obs, model):
     """
     Coefficient of determination
@@ -208,6 +252,7 @@ def calc_sample_rate(time, rel_error=1e-4):
 #        raise AssertionError
     return 1/deltas.mean()
 
+
 def findIntersection(fun1, fun2, x0):
     """
     Find Intersection points of two functions
@@ -237,6 +282,7 @@ def findIntersection(fun1, fun2, x0):
     """
     return sp.optimize.fsolve(lambda x : fun1(x) - fun2(x), x0)
 
+
 # TODO: replace this with some of the pyrain functions
 def find0(array, xi=0, yi=1, verbose=False, zerovalue=0.0):
     """
@@ -349,6 +395,7 @@ def find0(array, xi=0, yi=1, verbose=False, zerovalue=0.0):
 
     return y0, y0i
 
+
 def remove_items(list, value):
     """Remove items from list
     The given list wil be returned withouth the items equal to value.
@@ -368,6 +415,7 @@ def remove_items(list, value):
 
     return list
 
+
 class DictDB(object):
     """
     A dictionary based database class
@@ -475,6 +523,7 @@ class DictDB(object):
             if alltrue:
                 self.dict_sel[row] = self.dict_db[row]
 
+
 class DictDiff(object):
     """
     Calculate the difference between two dictionaries as:
@@ -517,6 +566,7 @@ class DictDiff(object):
         t=set(o for o in self.intersect if self.past_d[o] == self.current_d[o])
         return t
 
+
 def fit_exp(time, data, checkplot=True, method='linear', func=None, C0=0.0):
     """
     Note that all values in data have to be possitive for this method to work!
@@ -563,6 +613,7 @@ def fit_exp(time, data, checkplot=True, method='linear', func=None, C0=0.0):
 
     return fit, A, K, C
 
+
 def curve_fit_exp(time, data, checkplot=True, weights=None):
     """
     This code is based on a StackOverflow question/answer:
@@ -621,70 +672,47 @@ def curve_fit_exp(time, data, checkplot=True, weights=None):
 
     return
 
-def convert_to_utf8(filename):
-    # gather the encodings you think that the file may be
-    # encoded inside a tuple
-    encodings = ('windows-1253', 'iso-8859-7', 'macgreek')
 
-    # try to open the file and exit if some IOError occurs
-    try:
-        f = open(filename, 'r').read()
-    except Exception:
-        sys.exit(1)
+def readlines_try_encodings(fname):
+    """Read text file in binary and try to encode with a few common encodings.
+    Return the first succesful attempt.
+
+    Parameters
+    ----------
+    fname : str
+        Path to file to be read.
+
+    Returns
+    -------
+    lines : list of str
+        Returns same as open(fname).readlines().
+    """
+
+    # some potentially common encoding formats
+    encodings = ('utf-8', 'cp1252', 'iso-8859-1', 'latin-1', 'cp1250',
+                 'cp1251', 'cp1253', 'iso-8859-2', 'iso-8859-7', 'macgreek')
+
+    # open the file and read entire object in as bytes, we do this to avoid
+    # re-opening and reading the file for each encoding we try
+    with open(fname, 'rb') as fid:
+        fbytes = fid.read()
 
-    # now start iterating in our encodings tuple and try to
-    # decode the file
     for enc in encodings:
+        # convert to a file object
+        fid = io.BytesIO(fbytes)
         try:
-            # try to decode the file with the first encoding
-            # from the tuple.
-            # if it succeeds then it will reach break, so we
-            # will be out of the loop (something we want on
-            # success).
-            # the data variable will hold our decoded text
-            data = f.decode(enc)
+            # we want to mimic the standard readlines behaviour
+            text_wrapper = io.TextIOWrapper(fid, encoding=enc)
+            lines = text_wrapper.readlines()
+        except UnicodeDecodeError:
+            print('encoding with', enc, 'does not work for:', fname)
+            print('unicode error with %s , trying different encoding' % enc)
+            pass
+        else:
+            # print('opening the file with encoding:  %s ' % enc)
             break
-        except Exception:
-            # if the first encoding fail, then with the continue
-            # keyword will start again with the second encoding
-            # from the tuple an so on.... until it succeeds.
-            # if for some reason it reaches the last encoding of
-            # our tuple without success, then exit the program.
-            if enc == encodings[-1]:
-                sys.exit(1)
-            continue
-
-    # now get the absolute path of our filename and append .bak
-    # to the end of it (for our backup file)
-    fpath = os.path.abspath(filename)
-    newfilename = fpath + '.bak'
-    # and make our backup file with shutil
-    shutil.copy(filename, newfilename)
-
-    # and at last convert it to utf-8
-    f = open(filename, 'w')
-    try:
-        f.write(data.encode('utf-8'))
-    except Exception(e):
-        print(e)
-    finally:
-        f.close()
+    return lines
 
-def to_lower_case(proot):
-    """
-    Rename all the files in the subfolders of proot to lower case, and
-    also the subfolder name when it the folder name starts with DLC
-    """
-    # find all dlc defintions in the subfolders
-    for root, dirs, files in os.walk(proot):
-        for fname in files:
-            orig = os.path.join(root, fname)
-            rename = os.path.join(root, fname.lower())
-            os.rename(orig, rename)
-        base = root.split(os.path.sep)[-1]
-        if base[:3] == 'DLC':
-            new = root.replace(base, base.lower())
-            os.rename(root, new)
 
 def read_excel_files(proot, fext='xlsx', pignore=None, sheet=0,
                      pinclude=None, silent=False):
@@ -741,7 +769,7 @@ def read_excel_files(proot, fext='xlsx', pignore=None, sheet=0,
                 if fext == 'csv':
                     df = pd.read_csv(f_target)
                 else:
-                    df = pd.read_excel(f_target, sheetname=sheet)
+                    df = pd.read_excel(f_target, sheet_name=sheet)
                 df_list[f_target.replace('.'+fext, '')] = df
                 if not silent:
                     print(': sucesfully included %i case(s)' % len(df))
@@ -751,21 +779,6 @@ def read_excel_files(proot, fext='xlsx', pignore=None, sheet=0,
 
     return df_list
 
-def convert_xlsx2csv(fpath, sheet='Sheet1', fext='xlsx'):
-    """
-    Convert xlsx load case definitions to csv so we can track them with git
-    """
-
-    for root, dirs, files in os.walk(fpath):
-        for file_name in files:
-            if not file_name.split('.')[-1] == fext:
-                continue
-            fxlsx = os.path.join(root, file_name)
-            print(fxlsx)
-            xl = pd.ExcelFile(fxlsx)
-            df = xl.parse(sheet)
-            fcsv = fxlsx.replace(fext, 'csv')
-            df.to_csv(fcsv, sep=';')
 
 def check_df_dict(df_dict):
     """
@@ -1113,10 +1126,12 @@ def df_dict_check_datatypes(df_dict):
     # we can not pop/delete items from a dict while iterating over it
     df_dict2 = {}
     for colkey, col in df_dict.items():
+        if len(col)==0:
+            pass
         # if we have a list, convert to string
-        if type(col[0]).__name__ == 'list':
+        elif type(col[0]).__name__ == 'list':
             for ii, item in enumerate(col):
-                col[ii] = '**'.join(item)
+                col[ii] = '*;*'.join(item)
         # if we already have an array (statistics) or a list of numbers
         # do not try to cast into another data type, because downcasting
         # in that case will not raise any exception
@@ -1126,21 +1141,21 @@ def df_dict_check_datatypes(df_dict):
         # in case we have unicodes instead of strings, we need to convert
         # to strings otherwise the saved .h5 file will have pickled elements
         try:
-            df_dict2[str(colkey)] = np.array(col, dtype=np.int32)
+            df_dict2[str(colkey)] = np.array(col, dtype=int)
         except OverflowError:
             try:
-                df_dict2[str(colkey)] = np.array(col, dtype=np.int64)
+                df_dict2[str(colkey)] = np.array(col, dtype=int)
             except OverflowError:
-                df_dict2[str(colkey)] = np.array(col, dtype=np.float64)
+                df_dict2[str(colkey)] = np.array(col, dtype=float)
         except ValueError:
             try:
-                df_dict2[str(colkey)] = np.array(col, dtype=np.float64)
+                df_dict2[str(colkey)] = np.array(col, dtype=float)
             except ValueError:
-                df_dict2[str(colkey)] = np.array(col, dtype=np.str)
+                df_dict2[str(colkey)] = np.array(col, dtype=str)
         except TypeError:
             # in all other cases, make sure we have converted them to
             # strings and NOT unicode
-            df_dict2[str(colkey)] = np.array(col, dtype=np.str)
+            df_dict2[str(colkey)] = np.array(col, dtype=str)
         except Exception as e:
             print('failed to convert column %s to single data type' % colkey)
             raise(e)
@@ -1238,6 +1253,84 @@ class Tests(unittest.TestCase):
     def setUp(self):
         pass
 
+    def test_sanitize_wine_prefix(self):
+
+        wine_prefix = sanitize_wine_prefix('.wine32')
+        self.assertEqual(wine_prefix, '$HOME/.wine32')
+
+        wine_prefix = sanitize_wine_prefix('~/.wine32')
+        self.assertEqual(wine_prefix, '$HOME/.wine32')
+
+        wine_prefix = sanitize_wine_prefix('wine32')
+        self.assertEqual(wine_prefix, '$HOME/wine32')
+
+        notok = ['../wine', '/root/wine32', 'mnt/mimer/wine32', '.wine32/qq']
+        for wine_prefix in notok:
+            with self.assertRaises(ValueError):
+                sanitize_wine_prefix(wine_prefix)
+
+    def test_path_sanitize(self):
+
+        paths = ['./not/allowed/',
+                 '$not/$allowed/',
+                 ' forget/about/it',
+                 'forget/about/it ',
+                 '!are\you\n\\NUTS??',
+                 'don"t/dothis',
+                 'don\'t',
+                 'oh_look/a/ space'
+                 '/not/../at/all',
+                 '/not/at/all../',
+                 'you/think\\youneed\\backspaces',
+                 '../../for/db/dirs/this/is/notsafe',
+                 'auch/aa/wild/c*rd',
+                 'no/weird/-sub/dirs',
+                 'no/weird/sub-/dirs',
+                 'no/weird/sub./dirs',
+                 'no/weird/-sub-/dirs',
+                 '..nope/either8',
+                 '&',
+                 ]
+        for path in paths:
+            with self.assertRaises(ValueError):
+                path_sanitize(path)
+            with self.assertRaises(ValueError):
+                path_sanitize('/' + path)
+            with self.assertRaises(ValueError):
+                path_sanitize(path, allowdd=True)
+            with self.assertRaises(ValueError):
+                path_sanitize(path, allowabs=True)
+
+        # with allowdd active
+        pathsok = ['../for/db/dirs/this/is/safe',
+                   'what/aa/nice/path/',
+                   'whatanicepath',
+                   'ok/p-_/u_/in/this/case',
+                   'ev.en.this.__is/oo/kk---000/',
+                   ]
+        for path in pathsok:
+            path_sanitize(path, allowdd=True)
+            with self.assertRaises(ValueError):
+                path_sanitize('/' + path, allowdd=True)
+
+        # with allowabs active
+        pathsok = ['/for/db/dirs/this/is/safe',
+                   '/what/aa/nice/path/',
+                   '/whatanicepath',
+                   '/ok/p-_/u_/in/this/case',
+                   '/ev.en.this.__is/oo/kk---000/',
+                   ]
+        for path in pathsok:
+            path_sanitize(path, allowabs=True)
+            with self.assertRaises(ValueError):
+                path_sanitize('..' + path, allowabs=True)
+
+        with self.assertRaises(ValueError):
+            path_sanitize('../../for/db/dirs/this/is/notsafe', allowdd=True)
+
+        with self.assertRaises(ValueError):
+            path_sanitize('./for/db/dirs/this/is/notsafe', allowabs=True)
+
     def test_rebin1(self):
         hist = np.array([2,5,5,9,2,6])
         bins = np.arange(7)

wetb/prepost/mplutils.py

@@ -4,31 +4,15 @@ Created on Wed Nov 23 11:22:50 2011
 
 @author: dave
 """
-from __future__ import division
-from __future__ import print_function
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import range
-from builtins import int
-from builtins import dict
-from builtins import round
-from future import standard_library
-standard_library.install_aliases()
-
-
-
-
 # external libraries
 import numpy as np
+from scipy.signal import find_peaks
 
 import matplotlib as mpl
+from matplotlib.figure import Figure
 # use a headless backend
-from matplotlib.backends.backend_agg import FigureCanvasAgg as FigCanvas
-# wafo is an optional dependency only required for non default PSD peak marking
-try:
-    import wafo
-except ImportError:
-    pass
+# from matplotlib.backends.backend_agg import FigureCanvasAgg as FigCanvas
+
 
 
 def make_fig(nrows=1, ncols=1, figsize=(12,8), dpi=120):
@@ -50,7 +34,8 @@ def make_fig(nrows=1, ncols=1, figsize=(12,8), dpi=120):
     return subplots(nrows=nrows, ncols=ncols, figsize=figsize, dpi=dpi)
 
 
-def subplots(nrows=1, ncols=1, figsize=(12,8), dpi=120, num=0, subplot_kw={}):
+def subplots(nrows=1, ncols=1, figsize=(12,8), dpi=120, num=0, subplot_kw={},
+             gridspec_kw={}):
     """
 
     Equivalent function of pyplot.subplots(). The difference is that this one
@@ -69,15 +54,21 @@ def subplots(nrows=1, ncols=1, figsize=(12,8), dpi=120, num=0, subplot_kw={}):
 
     """
 
-    fig = mpl.figure.Figure(figsize=figsize, dpi=dpi)
-    canvas = FigCanvas(fig)
-    fig.set_canvas(canvas)
-    axes = np.ndarray((nrows, ncols), dtype=np.object)
-    plt_nr = 1
-    for row in range(nrows):
-        for col in range(ncols):
-            axes[row,col] = fig.add_subplot(nrows, ncols, plt_nr, **subplot_kw)
-            plt_nr += 1
+    fig = Figure(figsize=figsize, dpi=dpi)
+    axes = fig.subplots(nrows=nrows, ncols=ncols,
+                        subplot_kw=subplot_kw, gridspec_kw=gridspec_kw)
+    try:
+        axes = axes.reshape((nrows,ncols))
+    except AttributeError:
+        pass
+    # canvas = FigCanvas(fig)
+    # fig.set_canvas(canvas)
+    # axes = np.ndarray((nrows, ncols), dtype=object)
+    # plt_nr = 1
+    # for row in range(nrows):
+    #     for col in range(ncols):
+    #         axes[row,col] = fig.add_subplot(nrows, ncols, plt_nr, **subplot_kw)
+    #         plt_nr += 1
     return fig, axes
 
 
@@ -142,12 +133,12 @@ def one_legend(*args, **kwargs):
     return leg
 
 
-def p4psd(ax, rpm_mean, p_max=17, y_pos_rel=0.25, color='g', ls='--',
+def p4psd(ax, rpm_mean, fmax=10, y_pos_rel=0.25, color='g', ls='--', ps=None,
           col_text='w'):
     """
     Add the P's on a PSD
 
-    fn_max is the maximum value on the plot (ax.xlim). This only works when
+    fmax is the maximum value on the plot (ax.xlim). This only works when
     setting the xlim of the plot before calling p4psd.
 
     Parameters
@@ -157,22 +148,31 @@ def p4psd(ax, rpm_mean, p_max=17, y_pos_rel=0.25, color='g', ls='--',
 
     rpm_mean
 
-    p_max : int, default=17
+    fmax : int, default=17
+        stop plotting p's after fmax (when ps is None)
+
+    ps : iterable of ints, default=None
+        specify which p's to plot (ignores fmax)
 
     y_pos_rel : int or list, default=0.25
     """
+
+    if ps is None:
+        pmax = int(60*fmax/rpm_mean)
+        ps = list(range(1, pmax))
+    else:
+        pmax = len(ps)
+
     if isinstance(y_pos_rel, float) or isinstance(y_pos_rel, int):
-        y_pos_rel = [y_pos_rel]*p_max
+        y_pos_rel = [y_pos_rel]*pmax
 
     f_min = ax.get_xlim()[0]
     f_max = ax.get_xlim()[1]
 
     # add the P's
     bbox = dict(boxstyle="round", edgecolor=color, facecolor=color)
-    for i, p in enumerate(range(1, p_max)):
+    for i, p in enumerate(ps):
         p_freq = p * rpm_mean / 60.0
-        if p_freq > f_max:
-            break
         if p%3 == 0:
             alpha=0.5
             ax.axvline(x=p_freq, linewidth=1, color=color, alpha=0.7, ls=ls)
@@ -192,11 +192,25 @@ def p4psd(ax, rpm_mean, p_max=17, y_pos_rel=0.25, color='g', ls='--',
     return ax
 
 
-def peaks(ax, freqs, Pxx, fn_max, min_h, nr_peaks=15, col_line='k',
+def peaks(ax, freqs, Pxx, fn_max, min_h, nr_peaks=15, min_p=0, col_line='k',
           ypos_mean=0.14, col_text='w', ypos_delta=0.06, bbox_alpha=0.5,
-          verbose=False):
+          verbose=False, format_text='%2.2f', period=False):
     """
     indicate the peaks
+
+    Parameters
+    ----------
+
+    period : boolean, default=False
+        If True, the period instead of the frequency is given in the text
+
+    min_h : float
+        The threshold in the rainflowfilter (default 0.05*range(S(:))).
+        A zero value will return all the peaks of S.
+
+    min_p : float, 0..1
+        Only the peaks that are higher than min_p*max(max(S))
+        min_p*(the largest peak in S) are returned (default  0).
     """
     i_fn_max = np.abs(freqs - fn_max).argmin()
     # ignore everything above fn_max
@@ -204,20 +218,21 @@ def peaks(ax, freqs, Pxx, fn_max, min_h, nr_peaks=15, col_line='k',
     Pxx = Pxx[:i_fn_max]
     Pxx_log = 10.*np.log10(Pxx)
     try:
-        pi = wafo.misc.findpeaks(Pxx_log, n=len(Pxx), min_h=min_h)
+        # pi = wafo.misc.findpeaks(Pxx_log, n=len(Pxx), min_h=min_h, min_p=min_p)
+        pi, pi_props = find_peaks(Pxx_log, height=Pxx_log.min(), prominence=(20, None))
         if verbose:
             print('len Pxx', len(Pxx_log), 'nr of peaks:', len(pi))
     except Exception as e:
         if verbose:
             print('len Pxx', len(Pxx_log))
-            print('*** wafo.misc.findpeaks FAILED ***')
+            # print('*** wafo.misc.findpeaks FAILED ***')
             print(e)
         return ax
 
-    # only take the nr_peaks most significant heights
-    pi = pi[:nr_peaks]
-    # and sort them accoriding to frequency (or index)
-    pi.sort()
+    # # only take the nr_peaks most significant heights
+    # pi = pi[:nr_peaks]
+    # # and sort them accoriding to frequency (or index)
+    # pi.sort()
 
     # mark the peaks with a circle
 #    ax.plot(freqs[pi], Pxx[:xlim][pi], 'o')
@@ -229,8 +244,6 @@ def peaks(ax, freqs, Pxx, fn_max, min_h, nr_peaks=15, col_line='k',
 #        Pxx_peak = Pxx_log[ii]
         # take the average frequency and plot vertical line
         ax.axvline(x=freq_peak, linewidth=1, color=col_line, alpha=0.6)
-        # and the value in a text box
-        textbox = '%2.2f' % freq_peak
         if switch:
             # locate at the min value (down the plot), but a little
             # lower so it does not interfere with the plot itself
@@ -249,6 +262,10 @@ def peaks(ax, freqs, Pxx, fn_max, min_h, nr_peaks=15, col_line='k',
 #        print textbox
 #        print yrange_plot
         xrel = freq_peak/fn_max
+        # and the value in a text box
+        if period:
+            freq_peak = 1/freq_peak
+        textbox = format_text % freq_peak
         ax.text(xrel, text_ypos, textbox, fontsize=10, transform=ax.transAxes,
                  va='bottom', color=col_text, bbox=dict(boxstyle="round",
                  ec=col_line, fc=col_line, alpha=bbox_alpha,))
@@ -282,8 +299,8 @@ def match_yticks(ax1, ax2, nr_ticks_forced=None, extend=False):
 
 
 def psd(ax, time, sig, nfft=None, res_param=250, f0=0, f1=None, nr_peaks=10,
-        min_h=15, mark_peaks=False, col='r-', label=None, alpha=1.0,
-        ypos_peaks=0.9, ypos_peaks_delta=0.12):
+        min_h=15, min_p=0, mark_peaks=False, col='r-', label=None, alpha=1.0,
+        ypos_peaks=0.9, ypos_peaks_delta=0.12, format_text='%2.2f', period=False):
     """Only plot the psd on a given axis and optionally mark the peaks.
     """
 
@@ -297,7 +314,11 @@ def psd(ax, time, sig, nfft=None, res_param=250, f0=0, f1=None, nr_peaks=10,
         nfft = len(sig)
 
     # calculate the PSD
-    Pxx, freqs = mpl.mlab.psd(sig, NFFT=nfft, Fs=sps)
+
+    Pxx, freqs = mpl.mlab.psd(sig, NFFT=nfft, Fs=sps) # window=None
+    # for clean signals (steps, sinus) use another window
+    # for stochastic signal, the default window Hanning is good
+    # in Scipy you can provide strings for the different windows, not so in mlab
 
     i0 = np.abs(freqs - f0).argmin()
     i1 = np.abs(freqs - f1).argmin()
@@ -305,9 +326,9 @@ def psd(ax, time, sig, nfft=None, res_param=250, f0=0, f1=None, nr_peaks=10,
     # plotting psd, marking peaks
     ax.plot(freqs[i0:i1], Pxx[i0:i1], col, label=label, alpha=alpha)
     if mark_peaks:
-        ax = peaks(ax, freqs[i0:i1], Pxx[i0:i1], fn_max=f1,
-                   nr_peaks=nr_peaks, col_line=col[:1],
-                   ypos_delta=ypos_peaks_delta, bbox_alpha=0.5,
+        ax = peaks(ax, freqs[i0:i1], Pxx[i0:i1], fn_max=f1, min_p=min_p,
+                   nr_peaks=nr_peaks, col_line=col[:1], format_text=format_text,
+                   ypos_delta=ypos_peaks_delta, bbox_alpha=0.5, period=period,
                    ypos_mean=ypos_peaks, min_h=min_h, col_text='w')
 
     return ax
@@ -379,6 +400,40 @@ def time_psd(results, labels, axes, alphas=[1.0, 0.7], colors=['k-', 'r-'],
     return axes
 
 
+def get_list_colors(nr, cmap='magma'):
+    """Returns a list of color tuples
+
+    Paramters
+    ---------
+
+    nr : int
+        number of colors
+
+    cmap : str
+        a matplotlib color map name, for example: viridis, plasma, magma, hot,
+        cool, binary, see also https://matplotlib.org/users/colormaps.html
+
+    Returns
+    -------
+
+    clist : list
+        List continaing 'nr' color tuples (with 3 elements).
+
+    """
+    # define the number of positions you want to have the color for
+    # select a color map
+    cmap = mpl.cm.get_cmap(cmap, nr)
+    # other color maps: 	cubehelix, afmhot, hot
+    # convert to array
+    cmap_arr = cmap(np.arange(nr))
+    # and now you have each color as an RGB tuple as
+    clist = []
+    for i in cmap_arr:
+        clist.append(tuple(i[0:3]))
+
+    return clist
+
+
 if __name__ == '__main__':
 
     pass

wetb/prepost/parpost.py

+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 24 13:41:56 2018
+
+@author: dave
+"""
+
+import time
+import os
+from multiprocessing import Pool
+from glob import glob
+from argparse import ArgumentParser
+import json
+
+import numpy as np
+import pandas as pd
+# import scipy.io as sio
+import tables as tbl
+
+from tqdm import tqdm
+
+from wetb.utils.envelope import compute_envelope
+from wetb.prepost import windIO
+from wetb.prepost.Simulations import EnvelopeClass
+from wetb.prepost.simchunks import AppendDataFrames
+
+
+def logcheck(fname):
+    """Check the log file of a single HAWC2 simulation and save results to
+    textfile.
+    """
+    fsave = None
+    mode = 'w'
+
+    logf = windIO.LogFile()
+    logf.readlog(fname)
+    contents = logf._msglistlog2csv('')
+    if fsave is None:
+        fsave = fname.replace('.log', '.csv')
+    with open(fsave, mode) as f:
+        f.write(contents)
+
+
+def add_channels(res):
+    """Add channels in memory so they can be included in the statistics calculations
+    """
+
+    # EXAMPLE: tower bottom resultant moment
+    # chitx = res.ch_dict['tower-tower-node-001-momentvec-x']['chi']
+    # chity = res.ch_dict['tower-tower-node-001-momentvec-y']['chi']
+    # mx = res.sig[:,chitx]
+    # my = res.sig[:,chity]
+    # data = np.sqrt(mx*mx + my*my).reshape(len(res.sig),1)
+    # # add the channel
+    # res.add_channel(data, 'TB_res', 'kNm', 'Tower bottom resultant bending moment')
+
+    # blade resultant loads
+    chitx = res.ch_dict['blade1-blade1-node-004-forcevec-x']['chi']
+    chity = res.ch_dict['blade1-blade1-node-004-forcevec-y']['chi']
+    x = res.sig[:,chitx]
+    y = res.sig[:,chity]
+    data = np.sqrt(x*x + y*y).reshape(len(res.sig),1)
+    res.add_channel(data, 'BR_F_res', 'kN', 'Blade root resultant shear force')
+
+    chitx = res.ch_dict['blade1-blade1-node-004-momentvec-x']['chi']
+    chity = res.ch_dict['blade1-blade1-node-004-momentvec-y']['chi']
+    x = res.sig[:,chitx]
+    y = res.sig[:,chity]
+    data = np.sqrt(x*x + y*y).reshape(len(res.sig),1)
+    res.add_channel(data, 'BR_M_res', 'kNm', 'Blade root resultant bending moment')
+
+    return res
+
+
+def loads_at_extreme(dfres, chans_selmax, chans_atmax):
+    """For a range of channels defined in chans_list, list their value respective
+    values at the time when a max/min occurs.
+
+    It makes most sense to use chans_selmax == chans_atmax (BLADED style reports)
+
+    Other use is to track loads at for example the maximum yaw angle.
+    """
+    statparams = ['max', 'min']
+
+    # create multi-index: first is the selected channel, second is the stat param
+    chans_selmax_ = []
+    for k in chans_selmax: chans_selmax_.extend([k ,k])
+    index = pd.MultiIndex.from_arrays([chans_selmax_, statparams*len(chans_selmax)])
+    # initiale the table
+    dfextr = pd.DataFrame(np.zeros((len(chans_selmax_), (len(chans_atmax)))),
+                          index=index, columns=chans_atmax)
+
+    # vectorised slightly faster for 8x8 channels
+    # # 9.08 ms ± 72.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
+    # select the index at which the max/minimum occurs for the channels in chans_selmax
+    idxmax = dfres[chans_selmax].idxmax(axis=0).values
+    idxmin = dfres[chans_selmax].idxmin(axis=0).values
+    dfextr.loc[(chans_selmax, 'max'), chans_atmax] = dfres.loc[idxmax, chans_atmax].values
+    dfextr.loc[(chans_selmax, 'min'), chans_atmax] = dfres.loc[idxmin, chans_atmax].values
+
+    # # looping is slower for 8x8 channels
+    # # 11.1 ms ± 21.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
+    # for chan in chans_selmax:
+    #     idxmax = dfres[chan].idxmax(axis=0)
+    #     idxmin = dfres[chan].idxmin(axis=0)
+    #     for atchan in chans_atmax:
+    #         dfextr.loc[(chan, 'max'),atchan] = dfres.loc[idxmax, atchan]
+    #         dfextr.loc[(chan, 'min'),atchan] = dfres.loc[idxmin, atchan]
+
+    # dfextr.loc[pd.IndexSlice[:,'max'],chan] = dfres.iloc[imax, ichan]
+    # dfextr.loc[pd.IndexSlice[:,'min'],chan] = dfres.iloc[imin, ichan]
+    # dfextr.loc[pd.IndexSlice[:,:],:].values
+    # assert np.allclose(dfextr.loc[:,:].values, dfextr2.loc[:,:].values)
+
+    return dfextr
+
+
+def envelope(res, int_env=False, Nx=300):
+
+    # define which sensors to use by using the unique_ch_name feature
+    ch_list = []
+
+    # find all available outputs for tower, blades of type momentvec Mx
+    sel = ((  (res.ch_df['bodyname']=='tower')
+            | (res.ch_df['bodyname']=='blade1')
+            | (res.ch_df['bodyname']=='blade2')
+            | (res.ch_df['bodyname']=='blade3') )
+           & (res.ch_df['sensortype']=='momentvec') & (res.ch_df['component']=='x')
+          )
+    nodeselx = res.ch_df[sel]['unique_ch_name'].values.tolist()
+
+    # all momentvec My
+    sel = ((  (res.ch_df['bodyname']=='tower')
+            | (res.ch_df['bodyname']=='blade1')
+            | (res.ch_df['bodyname']=='blade2')
+            | (res.ch_df['bodyname']=='blade3') )
+           & (res.ch_df['sensortype']=='momentvec') & (res.ch_df['component']=='y')
+          )
+    nodesely = set(res.ch_df[sel]['unique_ch_name'].values.tolist())
+
+    for nodex in nodeselx:
+        nodey = nodex.replace('-x', '-y')
+        # make sure both Mx and My are defined
+        if nodey not in nodesely:
+            continue
+        ch_list.append([nodex, nodey])
+
+    fname = res.FileName + '_envelopes.h5'
+    filt = tbl.Filters(complevel=9)
+    with tbl.open_file(fname, mode="w", title=str(SIM_ID), filters=filt) as h5f:
+        groupname = str(os.path.basename(res.FileName))
+        groupname = groupname.replace('-', '_').replace('.', '_')
+        ctab = h5f.create_group("/", groupname)
+
+        # envelope = {}
+        for ch_names in ch_list:
+
+            ichans = []
+            for ch_name in ch_names:
+                ichans.append(res.ch_dict[ch_name]['chi'])
+            cloud = res.sig[:, ichans]
+            # Compute a Convex Hull, the vertices number varies according to
+            # the shape of the poligon
+            vertices = compute_envelope(cloud, int_env=int_env, Nx=Nx)
+            # envelope[ch_names[0]] = vertices
+
+            # # save as simple text file
+            # title = ch_names[0]
+            # fname = res.FileName.replace('res/', 'prepost-data/env_txt/') + f'_{title}.txt'
+            # os.makedirs(os.path.dirname(fname), exist_ok=True)
+            # np.savetxt(fname, vertices)
+
+            # save as HDF5 table
+            title = str(ch_names[0].replace('-', '_'))
+            csv_table = h5f.create_table(ctab, title,
+                                          EnvelopeClass.section,
+                                          title=title)
+            tablerow = csv_table.row
+            for row in vertices: #envelope[ch_names[0]]:
+                tablerow['Mx'] = float(row[0])
+                tablerow['My'] = float(row[1])
+                if len(row)>2:
+                    tablerow['Mz'] = float(row[2])
+                    if len(row)>3:
+                        tablerow['Fx'] = float(row[3])
+                        tablerow['Fy'] = float(row[4])
+                        tablerow['Fz'] = float(row[5])
+                    else:
+                        tablerow['Fx'] = 0.0
+                        tablerow['Fy'] = 0.0
+                        tablerow['Fz'] = 0.0
+                else:
+                    tablerow['Mz'] = 0.0
+                    tablerow['Fx'] = 0.0
+                    tablerow['Fy'] = 0.0
+                    tablerow['Fz'] = 0.0
+                tablerow.append()
+            csv_table.flush()
+    # h5f.close()
+
+
+def calc(fpath, do_envelope, no_bins, m, atmax):
+
+    t0 = time.time()
+    i0 = 0
+    i1 = None
+    ftype = '.csv'
+
+    # remove the extension
+    ext = fpath.split('.')[-1]
+    if ext in ('sel', 'dat', 'hdf5'):
+        fpath = '.'.join(fpath.split('.')[:-1])
+    else:
+        print('invalid file extension, ignored:', fpath)
+        return
+
+    fdir = os.path.dirname(fpath)
+    fname = os.path.basename(fpath)
+    try:
+        res = windIO.LoadResults(fdir, fname, usecols=None, readdata=True)
+    except Exception as e:
+        print(f'loading {fpath} failed:')
+        print(e)
+        return
+    neq = res.sig[-1,0] - res.sig[0,0]
+
+    # add channels in memory so they can be included in the statistics
+    # they are not saved back to disk
+    res = add_channels(res)
+
+    # convert to DataFrame
+    dfres = res.sig2df()
+
+    # save the envelope
+    if do_envelope:
+        envelope(res, int_env=False, Nx=300)
+
+    # extremes and corresponding loads at the same time
+    if atmax is not None:
+        for node, val in atmax.items():
+            dfextr = loads_at_extreme(dfres, val['chans_selmax'], val['chans_selmax'])
+            dfextr.to_excel(fpath + f'_{node}_loads_at_extreme.xlsx')
+
+    # statistics
+    df_stats = res.statsdel_df(i0=i0, i1=i1, statchans='all', neq=neq,
+                               no_bins=no_bins, m=m, delchans='all')
+
+    # add path and fname as columns
+    df_stats['case_id'] = fname
+    df_stats['res_dir'] = fdir
+
+    if ftype == '.csv':
+        df_stats.to_csv(fpath+ftype)
+    elif ftype == '.h5':
+        df_stats.to_hdf(fpath+ftype, 'table', complib='zlib', complevel=9)
+    print('% 7.03f  ' % (time.time() - t0), fname)
+
+
+def par(logdir, resdir, cpus=30, chunksize=30, nostats=False, nolog=False,
+        do_envelope=True, no_bins=46, m=[3,4,6,8,9,10,12], atmax=None):
+    """Sophia has 32 CPUs per node.
+    """
+
+    # log file analysis
+    if not nolog:
+        files = glob(os.path.join(logdir, '**', '*.log'), recursive=True)
+        print(f'start processing {len(files)} logfiles from dir: {logdir}')
+        with Pool(processes=cpus) as pool:
+            # chunksize = 10 #this may take some guessing ...
+            for ind, res in enumerate(pool.imap(logcheck, files), chunksize):
+                pass
+    print()
+
+    # consider both hdf5 and sel outputs
+    if not nostats:
+        files = glob(os.path.join(resdir, '**', '*.hdf5'), recursive=True)
+        files.extend(glob(os.path.join(resdir, '**', '*.sel'), recursive=True))
+        print(f'start processing {len(files)} resfiles from dir: {resdir}')
+        # prepare the other arguments
+        combis = [(k, do_envelope, no_bins, m, atmax) for k in files]
+
+        # sequential loop
+        # for combi in combis:
+        #     calc(*combi)
+
+        # start the parallal for-loop using X number of cpus
+        with Pool(processes=cpus) as pool:
+            # This method chops the iterable into a number of chunks which it submits
+            # to the process pool as separate tasks. The (approximate) size of these
+            # chunks can be specified by setting chunksize to a positive integer.
+            for ind, res in enumerate(pool.starmap(calc, combis), chunksize):
+                pass
+    print()
+
+
+def merge_stats(post_dir='prepost-data'):
+    """Merge stats for each case into one stats table
+    """
+
+    print('start merging statistics into single table...')
+
+    os.makedirs(post_dir, exist_ok=True)
+    # -------------------------------------------------------------------------
+    # MERGE POSTPRO ON NODE APPROACH INTO ONE DataFrame
+    # -------------------------------------------------------------------------
+    path_pattern = os.path.join('res', '**', '*.csv')
+    csv_fname = '%s_statistics.csv' % SIM_ID
+    # if zipchunks:
+    #     path_pattern = os.path.join(post_dir, 'statsdel_chnk*.tar.xz')
+    fcsv = os.path.join(post_dir, csv_fname)
+    mdf = AppendDataFrames(tqdm=tqdm)
+    cols = mdf.txt2txt(fcsv, path_pattern, tarmode='r:xz', header=0, sep=',',
+                       header_fjoined=None, recursive=True)#, fname_col='[case_id]')
+    # and convert to df: takes 2 minutes
+    fdf = fcsv.replace('.csv', '.h5')
+    store = pd.HDFStore(fdf, mode='w', complevel=9, complib='zlib')
+    colnames = cols.split(',')
+    # the first column is the channel name but the header is emtpy
+    assert colnames[0] == ''
+    colnames[0] = 'channel'
+    dtypes = {col:np.float64 for col in colnames}
+    dtypes['channel'] = str
+    dtypes['case_id'] = str
+    dtypes['res_dir'] = str
+    # when using min_itemsize the column names should be valid variable names
+    # mitemsize = {'channel':60, '[case_id]':60}
+    mdf.csv2df_chunks(store, fcsv, chunksize=1000000, min_itemsize={}, sep=',',
+                      colnames=colnames, dtypes=dtypes, header=0)
+    store.close()
+    print(f'\nsaved at: {fcsv}')
+
+
+def merge_logs(post_dir='prepost-data'):
+
+    print('start merging log analysis into single table...')
+
+    os.makedirs(post_dir, exist_ok=True)
+    # -------------------------------------------------------------------------
+    # MERGE POSTPRO ON NODE APPROACH INTO ONE DataFrame
+    # -------------------------------------------------------------------------
+    lf = windIO.LogFile()
+    path_pattern = os.path.join('logfiles', '**', '*.csv')
+    # if zipchunks:
+    #     path_pattern = os.path.join(POST_DIR, 'loganalysis_chnk*.tar.xz')
+    csv_fname = '%s_ErrorLogs.csv' % SIM_ID
+    fcsv = os.path.join(post_dir, csv_fname)
+    mdf = AppendDataFrames(tqdm=tqdm)
+    # individual log file analysis does not have header, make sure to include
+    # a line for the header
+    cols = mdf.txt2txt(fcsv, path_pattern, tarmode='r:xz', header=None,
+                       header_fjoined=lf._header(), recursive=True)
+
+    # FIXME: this is due to bug in log file analysis. What is going on here??
+    # fix that some cases do not have enough columns
+    with open(fcsv.replace('.csv', '2.csv'), 'w') as f1:
+        with open(fcsv) as f2:
+            for line in f2.readlines():
+                if len(line.split(';'))==96:
+                    line = line.replace(';0.00000000000;nan;-0.0000;',
+                                        '0.00000000000;nan;-0.0000;')
+                f1.write(line)
+
+    # convert from CSV to DataFrame and save as Excel
+    df = lf.csv2df(fcsv.replace('.csv', '2.csv'))
+    # since this is mainly text it doesn't make any sense to convert to hdf5
+    # for one DLB example, csv is 350 kB, hdf 2 MB.
+    # df.to_hdf(fcsv.replace('.csv', '.h5'), 'table')
+    df.to_excel(fcsv.replace('.csv', '.xlsx'))
+    print(f"\nsaved at: {fcsv.replace('.csv', '2.csv')}")
+
+
+def save_unique_chan_names(post_dir='prepost-data'):
+
+    fdf = os.path.join(post_dir, '%s_statistics.h5' % SIM_ID)
+    df_stats = pd.read_hdf(fdf, 'table')
+    chans = df_stats['channel'].unique()
+    # TODO: print table with HTC channel names on one side, and the unique
+    # channel names on the other. Problem: one HTC line can have multiple channels
+    # do not sort, leave in same order as in the sel file
+    # chans.sort()
+    fname = os.path.join(post_dir, '%s_unique-channel-names.csv' % SIM_ID)
+    pd.DataFrame(chans).to_csv(fname)
+
+
+if __name__ == '__main__':
+
+    parser = ArgumentParser(description = "Parallel post-processing with the "
+                            "Python build-in multiprocessing module.")
+    parser.add_argument('-n', type=int, default=2, action='store',
+                        dest='cpus', help='Number of CPUs to use on the node')
+    parser.add_argument('-c', type=int, default=80, action='store',
+                        dest='chunksize', help='Chunksize Pool.')
+
+    parser.add_argument('--res', type=str, default='res', action='store',
+                        dest='resdir', help='Directory containing result files')
+    parser.add_argument('--log', type=str, default='log', action='store',
+                        dest='logdir', help='Directory containing HAWC2 log files')
+
+    parser.add_argument('--nologs', default=False, action='store_true',
+                        dest='nologs', help='Do not perform the logfile analysis.')
+    parser.add_argument('--nostats', default=False, action='store_true',
+                        dest='nostats', help='Do not calculate statistics.')
+    parser.add_argument('--envelope', default=False, action='store_true',
+                        dest='do_envelope', help='Calculate envelopes.')
+    parser.add_argument('--no_bins', default=46, action='store',
+                        dest='no_bins', help='Number of bins for the rainflow counting.')
+    parser.add_argument('--atmax', default=False, action='store', type=str,
+                        dest='atmax', help='File name that holds the channels to create '
+                        'the table at which maxima and corresponding values are selected.')
+    opt = parser.parse_args()
+
+    if opt.atmax:
+        with open(opt.atmax) as fio:
+            atmax_js = json.load(fio)
+
+    SIM_ID = os.path.basename(os.getcwd())
+    par(opt.logdir, opt.resdir, cpus=opt.cpus, chunksize=opt.chunksize,
+        nolog=opt.nologs, nostats=opt.nostats, do_envelope=opt.do_envelope,
+        no_bins=opt.no_bins, m=[3,4,6,8,9,10,12], atmax=atmax_js)
+    if not opt.nostats:
+        merge_stats(post_dir='prepost-data')
+        save_unique_chan_names(post_dir='prepost-data')
+    if not opt.nologs:
+        merge_logs(post_dir='prepost-data')