Skip to content

Models

The models in SAMBA are two expansions of a toy model, where the full toy model is given by

$$ F(g) = \int_{-\infty}^{\infty} dx~ e^{-\frac{x^{2}}{2} - g^{2} x^{4}} = \frac{e^{\frac{1}{32 g^{2}}}}{2 \sqrt{2}g} K_{\frac{1}{4}}\left(\frac{1}{32 g^{2}} \right). $$

The two expansions are limits taken at $g = 0$ and $g = \infty$:

$$ F_{s}^{N_s}(g) = \sum_{k=0}^{N_{s}} s_{k} g^{k}, $$

and

$$ F_{l}^{N_{l}}(g) = \frac{1}{\sqrt{g}} \sum_{k=0}^{N_{l}} l_{k} g^{-k}, $$

with coefficients given as:

$$ s_{2k} = \frac{\sqrt{2} \Gamma{(2k + 1/2)}}{k!} (-4)^{k},~~~~~s_{2k + 1} = 0 $$

and

$$ l_{k} = \frac{\Gamma{\left(\frac{k}{2} + \frac{1}{4}\right)}}{2k!} \left(-\frac{1}{2}\right)^{k}. $$

We also include models for the uncertainties of each expansion, given in the uninformative limit, for the small-$g$ expansion, as

$$ \sigma_{N_s}(g)= \Gamma(N_s+3) g^{N_s + 2} \bar{c}, $$

if $N_s$ is even, and

$$ \sigma_{N_s}(g)= \Gamma(N_s+2) g^{N_s+1} \bar{c}, $$

if $N_s$ is odd. For the large-$g$ limit,

$$ \sigma_{N_l}(g)=\frac{1}{\Gamma(N_l+2)} \frac{1}{g^{N_l+3/2}} \bar{d}. $$

We also devise expressions for the informative limit, for the small-$g$ expansion, as

$$ \sigma_{N_s}(g)= \Gamma(N_s/2+1) (4g)^{N_s + 2} \bar{c}, $$

if $N_s$ is even, and

$$ \sigma_{N_s}(g)= \Gamma(N_s/2+1/2) (4g)^{N_s+1} \bar{c}, $$

if $N_s$ is odd. For the large-$g$ limit,

$$ \sigma_{N_l}(g)=\left(\frac{1}{4g}\right)^{N_l + 3/2} \frac{1}{\Gamma(N_l/2+3/2)} \bar{d}. $$

Models

Source code in samba/models.py 
 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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
