HistogramFreeMUCA/dataSet_methods.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
import numpy as np
import sympy
from plotMethods import *


def construct_basisSet(maxBasisTerms, ymax, ymin):
    '''
    Constructs symbolic expressions and lambda functions of the basis set and corresponding derivative.

    Parameters
    ----------
    maxBasisTerms :
        - maximum number of terms allowed in basis set expansion

    ymax, ymin :
        - min. and max. of f(x) within the integration domain

    Returns
    -------
    lncBasis_sym :
        - a list that stores the basis set for ln(c) in SymPy symbolic expressions
    lncBasis :
        - a list that stores the basis set for ln(c) in SymPy lambda functions
    lngBasis_sym :
        - a list that stores the basis set for ln(g) in SymPy symbolic expressions
    lngBasis :
        - a list that stores the basis set for ln(g) in SymPy lambda functions
    '''

    i  = 0
    q  = sympy.Symbol('q')                 # Define 'q' as a symbol, not a parameter
    lncBasis  = []
    lngBasis = []

    while (i < maxBasisTerms):
        lncBasis.append("sin(%f*%f*((q-%f)/(%f-%f)))" %(i, np.pi, ymin, ymax, ymin))
        i = i + 1

    lncBasis_sym = sympy.sympify(lncBasis)
    lncBasis = sympy.lambdify(q,lncBasis_sym)    # Applies expression 'lncBasis_sym' onto argument 'q'
    # ln(g) basis set - derivative of ln(c) basis
    i = 0
    while (i < maxBasisTerms):
        lngBasis.append(str(sympy.diff(lncBasis_sym[i], q)))
        i = i + 1
    lngBasis_sym = sympy.sympify(lngBasis)
    lngBasis = sympy.lambdify(q, lngBasis_sym)

    return lncBasis_sym, lncBasis, lngBasis_sym, lngBasis


def fill_dataSet(function, lngCoeff, lngBasis, points_per_dataSet, xmin, xmax, ymin, ymax, plotCDF):
    '''
    Generates data set; calculates the cumulative distribution function (CDF) and reminder.

    Parameters
    ----------
    function :
        - symbolic expression of the integrand, f(x)
    lngCoeff :
        - coefficients of the basis terms for the natural log of the instantaneous density of states, g
    lngBasis :
        - natural log of the instantaneous density of states, g
        - it acts as the sampling weight for determining acceptance
    points_per_dataSet :
        - number of data points in the data set
        - determines the upper bound of the data array
    xmin, xmax :
        - integration domain
    ymin, ymax :
        - min. and max. of f(x) within the integration domain
    plotCDF :
        - a flag to plot a graph of the CDF; 0 = off, 1 = on

    Returns
    -------
    data :
        - 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)
    '''

    data = np.zeros((points_per_dataSet, 3))

    # Create Lambda function from input function
    x = sympy.Symbol("x")
    f_sym = sympy.sympify(function)
    f = sympy.lambdify(x, f_sym)

    # Fill data set: Accept or reject data points based on acceptance criterion
    j = 0

    rand_xinit = np.random.uniform(xmin, xmax)
    yold = f(rand_xinit)
    while (j < points_per_dataSet):
        rand_x = np.random.uniform(xmin, xmax)
        y = f(rand_x)
        g_new = np.exp(np.dot(np.asarray(lngCoeff),np.asarray(lngBasis(y))))
        g_old = np.exp(np.dot(np.asarray(lngCoeff),np.asarray(lngBasis(yold))))

        if (g_new < g_old):
            data[j,0] = y
            yold = y
        elif ((g_old/g_new) > np.random.uniform(0,1)):
            data[j,0] = y
            yold = y
        else:
            data[j,0] = yold
        j = j + 1

    # Sort data set
    data[:,0] = np.sort(data[:,0])

    # Compute CDF of data set
    j = 0
    while (j < points_per_dataSet):
            data[j,1] = float ((j+1.0)/points_per_dataSet)
            j = j + 1

    # Compute remainder
    j = 0
    temp = np.zeros((points_per_dataSet))
    while (j < points_per_dataSet):
        data[j,2] = (data[j,1] - ((data[j,0] - data[0,0]) /
                                      (data[points_per_dataSet - 1,0] - data[0,0])))
        temp[j] = ((data[j,0] - data[0,0]) /
                                      (data[points_per_dataSet - 1,0] - data[0,0]))
        j = j + 1

    # Plot ECDF, F_0, and empirical remainder
    if (plotCDF == 1):
        show_CDF_Plot(data,temp)

    return data