<commit_before>import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from dynaphopy.analysis.fitting import fitting_functions

#from scipy.optimize import curve_fit, minimize_scalar, root

h_planck = 4.135667662e-3  # eV/ps
h_planck_bar = 6.58211951e-4  # eV/ps
kb_boltzmann = 8.6173324e-5  # eV/K


def degenerate_sets(freqs, cutoff=1e-4):
    indices = []
    done = []
    for i in range(len(freqs)):
        if i in done:
            continue
        else:
            f_set = [i]
            done.append(i)
        for j in range(i + 1, len(freqs)):
            if (np.abs(freqs[f_set] - freqs[j]) < cutoff).any():
                f_set.append(j)
                done.append(j)
        indices.append(f_set[:])

    return indices


def average_phonon(index, data, degeneracy):
    for i, set in enumerate(degeneracy):
        if index in set:
            return np.average([data[:, j] for j in degeneracy[i]], axis=0)


def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequencies=None,
                            fitting_function_type=0,
                            show_plots=True,
                            use_degeneracy=True,
                            show_occupancy=True):

    widths = []
    positions = []
    errors = []
    dt_Q2_s = []

    for i in range(original.shape[1]):

        if use_degeneracy:
            degeneracy = degenerate_sets(harmonic_frequencies)
            power_spectrum = average_phonon(i, original, degeneracy)

        else:
            power_spectrum = original[:, i]

        guess_height = np.max(power_spectrum)
        guess_position = test_frequencies_range[np.argmax(power_spectrum)]

        Fitting_function_class = fitting_functions.fitting_functions[fitting_function_type]
        fitting_function = Fitting_function_class(test_frequencies_range,
                                                  power_spectrum,
                                                  guess_height=guess_height,
                                                  guess_position=guess_position)

        fitting_parameters = fitting_function.get_fitting()

        if not fitting_parameters['all_good']:
            positions.append(0)
            widths.append(0)
            errors.append(0)
            dt_Q2_s.append(0)
            continue

        frequency = fitting_parameters['peak_position']
        area = fitting_parameters['area']
        width = fitting_parameters['width']
        base_line = fitting_parameters['base_line']
        maximum = fitting_parameters['maximum']
        error = fitting_parameters['average_relative_error']

        total_integral = simps(power_spectrum, x=test_frequencies_range)/ (2 * np.pi)

        # Calculated properties
        dt_Q2_lor = 2 * 2 * area
        dt_Q2_tot = 2 * 2 * total_integral

        # Only within harmonic approximation
        Q2_lor = dt_Q2_lor / pow(frequency * 2 * np.pi, 2)
        Q2_tot = dt_Q2_tot / pow(frequency * 2 * np.pi,2)

        occupancy_lor = dt_Q2_lor / (frequency * h_planck_bar) - 0.5
        occupancy_tot = dt_Q2_tot / (frequency * h_planck_bar) - 0.5

    #    fit_temperature = dt_Q2_lor / kb_boltzmann # High temperature limit
        fit_temperature = h_planck_bar * frequency / (kb_boltzmann * np.log((1.0 / occupancy_lor + 1.0)))

        #Print section
        print ('\nPeak # {0}'.format(i+1))
        print ('----------------------------------------------')
        print ('Width (FWHM)               {0:15.6f} THz'.format(width))
        print ('Position                   {0:15.6f} THz'.format(frequency))
        print ('Area (1/2<K>) ({0:.10s}) {1:15.6f} eV'.format(fitting_function.curve_name, area))      # 1/2 Kinetic energy
        print ('Area (1/2<K>) (Total)      {0:15.6f} eV'.format(total_integral))   # 1/2 Kinetic energy
        print ('<|dQ/dt|^2>                {0:15.6f} eV'.format(dt_Q2_lor))        # Kinetic energy
 #       print '<|dQ/dt|^2> (tot):        ', dt_Q2_tot, 'eV'        # Kinetic energy
 #       print '<|Q|^2> (lor):          ', Q2_lor, 'u * Angstrom^2'
 #       print '<|Q|^2> (tot):          ', Q2_tot, 'u * Angstrom^2'
        if show_occupancy:
            print ('Occupation number          {0:15.6f}'.format(occupancy_lor))
     #       print 'Occupation number(tot): ', occupancy_tot
            print ('Fit temperature            {0:15.6f} K'.format(fit_temperature))
     #       print 'Fit temperature (tot)   ', dt_Q2_tot / kb_bolzman, 'K'
        print ('Base line                  {0:15.6f} eV * ps'.format(base_line))
        print ('Maximum height             {0:15.6f} eV * ps'.format(maximum))
        print ('Fit relative error         {0:15.6f}'.format(error))

        if 'asymmetry' in fitting_parameters:<commit_msg>Added standard errors in fitting log<commit_after>import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from dynaphopy.analysis.fitting import fitting_functions

#from scipy.optimize import curve_fit, minimize_scalar, root