class Models():

    def __init__(self, loworder, highorder):

        r'''
        The class containing the expansion models from Honda's paper
        and the means to plot them. 

        Example:
            Models(loworder=np.array([2]), highorder=np.array([5]))

        Parameters:
            loworder (numpy.ndarray, int, float): The value of N_s to be used 
                to truncate the small-g expansion. Can be an array of multiple 
                values or one. 

            highorder (numpy.ndarray, int, float): The value of N_l to be used 
                to truncate the large-g expansion. Can be an array of multiple 
                values or one. 

        Returns:
            None.
        '''

        #check type and assign to class variable
        if isinstance(loworder, float) == True or isinstance(loworder, int) == True:
            loworder = np.array([loworder])

        #check type and assign to class variable
        if isinstance(highorder, float) == True or isinstance(highorder, int) == True:
            highorder = np.array([highorder])

        self.loworder = loworder
        self.highorder = highorder 

        return None 


    def low_g(self, g):

        '''
        A function to calculate the small-g divergent asymptotic expansion for a 
        given range in the coupling constant, g.

        Example:
            Models.low_g(g=np.linspace(0.0, 0.5, 20))

        Parameters:
            g (linspace): The linspace of the coupling constant for 
                this calculation. 

        Returns:
            output (numpy.ndarray): The array of values of the expansion in 
                small-g at each point in g_true space, for each value of loworder 
                (highest power the expansion reaches).
        '''

        output = []

        for order in self.loworder:
            low_c = np.empty([int(order)+1])
            low_terms = np.empty([int(order) + 1])

            #if g is an array, execute here
            try:
                value = np.empty([len(g)])

                #loop over array in g
                for i in range(len(g)):      

                    #loop over orders
                    for k in range(int(order)+1):

                        if k % 2 == 0:
                            low_c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2))
                        else:
                            low_c[k] = 0

                        low_terms[k] = low_c[k] * g[i]**(k) 

                    value[i] = np.sum(low_terms)

                output.append(value)
                data = np.array(output, dtype = np.float64)

            #if g is a single value, execute here
            except:
                value = 0.0
                for k in range(int(order)+1):

                    if k % 2 == 0:
                        low_c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2))
                    else:
                        low_c[k] = 0

                    low_terms[k] = low_c[k] * g**(k) 

                value = np.sum(low_terms)
                data = value

        return data


    def high_g(self, g):

        r'''
        A function to calculate the large-g convergent Taylor expansion for a given range in the coupling 
        constant, g.

        Example:
            Models.high_g(g=np.linspace(1e-6,1.0,100))

        Parameters:
            g (linspace): The linspace of the coupling constant for this calculation.

        Returns:
            output (numpy.ndarray): The array of values of the expansion at large-g at 
                each point in g_true space, for each value of highorder (highest power 
                the expansion reaches).
        '''

        output = []

        for order in self.highorder:
            high_c = np.empty([int(order) + 1])
            high_terms = np.empty([int(order) + 1])

            #if g is an array, execute here
            try:
                value = np.empty([len(g)])

                #loop over array in g
                for i in range(len(g)):

                    #loop over orders
                    for k in range(int(order)+1):

                        high_c[k] = special.gamma(k/2.0 + 0.25) * (-0.5)**k / (2.0 * math.factorial(k))

                        high_terms[k] = (high_c[k] * g[i]**(-k)) / np.sqrt(g[i])

                    #sum the terms for each value of g
                    value[i] = np.sum(high_terms)

                output.append(value)

                data = np.array(output, dtype = np.float64)

            #if g is a single value, execute here           
            except:
                value = 0.0

                #loop over orders
                for k in range(int(order)+1):

                    high_c[k] = special.gamma(k/2.0 + 0.25) * (-0.5)**k / (2.0 * math.factorial(k))

                    high_terms[k] = (high_c[k] * g**(-k)) / np.sqrt(g) 

                #sum the terms for each value of g
                value = np.sum(high_terms)
                data = value

        return data 


    def true_model(self, g):

        r'''
        The true model of the zero-dimensional phi^4 theory partition function using an input 
        linspace.

        Example:
            Models.true_model(g=np.linspace(0.0, 0.5, 100))

        Parameters:
            g (linspace): The linspace for g desired to calculate the true model. This can be 
                the g_true linspace, g_data linspace, or another linspace of the user's choosing. 

        Returns:
            model (numpy.ndarray): The model calculated at each point in g space. 
        '''

        #define a function for the integrand
        def function(x,g):
            return np.exp(-(x**2.0)/2.0 - (g**2.0 * x**4.0)) 

        #initialization
        self.model = np.zeros([len(g)])

        #perform the integral for each g
        for i in range(len(g)):

            self.model[i], self.err = integrate.quad(function, -np.inf, np.inf, args=(g[i],))

        return self.model 


    def plot_models(self, g, only_true=False):

        r'''
        A plotting function to produce a figure of the model expansions calculated in Models.low_g and Models.high_g, 
        and including the true model calculated using Mixing.true_model.

        Example:
            Mixing.plot_models(g=np.linspace(0.0, 0.5, 100))

        Parameters:
            g (linspace): The linspace in on which the models will be plotted here. 

        Returns:
            None.
        '''

        #uncertainties
        unc = Uncertainties()

        #set up the plot
        fig = plt.figure(figsize=(8,6), dpi=600)
        ax = plt.axes()
        ax.tick_params(axis='x', labelsize=18)
        ax.tick_params(axis='y', labelsize=18)
        ax.locator_params(nbins=8)
        ax.xaxis.set_minor_locator(AutoMinorLocator())
        ax.yaxis.set_minor_locator(AutoMinorLocator())
        ax.set_xlabel('g', fontsize=22)
        ax.set_ylabel('F(g)', fontsize=22)
        ax.set_facecolor("white")

        #x and y limits
        if min(g) == 1e-6:
            ax.set_xlim(0.0, max(g))
        else:
            ax.set_xlim(min(g), max(g))
        ax.set_ylim(1.0,3.0)
      #  ax.set_ylim(1.2,3.2)
      #  ax.set_yticks([1.2, 1.6, 2.0, 2.4, 2.8, 3.2])

        #title option
        ax.set_title('F(g): expansions and true model', fontsize=22)

        #plot the true model 
        ax.plot(g, self.true_model(g), 'k', label='True model')

        if only_true is False:

            #add linestyle cycler
            linestyle_cycler = cycler(linestyle=['dashed'])
            ax.set_prop_cycle(linestyle_cycler)

            #for each large-g order, calculate and plot
            for i,j in zip(range(len(self.loworder)), self.loworder):
                ax.plot(g, self.low_g(g)[i,:], color='r', label=r'$f_s$ ($N_s$ = {})'.format(j))

            #for each large-g order, calculate and plot
            for i,j in zip(range(len(self.highorder)), self.highorder):
                ax.plot(g, self.high_g(g)[i,:], color='b', label=r'$f_l$ ($N_l$ = {})'.format(j))

            #calculate uncertainty bands
            std_low = []
            std_high = []

            for i in (self.loworder):
                std_low.append(np.sqrt(unc.variance_low(g,i)))

            for i in (self.highorder):
                std_high.append(np.sqrt(unc.variance_high(g,i)))

            #plot uncertainties as fill-between
            for i,j in zip(range(len(self.loworder)), self.loworder):
                ax.plot(g, self.low_g(g)[i,:]-std_low[i], color='red', linestyle='dotted')
                ax.plot(g, self.low_g(g)[i,:]+std_low[i], color='red', linestyle='dotted')
