<commit_before>'''
Created on July 10, 2013

@author: alfoa
'''
from __future__ import division, print_function, unicode_literals, absolute_import
import warnings
warnings.simplefilter('default',DeprecationWarning)

#External Modules------------------------------------------------------------------------------------
import sys
import numpy as np
from sklearn import tree
from scipy import spatial
from scipy import interpolate
from scipy import integrate
import os
from glob import glob
import copy
import Datas
import math
#External Modules End--------------------------------------------------------------------------------

#Internal Modules------------------------------------------------------------------------------------
from utils import toString, toBytes, first, returnPrintTag, returnPrintPostTag
from BaseClasses import Assembler
import SupervisedLearning
#Internal Modules End--------------------------------------------------------------------------------

#def error(*objs):
#  print("ERROR: ", *objs, file=sys.stderr)

'''
  ***************************************
  *  SPECIALIZED PostProcessor CLASSES  *
  ***************************************
'''

class BasePostProcessor(Assembler):
  '''This is the base class for postprocessors'''
  def __init__(self):
    self.type              = self.__class__.__name__  # pp type
    self.name              = self.__class__.__name__  # pp name
    self.assemblerObjects  = {}                       # {MainClassName(e.g.Distributions):[class(e.g.Models),type(e.g.ROM),objectName]}
    self.requiredAssObject = (False,([],[]))          # tuple. first entry boolean flag. True if the XML parser must look for assembler objects;
                                                      # second entry tuple.first entry list of object can be retrieved, second entry multiplicity (-1,-2,-n means optional (max 1 object,2 object, no number limit))
    self.debug             = False


  def whatDoINeed(self):
    '''
    This method is used mainly by the Simulation class at the Step construction stage.
    It is used for inquiring the class, which is implementing the method, about the kind of objects the class needs to
    be initialize. It is an abstract method -> It must be implemented in the derived class!
    NB. In this implementation, the method only calls the self.interface.whatDoINeed() method
    @ In , None, None
    @ Out, needDict, dictionary of objects needed (class:tuple(object type{if None, Simulation does not check the type}, object name))
    '''
    needDict = self._localWhatDoINeed()
    for val in self.assemblerObjects.values():
      for value in val:
        if value[0] not in needDict.keys(): needDict[value[0]] = []
        needDict[value[0]].append((value[1],value[2]))
    return needDict

  def _localWhatDoINeed(self):
    '''
    local whatDoINeed method.
    In here there is the common implementation if the  self.assemblerObjects dictionary has the form:
    {MainClassName(e.g.Distributions):[class(e.g.Models),type(e.g.ROM),objectName]}
    '''
    return {}

  def generateAssembler(self,initDict):
    '''
    This method is used mainly by the Simulation class at the Step construction stage.
    It is used for sending to the instanciated class, which is implementing the method, the objects that have been requested through "whatDoINeed" method
    It is an abstract method -> It must be implemented in the derived class!
    NB. In this implementation, the method only calls the self.interface.generateAssembler(initDict) method
    @ In , initDict, dictionary ({'mainClassName(e.g., DataBases):{specializedObjectName(e.g.,DataBaseForSystemCodeNamedWolf):ObjectInstance}'})
    @ Out, None, None
    '''
    self._localGenerateAssembler(initDict)

  def _localGenerateAssembler(self,initDict):
    ''' see generateAssembler method '''
    pass

  def initialize(self, runInfo, inputs, initDict) :
    #if 'externalFunction' in initDict.keys(): self.externalFunction = initDict['externalFunction']
    self.inputs           = inputs

  def _readMoreXML(self,xmlNode):
    self.type = xmlNode.tag
    self.name = xmlNode.attrib['name']
    self.printTag = self.type.ljust(25)
    if 'debug' in xmlNode.attrib.keys():self.debug = bool(xmlNode.attrib['debug'])
    if self.requiredAssObject[0]:
        testObjects = {}
        for token in self.requiredAssObject[1][0]:
            testObjects[token] = 0
        found = False
        for subNode in xmlNode:
            for token in self.requiredAssObject[1][0]:
                if subNode.tag in token:
                    found = True
                    if 'class' not in subNode.attrib.keys(): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> In '+self.type+' PostProcessor ' + self.name+ ', block ' + subNode.tag + ' does not have the attribute class!!')
                    if  subNode.tag not in self.assemblerObjects.keys(): self.assemblerObjects[subNode.tag] = []
                    self.assemblerObjects[subNode.tag].append([subNode.attrib['class'],subNode.attrib['type'],subNode.text])
                    testObjects[token] += 1
        if not found:
            for tofto in self.requiredAssObject[1][0]:
                if not str(self.requiredAssObject[1][1][0]).strip().startswith('-'):
                    raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> the required object ' +tofto+ ' is missed in the definition of the '+self.type+' PostProcessor!')
        # test the objects found
        else:
            for cnt,tofto in enumerate(self.requiredAssObject[1][0]):
                numerosity = str(self.requiredAssObject[1][1][cnt])
                if numerosity.strip().startswith('-'):
                # optional
                    if tofto in testObjects.keys():
                        numerosity = numerosity.replace('-', '').replace('n',str(testObjects[tofto]))
                        if testObjects[tofto] != int(numerosity): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> Only '+numerosity+' '+tofto+' object/s is/are optionally required. PostProcessor '+self.name + ' got '+str(testObjects[tofto]) + '!')
                else:
                # required
                    if tofto not in testObjects.keys(): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> Required object/s "'+tofto+'" not found. PostProcessor '+self.name + '!')
                    else:
                        numerosity = numerosity.replace('n',str(testObjects[tofto]))
                        if testObjects[tofto] != int(numerosity): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> Only '+numerosity+' '+tofto+' object/s is/are required. PostProcessor '+self.name + ' got '+str(testObjects[tofto]) + '!')
    self._localReadMoreXML(xmlNode)
  def inputToInternal(self,currentInput): return [(copy.deepcopy(currentInput))]

  def run(self, Input): pass