h_planck = 4.135667662e-3  # eV/ps
h_planck_bar = 6.58211951e-4  # eV/ps
kb_boltzmann = 8.6173324e-5  # eV/K


def degenerate_sets(freqs, cutoff=1e-4):
    indices = []
    done = []
    for i in range(len(freqs)):
        if i in done:
            continue
        else:
            f_set = [i]
            done.append(i)
        for j in range(i + 1, len(freqs)):
            if (np.abs(freqs[f_set] - freqs[j]) < cutoff).any():
                f_set.append(j)
                done.append(j)
        indices.append(f_set[:])

    return indices


def average_phonon(index, data, degeneracy):
    for i, set in enumerate(degeneracy):
        if index in set:
            return np.average([data[:, j] for j in degeneracy[i]], axis=0)


def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequencies=None,
                            fitting_function_type=0,
                            show_plots=True,
                            use_degeneracy=True,
                            show_occupancy=True):

    widths = []
    positions = []
    errors = []
    dt_Q2_s = []

    for i in range(original.shape[1]):

        if use_degeneracy:
            degeneracy = degenerate_sets(harmonic_frequencies)
            power_spectrum = average_phonon(i, original, degeneracy)

        else:
            power_spectrum = original[:, i]

        guess_height = np.max(power_spectrum)
        guess_position = test_frequencies_range[np.argmax(power_spectrum)]

        Fitting_function_class = fitting_functions.fitting_functions[fitting_function_type]
        fitting_function = Fitting_function_class(test_frequencies_range,
                                                  power_spectrum,
                                                  guess_height=guess_height,
                                                  guess_position=guess_position)

        fitting_parameters = fitting_function.get_fitting()

        if not fitting_parameters['all_good']:
            positions.append(0)
            widths.append(0)
            errors.append(0)
            dt_Q2_s.append(0)
            continue

        frequency = fitting_parameters['peak_position']
        area = fitting_parameters['area']
        width = fitting_parameters['width']
        base_line = fitting_parameters['base_line']
        maximum = fitting_parameters['maximum']
        error = fitting_parameters['average_relative_error']
        standard_errors = fitting_parameters['standard_errors']
        total_integral = simps(power_spectrum, x=test_frequencies_range)/ (2 * np.pi)

        # Calculated properties
        dt_Q2_lor = 2 * 2 * area
        dt_Q2_tot = 2 * 2 * total_integral

        # Only within harmonic approximation
        Q2_lor = dt_Q2_lor / pow(frequency * 2 * np.pi, 2)
        Q2_tot = dt_Q2_tot / pow(frequency * 2 * np.pi,2)

        occupancy_lor = dt_Q2_lor / (frequency * h_planck_bar) - 0.5
        occupancy_tot = dt_Q2_tot / (frequency * h_planck_bar) - 0.5

    #    fit_temperature = dt_Q2_lor / kb_boltzmann # High temperature limit
        fit_temperature = h_planck_bar * frequency / (kb_boltzmann * np.log((1.0 / occupancy_lor + 1.0)))

        #Print section
        print ('\nPeak # {0}'.format(i+1))
        print ('----------------------------------------------')
        print ('Width (FWHM)               {0:15.6f} +- {1:8.3e} THz'.format(width, standard_errors[0]))
        print ('Position                   {0:15.6f} +- {1:8.3e} THz'.format(frequency, standard_errors[1]))
        print ('Area (1/2<K>) ({0:.10s}) {1:15.6f} eV'.format(fitting_function.curve_name, area))      # 1/2 Kinetic energy
        print ('Area (1/2<K>) (Total)      {0:15.6f} eV'.format(total_integral))   # 1/2 Kinetic energy
        print ('<|dQ/dt|^2>                {0:15.6f} eV'.format(dt_Q2_lor))        # Kinetic energy
 #       print '<|dQ/dt|^2> (tot):        ', dt_Q2_tot, 'eV'        # Kinetic energy
 #       print '<|Q|^2> (lor):          ', Q2_lor, 'u * Angstrom^2'
 #       print '<|Q|^2> (tot):          ', Q2_tot, 'u * Angstrom^2'
        if show_occupancy:
            print ('Occupation number          {0:15.6f}'.format(occupancy_lor))
     #       print 'Occupation number(tot): ', occupancy_tot
            print ('Fit temperature            {0:15.6f} K'.format(fit_temperature))
     #       print 'Fit temperature (tot)   ', dt_Q2_tot / kb_bolzman, 'K'
        print ('Base line                  {0:15.6f} eV * ps'.format(base_line))
        print ('Maximum height             {0:15.6f} eV * ps'.format(maximum))
        print ('Fit relative error         {0:15.6f}'.format(error))

        if 'asymmetry' in fitting_parameters:
{'fraction_non_alphanumeric': 0.12112763320941759, 'fraction_numerical': 0.02808756712102437, 'mean_word_length': 1.9928924598269468, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '960188', 'n_tokens_mistral': 3225, 'n_tokens_neox': 2961, 'n_words': 612}