#                 ax.fill_between(g, self.low_g(g)[i,:]-std_low[i], self.low_g(g)[i,:]+std_low[i],
#                             zorder=i-5, facecolor='red', alpha=0.1, lw=0.6, label=r'$f_s$ ($N_s$ = {}) 68\% CI'.format(j))

            for i,j in zip(range(len(self.highorder)), self.highorder):
                ax.plot(g, self.high_g(g)[i,:]-std_high[i], color='blue', linestyle='dotted')
                ax.plot(g, self.high_g(g)[i,:]+std_high[i], color='blue', linestyle='dotted')
#                 ax.fill_between(g, self.high_g(g)[i,:]-std_high[i], self.high_g(g)[i,:]+std_high[i],
#                                 zorder=i-5, facecolor='blue', alpha=0.1, lw=0.6, label=r'$f_l$ ($N_l$ = {}) 68\% CI'.format(j))

        ax.legend(fontsize=16, loc='upper right')
        plt.show()

        return fig


    def residuals(self):

        r'''
        A calculation and plot of the residuals of the model expansions vs the true model values at each point in g.
        g is set internally for this plot, as the plot must be shown in loglog format to see the power law of the
        residuals. 

        Example:
            Mixing.residuals()

        Parameters:
            None.

        Returns:
            None. 
        '''

        #set up the plot
        fig = plt.figure(figsize=(8,6), dpi=600)
        ax = plt.axes()
        ax.tick_params(axis='x', labelsize=18)
        ax.tick_params(axis='y', labelsize=18)
        ax.xaxis.set_minor_locator(AutoMinorLocator())
        ax.yaxis.set_minor_locator(AutoMinorLocator())
        ax.set_xlabel('g', fontsize=22)
        ax.set_ylabel('Residual', fontsize=22)
        ax.set_title('F(g): residuals', fontsize=22)
        ax.set_xlim(1e-2, 10.)
        ax.set_ylim(1e-6, 1e17)

        #set range for g
        g_ext = np.logspace(-6., 6., 800)

        #set up marker cycler
        marker_cycler = cycler(marker=['.', '*', '+', '.', '*', '+'])
        ax.set_prop_cycle(marker_cycler)

        #calculate true model
        value_true = self.true_model(g_ext)

        #for each small-g order, plot
        valuelow = np.zeros([len(self.loworder), len(g_ext)])
        residlow = np.zeros([len(self.loworder), len(g_ext)])

        for i,j in zip(range(len(self.loworder)), self.loworder):
            valuelow[i,:] = self.low_g(g_ext)[i]
            residlow[i,:] = (valuelow[i,:] - value_true)/value_true
            ax.loglog(g_ext, abs(residlow[i,:]), 'r', linestyle="None", label=r"$F_s({})$".format(j))

        #for each large-g order, plot
        valuehi = np.zeros([len(self.highorder), len(g_ext)])
        residhi = np.zeros([len(self.highorder), len(g_ext)])

        for i,j in zip(range(len(self.highorder)), self.highorder):
            valuehi[i,:] = self.high_g(g_ext)[i]
            residhi[i,:] = (valuehi[i,:] - value_true)/value_true
            ax.loglog(g_ext, abs(residhi[i,:]), 'b', linestyle="None", label=r"$F_l({})$".format(j))

        ax.legend(fontsize=18)
        plt.show()

__init__(loworder, highorder)

The class containing the expansion models from Honda's paper and the means to plot them.

Example

Models(loworder=np.array([2]), highorder=np.array([5]))

Parameters:

Name Type Description Default
loworder (ndarray, int, float)

The value of N_s to be used to truncate the small-g expansion. Can be an array of multiple values or one.

required
highorder (ndarray, int, float)

The value of N_l to be used to truncate the large-g expansion. Can be an array of multiple values or one.

required

Returns:

Type Description

None.
Source code in samba/models.py 
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
def __init__(self, loworder, highorder):

    r'''
    The class containing the expansion models from Honda's paper
    and the means to plot them. 

    Example:
        Models(loworder=np.array([2]), highorder=np.array([5]))

    Parameters:
        loworder (numpy.ndarray, int, float): The value of N_s to be used 
            to truncate the small-g expansion. Can be an array of multiple 
            values or one. 

        highorder (numpy.ndarray, int, float): The value of N_l to be used 
            to truncate the large-g expansion. Can be an array of multiple 
            values or one. 

    Returns:
        None.
    '''

    #check type and assign to class variable
    if isinstance(loworder, float) == True or isinstance(loworder, int) == True:
        loworder = np.array([loworder])

    #check type and assign to class variable
    if isinstance(highorder, float) == True or isinstance(highorder, int) == True:
        highorder = np.array([highorder])

    self.loworder = loworder
    self.highorder = highorder 

    return None