class SafestPoint(BasePostProcessor):
  '''
  It searches for the probability-weighted safest point inside the space of the system controllable variables
  '''
  def __init__(self):
    BasePostProcessor.__init__(self)
    self.controllableDist = {}                                    #dictionary created upon the .xml input file reading. It stores the distributions for each controllale variable.
    self.nonControllableDist = {}                                 #dictionary created upon the .xml input file reading. It stores the distributions for each non-controllale variable.
    self.controllableGrid = {}                                    #dictionary created upon the .xml input file reading. It stores the grid type ('value' or 'CDF'), the number of steps and the step length for each controllale variable.
    self.nonControllableGrid = {}                                 #dictionary created upon the .xml input file reading. It stores the grid type ('value' or 'CDF'), the number of steps and the step length for each non-controllale variable.
    self.gridInfo = {}                                            #dictionary contaning the grid type ('value' or 'CDF'), the grid construction type ('equal', set by default) and the list of sampled points for each variable.
    self.controllableOrd = []                                     #list contaning the controllable variables' names in the same order as they appear inside the controllable space (self.controllableSpace)
    self.nonControllableOrd = []                                  #list contaning the controllable variables' names in the same order as they appear inside the non-controllable space (self.nonControllableSpace)
    self.surfPointsMatrix = None                                  #2D-matrix containing the coordinates of the points belonging to the failure boundary (coordinates are derived from both the controllable and non-controllable space)
    self.stat = returnInstance('BasicStatistics')                 #instantiation of the 'BasicStatistics' processor, which is used to compute the expected value of the safest point through the coordinates and probability values collected in the 'run' function
    self.stat.what = ['expectedValue']
    self.requiredAssObject = (True,(['Distribution'],['n']))
    self.printTag = returnPrintTag('POSTPROCESSOR SAFESTPOINT')

  def _localGenerateAssembler(self,initDict):
    ''' see generateAssembler method '''
    for varName, distName in self.controllableDist.items():
      if distName not in initDict['Distributions'].keys():
        raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> distribution ' +distName+ ' not found.')
      self.controllableDist[varName] = initDict['Distributions'][distName]
    for varName, distName in self.nonControllableDist.items():
      if distName not in initDict['Distributions'].keys():
        raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> distribution ' +distName+ ' not found.')
      self.nonControllableDist[varName] = initDict['Distributions'][distName]

  def _localReadMoreXML(self,xmlNode):
    for child in xmlNode:
      if child.tag == 'controllable':
        for childChild in child:
          if childChild.tag == 'variable':
            varName = childChild.attrib['name']
            for childChildChild in childChild:
              if childChildChild.tag == 'distribution':
                self.controllableDist[varName] = childChildChild.text
              elif childChildChild.tag == 'grid':
                if 'type' in childChildChild.attrib.keys():
                  if 'steps' in childChildChild.attrib.keys():
                    self.controllableGrid[varName] = (childChildChild.attrib['type'], int(childChildChild.attrib['steps']), float(childChildChild.text))
                  else:
                    raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> number of steps missing after the grid call.')
                else:
                  raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> grid type missing after the grid call.')
              else:
                raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid labels after the variable call. Only "distribution" and "grid" are accepted.')
          else:
            raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid or missing labels after the controllable variables call. Only "variable" is accepted.')
      elif child.tag == 'non-controllable':
        for childChild in child:
          if childChild.tag == 'variable':
            varName = childChild.attrib['name']
            for childChildChild in childChild:
              if childChildChild.tag == 'distribution':
                self.nonControllableDist[varName] = childChildChild.text
              elif childChildChild.tag == 'grid':
                if 'type' in childChildChild.attrib.keys():
                  if 'steps' in childChildChild.attrib.keys():
                    self.nonControllableGrid[varName] = (childChildChild.attrib['type'], int(childChildChild.attrib['steps']), float(childChildChild.text))
                  else:
                    raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> number of steps missing after the grid call.')
                else:
                  raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> grid type missing after the grid call.')
              else:
                raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid labels after the variable call. Only "distribution" and "grid" are accepted.')
          else:
            raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid or missing labels after the controllable variables call. Only "variable" is accepted.')
      #else:
      #  if child.tag != 'Assembler': raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid or missing labels after the post-processor call. Only "controllable", "non-controllable" and "Assembler" are accepted.')
    if self.debug:
      print('CONTROLLABLE DISTRIBUTIONS:')
      print(self.controllableDist)
      print('CONTROLLABLE GRID:')
      print(self.controllableGrid)
      print('NON-CONTROLLABLE DISTRIBUTIONS:')
      print(self.nonControllableDist)
      print('NON-CONTROLLABLE GRID:')
      print(self.nonControllableGrid)

  def initialize(self,runInfo,inputs,initDict):
    self.__gridSetting__()
    self.__gridGeneration__()
    self.inputToInternal(inputs)
    self.stat.parameters['targets'] = self.controllableOrd
    self.stat.initialize(runInfo,inputs,initDict)
    if self.debug:
      print('GRID INFO:')
      print(self.gridInfo)
      print('N-DIMENSIONAL CONTROLLABLE SPACE:')
      print(self.controllableSpace)
      print('N-DIMENSIONAL NON-CONTROLLABLE SPACE:')
      print(self.nonControllableSpace)
      print('CONTROLLABLE VARIABLES ORDER:')
      print(self.controllableOrd)
      print('NON-CONTROLLABLE VARIABLES ORDER:')
      print(self.nonControllableOrd)
      print('SURFACE POINTS MATRIX:')
      print(self.surfPointsMatrix)

  def __gridSetting__(self,constrType='equal'):
    for varName in self.controllableGrid.keys():
      if self.controllableGrid[varName][0] == 'value':
        self.__stepError__(float(self.controllableDist[varName].lowerBound),float(self.controllableDist[varName].upperBound),self.controllableGrid[varName][1],self.controllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.controllableGrid[varName][0], constrType, [float(self.controllableDist[varName].lowerBound)+self.controllableGrid[varName][2]*i for i in range(self.controllableGrid[varName][1]+1)])
      elif self.controllableGrid[varName][0] == 'CDF':
        self.__stepError__(0,1,self.controllableGrid[varName][1],self.controllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.controllableGrid[varName][0], constrType, [self.controllableGrid[varName][2]*i for i in range(self.controllableGrid[varName][1]+1)])
      else:
        raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> inserted invalid grid type. Only "value" and "CDF" are accepted.')
    for varName in self.nonControllableGrid.keys():
      if self.nonControllableGrid[varName][0] == 'value':
        self.__stepError__(float(self.nonControllableDist[varName].lowerBound),float(self.nonControllableDist[varName].upperBound),self.nonControllableGrid[varName][1],self.nonControllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.nonControllableGrid[varName][0], constrType, [float(self.nonControllableDist[varName].lowerBound)+self.nonControllableGrid[varName][2]*i for i in range(self.nonControllableGrid[varName][1]+1)])
      elif self.nonControllableGrid[varName][0] == 'CDF':
        self.__stepError__(0,1,self.nonControllableGrid[varName][1],self.nonControllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.nonControllableGrid[varName][0], constrType, [self.nonControllableGrid[varName][2]*i for i in range(self.nonControllableGrid[varName][1]+1)])
      else:
        raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> inserted invalid grid type. Only "value" and "CDF" are accepted.')

  def __stepError__(self,lowerBound,upperBound,steps,tol,varName):
    if upperBound-lowerBound<steps*tol:
      raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> inserted number of steps or tolerance for variable ' +varName+ ' exceeds its limit.')

  def __gridGeneration__(self):
    NotchesByVar = [None]*len(self.controllableGrid.keys())
    controllableSpaceSize = None
    for varId, varName in enumerate(self.controllableGrid.keys()):
      NotchesByVar[varId] = self.controllableGrid[varName][1]+1
      self.controllableOrd.append(varName)
    controllableSpaceSize = tuple(NotchesByVar+[len(self.controllableGrid.keys())])
    self.controllableSpace = np.zeros(controllableSpaceSize)
    iterIndex = np.nditer(self.controllableSpace,flags=['multi_index'])
    while not iterIndex.finished:
      coordIndex = iterIndex.multi_index[-1]
      varName = self.controllableGrid.keys()[coordIndex]
      notchPos = iterIndex.multi_index[coordIndex]
      if self.gridInfo[varName][0] == 'CDF':
        valList = []
        for probVal in self.gridInfo[varName][2]:
          valList.append(self.controllableDist[varName].cdf(probVal))
        self.controllableSpace[iterIndex.multi_index] = valList[notchPos]
      else:
        self.controllableSpace[iterIndex.multi_index] = self.gridInfo[varName][2][notchPos]
      iterIndex.iternext()
    NotchesByVar = [None]*len(self.nonControllableGrid.keys())
    nonControllableSpaceSize = None
    for varId, varName in enumerate(self.nonControllableGrid.keys()):
      NotchesByVar[varId] = self.nonControllableGrid[varName][1]+1
      self.nonControllableOrd.append(varName)
    nonControllableSpaceSize = tuple(NotchesByVar+[len(self.nonControllableGrid.keys())])
    self.nonControllableSpace = np.zeros(nonControllableSpaceSize)
    iterIndex = np.nditer(self.nonControllableSpace,flags=['multi_index'])
    while not iterIndex.finished:
      coordIndex = iterIndex.multi_index[-1]
      varName = self.nonControllableGrid.keys()[coordIndex]
      notchPos = iterIndex.multi_index[coordIndex]
      if self.gridInfo[varName][0] == 'CDF':
        valList = []
        for probVal in self.gridInfo[varName][2]:
          valList.append(self.nonControllableDist[varName].cdf(probVal))
        self.nonControllableSpace[iterIndex.multi_index] = valList[notchPos]
      else:
        self.nonControllableSpace[iterIndex.multi_index] = self.gridInfo[varName][2][notchPos]
      iterIndex.iternext()

  def inputToInternal(self,currentInput):
    for item in currentInput:
      if item.type == 'TimePointSet':
        self.surfPointsMatrix = np.zeros((len(item.getParam('output',item.getParaKeys('outputs')[-1])),len(self.gridInfo.keys())+1))
        k=0
        for varName in self.controllableOrd:
          self.surfPointsMatrix[:,k] = item.getParam('input',varName)
          k+=1
        for varName in self.nonControllableOrd:
          self.surfPointsMatrix[:,k] = item.getParam('input',varName)
          k+=1
        self.surfPointsMatrix[:,k] = item.getParam('output',item.getParaKeys('outputs')[-1])

  def run(self,Input):
    nearestPointsInd = []
    dataCollector = Datas.returnInstance('TimePointSet')
    dataCollector.type = 'TimePointSet'
    surfTree = spatial.KDTree(copy.copy(self.surfPointsMatrix[:,0:self.surfPointsMatrix.shape[-1]-1]))
    self.controllableSpace.shape = (np.prod(self.controllableSpace.shape[0:len(self.controllableSpace.shape)-1]),self.controllableSpace.shape[-1])
    self.nonControllableSpace.shape = (np.prod(self.nonControllableSpace.shape[0:len(self.nonControllableSpace.shape)-1]),self.nonControllableSpace.shape[-1])
    if self.debug:
      print('RESHAPED CONTROLLABLE SPACE:')
      print(self.controllableSpace)
      print('RESHAPED NON-CONTROLLABLE SPACE:')
      print(self.nonControllableSpace)
    for ncLine in range(self.nonControllableSpace.shape[0]):
      queryPointsMatrix = np.append(self.controllableSpace,np.tile(self.nonControllableSpace[ncLine,:],(self.controllableSpace.shape[0],1)),axis=1)
      print('QUERIED POINTS MATRIX:')
      print(queryPointsMatrix)
      nearestPointsInd = surfTree.query(queryPointsMatrix)[-1]
      distList = []
      indexList = []
      probList = []
      for index in range(len(nearestPointsInd)):
        if self.surfPointsMatrix[np.where(np.prod(surfTree.data[nearestPointsInd[index],0:self.surfPointsMatrix.shape[-1]-1] == self.surfPointsMatrix[:,0:self.surfPointsMatrix.shape[-1]-1],axis=1))[0][0],-1] == 1:
          distList.append(np.sqrt(np.sum(np.power(queryPointsMatrix[index,0:self.controllableSpace.shape[-1]]-surfTree.data[nearestPointsInd[index],0:self.controllableSpace.shape[-1]],2))))
          indexList.append(index)
      if distList == []:
        raise Exception(self.printTag+': ' +returnPrintPostTag("ERROR") + '-> no safest point found for the current set of non-controllable variables: ' +str(self.nonControllableSpace[ncLine,:])+ '.')
      else:
        for cVarIndex in range(len(self.controllableOrd)):
          dataCollector.updateInputValue(self.controllableOrd[cVarIndex],copy.copy(queryPointsMatrix[indexList[distList.index(max(distList))],cVarIndex]))
        for ncVarIndex in range(len(self.nonControllableOrd)):
          dataCollector.updateInputValue(self.nonControllableOrd[ncVarIndex],copy.copy(queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex]))
          if queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex] == self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].lowerBound:
            if self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].type == 'Bernoulli':
              prob = 1-self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].p
            else:
              if self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][0] == 'CDF':
                prob = self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2)
              else:
                prob = self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].lowerBound+self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))
          elif queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex] == self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].upperBound:
            if self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].type == 'Bernoulli':
              prob = self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].p
            else:
              if self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][0] == 'CDF':
                prob = self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2)
              else:
                prob = 1-self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].upperBound-self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))
          else:
            if self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][0] == 'CDF':
              prob = self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]
            else:
              prob = self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex]+self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))-self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex]-self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))
          probList.append(prob)
      dataCollector.updateOutputValue('Probability',np.prod(probList))
      dataCollector.updateMetadata('ProbabilityWeight',np.prod(probList))
    dataCollector.updateMetadata('ExpectedSafestPointCoordinates',self.stat.run(dataCollector)['expectedValue'])
    if self.debug:
      print(dataCollector.getParametersValues('input'))
      print(dataCollector.getParametersValues('output'))
      print(dataCollector.getMetadata('ExpectedSafestPointCoordinates'))
    return dataCollector

  def collectOutput(self,finishedjob,output):
    if finishedjob.returnEvaluation() == -1:
      raise Exception(self.printTag+': ' +returnPrintPostTag("ERROR") + '-> no available output to collect (the run is likely not over yet).')
    else:
      dataCollector = finishedjob.returnEvaluation()[1]
      if output.type != 'TimePointSet':
        raise Exception(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> output item type must be "TimePointSet".')
      else:
        if not output.isItEmpty():
          raise Exception(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> output item must be empty.')
        else:
          for key,value in dataCollector.getParametersValues('input').items():
            for val in value: output.updateInputValue(key, val)
          for key,value in dataCollector.getParametersValues('output').items():
            for val in value: output.updateOutputValue(key,val)
          for key,value in dataCollector.getAllMetadata().items(): output.updateMetadata(key,value)

class ComparisonStatistics(BasePostProcessor):
  '''
  ComparisonStatistics is to calculate statistics that compare
  two different codes or code to experimental data.
  '''

  class CompareGroup:
    def __init__(self):
      self.dataPulls = []
      self.referenceData = {}

  def __init__(self):
    BasePostProcessor.__init__(self)
    self.dataDict = {} #Dictionary of all the input data, keyed by the name
    self.compare_groups = [] #List of each of the groups that will be compared
    #self.dataPulls = [] #List of data references that will be used
    #self.referenceData = [] #List of reference (experimental) data
    self.methodInfo = {} #Information on what stuff to do.

  def inputToInternal(self,currentInput):
    return [(currentInput)]

  def initialize(self, runInfo, inputs, initDict):
    BasePostProcessor.initialize(self, runInfo, inputs, initDict)
    #print("runInfo",runInfo,"inputs",inputs,"initDict",initDict)

  def _localReadMoreXML(self,xmlNode):
    for outer in xmlNode:
      if outer.tag == 'compare':
        compare_group = ComparisonStatistics.CompareGroup()
        for child in outer:
          if child.tag == 'data':
            dataName = child.text
            splitName = dataName.split("|")
            name, kind = splitName[:2]
            rest = splitName[2:]
            compare_group.dataPulls.append([name, kind, rest])
            #print("xml dataName",dataName,self.dataPulls[-1])
          elif child.tag == 'reference':
            compare_group.referenceData = dict(child.attrib)
        self.compare_groups.append(compare_group)
      if outer.tag == 'kind':
        self.methodInfo['kind'] = outer.text
        if 'num_bins' in outer.attrib:
          self.methodInfo['num_bins'] = int(outer.attrib['num_bins'])
        if 'bin_method' in outer.attrib:
          self.methodInfo['bin_method'] = outer.attrib['bin_method'].lower()


  def run(self, Input): # inObj,workingDir=None):
    '''
     Function to finalize the filter => execute the filtering
     @ Out, None      : Print of the CSV file
    '''
    for aInput in Input:
      self.dataDict[aInput.name] = aInput
    #print("input",Input,"input name",Input.name,"input input",Input.getParametersValues('inputs'),
    #      "input output",Input.getParametersValues('outputs'))

  def collectOutput(self,finishedjob,output):
    print("finishedjob",finishedjob,"output",output)
    dataToProcess = []
    for compare_group in self.compare_groups:
      dataPulls = compare_group.dataPulls
      reference = compare_group.referenceData
      found_datas = []
      for name, kind, rest in dataPulls:
        data = self.dataDict[name].getParametersValues(kind)
        #print("dataPull",dataPull) #("result",self.dataDict[name].getParametersValues(kind))
        if len(rest) == 1:
          #print("dataPart",data[rest[0]])
          found_datas.append(data[rest[0]])
      dataToProcess.append((dataPulls,found_datas,reference))
    #print("dataToProcess",dataToProcess)
    csv = open(output,"w")
    def print_csv(*args):
      print(*args,file=csv,sep=',')
    for dataPulls, datas, reference in dataToProcess:
      for dataPull, data in zip(dataPulls,datas):
        data_stats = process_data(dataPull, data, self.methodInfo)
        data_keys = set(data_stats.keys())
        print_csv('"'+str(dataPull)+'"')
        print_csv('"num_bins"',data_stats['num_bins'])
        counts = data_stats['counts']
        bins = data_stats['bins']
        count_sum = sum(counts)
        bin_boundaries = [data_stats['low']]+bins+[data_stats['high']]
        print_csv('"bin_boundary"','"bin_midpoint"','"bin_count"','"normalized_bin_count"','"f_prime"','"cdf"')
        cdf = [0.0]*len(counts)
        midpoints = [0.0]*len(counts)
        cdf_sum = 0.0
        for i in range(len(counts)):
          f_0 = counts[i]/count_sum
          cdf_sum += f_0
          cdf[i] = cdf_sum
          midpoints[i] = (bin_boundaries[i]+bin_boundaries[i+1])/2.0
        cdf_func = create_interp(midpoints,cdf,0.0,1.0,'quadratic')
        f_prime_data = [0.0]*len(counts)
        for i in range(len(counts)):
          h = bin_boundaries[i+1] - bin_boundaries[i]
          n_count = counts[i]/count_sum #normalized count
          f_0 = cdf[i]
          if i + 1 < len(counts):
            f_1 = cdf[i+1]
          else:
            f_1 = 1.0
          if i + 2 < len(counts):
            f_2 = cdf[i+2]
          else:
            f_2 = 1.0
          #f_prime = (f_1 - f_0)/h
          #print(f_0,f_1,f_2,h,f_prime)
          f_prime = (-1.5*f_0 + 2.0*f_1 + -0.5*f_2)/h
          f_prime_data[i] = f_prime
          print_csv(bin_boundaries[i+1],midpoints[i],counts[i],n_count,f_prime,cdf[i])
        pdf_func = create_interp(midpoints,f_prime_data,0.0,0.0,'quadratic')
        data_keys -= set({'num_bins','counts','bins'})
        for key in data_keys:
          print_csv('"'+key+'"',data_stats[key])
        print("data_stats",data_stats)
        print_graphs(csv, reference, data_stats, cdf_func,pdf_func)
        print_csv()

def normal(x,mu=0.0,sigma=1.0):
  return (1.0/(sigma*math.sqrt(2*math.pi)))*math.exp(-(x - mu)**2/(2.0*sigma**2))

def normal_cdf(x,mu=0.0,sigma=1.0):
  return 0.5*(1.0+math.erf((x-mu)/(sigma*math.sqrt(2.0))))

def skew_normal(x,alpha,xi,omega):
  def phi(x):
    return (1.0/math.sqrt(2*math.pi))*math.exp(-(x**2)/2.0)

  def Phi(x):
    return 0.5*(1+math.erf(x/math.sqrt(2)))

  return (2.0/omega)*phi((x-xi)/omega)*Phi(alpha*(x-xi)/omega)

def create_interp(x, y, low_fill, high_fill, kind='linear'):
  interp = interpolate.interp1d(x, y, kind)
  low = x[0]
  high = x[-1]
  def my_interp(x):
    try:
      return interp(x)+0.0
    except ValueError:
      if x <= low:
        return low_fill
      else:
        return high_fill
  return my_interp

def simpson(f, a, b, n):
  h = (b - a) / float(n)
  sum = f(a) + f(b)
  for i in range(1,n, 2):
    sum += 4*f(a + i*h)
  for i in range(2, n-1, 2):
    sum += 2*f(a + i*h)

  return sum * h / 3.0

def print_graphs(csv, reference, data_stats, cdf_func,pdf_func):
  if "mean" not in reference:
    return
  ref_mean = float(reference["mean"])
  ref_stddev = float(reference["sigma"])
  calc_mean = data_stats["mean"]
  calc_stddev = data_stats["stdev"]
  calc_omega = data_stats["omega"]
  calc_xi = data_stats["xi"]
  calc_alpha = data_stats["alpha"]
  low = min(ref_mean - 3.0*ref_stddev,calc_mean - 3.0*calc_stddev)
  high = max(ref_mean + 3.0*ref_stddev,calc_mean + 3.0*calc_stddev)
  low_low = min(ref_mean - 5.0*ref_stddev,calc_mean - 5.0*calc_stddev)
  high_high = max(ref_mean + 5.0*ref_stddev,calc_mean + 5.0*calc_stddev)
  print("Graph from ",low,"to",high)
  n = int(math.ceil((high-low)/data_stats['min_bin_size']))
  interval = (high - low)/n
  def print_csv(*args):
    print(*args,file=csv,sep=',')

  print_csv('"x"','"reference"','"reference_cdf"','"calculated"','"interpolated_cdf"','"interpolated_pdf"')

  for i in range(n):
    x = low+interval*i
    print_csv(x,normal(x,ref_mean,ref_stddev),normal_cdf(x,ref_mean,ref_stddev),skew_normal(x,calc_alpha,calc_xi,calc_omega),cdf_func(x),pdf_func(x))

  def f_z(z):
    return simpson(lambda x: normal(x,ref_mean,ref_stddev)*skew_normal(x-z,calc_alpha,calc_xi,calc_omega), low_low, high_high, 1000)

  mid_z = ref_mean-calc_mean
  low_z = mid_z - 3.0*max(ref_stddev,calc_stddev)
  high_z = mid_z + 3.0*max(ref_stddev,calc_stddev)
  print_csv('"z"','"f_z(z)"')
  z_n = 20
  interval_z = (high_z - low_z)/z_n
  for i in range(z_n):
    z = low_z + interval_z*i
    print_csv(z,f_z(z))
  cdf_area_difference = simpson(lambda x:abs(cdf_func(x)-normal_cdf(x,ref_mean,ref_stddev)),low_low,high_high,100000)

  def first_moment_simpson(f, a, b, n):
    return simpson(lambda x:x*f(x), a, b, n)

  pdf_common_area = simpson(lambda x:min(normal(x,ref_mean,ref_stddev),
                                         skew_normal(x,calc_alpha,calc_xi,calc_omega)),
                            low_low,high_high,100000)<commit_msg>Switching to generic calculating of the cdf area difference and others.<commit_after>'''
Created on July 10, 2013

@author: alfoa
'''
from __future__ import division, print_function, unicode_literals, absolute_import
import warnings
warnings.simplefilter('default',DeprecationWarning)

#External Modules------------------------------------------------------------------------------------
import sys
import numpy as np
from sklearn import tree
from scipy import spatial
from scipy import interpolate
from scipy import integrate
import os
from glob import glob
import copy
import Datas
import math
#External Modules End--------------------------------------------------------------------------------

#Internal Modules------------------------------------------------------------------------------------
from utils import toString, toBytes, first, returnPrintTag, returnPrintPostTag
from BaseClasses import Assembler
import SupervisedLearning
#Internal Modules End--------------------------------------------------------------------------------

#def error(*objs):
#  print("ERROR: ", *objs, file=sys.stderr)

'''
  ***************************************
  *  SPECIALIZED PostProcessor CLASSES  *
  ***************************************
'''

class BasePostProcessor(Assembler):
  '''This is the base class for postprocessors'''
  def __init__(self):
    self.type              = self.__class__.__name__  # pp type
    self.name              = self.__class__.__name__  # pp name
    self.assemblerObjects  = {}                       # {MainClassName(e.g.Distributions):[class(e.g.Models),type(e.g.ROM),objectName]}
    self.requiredAssObject = (False,([],[]))          # tuple. first entry boolean flag. True if the XML parser must look for assembler objects;
                                                      # second entry tuple.first entry list of object can be retrieved, second entry multiplicity (-1,-2,-n means optional (max 1 object,2 object, no number limit))
    self.debug             = False


  def whatDoINeed(self):
    '''
    This method is used mainly by the Simulation class at the Step construction stage.
    It is used for inquiring the class, which is implementing the method, about the kind of objects the class needs to
    be initialize. It is an abstract method -> It must be implemented in the derived class!
    NB. In this implementation, the method only calls the self.interface.whatDoINeed() method
    @ In , None, None
    @ Out, needDict, dictionary of objects needed (class:tuple(object type{if None, Simulation does not check the type}, object name))
    '''
    needDict = self._localWhatDoINeed()
    for val in self.assemblerObjects.values():
      for value in val:
        if value[0] not in needDict.keys(): needDict[value[0]] = []
        needDict[value[0]].append((value[1],value[2]))
    return needDict

  def _localWhatDoINeed(self):
    '''
    local whatDoINeed method.
    In here there is the common implementation if the  self.assemblerObjects dictionary has the form:
    {MainClassName(e.g.Distributions):[class(e.g.Models),type(e.g.ROM),objectName]}
    '''
    return {}

  def generateAssembler(self,initDict):
    '''
    This method is used mainly by the Simulation class at the Step construction stage.
    It is used for sending to the instanciated class, which is implementing the method, the objects that have been requested through "whatDoINeed" method
    It is an abstract method -> It must be implemented in the derived class!
    NB. In this implementation, the method only calls the self.interface.generateAssembler(initDict) method
    @ In , initDict, dictionary ({'mainClassName(e.g., DataBases):{specializedObjectName(e.g.,DataBaseForSystemCodeNamedWolf):ObjectInstance}'})
    @ Out, None, None
    '''
    self._localGenerateAssembler(initDict)

  def _localGenerateAssembler(self,initDict):
    ''' see generateAssembler method '''
    pass

  def initialize(self, runInfo, inputs, initDict) :
    #if 'externalFunction' in initDict.keys(): self.externalFunction = initDict['externalFunction']
    self.inputs           = inputs

  def _readMoreXML(self,xmlNode):
    self.type = xmlNode.tag
    self.name = xmlNode.attrib['name']
    self.printTag = self.type.ljust(25)
    if 'debug' in xmlNode.attrib.keys():self.debug = bool(xmlNode.attrib['debug'])
    if self.requiredAssObject[0]:
        testObjects = {}
        for token in self.requiredAssObject[1][0]:
            testObjects[token] = 0
        found = False
        for subNode in xmlNode:
            for token in self.requiredAssObject[1][0]:
                if subNode.tag in token:
                    found = True
                    if 'class' not in subNode.attrib.keys(): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> In '+self.type+' PostProcessor ' + self.name+ ', block ' + subNode.tag + ' does not have the attribute class!!')
                    if  subNode.tag not in self.assemblerObjects.keys(): self.assemblerObjects[subNode.tag] = []
                    self.assemblerObjects[subNode.tag].append([subNode.attrib['class'],subNode.attrib['type'],subNode.text])
                    testObjects[token] += 1
        if not found:
            for tofto in self.requiredAssObject[1][0]:
                if not str(self.requiredAssObject[1][1][0]).strip().startswith('-'):
                    raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> the required object ' +tofto+ ' is missed in the definition of the '+self.type+' PostProcessor!')
        # test the objects found
        else:
            for cnt,tofto in enumerate(self.requiredAssObject[1][0]):
                numerosity = str(self.requiredAssObject[1][1][cnt])
                if numerosity.strip().startswith('-'):
                # optional
                    if tofto in testObjects.keys():
                        numerosity = numerosity.replace('-', '').replace('n',str(testObjects[tofto]))
                        if testObjects[tofto] != int(numerosity): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> Only '+numerosity+' '+tofto+' object/s is/are optionally required. PostProcessor '+self.name + ' got '+str(testObjects[tofto]) + '!')
                else:
                # required
                    if tofto not in testObjects.keys(): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> Required object/s "'+tofto+'" not found. PostProcessor '+self.name + '!')
                    else:
                        numerosity = numerosity.replace('n',str(testObjects[tofto]))
                        if testObjects[tofto] != int(numerosity): raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> Only '+numerosity+' '+tofto+' object/s is/are required. PostProcessor '+self.name + ' got '+str(testObjects[tofto]) + '!')
    self._localReadMoreXML(xmlNode)
  def inputToInternal(self,currentInput): return [(copy.deepcopy(currentInput))]

  def run(self, Input): pass

class SafestPoint(BasePostProcessor):
  '''
  It searches for the probability-weighted safest point inside the space of the system controllable variables
  '''
  def __init__(self):
    BasePostProcessor.__init__(self)
    self.controllableDist = {}                                    #dictionary created upon the .xml input file reading. It stores the distributions for each controllale variable.
    self.nonControllableDist = {}                                 #dictionary created upon the .xml input file reading. It stores the distributions for each non-controllale variable.
    self.controllableGrid = {}                                    #dictionary created upon the .xml input file reading. It stores the grid type ('value' or 'CDF'), the number of steps and the step length for each controllale variable.
    self.nonControllableGrid = {}                                 #dictionary created upon the .xml input file reading. It stores the grid type ('value' or 'CDF'), the number of steps and the step length for each non-controllale variable.
    self.gridInfo = {}                                            #dictionary contaning the grid type ('value' or 'CDF'), the grid construction type ('equal', set by default) and the list of sampled points for each variable.
    self.controllableOrd = []                                     #list contaning the controllable variables' names in the same order as they appear inside the controllable space (self.controllableSpace)
    self.nonControllableOrd = []                                  #list contaning the controllable variables' names in the same order as they appear inside the non-controllable space (self.nonControllableSpace)
    self.surfPointsMatrix = None                                  #2D-matrix containing the coordinates of the points belonging to the failure boundary (coordinates are derived from both the controllable and non-controllable space)
    self.stat = returnInstance('BasicStatistics')                 #instantiation of the 'BasicStatistics' processor, which is used to compute the expected value of the safest point through the coordinates and probability values collected in the 'run' function
    self.stat.what = ['expectedValue']
    self.requiredAssObject = (True,(['Distribution'],['n']))
    self.printTag = returnPrintTag('POSTPROCESSOR SAFESTPOINT')

  def _localGenerateAssembler(self,initDict):
    ''' see generateAssembler method '''
    for varName, distName in self.controllableDist.items():
      if distName not in initDict['Distributions'].keys():
        raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> distribution ' +distName+ ' not found.')
      self.controllableDist[varName] = initDict['Distributions'][distName]
    for varName, distName in self.nonControllableDist.items():
      if distName not in initDict['Distributions'].keys():
        raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> distribution ' +distName+ ' not found.')
      self.nonControllableDist[varName] = initDict['Distributions'][distName]

  def _localReadMoreXML(self,xmlNode):
    for child in xmlNode:
      if child.tag == 'controllable':
        for childChild in child:
          if childChild.tag == 'variable':
            varName = childChild.attrib['name']
            for childChildChild in childChild:
              if childChildChild.tag == 'distribution':
                self.controllableDist[varName] = childChildChild.text
              elif childChildChild.tag == 'grid':
                if 'type' in childChildChild.attrib.keys():
                  if 'steps' in childChildChild.attrib.keys():
                    self.controllableGrid[varName] = (childChildChild.attrib['type'], int(childChildChild.attrib['steps']), float(childChildChild.text))
                  else:
                    raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> number of steps missing after the grid call.')
                else:
                  raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> grid type missing after the grid call.')
              else:
                raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid labels after the variable call. Only "distribution" and "grid" are accepted.')
          else:
            raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid or missing labels after the controllable variables call. Only "variable" is accepted.')
      elif child.tag == 'non-controllable':
        for childChild in child:
          if childChild.tag == 'variable':
            varName = childChild.attrib['name']
            for childChildChild in childChild:
              if childChildChild.tag == 'distribution':
                self.nonControllableDist[varName] = childChildChild.text
              elif childChildChild.tag == 'grid':
                if 'type' in childChildChild.attrib.keys():
                  if 'steps' in childChildChild.attrib.keys():
                    self.nonControllableGrid[varName] = (childChildChild.attrib['type'], int(childChildChild.attrib['steps']), float(childChildChild.text))
                  else:
                    raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> number of steps missing after the grid call.')
                else:
                  raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> grid type missing after the grid call.')
              else:
                raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid labels after the variable call. Only "distribution" and "grid" are accepted.')
          else:
            raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid or missing labels after the controllable variables call. Only "variable" is accepted.')
      #else:
      #  if child.tag != 'Assembler': raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> invalid or missing labels after the post-processor call. Only "controllable", "non-controllable" and "Assembler" are accepted.')
    if self.debug:
      print('CONTROLLABLE DISTRIBUTIONS:')
      print(self.controllableDist)
      print('CONTROLLABLE GRID:')
      print(self.controllableGrid)
      print('NON-CONTROLLABLE DISTRIBUTIONS:')
      print(self.nonControllableDist)
      print('NON-CONTROLLABLE GRID:')
      print(self.nonControllableGrid)

  def initialize(self,runInfo,inputs,initDict):
    self.__gridSetting__()
    self.__gridGeneration__()
    self.inputToInternal(inputs)
    self.stat.parameters['targets'] = self.controllableOrd
    self.stat.initialize(runInfo,inputs,initDict)
    if self.debug:
      print('GRID INFO:')
      print(self.gridInfo)
      print('N-DIMENSIONAL CONTROLLABLE SPACE:')
      print(self.controllableSpace)
      print('N-DIMENSIONAL NON-CONTROLLABLE SPACE:')
      print(self.nonControllableSpace)
      print('CONTROLLABLE VARIABLES ORDER:')
      print(self.controllableOrd)
      print('NON-CONTROLLABLE VARIABLES ORDER:')
      print(self.nonControllableOrd)
      print('SURFACE POINTS MATRIX:')
      print(self.surfPointsMatrix)

  def __gridSetting__(self,constrType='equal'):
    for varName in self.controllableGrid.keys():
      if self.controllableGrid[varName][0] == 'value':
        self.__stepError__(float(self.controllableDist[varName].lowerBound),float(self.controllableDist[varName].upperBound),self.controllableGrid[varName][1],self.controllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.controllableGrid[varName][0], constrType, [float(self.controllableDist[varName].lowerBound)+self.controllableGrid[varName][2]*i for i in range(self.controllableGrid[varName][1]+1)])
      elif self.controllableGrid[varName][0] == 'CDF':
        self.__stepError__(0,1,self.controllableGrid[varName][1],self.controllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.controllableGrid[varName][0], constrType, [self.controllableGrid[varName][2]*i for i in range(self.controllableGrid[varName][1]+1)])
      else:
        raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> inserted invalid grid type. Only "value" and "CDF" are accepted.')
    for varName in self.nonControllableGrid.keys():
      if self.nonControllableGrid[varName][0] == 'value':
        self.__stepError__(float(self.nonControllableDist[varName].lowerBound),float(self.nonControllableDist[varName].upperBound),self.nonControllableGrid[varName][1],self.nonControllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.nonControllableGrid[varName][0], constrType, [float(self.nonControllableDist[varName].lowerBound)+self.nonControllableGrid[varName][2]*i for i in range(self.nonControllableGrid[varName][1]+1)])
      elif self.nonControllableGrid[varName][0] == 'CDF':
        self.__stepError__(0,1,self.nonControllableGrid[varName][1],self.nonControllableGrid[varName][2],varName)
        self.gridInfo[varName] = (self.nonControllableGrid[varName][0], constrType, [self.nonControllableGrid[varName][2]*i for i in range(self.nonControllableGrid[varName][1]+1)])
      else:
        raise NameError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> inserted invalid grid type. Only "value" and "CDF" are accepted.')

  def __stepError__(self,lowerBound,upperBound,steps,tol,varName):
    if upperBound-lowerBound<steps*tol:
      raise IOError(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> inserted number of steps or tolerance for variable ' +varName+ ' exceeds its limit.')

  def __gridGeneration__(self):
    NotchesByVar = [None]*len(self.controllableGrid.keys())
    controllableSpaceSize = None
    for varId, varName in enumerate(self.controllableGrid.keys()):
      NotchesByVar[varId] = self.controllableGrid[varName][1]+1
      self.controllableOrd.append(varName)
    controllableSpaceSize = tuple(NotchesByVar+[len(self.controllableGrid.keys())])
    self.controllableSpace = np.zeros(controllableSpaceSize)
    iterIndex = np.nditer(self.controllableSpace,flags=['multi_index'])
    while not iterIndex.finished:
      coordIndex = iterIndex.multi_index[-1]
      varName = self.controllableGrid.keys()[coordIndex]
      notchPos = iterIndex.multi_index[coordIndex]
      if self.gridInfo[varName][0] == 'CDF':
        valList = []
        for probVal in self.gridInfo[varName][2]:
          valList.append(self.controllableDist[varName].cdf(probVal))
        self.controllableSpace[iterIndex.multi_index] = valList[notchPos]
      else:
        self.controllableSpace[iterIndex.multi_index] = self.gridInfo[varName][2][notchPos]
      iterIndex.iternext()
    NotchesByVar = [None]*len(self.nonControllableGrid.keys())
    nonControllableSpaceSize = None
    for varId, varName in enumerate(self.nonControllableGrid.keys()):
      NotchesByVar[varId] = self.nonControllableGrid[varName][1]+1
      self.nonControllableOrd.append(varName)
    nonControllableSpaceSize = tuple(NotchesByVar+[len(self.nonControllableGrid.keys())])
    self.nonControllableSpace = np.zeros(nonControllableSpaceSize)
    iterIndex = np.nditer(self.nonControllableSpace,flags=['multi_index'])
    while not iterIndex.finished:
      coordIndex = iterIndex.multi_index[-1]
      varName = self.nonControllableGrid.keys()[coordIndex]
      notchPos = iterIndex.multi_index[coordIndex]
      if self.gridInfo[varName][0] == 'CDF':
        valList = []
        for probVal in self.gridInfo[varName][2]:
          valList.append(self.nonControllableDist[varName].cdf(probVal))
        self.nonControllableSpace[iterIndex.multi_index] = valList[notchPos]
      else:
        self.nonControllableSpace[iterIndex.multi_index] = self.gridInfo[varName][2][notchPos]
      iterIndex.iternext()

  def inputToInternal(self,currentInput):
    for item in currentInput:
      if item.type == 'TimePointSet':
        self.surfPointsMatrix = np.zeros((len(item.getParam('output',item.getParaKeys('outputs')[-1])),len(self.gridInfo.keys())+1))
        k=0
        for varName in self.controllableOrd:
          self.surfPointsMatrix[:,k] = item.getParam('input',varName)
          k+=1
        for varName in self.nonControllableOrd:
          self.surfPointsMatrix[:,k] = item.getParam('input',varName)
          k+=1
        self.surfPointsMatrix[:,k] = item.getParam('output',item.getParaKeys('outputs')[-1])

  def run(self,Input):
    nearestPointsInd = []
    dataCollector = Datas.returnInstance('TimePointSet')
    dataCollector.type = 'TimePointSet'
    surfTree = spatial.KDTree(copy.copy(self.surfPointsMatrix[:,0:self.surfPointsMatrix.shape[-1]-1]))
    self.controllableSpace.shape = (np.prod(self.controllableSpace.shape[0:len(self.controllableSpace.shape)-1]),self.controllableSpace.shape[-1])
    self.nonControllableSpace.shape = (np.prod(self.nonControllableSpace.shape[0:len(self.nonControllableSpace.shape)-1]),self.nonControllableSpace.shape[-1])
    if self.debug:
      print('RESHAPED CONTROLLABLE SPACE:')
      print(self.controllableSpace)
      print('RESHAPED NON-CONTROLLABLE SPACE:')
      print(self.nonControllableSpace)
    for ncLine in range(self.nonControllableSpace.shape[0]):
      queryPointsMatrix = np.append(self.controllableSpace,np.tile(self.nonControllableSpace[ncLine,:],(self.controllableSpace.shape[0],1)),axis=1)
      print('QUERIED POINTS MATRIX:')
      print(queryPointsMatrix)
      nearestPointsInd = surfTree.query(queryPointsMatrix)[-1]
      distList = []
      indexList = []
      probList = []
      for index in range(len(nearestPointsInd)):
        if self.surfPointsMatrix[np.where(np.prod(surfTree.data[nearestPointsInd[index],0:self.surfPointsMatrix.shape[-1]-1] == self.surfPointsMatrix[:,0:self.surfPointsMatrix.shape[-1]-1],axis=1))[0][0],-1] == 1:
          distList.append(np.sqrt(np.sum(np.power(queryPointsMatrix[index,0:self.controllableSpace.shape[-1]]-surfTree.data[nearestPointsInd[index],0:self.controllableSpace.shape[-1]],2))))
          indexList.append(index)
      if distList == []:
        raise Exception(self.printTag+': ' +returnPrintPostTag("ERROR") + '-> no safest point found for the current set of non-controllable variables: ' +str(self.nonControllableSpace[ncLine,:])+ '.')
      else:
        for cVarIndex in range(len(self.controllableOrd)):
          dataCollector.updateInputValue(self.controllableOrd[cVarIndex],copy.copy(queryPointsMatrix[indexList[distList.index(max(distList))],cVarIndex]))
        for ncVarIndex in range(len(self.nonControllableOrd)):
          dataCollector.updateInputValue(self.nonControllableOrd[ncVarIndex],copy.copy(queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex]))
          if queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex] == self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].lowerBound:
            if self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].type == 'Bernoulli':
              prob = 1-self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].p
            else:
              if self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][0] == 'CDF':
                prob = self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2)
              else:
                prob = self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].lowerBound+self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))
          elif queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex] == self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].upperBound:
            if self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].type == 'Bernoulli':
              prob = self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].p
            else:
              if self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][0] == 'CDF':
                prob = self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2)
              else:
                prob = 1-self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].upperBound-self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))
          else:
            if self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][0] == 'CDF':
              prob = self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]
            else:
              prob = self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex]+self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))-self.nonControllableDist[self.nonControllableOrd[ncVarIndex]].cdf(queryPointsMatrix[indexList[distList.index(max(distList))],len(self.controllableOrd)+ncVarIndex]-self.nonControllableGrid[self.nonControllableOrd[ncVarIndex]][2]/float(2))
          probList.append(prob)
      dataCollector.updateOutputValue('Probability',np.prod(probList))
      dataCollector.updateMetadata('ProbabilityWeight',np.prod(probList))
    dataCollector.updateMetadata('ExpectedSafestPointCoordinates',self.stat.run(dataCollector)['expectedValue'])
    if self.debug:
      print(dataCollector.getParametersValues('input'))
      print(dataCollector.getParametersValues('output'))
      print(dataCollector.getMetadata('ExpectedSafestPointCoordinates'))
    return dataCollector

  def collectOutput(self,finishedjob,output):
    if finishedjob.returnEvaluation() == -1:
      raise Exception(self.printTag+': ' +returnPrintPostTag("ERROR") + '-> no available output to collect (the run is likely not over yet).')
    else:
      dataCollector = finishedjob.returnEvaluation()[1]
      if output.type != 'TimePointSet':
        raise Exception(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> output item type must be "TimePointSet".')
      else:
        if not output.isItEmpty():
          raise Exception(self.printTag+': ' +returnPrintPostTag('ERROR') + '-> output item must be empty.')
        else:
          for key,value in dataCollector.getParametersValues('input').items():
            for val in value: output.updateInputValue(key, val)
          for key,value in dataCollector.getParametersValues('output').items():
            for val in value: output.updateOutputValue(key,val)
          for key,value in dataCollector.getAllMetadata().items(): output.updateMetadata(key,value)

class ComparisonStatistics(BasePostProcessor):
  '''
  ComparisonStatistics is to calculate statistics that compare
  two different codes or code to experimental data.
  '''

  class CompareGroup:
    def __init__(self):
      self.dataPulls = []
      self.referenceData = {}

  def __init__(self):
    BasePostProcessor.__init__(self)
    self.dataDict = {} #Dictionary of all the input data, keyed by the name
    self.compare_groups = [] #List of each of the groups that will be compared
    #self.dataPulls = [] #List of data references that will be used
    #self.referenceData = [] #List of reference (experimental) data
    self.methodInfo = {} #Information on what stuff to do.

  def inputToInternal(self,currentInput):
    return [(currentInput)]

  def initialize(self, runInfo, inputs, initDict):
    BasePostProcessor.initialize(self, runInfo, inputs, initDict)
    #print("runInfo",runInfo,"inputs",inputs,"initDict",initDict)

  def _localReadMoreXML(self,xmlNode):
    for outer in xmlNode:
      if outer.tag == 'compare':
        compare_group = ComparisonStatistics.CompareGroup()
        for child in outer:
          if child.tag == 'data':
            dataName = child.text
            splitName = dataName.split("|")
            name, kind = splitName[:2]
            rest = splitName[2:]
            compare_group.dataPulls.append([name, kind, rest])
            #print("xml dataName",dataName,self.dataPulls[-1])
          elif child.tag == 'reference':
            compare_group.referenceData = dict(child.attrib)
        self.compare_groups.append(compare_group)
      if outer.tag == 'kind':
        self.methodInfo['kind'] = outer.text
        if 'num_bins' in outer.attrib:
          self.methodInfo['num_bins'] = int(outer.attrib['num_bins'])
        if 'bin_method' in outer.attrib:
          self.methodInfo['bin_method'] = outer.attrib['bin_method'].lower()


  def run(self, Input): # inObj,workingDir=None):
    '''
     Function to finalize the filter => execute the filtering
     @ Out, None      : Print of the CSV file
    '''
    for aInput in Input:
      self.dataDict[aInput.name] = aInput
    #print("input",Input,"input name",Input.name,"input input",Input.getParametersValues('inputs'),
    #      "input output",Input.getParametersValues('outputs'))

  def collectOutput(self,finishedjob,output):
    print("finishedjob",finishedjob,"output",output)
    dataToProcess = []
    for compare_group in self.compare_groups:
      dataPulls = compare_group.dataPulls
      reference = compare_group.referenceData
      found_datas = []
      for name, kind, rest in dataPulls:
        data = self.dataDict[name].getParametersValues(kind)
        #print("dataPull",dataPull) #("result",self.dataDict[name].getParametersValues(kind))
        if len(rest) == 1:
          #print("dataPart",data[rest[0]])
          found_datas.append(data[rest[0]])
      dataToProcess.append((dataPulls,found_datas,reference))
    #print("dataToProcess",dataToProcess)
    csv = open(output,"w")
    def print_csv(*args):
      print(*args,file=csv,sep=',')
    for dataPulls, datas, reference in dataToProcess:
      for dataPull, data in zip(dataPulls,datas):
        data_stats = process_data(dataPull, data, self.methodInfo)
        data_keys = set(data_stats.keys())
        print_csv('"'+str(dataPull)+'"')
        print_csv('"num_bins"',data_stats['num_bins'])
        counts = data_stats['counts']
        bins = data_stats['bins']
        count_sum = sum(counts)
        bin_boundaries = [data_stats['low']]+bins+[data_stats['high']]
        print_csv('"bin_boundary"','"bin_midpoint"','"bin_count"','"normalized_bin_count"','"f_prime"','"cdf"')
        cdf = [0.0]*len(counts)
        midpoints = [0.0]*len(counts)
        cdf_sum = 0.0
        for i in range(len(counts)):
          f_0 = counts[i]/count_sum
          cdf_sum += f_0
          cdf[i] = cdf_sum
          midpoints[i] = (bin_boundaries[i]+bin_boundaries[i+1])/2.0
        cdf_func = create_interp(midpoints,cdf,0.0,1.0,'quadratic')
        f_prime_data = [0.0]*len(counts)
        for i in range(len(counts)):
          h = bin_boundaries[i+1] - bin_boundaries[i]
          n_count = counts[i]/count_sum #normalized count
          f_0 = cdf[i]
          if i + 1 < len(counts):
            f_1 = cdf[i+1]
          else:
            f_1 = 1.0
          if i + 2 < len(counts):
            f_2 = cdf[i+2]
          else:
            f_2 = 1.0
          #f_prime = (f_1 - f_0)/h
          #print(f_0,f_1,f_2,h,f_prime)
          f_prime = (-1.5*f_0 + 2.0*f_1 + -0.5*f_2)/h
          f_prime_data[i] = f_prime
          print_csv(bin_boundaries[i+1],midpoints[i],counts[i],n_count,f_prime,cdf[i])
        pdf_func = create_interp(midpoints,f_prime_data,0.0,0.0,'quadratic')
        data_keys -= set({'num_bins','counts','bins'})
        for key in data_keys:
          print_csv('"'+key+'"',data_stats[key])
        print("data_stats",data_stats)
        if "mean" in reference:
          ref_data_stats = {"mean":float(reference["mean"]),
                            "stdev":float(reference["sigma"]),
                            "min_bin_size":float(reference["sigma"])/2.0}
          ref_pdf = lambda x:normal(x,ref_data_stats["mean"],ref_data_stats["stdev"])
          ref_cdf = lambda x:normal_cdf(x,ref_data_stats["mean"],ref_data_stats["stdev"])
          print_graphs(csv, [(ref_data_stats, ref_cdf, ref_pdf,"ref"),
                             (data_stats, cdf_func, pdf_func,"calculated")])
        print_csv()

def normal(x,mu=0.0,sigma=1.0):
  return (1.0/(sigma*math.sqrt(2*math.pi)))*math.exp(-(x - mu)**2/(2.0*sigma**2))

def normal_cdf(x,mu=0.0,sigma=1.0):
  return 0.5*(1.0+math.erf((x-mu)/(sigma*math.sqrt(2.0))))

def skew_normal(x,alpha,xi,omega):
  def phi(x):
    return (1.0/math.sqrt(2*math.pi))*math.exp(-(x**2)/2.0)

  def Phi(x):
    return 0.5*(1+math.erf(x/math.sqrt(2)))

  return (2.0/omega)*phi((x-xi)/omega)*Phi(alpha*(x-xi)/omega)

def create_interp(x, y, low_fill, high_fill, kind='linear'):
  interp = interpolate.interp1d(x, y, kind)
  low = x[0]
  high = x[-1]
  def my_interp(x):
    try:
      return interp(x)+0.0
    except ValueError:
      if x <= low:
        return low_fill
      else:
        return high_fill
  return my_interp

def simpson(f, a, b, n):
  h = (b - a) / float(n)
  sum = f(a) + f(b)
  for i in range(1,n, 2):
    sum += 4*f(a + i*h)
  for i in range(2, n-1, 2):
    sum += 2*f(a + i*h)

  return sum * h / 3.0

def print_graphs(csv, functions):
  """prints graphs of the functions.
  The functions are a list of (data_stats_dict, cdf_function, pdf_function,name)
  """

  means = [x[0]["mean"] for x in functions]
  stddevs = [x[0]["stdev"] for x in functions]
  cdfs = [x[1] for x in functions]
  pdfs = [x[2] for x in functions]
  names = [x[3] for x in functions]
  low = min([m - 3.0*s for m,s in zip(means,stddevs)])
  high = max([m + 3.0*s for m,s in zip(means,stddevs)])
  low_low = min([m - 5.0*s for m,s in zip(means,stddevs)])
  high_high = max([m + 5.0*s for m,s in zip(means,stddevs)])
  min_bin_size = min([x[0]["min_bin_size"] for x in functions])
  print("Graph from ",low,"to",high)
  n = int(math.ceil((high-low)/min_bin_size))
  interval = (high - low)/n
  def print_csv(*args):
    print(*args,file=csv,sep=',')

  def print_csv_part(*args):
    print(*args,file=csv,sep=',',end=',')

  print_csv_part('"x"')
  for name in names:
    print_csv_part('"'+name+'_cdf"','"'+name+'_pdf"')
  print_csv()

  for i in range(n):
    x = low+interval*i
    print_csv_part(x)
    for stats, cdf, pdf, name in functions:
      print_csv_part(cdf(x),pdf(x))
    print_csv()

  def f_z(z):
    return simpson(lambda x: pdfs[0](x)*pdfs[1](x-z), low_low, high_high, 1000)

  mid_z = means[0]-means[1]
  low_z = mid_z - 3.0*max(stddevs[0],stddevs[1])
  high_z = mid_z + 3.0*max(stddevs[0],stddevs[1])
  print_csv('"z"','"f_z(z)"')
  z_n = 20
  interval_z = (high_z - low_z)/z_n
  for i in range(z_n):
    z = low_z + interval_z*i
    print_csv(z,f_z(z))
  cdf_area_difference = simpson(lambda x:abs(cdfs[1](x)-cdfs[0](x)),low_low,high_high,100000)

  def first_moment_simpson(f, a, b, n):
    return simpson(lambda x:x*f(x), a, b, n)

  pdf_common_area = simpson(lambda x:min(pdfs[0](x),pdfs[1](x)),
                            low_low,high_high,100000)
{'fraction_non_alphanumeric': 0.13876498569967843, 'fraction_numerical': 0.0083282701797543, 'mean_word_length': 3.341912237807232, 'pattern_counts': {'":': 3, '<': 11, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 32}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7559626', 'n_tokens_mistral': 21094, 'n_tokens_neox': 20477, 'n_words': 4089}