HistogramFreeMUCA/correction.py at master · yingwaili/HistogramFreeMUCA · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
import numpy as np
import sys
from KS import *
from plotMethods import *

def correction(dataSet, ymin, ymax, KS_cutoff, numTerms, lncBasis, plotFT):
    '''
    Calculates the coefficents that make up the log of the correction expression.

    Parameters
    ----------
    dataSet :
        - the data set generated by the function "fill_dataSets"
        - an array with the following indexing scheme: [data point number, data values]
          data values as follows:  [x,0]: sorted data set
                                   [x,1]: CDF
                                   [x,2]: remainder (CDF - straight line)
    ymin, ymax :
        - min. and max. of f(x) within the integration domain
    KS_cutoff :
        - the cutoff value for the Kolmogorov-Smirnov (KS) test, above which the data are considered to be well fitted
    numTerms :
        - number of terms in the basis expansion
    lncBasis :
        - a list that stores the basis set for ln(c) in SymPy lambda functions
    plotFT :
        - a flag to plot and compare the empirical remainder and Rnalytic remainder

    Returns
    -------
    lncCoeff :
        - an array that stores the Fourier coefficients for the data set
          (up to the term where KS test reaches KS_cutoff; afterwards the coefficients are zero)
    '''

    i = 0

    lncCoeff = np.zeros(numTerms)

    # track convergence of KS test (0: not converged; 1: converged)
    converged = 0

    # calculate Fourier coefficients for the log of correction (lncCoeff) up until KS Test >= 0.5
    while (i < numTerms):
        if (converged == 0):
            lncCoeff[i] = fourierCoeff(dataSet, i, ymin, ymax)
        stat = two_sided_KS_test(dataSet, lncBasis, lncCoeff, ymin, ymax)

        if (stat > KS_cutoff):
            converged = 1
            # plot remainder and empirical remainder
            if (plotFT == 1):
                show_FT_Plot(dataSet, lncCoeff, lncBasis, ymin, ymax)
            i = numTerms # Exit loop

        if ((i == (numTerms - 1)) and (converged == 0)):
            print("Error: Kolmogorov-Smirnov Probability Not Satisfied Given Max. Number of Terms")
            sys.exit()
        i = i + 1

    return lncCoeff


def correction_improved(dataSet, ymin, ymax, KS_cutoff, numTerms, lncBasis, s, plotFT):

    '''
    - Implements the improved scheme for constructing the log of the correction expression
    - Same parameters as in the function "correction"
    - Returns coefficients (multiplied by a small factor s)
    '''

    lncCoeff = np.zeros(numTerms)

    # convergence of KS test
    converged = 0

    # randomly pick starting coefficient for Fourier transform
    i_init = np.random.random_integers(0, (numTerms-1))
    i = i_init

    # calculate Fourier coefficients for the log of correction (lncCoeff) up until KS Test >= 0.5
    while (i < (numTerms + i_init)):
        # "Periodic Boundary" over array
        term = ((i + numTerms) % numTerms)

        if (converged == 0):
            lncCoeff[term] = fourierCoeff(dataSet, term, ymin, ymax)
        stat = two_sided_KS_test(dataSet, lncBasis, lncCoeff, ymin, ymax)

        if (stat > KS_cutoff):
            converged = 1
            i = numTerms + i_init # Exit loop
            # plot remainder and empirical remainder
            if (plotFT == 1):
                show_FT_Plot(dataSet, lncCoeff, lncBasis, ymin, ymax)
        i = i + 1

    return (s * lncCoeff)


def fourierCoeff(dataSet, termIndex, ymin, ymax):

    ''' Calculate the Fourier coefficients '''

    N = np.size(dataSet[:,0])
    i = 0
    C = 0.0
    yRange = ymax - ymin
    if (termIndex == 0):
        C = 0.0
    else:
        while (i < N-1):
            y1 = (dataSet[i,0] - ymin) / (yRange)
            y2 = (dataSet[i+1,0] - ymin) / (yRange)
            C += (dataSet[i,1] - y1) * np.cos(termIndex * np.pi * y1) / (termIndex * np.pi)
            C -= (dataSet[i,1] - y2) * np.cos(termIndex * np.pi * y2) / (termIndex * np.pi)
            C += np.sin(termIndex * np.pi * y1) / (termIndex * termIndex * np.pi * np.pi)
            C -= np.sin(termIndex * np.pi * y2) / (termIndex * termIndex * np.pi * np.pi)
            i = i + 1
        C = C * 2.0
    return (C)