high_g(g)

A function to calculate the large-g convergent Taylor expansion for a given range in the coupling constant, g.

Example

Models.high_g(g=np.linspace(1e-6,1.0,100))

Parameters:

Name Type Description Default
g linspace

The linspace of the coupling constant for this calculation.

 required

Returns:

Name Type Description
output ndarray

The array of values of the expansion at large-g at each point in g_true space, for each value of highorder (highest power the expansion reaches).
Source code in samba/models.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
def high_g(self, g):

    r'''
    A function to calculate the large-g convergent Taylor expansion for a given range in the coupling 
    constant, g.

    Example:
        Models.high_g(g=np.linspace(1e-6,1.0,100))

    Parameters:
        g (linspace): The linspace of the coupling constant for this calculation.

    Returns:
        output (numpy.ndarray): The array of values of the expansion at large-g at 
            each point in g_true space, for each value of highorder (highest power 
            the expansion reaches).
    '''

    output = []

    for order in self.highorder:
        high_c = np.empty([int(order) + 1])
        high_terms = np.empty([int(order) + 1])

        #if g is an array, execute here
        try:
            value = np.empty([len(g)])

            #loop over array in g
            for i in range(len(g)):

                #loop over orders
                for k in range(int(order)+1):

                    high_c[k] = special.gamma(k/2.0 + 0.25) * (-0.5)**k / (2.0 * math.factorial(k))

                    high_terms[k] = (high_c[k] * g[i]**(-k)) / np.sqrt(g[i])

                #sum the terms for each value of g
                value[i] = np.sum(high_terms)

            output.append(value)

            data = np.array(output, dtype = np.float64)

        #if g is a single value, execute here           
        except:
            value = 0.0

            #loop over orders
            for k in range(int(order)+1):

                high_c[k] = special.gamma(k/2.0 + 0.25) * (-0.5)**k / (2.0 * math.factorial(k))

                high_terms[k] = (high_c[k] * g**(-k)) / np.sqrt(g) 

            #sum the terms for each value of g
            value = np.sum(high_terms)
            data = value

    return data

low_g(g)

A function to calculate the small-g divergent asymptotic expansion for a given range in the coupling constant, g.

Example

Models.low_g(g=np.linspace(0.0, 0.5, 20))

Parameters:

Name Type Description Default
g linspace

The linspace of the coupling constant for this calculation.

required

Returns:

Name Type Description
output ndarray

The array of values of the expansion in small-g at each point in g_true space, for each value of loworder (highest power the expansion reaches).
Source code in samba/models.py 
 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
def low_g(self, g):

    '''
    A function to calculate the small-g divergent asymptotic expansion for a 
    given range in the coupling constant, g.

    Example:
        Models.low_g(g=np.linspace(0.0, 0.5, 20))

    Parameters:
        g (linspace): The linspace of the coupling constant for 
            this calculation. 

    Returns:
        output (numpy.ndarray): The array of values of the expansion in 
            small-g at each point in g_true space, for each value of loworder 
            (highest power the expansion reaches).
    '''

    output = []

    for order in self.loworder:
        low_c = np.empty([int(order)+1])
        low_terms = np.empty([int(order) + 1])

        #if g is an array, execute here
        try:
            value = np.empty([len(g)])

            #loop over array in g
            for i in range(len(g)):      

                #loop over orders
                for k in range(int(order)+1):

                    if k % 2 == 0:
                        low_c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2))
                    else:
                        low_c[k] = 0

                    low_terms[k] = low_c[k] * g[i]**(k) 

                value[i] = np.sum(low_terms)

            output.append(value)
            data = np.array(output, dtype = np.float64)

        #if g is a single value, execute here
        except:
            value = 0.0
            for k in range(int(order)+1):

                if k % 2 == 0:
                    low_c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2))
                else:
                    low_c[k] = 0

                low_terms[k] = low_c[k] * g**(k) 

            value = np.sum(low_terms)
            data = value

    return data

plot_models(g, only_true=False)

A plotting function to produce a figure of the model expansions calculated in Models.low_g and Models.high_g, and including the true model calculated using Mixing.true_model.

Example

Mixing.plot_models(g=np.linspace(0.0, 0.5, 100))

Parameters:

Name Type Description Default
g linspace

The linspace in on which the models will be plotted here.

required

Returns:

Type Description

