text: import shutil
  from pathlib import Path
  import itertools
  import math
  import numpy as np
  import pandas as pd
  from statistics import mean
  from scipy.optimize import minimize_scalar
  from matplotlib import pyplot as plt
  
  def removeChars(s):
      for c in [' ', '\\', '/', '^']:
          s = s.replace(c, '')
      return s
  
  
  def rchop(s, suffix):
      if suffix and s.endswith(suffix):
          return s[:-len(suffix)]
      return s
  
  def GetAggGroup(row, word_group):
    comps = row['Filename'].split('_')
    comps[4] = 'k'
    key = '_'.join(comps)
    return '&'.join(word_group[key])
  
  def GetMean(row, formant):
      colIndex = range(2, 21)
      sum = 0
      for i in colIndex:
        sum += row[formant + '_' + str(i) + '_pp']
      return sum/19.0
  
  def GetPerson(row):
    comps = row['Filename'].split('_')
    return int(comps[2])
  
  def GetPersonPos(row):
    comps = row['Filename'].split('_')
    return 'S_a_' + comps[2] + '_' + comps[4]
  
  def GetPos(row):
    comps = row['Filename'].split('_')
    return comps[4]
  
  def LoadFormantData():
      all_data = []
      for input in sorted(input_base_dir.glob('S*.CSV')):
          print(input)
          single_df = pd.read_csv(input, converters={
              'Annotation': removeChars}, na_values=['--undefined--', 'null'], skipinitialspace=True, sep=r"\s*[,]\s*", engine='python')
          single_df.drop(single_df.filter(regex="Unname"), axis=1, inplace=True)
          clean_df = single_df.dropna(subset=['Annotation'] + kCols)
          clean_df = clean_df.copy()
          clean_df['Table'] = input
          num_nan = len(single_df) - len(clean_df)
          if num_nan > 0:
              print(input, 'Dropped', num_nan)
          all_data.append(clean_df)
      df = pd.concat(all_data, ignore_index=True)
      df = df[df.Annotation != 'null']
      # df.drop_duplicates(subset=['Filename', 'Annotation'], keep='first', inplace=True)
      print('Num files', len(all_data))
      print('Final', df.shape)
      output_df = pd.concat(
          [df[['Filename']],
           df[['Annotation']],
           df[['F1_mean']],
           df[['F1_stdev']],
           df[['F2_mean']],
           df[['F2_stdev']],
           df[kCols],
           ], axis=1)
      return output_df
  
  def AddDataPoint(F, Z, row, cols):
    for col in cols:
      F.append(float(row[col]))
      Z.append(float(row[col+'_Z']))
  
  
  def NormalizeColumn(df, formant):
      colIndex = range(1, 22)
      cols = []
      for i in colIndex:
        col = formant + '_' + str(i)
        cols.append(col)
        df[col + '_Z'] = (df[col] - df[formant+ '_mean']) / df[formant+ '_stdev']
      print(cols)
  
      F_a = []
      F_i = []
      F_u = []
      F_c = []
      F_e = []
      F_ic = []
      F_ue = []
      Z_a = []
      Z_i = []
      Z_u = []
      Z_c = []
      Z_e = []
      Z_ic = []
      Z_ue = []
      for _, row in df.iterrows():
          comps = row['Filename'].split('_')
          age_gender = int(row['Filename'].split('_')[2])
          if not age_gender % 2 == 1:  # not male
              continue
          if not comps[0].startswith('norm'):
              continue
          if row['Annotation'] == 'a':
            AddDataPoint(F_a, Z_a, row, cols)
            continue
          if row['Annotation'] == 'i':
            AddDataPoint(F_i, Z_i, row, cols)
            continue
          if row['Annotation'] == 'u':
            AddDataPoint(F_u, Z_u, row, cols)
            continue
          if row['Annotation'] == 'c':
            AddDataPoint(F_c, Z_c, row, cols)
            continue
          if row['Annotation'] == 'e':
            AddDataPoint(F_e, Z_e, row, cols)
            continue
          if row['Annotation'] == 'ic':
            AddDataPoint(F_ic, Z_ic, row, cols)
            continue
          if row['Annotation'] == 'ue':
            AddDataPoint(F_ue, Z_ue, row, cols)
      F_bar_a = mean(F_a)
      F_bar_i = mean(F_i)
      F_bar_u = mean(F_u)
      F_bar_c = mean(F_c)
      F_bar_e = mean(F_e)
      F_bar_ic = mean(F_ic)
      F_bar_ue = mean(F_ue)
      F_bar = mean(F_a + F_i + F_u + F_c + F_e + F_ic + F_ue)
      Z_bar_a = mean(Z_a)
      Z_bar_i = mean(Z_i)
      Z_bar_u = mean(Z_u)
      Z_bar_c = mean(Z_c)
      Z_bar_e = mean(Z_e)
      Z_bar_ic = mean(Z_ic)
      Z_bar_ue = mean(Z_ue)
      Z_bar = mean(Z_a + Z_i + Z_u + Z_c + Z_e + Z_ic + Z_ue)
  
      numerator = 0
      numerator += abs(F_bar_a - F_bar)
      numerator += abs(F_bar_i - F_bar)
      numerator += abs(F_bar_u - F_bar)
      numerator += abs(F_bar_c - F_bar)
      numerator += abs(F_bar_e - F_bar)
      numerator += abs(F_bar_ic - F_bar)
      numerator += abs(F_bar_ue - F_bar)
      denominator = 0
      denominator = abs(Z_bar_a - Z_bar)
      denominator += abs(Z_bar_i - Z_bar)
      denominator += abs(Z_bar_u - Z_bar)
      denominator += abs(Z_bar_c - Z_bar)
      denominator += abs(Z_bar_e - Z_bar)
      denominator += abs(Z_bar_ic - Z_bar)
      denominator += abs(Z_bar_ue - Z_bar)
      R = numerator / denominator
      df[formant+'_R'] = R
      for col in cols:
        df[col + '_p'] = df[col + '_Z'] * df[formant + '_R']
      F_p = []
      for _, row in df.iterrows():
          comps = row['Filename'].split('_')
          age_gender = int(row['Filename'].split('_')[2])
          if not age_gender % 2 == 1:  # not male
              continue
          if not comps[0].startswith('norm'):
              continue
          for col in cols:
            F_p.append(float(row[col+'_p']))
      F_p_bar = mean(F_p)
      df[formant+'_bar'] = F_bar
      df[formant+'_p_bar'] = F_p_bar
      for col in cols:
        df[col + '_pp'] = df[col + '_p'] - \
          (df[formant + '_p_bar'] - df[formant + '_bar'])
        df[col + '_pp_bark'] = df[col+'_pp'].apply(lambda x:7 * math.asinh(x/650))
      return df
  
  
  def ComputeDelta(df):
      df['deltaF1'] = df['F1_16_pp'] - df['F1_6_pp']
      df['deltaF2'] = df['F2_16_pp'] - df['F2_6_pp']
      df['mean_F1'] = df.apply(GetMean, formant = 'F1', axis=1)
      df['mean_F2'] = df.apply(GetMean, formant = 'F2', axis=1)
      df['bark_mean_F1'] = df['mean_F1'].apply(lambda x:7 * math.asinh(x/650))
      df['bark_mean_F2'] = df['mean_F2'].apply(lambda x:7 * math.asinh(x/650))
      # df['delta_barkF1'] = df['F1_16_pp_bark'] - df['F1_6_pp_bark']
      # df['delta_barkF2'] = df['F2_16_pp_bark'] - df['F2_6_pp_bark']
      # df['BD_1'] = df['F1_1_pp_bark'] - df['F1_21_pp_bark']
      # df['BD_2'] = df['F2_1_pp_bark'] - df['F2_21_pp_bark']
      # df['BD'] = (df['BD_1']**2 + df['BD_2']**2)**(1/2)
      return df
  
  
  def GetBreak(row, formant):
      x = np.arange(0, 19)
      y = row[formant+'_2': formant+'_20'].to_numpy(dtype='float')
      coeff = np.polyfit(x, y, 4)
      line = np.poly1d(coeff)
      linedd = np.polyder(line, 2)
      linedd_max = minimize_scalar(-linedd, bounds=(0, 19), method='bounded')
      return linedd_max.x / 19.0
  
  
  def ComputeBreak(df):
      df['breakF1'] = df.apply(lambda row: GetBreak(row, 'F1'), axis=1)
      df['breakF2'] = df.apply(lambda row: GetBreak(row, 'F2'), axis=1)
      return df
  
  def ComputeAggGroup(df):
    word_group = dict()
    for _, row in df.iterrows():
      comps = row['Filename'].split('_')
      pos = comps[4]
      comps[4] = 'k'
      key = '_'.join(comps)
      if key in word_group:
        word_group[key].append('S'+pos+'='+row['Annotation'])
      else:
        word_group[key] = ['S'+pos+'='+row['Annotation']]
    df['AggregationGroup'] = df.apply(GetAggGroup, word_group = word_group, axis=1)
    df['PersonPos'] = df.apply(GetPersonPos, axis=1)
    df['Pos'] = df.apply(GetPos, axis=1)
    return df
  
  def GetBDMaxIndex(row):
    for i in range(1,21):
      if row['BD_'+str(i)] == row['BD_max']:
        return i
    raise Exception("Should not reach here")
  
  def GetBDAvg(row):
    bd_sum = 0
    for i in range(1,21):
      bd_sum += row['BD_'+str(i)]
    return bd_sum / 20
  
  def ComputeBD(df):
    df['BD'] = ((df['F1_1_pp_bark'] - df['F1_21_pp_bark'])**2 + 
                (df['F2_1_pp_bark'] - df['F2_21_pp_bark'])**2)**(0.5)
    for i in range(1,21):
      si = str(i)
      sj = str(i+1)
      df['BD_'+si] = ((df['F1_'+si+'_pp_bark'] - df['F1_'+sj+'_pp_bark'])**2 + 
                      (df['F2_'+si+'_pp_bark'] - df['F2_'+sj+'_pp_bark'])**2)**(0.5)
    colBD = ['BD_'+str(i) for i in range(1,21)]
    df['BD_max'] = df[colBD].max(axis=1)
    df['BD_max_index'] = df.apply(lambda row: GetBDMaxIndex(row), axis=1)
    df['BD_avg'] = df.apply(lambda row: GetBDAvg(row), axis=1)
    return df
  
  
  input_base_dir = Path('./analysis/item_ae/normalized/')
  output_base_dir = Path('./analysis/output_ae/')
  shutil.rmtree(output_base_dir, ignore_errors=True)
  output_base_dir.mkdir(parents=True, exist_ok=True)
  
  cols1 = ['F1_' + str(i) for i in range(1, 22)]
  cols2 = ['F2_' + str(i) for i in range(1, 22)]
  kCols = cols1 + cols2
  df_formant = LoadFormantData()
  df_formant = ComputeBreak(df_formant)
  df_formant = NormalizeColumn(df_formant, 'F1')
  df_formant = NormalizeColumn(df_formant, 'F2')
  df_formant = ComputeDelta(df_formant)
  df_formant = ComputeAggGroup(df_formant)
  df_formant = ComputeBD(df_formant)
  # df_formant['Person'] = df_formant.apply(lambda row: GetPerson(row), axis=1)
  output_df = pd.concat(
      [df_formant[['Filename']],
       df_formant[['Annotation']],
       df_formant[['PersonPos']],
       df_formant[['AggregationGroup']],
       df_formant[['breakF1']],
       df_formant[['breakF2']],
       df_formant[['deltaF1']],
       df_formant[['deltaF2']],
       df_formant[['mean_F1']],
       df_formant[['mean_F2']],
       df_formant[['bark_mean_F1']],
       df_formant[['bark_mean_F2']],
       df_formant[['BD']],
       df_formant[['BD_avg']],
       df_formant[['BD_max']],
       df_formant[['BD_max_index']],
      # df_formant.loc[:, df_formant.columns.str.startswith("BD")],
      #  df_formant.loc[:, df_formant.columns.str.startswith("F1_")],
      #  df_formant.loc[:, df_formant.columns.str.startswith("F2_")],
       ], axis=1)
  output_df.to_csv(output_base_dir / 'S_normalized.csv', index=False)
  output_bd_df = pd.concat(
      [df_formant[['Filename']],
       df_formant[['Annotation']],
       df_formant[['PersonPos']],
       df_formant[['AggregationGroup']],
       df_formant.loc[:, df_formant.columns.str.startswith("BD")],
      #  df_formant.loc[:, df_formant.columns.str.startswith("F1_")],
      #  df_formant.loc[:, df_formant.columns.str.startswith("F2_")],
       ], axis=1)
  output_bd_df.to_csv(output_base_dir / 'S_normalized_BD.csv', index=False)