text: # -*- coding: utf-8 -*-
  """
  Created on Wed May 27 14:18:37 2015
  Updated on July 30 14:18:37 2015
  @author: ke291
  
  Equivalent and compact port of DP4.jar to python. The results
  produced are essentially equivalent, but not identical due to different
  floating point precision used in the Python (53 bits) and Java (32 bits)
  implementation.
  """
  from scipy import stats
  import bisect
  import os
  import numpy as np
  
  # Standard DP4 parameters
  meanC = 0.0
  meanH = 0.0
  stdevC = 2.269372270818724
  stdevH = 0.18731058105269952
  
  
  class DP4data:
      def __init__(self):
          self.Cshifts = []  # Carbon shifts used in DP4 calculation
          self.Cexp = []  # Carbon experimental shifts used in DP4 calculation
          self.Clabels = []  # Carbon atom labels
          self.Hshifts = []  # Proton shifts used in DP4 calculation
          self.Hexp = []  # Proton experimental shifts used in DP4 calculation
          self.Hlabels = []  # Proton atom labels
          self.Cscaled = []  # Internally scaled carbon shifts
          self.Hscaled = []  # Internally scaled proton shifts
          self.Cerrors = []  # Scaled Carbon prediction errors
          self.Herrors = []  # Scaled Proton prediction errors
          self.Cprobs = []  # Scaled carbon prediction error probabilities
          self.Hprobs = []  # Scaled proton prediction error probabilities
          self.CDP4probs = []  # Carbon DP4 probabilities
          self.HDP4probs = []  # Proton DP4 probabilities
          self.DP4probs = []  # combined Carbon and Proton DP4 probabilities
          self.output = str()  # final DP4 output
  
  
  def ProcessIsomers(DP4data, Isomers):
      # extract calculated and experimental shifts and add to DP4data instance
  
      # Carbon
  
      # make sure any shifts with missing peaks are removed from all isomers
  
      removedC = []
  
      removedH = []
  
      for iso in Isomers:
  
          DP4data.Cshifts.append([])
          DP4data.Cexp.append([])
          DP4data.Clabels.append([])
  
          for shift, exp, label in zip(iso.Cshifts, iso.Cexp, iso.Clabels):
  
              if exp != '':
                  DP4data.Cshifts[-1].append(shift)
                  DP4data.Cexp[-1].append(exp)
                  DP4data.Clabels[-1].append(label)
  
              elif label not in removedC:
  
                  removedC.append(label)
  
      for l in removedC:
  
          for j, Clabel in enumerate(DP4data.Clabels):
  
              if l in Clabel:
                  i = Clabel.index(l)
  
                  DP4data.Cshifts[j].pop(i)
  
                  DP4data.Cexp[j].pop(i)
  
                  DP4data.Clabels[j].pop(i)
  
      # proton
      for iso in Isomers:
  
          DP4data.Hshifts.append([])
          DP4data.Hexp.append([])
          DP4data.Hlabels.append([])
  
          for shift, exp, label in zip(iso.Hshifts, iso.Hexp, iso.Hlabels):
  
              if exp != '':
                  DP4data.Hshifts[-1].append(shift)
                  DP4data.Hexp[-1].append(exp)
                  DP4data.Hlabels[-1].append(label)
  
              elif label not in removedH:
  
                  removedH.append(label)
  
      for l in removedH:
  
          for j, Hlabel in enumerate(DP4data.Hlabels):
  
              if l in Hlabel:
                  i = Hlabel.index(l)
  
                  DP4data.Hshifts[j].pop(i)
  
                  DP4data.Hexp[j].pop(i)
  
                  DP4data.Hlabels[j].pop(i)
  
      return DP4data
  
  
  def InternalScaling(DP4data):
      # perform internal scaling process
  
      # calculate prediction errors
  
      if len(DP4data.Cexp[0]) > 0:
  
          for Cshifts, Cexp in zip(DP4data.Cshifts, DP4data.Cexp):
              DP4data.Cscaled.append(ScaleNMR(Cshifts, Cexp))
  
          for Cscaled, Cexp in zip(DP4data.Cscaled, DP4data.Cexp):
              DP4data.Cerrors.append([Cscaled[i] - Cexp[i] for i in range(0, len(Cscaled))])
  
      if len(DP4data.Hexp[0]) > 0:
  
          for Hshifts, Hexp in zip(DP4data.Hshifts, DP4data.Hexp):
              DP4data.Hscaled.append(ScaleNMR(Hshifts, Hexp))
  
          for Hscaled, Hexp in zip(DP4data.Hscaled, DP4data.Hexp):
              DP4data.Herrors.append([Hscaled[i] - Hexp[i] for i in range(0, len(Hscaled))])
  
      return DP4data
  
  
  def ScaleNMR(calcShifts, expShifts):
      slope, intercept, r_value, p_value, std_err = stats.linregress(expShifts,
                                                                     calcShifts)
      scaled = [(x - intercept) / slope for x in calcShifts]
  
      return scaled
  
  
  def CalcProbs(DP4data, Settings):
      # calculates probability values for each scaled prediction error value using the chosen statistical model
  
      if Settings.StatsModel == 'g' or 'm':
  
          print(Settings.StatsParamFile)
  
          if Settings.StatsParamFile == 'none':
  
              print('No stats model provided, using default')
  
              for errors in DP4data.Cerrors:
                  DP4data.Cprobs.append([SingleGausProbability(e, meanC, stdevC) for e in errors])
  
              for errors in DP4data.Herrors:
                  DP4data.Hprobs.append([SingleGausProbability(e, meanH, stdevH) for e in errors])
  
          else:
  
              print('Using stats model provided')
  
              Cmeans, Cstdevs, Hmeans, Hstdevs = ReadParamFile(Settings.StatsParamFile, Settings.StatsModel)
  
              for errors in DP4data.Cerrors:
                  DP4data.Cprobs.append([MultiGausProbability(e, Cmeans, Cstdevs) for e in errors])
  
              for errors in DP4data.Herrors:
                  DP4data.Hprobs.append([MultiGausProbability(e, Hmeans, Hstdevs) for e in errors])
  
      return DP4data
  
  
  def SingleGausProbability(error, mean, stdev):
      z = abs((error - mean) / stdev)
      cdp4 = 2 * stats.norm.cdf(-z)
  
      return cdp4
  
  
  def MultiGausProbability(error, means, stdevs):
      res = 0
  
      for mean, stdev in zip(means, stdevs):
          res += stats.norm(mean, stdev).pdf(error)
  
      return res / len(means)
  
  
  def ReadParamFile(f, t):
      infile = open(f, 'r')
      inp = infile.readlines()
      infile.close()
  
      if t not in inp[0]:
          print("Wrong parameter file type, exiting...")
          quit()
  
      if t == 'm':
          Cmeans = [float(x) for x in inp[1].split(',')]
          Cstdevs = [float(x) for x in inp[2].split(',')]
          Hmeans = [float(x) for x in inp[3].split(',')]
          Hstdevs = [float(x) for x in inp[4].split(',')]
  
          return Cmeans, Cstdevs, Hmeans, Hstdevs
  
  
  def CalcDP4(DP4data):
      # Calculate Carbon DP4 probabilities
  
      for probs in DP4data.Cprobs:
  
          DP4data.CDP4probs.append(1)
  
          for p in probs:
              DP4data.CDP4probs[-1] *= p
  
      # Calculate Proton DP4 probabilities
  
      for probs in DP4data.Hprobs:
  
          DP4data.HDP4probs.append(1)
  
          for p in probs:
              DP4data.HDP4probs[-1] *= p
  
      # Calculate Combined DP4 probabilities
  
      for Hp, Cp in zip(DP4data.HDP4probs, DP4data.CDP4probs):
          DP4data.DP4probs.append(Hp * Cp)
  
      Cs = sum(DP4data.CDP4probs)
  
      Hs = sum(DP4data.HDP4probs)
  
      Ts = sum(DP4data.DP4probs)
  
      DP4data.CDP4probs = [i / Cs for i in DP4data.CDP4probs]
  
      DP4data.HDP4probs = [i / Hs for i in DP4data.HDP4probs]
  
      DP4data.DP4probs = [i / Ts for i in DP4data.DP4probs]
  
      return DP4data
  
  
  def PrintAssignment(DP4Data):
      isomer = 0
  
      for Clabels, Cshifts, Cexp, Cscaled in zip(DP4Data.Clabels, DP4Data.Cshifts, DP4Data.Cexp, DP4Data.Cscaled):
          DP4Data.output += ("\n\nAssigned C shifts for isomer " + str(isomer + 1) + ": ")
  
          PrintNMR(Clabels, Cshifts, Cscaled, Cexp, DP4Data)
  
          isomer += 1
  
      isomer = 0
  
      for Hlabels, Hshifts, Hexp, Hscaled in zip(DP4Data.Hlabels, DP4Data.Hshifts, DP4Data.Hexp, DP4Data.Hscaled):
          DP4Data.output += ("\n\nAssigned H shifts for isomer " + str(isomer + 1) + ": ")
  
          PrintNMR(Hlabels, Hshifts, Hscaled, Hexp, DP4Data)
  
          isomer += 1
  
  
  def PrintNMR(labels, values, scaled, exp, DP4Data):
      s = np.argsort(values)
  
      svalues = np.array(values)[s]
  
      slabels = np.array(labels)[s]
  
      sscaled = np.array(scaled)[s]
  
      sexp = np.array(exp)[s]
  
      DP4Data.output += ("\nlabel, calc, corrected, exp, error")
  
      for i in range(len(labels)):
          DP4Data.output += ("\n" + format(slabels[i], "6s") + ' ' + format(svalues[i], "6.2f") + ' '
                             + format(sscaled[i], "6.2f") + ' ' + format(sexp[i], "6.2f") + ' ' +
                             format(sexp[i] - sscaled[i], "6.2f"))
  
  
  def MakeOutput(DP4Data, Isomers, Settings):
      # add some info about the calculation
  
      DP4Data.output += Settings.InputFiles[0] + "\n"
  
      DP4Data.output += "\n" + "Solvent = " + Settings.Solvent
  
      DP4Data.output += "\n" + "Force Field = " + Settings.ForceField + "\n"
  
      if 'o' in Settings.Workflow:
          DP4Data.output += "\n" + "DFT optimisation Functional = " + Settings.oFunctional
          DP4Data.output += "\n" + "DFT optimisation Basis = " + Settings.oBasisSet
  
      if 'e' in Settings.Workflow:
          DP4Data.output += "\n" + "DFT energy Functional = " + Settings.eFunctional
          DP4Data.output += "\n" + "DFT energy Basis = " + Settings.eBasisSet
  
      if 'n' in Settings.Workflow:
          DP4Data.output += "\n" + "DFT NMR Functional = " + Settings.nFunctional
          DP4Data.output += "\n" + "DFT NMR Basis = " + Settings.nBasisSet
  
      if Settings.StatsParamFile != "none":
          DP4Data.output += "\n\nStats model = " + Settings.StatsParamFile
  
      DP4Data.output += "\n\nNumber of isomers = " + str(len(Isomers))
  
      c = 1
  
      for i in Isomers:
          DP4Data.output += "\nNumber of conformers for isomer " + str(c) + " = " + str(len(i.Conformers))
  
          c += 1
  
      PrintAssignment(DP4Data)
  
      DP4Data.output += ("\n\nResults of DP4 using Proton: ")
  
      for i, p in enumerate(DP4Data.HDP4probs):
          DP4Data.output += ("\nIsomer " + str(i + 1) + ": " + format(p * 100, "4.1f") + "%")
  
      DP4Data.output += ("\n\nResults of DP4 using Carbon: ")
  
      for i, p in enumerate(DP4Data.CDP4probs):
          DP4Data.output += ("\nIsomer " + str(i + 1) + ": " + format(p * 100, "4.1f") + "%")
  
      DP4Data.output += ("\n\nResults of DP4: ")
  
      for i, p in enumerate(DP4Data.DP4probs):
          DP4Data.output += ("\nIsomer " + str(i + 1) + ": " + format(p * 100, "4.1f") + "%")
  
      print("number of c protons = " + str(len(Isomers[0].Hlabels)))
      print("number of c carbons = " + str(len(Isomers[0].Clabels)))
  
      print("number of e protons = " + str(len(DP4Data.Hexp[0])))
      print("number of e carbons = " + str(len(DP4Data.Cexp[0])))
  
      print(DP4Data.output)
  
      if Settings.OutputFolder == '':
  
          out = open(str(os.getcwd()) + "/" + str(Settings.InputFiles[0] + "NMR.dp4"), "w+")
  
      else:
  
          out = open(Settings.OutputFolder + "/" + str(Settings.InputFiles[0] + "NMR.dp4"), "w+")
  
      out.write(DP4Data.output)
  
      out.close()
  
      return DP4Data