None.
Source code in samba/models.py 
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
    def plot_models(self, g, only_true=False):

        r'''
        A plotting function to produce a figure of the model expansions calculated in Models.low_g and Models.high_g, 
        and including the true model calculated using Mixing.true_model.

        Example:
            Mixing.plot_models(g=np.linspace(0.0, 0.5, 100))

        Parameters:
            g (linspace): The linspace in on which the models will be plotted here. 

        Returns:
            None.
        '''

        #uncertainties
        unc = Uncertainties()

        #set up the plot
        fig = plt.figure(figsize=(8,6), dpi=600)
        ax = plt.axes()
        ax.tick_params(axis='x', labelsize=18)
        ax.tick_params(axis='y', labelsize=18)
        ax.locator_params(nbins=8)
        ax.xaxis.set_minor_locator(AutoMinorLocator())
        ax.yaxis.set_minor_locator(AutoMinorLocator())
        ax.set_xlabel('g', fontsize=22)
        ax.set_ylabel('F(g)', fontsize=22)
        ax.set_facecolor("white")

        #x and y limits
        if min(g) == 1e-6:
            ax.set_xlim(0.0, max(g))
        else:
            ax.set_xlim(min(g), max(g))
        ax.set_ylim(1.0,3.0)
      #  ax.set_ylim(1.2,3.2)
      #  ax.set_yticks([1.2, 1.6, 2.0, 2.4, 2.8, 3.2])

        #title option
        ax.set_title('F(g): expansions and true model', fontsize=22)

        #plot the true model 
        ax.plot(g, self.true_model(g), 'k', label='True model')

        if only_true is False:

            #add linestyle cycler
            linestyle_cycler = cycler(linestyle=['dashed'])
            ax.set_prop_cycle(linestyle_cycler)

            #for each large-g order, calculate and plot
            for i,j in zip(range(len(self.loworder)), self.loworder):
                ax.plot(g, self.low_g(g)[i,:], color='r', label=r'$f_s$ ($N_s$ = {})'.format(j))

            #for each large-g order, calculate and plot
            for i,j in zip(range(len(self.highorder)), self.highorder):
                ax.plot(g, self.high_g(g)[i,:], color='b', label=r'$f_l$ ($N_l$ = {})'.format(j))

            #calculate uncertainty bands
            std_low = []
            std_high = []

            for i in (self.loworder):
                std_low.append(np.sqrt(unc.variance_low(g,i)))

            for i in (self.highorder):
                std_high.append(np.sqrt(unc.variance_high(g,i)))

            #plot uncertainties as fill-between
            for i,j in zip(range(len(self.loworder)), self.loworder):
                ax.plot(g, self.low_g(g)[i,:]-std_low[i], color='red', linestyle='dotted')
                ax.plot(g, self.low_g(g)[i,:]+std_low[i], color='red', linestyle='dotted')
#                 ax.fill_between(g, self.low_g(g)[i,:]-std_low[i], self.low_g(g)[i,:]+std_low[i],
#                             zorder=i-5, facecolor='red', alpha=0.1, lw=0.6, label=r'$f_s$ ($N_s$ = {}) 68\% CI'.format(j))

            for i,j in zip(range(len(self.highorder)), self.highorder):
                ax.plot(g, self.high_g(g)[i,:]-std_high[i], color='blue', linestyle='dotted')
                ax.plot(g, self.high_g(g)[i,:]+std_high[i], color='blue', linestyle='dotted')
#                 ax.fill_between(g, self.high_g(g)[i,:]-std_high[i], self.high_g(g)[i,:]+std_high[i],
#                                 zorder=i-5, facecolor='blue', alpha=0.1, lw=0.6, label=r'$f_l$ ($N_l$ = {}) 68\% CI'.format(j))

        ax.legend(fontsize=16, loc='upper right')
        plt.show()

        return fig

residuals()

A calculation and plot of the residuals of the model expansions vs the true model values at each point in g. g is set internally for this plot, as the plot must be shown in loglog format to see the power law of the residuals.

Example

Mixing.residuals()

Returns:

Type Description

None.
Source code in samba/models.py 
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
def residuals(self):

    r'''
    A calculation and plot of the residuals of the model expansions vs the true model values at each point in g.
    g is set internally for this plot, as the plot must be shown in loglog format to see the power law of the
    residuals. 

    Example:
        Mixing.residuals()

    Parameters:
        None.

    Returns:
        None. 
    '''

    #set up the plot
    fig = plt.figure(figsize=(8,6), dpi=600)
    ax = plt.axes()
    ax.tick_params(axis='x', labelsize=18)
    ax.tick_params(axis='y', labelsize=18)
    ax.xaxis.set_minor_locator(AutoMinorLocator())
    ax.yaxis.set_minor_locator(AutoMinorLocator())
    ax.set_xlabel('g', fontsize=22)
    ax.set_ylabel('Residual', fontsize=22)
    ax.set_title('F(g): residuals', fontsize=22)
    ax.set_xlim(1e-2, 10.)
    ax.set_ylim(1e-6, 1e17)

    #set range for g
    g_ext = np.logspace(-6., 6., 800)

    #set up marker cycler
    marker_cycler = cycler(marker=['.', '*', '+', '.', '*', '+'])
    ax.set_prop_cycle(marker_cycler)

    #calculate true model
    value_true = self.true_model(g_ext)

    #for each small-g order, plot
    valuelow = np.zeros([len(self.loworder), len(g_ext)])
    residlow = np.zeros([len(self.loworder), len(g_ext)])

    for i,j in zip(range(len(self.loworder)), self.loworder):
        valuelow[i,:] = self.low_g(g_ext)[i]
        residlow[i,:] = (valuelow[i,:] - value_true)/value_true
        ax.loglog(g_ext, abs(residlow[i,:]), 'r', linestyle="None", label=r"$F_s({})$".format(j))

    #for each large-g order, plot
    valuehi = np.zeros([len(self.highorder), len(g_ext)])
    residhi = np.zeros([len(self.highorder), len(g_ext)])

    for i,j in zip(range(len(self.highorder)), self.highorder):
        valuehi[i,:] = self.high_g(g_ext)[i]
        residhi[i,:] = (valuehi[i,:] - value_true)/value_true
        ax.loglog(g_ext, abs(residhi[i,:]), 'b', linestyle="None", label=r"$F_l({})$".format(j))

    ax.legend(fontsize=18)
    plt.show()

true_model(g)

The true model of the zero-dimensional phi^4 theory partition function using an input linspace.

Example

Models.true_model(g=np.linspace(0.0, 0.5, 100))

Parameters:

Name Type Description Default
g linspace

The linspace for g desired to calculate the true model. This can be the g_true linspace, g_data linspace, or another linspace of the user's choosing.

required

Returns:

Name Type Description
model ndarray

The model calculated at each point in g space.
Source code in samba/models.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
def true_model(self, g):

    r'''
    The true model of the zero-dimensional phi^4 theory partition function using an input 
    linspace.

    Example:
        Models.true_model(g=np.linspace(0.0, 0.5, 100))

    Parameters:
        g (linspace): The linspace for g desired to calculate the true model. This can be 
            the g_true linspace, g_data linspace, or another linspace of the user's choosing. 

    Returns:
        model (numpy.ndarray): The model calculated at each point in g space. 
    '''

    #define a function for the integrand
    def function(x,g):
        return np.exp(-(x**2.0)/2.0 - (g**2.0 * x**4.0)) 

    #initialization
    self.model = np.zeros([len(g)])

    #perform the integral for each g
    for i in range(len(g)):

        self.model[i], self.err = integrate.quad(function, -np.inf, np.inf, args=(g[i],))

    return self.model

Uncertainties

Source code in samba/models.py 
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
class Uncertainties:


    def __init__(self, error_model='informative'):

        r'''
        An accompanying class to Models() that possesses the truncation error models
        that are included as variances with the small-g and large-g expansions. 

        Example:
            Uncertainties()

        Parameters:
            error_model (str): The name of the error model to use in the calculation. 
                Options are 'uninformative' and 'informative'. Default is 'informative'.

        Returns:
            None.
        '''

        #assign error model 
        if error_model == 'uninformative':
           self.error_model = 1

        elif error_model == 'informative':
            self.error_model = 2

        else:
            raise ValueError("Please choose 'uninformative' or 'informative'.")

        return None


    def variance_low(self, g, loworder):


        r'''
        A function to calculate the variance corresponding to the small-g expansion model.

        Example:
            Bivariate.variance_low(g=np.linspace(1e-6, 0.5, 100), loworder=5)

        Parameters:
            g (numpy.linspace): The linspace over which this calculation is performed.

            loworder (int): The order to which we know our expansion model. Must be 
                passed one at a time if more than one model is to be calculated.

        Returns:
            var1 (numpy.ndarray): The array of variance values corresponding to each 
                value in the linspace of g. 
        '''

        #even order 
        if loworder % 2 == 0:

            #find coefficients
            c = np.empty([int(loworder + 2)])

            #model 1 for even orders
            if self.error_model == 1:

                for k in range(int(loworder + 2)):

                    if k % 2 == 0:
                        c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k) * math.factorial(k//2))
                    else:
                        c[k] = 0.0

                #rms value
                cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

                #variance 
                var1 = (cbar)**2.0 * (math.factorial(loworder + 2))**2.0 * g**(2.0*(loworder + 2))

            #model 2 for even orders
            elif self.error_model == 2:

                for k in range(int(loworder + 2)):

                    if k % 2 == 0:

                        #skip first coefficient
                        if k == 0:
                            c[k] = 0.0
                        else:
                            c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2) \
                                   * math.factorial(k//2 - 1) * 4.0**(k))
                    else:
                        c[k] = 0.0

                #rms value
                cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

                #variance
                var1 = (cbar)**2.0 * (math.factorial(loworder//2))**2.0 * (4.0 * g)**(2.0*(loworder + 2))

        #odd order
        else:

            #find coefficients
            c = np.empty([int(loworder + 1)])

            #model 1 for odd orders
            if self.error_model == 1:

                for k in range(int(loworder + 1)):

                    if k % 2 == 0:
                        c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k) * math.factorial(k//2))
                    else:
                        c[k] = 0.0

                #rms value
                cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

                #variance
                var1 = (cbar)**2.0 * (math.factorial(loworder + 1))**2.0 * g**(2.0*(loworder + 1))

            #model 2 for odd orders
            elif self.error_model == 2:

                for k in range(int(loworder + 1)):

                    if k % 2 == 0:

                        #skip first coefficient
                        if k == 0:
                            c[k] = 0.0
                        else:
                            c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2) \
                                    * math.factorial(k//2 - 1) * 4.0**(k))
                    else:
                        c[k] = 0.0

                #rms value
                cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

                #variance
                var1 = (cbar)**2.0 * (math.factorial((loworder-1)//2))**2.0 * (4.0 * g)**(2.0*(loworder + 1))

        return var1


    def variance_high(self, g, highorder):

        r'''
        A function to calculate the variance corresponding to the large-g expansion model.

        Example:
            Bivariate.variance_low(g=np.linspace(1e-6, 0.5, 100), highorder=23)

        Parameters:
            g (numpy.linspace): The linspace over which this calculation is performed.

            highorder (int): The order to which we know our expansion model. This must be 
                a single value.

        Returns:
            var2 (numpy.ndarray): The array of variance values corresponding to each value 
                in the linspace of g. 
        '''

        #find coefficients
        d = np.zeros([int(highorder) + 1])

        #model 1
        if self.error_model == 1:

            for k in range(int(highorder) + 1):

                d[k] = special.gamma(k/2.0 + 0.25) * (-0.5)**k * (math.factorial(k)) / (2.0 * math.factorial(k))

            #rms value (ignore first two coefficients in this model)
            dbar = np.sqrt(np.sum((np.asarray(d)[2:])**2.0) / (highorder-1))

            #variance
            var2 = (dbar)**2.0 * (g)**(-1.0) * (math.factorial(highorder + 1))**(-2.0) \
                    * g**(-2.0*highorder - 2)

        #model 2
        elif self.error_model == 2:

            for k in range(int(highorder) + 1):

                d[k] = special.gamma(k/2.0 + 0.25) * special.gamma(k/2.0 + 1.0) * 4.0**(k) \
                       * (-0.5)**k / (2.0 * math.factorial(k))

            #rms value
            dbar = np.sqrt(np.sum((np.asarray(d))**2.0) / (highorder + 1))

            #variance
            var2 = (dbar)**2.0 * g**(-1.0) * (special.gamma((highorder + 3)/2.0))**(-2.0) \
                    * (4.0 * g)**(-2.0*highorder - 2.0)

        return var2

__init__(error_model='informative')

An accompanying class to Models() that possesses the truncation error models that are included as variances with the small-g and large-g expansions.

Example

Uncertainties()

Parameters:

Name Type Description Default
error_model str

The name of the error model to use in the calculation. Options are 'uninformative' and 'informative'. Default is 'informative'.

 'informative'

Returns:

Type Description

None.
Source code in samba/models.py 
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
def __init__(self, error_model='informative'):

    r'''
    An accompanying class to Models() that possesses the truncation error models
    that are included as variances with the small-g and large-g expansions. 

    Example:
        Uncertainties()

    Parameters:
        error_model (str): The name of the error model to use in the calculation. 
            Options are 'uninformative' and 'informative'. Default is 'informative'.

    Returns:
        None.
    '''

    #assign error model 
    if error_model == 'uninformative':
       self.error_model = 1

    elif error_model == 'informative':
        self.error_model = 2

    else:
        raise ValueError("Please choose 'uninformative' or 'informative'.")

    return None

variance_high(g, highorder)

A function to calculate the variance corresponding to the large-g expansion model.

Example

Bivariate.variance_low(g=np.linspace(1e-6, 0.5, 100), highorder=23)

Parameters:

Name Type Description Default
g linspace

The linspace over which this calculation is performed.

 required
highorder int

The order to which we know our expansion model. This must be a single value.

 required

Returns:

Name Type Description
var2 ndarray

The array of variance values corresponding to each value in the linspace of g.
Source code in samba/models.py 
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
def variance_high(self, g, highorder):

    r'''
    A function to calculate the variance corresponding to the large-g expansion model.

    Example:
        Bivariate.variance_low(g=np.linspace(1e-6, 0.5, 100), highorder=23)

    Parameters:
        g (numpy.linspace): The linspace over which this calculation is performed.

        highorder (int): The order to which we know our expansion model. This must be 
            a single value.

    Returns:
        var2 (numpy.ndarray): The array of variance values corresponding to each value 
            in the linspace of g. 
    '''

    #find coefficients
    d = np.zeros([int(highorder) + 1])

    #model 1
    if self.error_model == 1:

        for k in range(int(highorder) + 1):

            d[k] = special.gamma(k/2.0 + 0.25) * (-0.5)**k * (math.factorial(k)) / (2.0 * math.factorial(k))

        #rms value (ignore first two coefficients in this model)
        dbar = np.sqrt(np.sum((np.asarray(d)[2:])**2.0) / (highorder-1))

        #variance
        var2 = (dbar)**2.0 * (g)**(-1.0) * (math.factorial(highorder + 1))**(-2.0) \
                * g**(-2.0*highorder - 2)

    #model 2
    elif self.error_model == 2:

        for k in range(int(highorder) + 1):

            d[k] = special.gamma(k/2.0 + 0.25) * special.gamma(k/2.0 + 1.0) * 4.0**(k) \
                   * (-0.5)**k / (2.0 * math.factorial(k))

        #rms value
        dbar = np.sqrt(np.sum((np.asarray(d))**2.0) / (highorder + 1))

        #variance
        var2 = (dbar)**2.0 * g**(-1.0) * (special.gamma((highorder + 3)/2.0))**(-2.0) \
                * (4.0 * g)**(-2.0*highorder - 2.0)

    return var2

variance_low(g, loworder)

A function to calculate the variance corresponding to the small-g expansion model.

Example

Bivariate.variance_low(g=np.linspace(1e-6, 0.5, 100), loworder=5)

Parameters:

Name Type Description Default
g linspace

The linspace over which this calculation is performed.

 required
loworder int

The order to which we know our expansion model. Must be passed one at a time if more than one model is to be calculated.

 required

Returns:

Name Type Description
var1 ndarray

The array of variance values corresponding to each value in the linspace of g.
Source code in samba/models.py 
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
def variance_low(self, g, loworder):


    r'''
    A function to calculate the variance corresponding to the small-g expansion model.

    Example:
        Bivariate.variance_low(g=np.linspace(1e-6, 0.5, 100), loworder=5)

    Parameters:
        g (numpy.linspace): The linspace over which this calculation is performed.

        loworder (int): The order to which we know our expansion model. Must be 
            passed one at a time if more than one model is to be calculated.

    Returns:
        var1 (numpy.ndarray): The array of variance values corresponding to each 
            value in the linspace of g. 
    '''

    #even order 
    if loworder % 2 == 0:

        #find coefficients
        c = np.empty([int(loworder + 2)])

        #model 1 for even orders
        if self.error_model == 1:

            for k in range(int(loworder + 2)):

                if k % 2 == 0:
                    c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k) * math.factorial(k//2))
                else:
                    c[k] = 0.0

            #rms value
            cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

            #variance 
            var1 = (cbar)**2.0 * (math.factorial(loworder + 2))**2.0 * g**(2.0*(loworder + 2))

        #model 2 for even orders
        elif self.error_model == 2:

            for k in range(int(loworder + 2)):

                if k % 2 == 0:

                    #skip first coefficient
                    if k == 0:
                        c[k] = 0.0
                    else:
                        c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2) \
                               * math.factorial(k//2 - 1) * 4.0**(k))
                else:
                    c[k] = 0.0

            #rms value
            cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

            #variance
            var1 = (cbar)**2.0 * (math.factorial(loworder//2))**2.0 * (4.0 * g)**(2.0*(loworder + 2))

    #odd order
    else:

        #find coefficients
        c = np.empty([int(loworder + 1)])

        #model 1 for odd orders
        if self.error_model == 1:

            for k in range(int(loworder + 1)):

                if k % 2 == 0:
                    c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k) * math.factorial(k//2))
                else:
                    c[k] = 0.0

            #rms value
            cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

            #variance
            var1 = (cbar)**2.0 * (math.factorial(loworder + 1))**2.0 * g**(2.0*(loworder + 1))

        #model 2 for odd orders
        elif self.error_model == 2:

            for k in range(int(loworder + 1)):

                if k % 2 == 0:

                    #skip first coefficient
                    if k == 0:
                        c[k] = 0.0
                    else:
                        c[k] = np.sqrt(2.0) * special.gamma(k + 0.5) * (-4.0)**(k//2) / (math.factorial(k//2) \
                                * math.factorial(k//2 - 1) * 4.0**(k))
                else:
                    c[k] = 0.0

            #rms value
            cbar = np.sqrt(np.sum((np.asarray(c))**2.0) / (loworder//2 + 1))

            #variance
            var1 = (cbar)**2.0 * (math.factorial((loworder-1)//2))**2.0 * (4.0 * g)**(2.0*(loworder + 1))

    return var1