text:
# Displaying details of the included models & equations

---

[Author] ERGM

---

A number of functions can be used to display the equations for the various models already encoded in Pyfolding.

1. We can list all of the models available by calling the ```list_models``` function.

2. If we want to know any further information, we can print individual equations as each function has a pre-formatted $\LaTeX$ equation. This can be displayed using the ```print_equation()``` function.

3. If we would like to get further details, for example the arguments of the equation, we can get those easily by request the fit function arguments using the ```fit_func_args``` function.
---

#### First off lets load pyfolding into this ipython notebook:

```python
# The following two command lines are always required for the correct running of PyFolding

# this command causes figures to be plotted inline with the text
%matplotlib inline

# this command imports all pyfolding models and pyfolding ising models
import pyfolding
from pyfolding import *

```

    <IPython.core.display.Javascript object>

    PyFolding: Jupyter autoscrolling has been disabled

---
#### Now lets list all the models available in pyfolding {does not list the Ising model}

```python
# Command imports pyfolding models
from pyfolding.models import *

# command lists models
list_models()

# After the model name:
#'Verified: True' signifies that we have tested the model and it functions as expected.
#'Verified:False' signifies that the model has not been tested.
```

    [('ChevronPolynomialFit', 'Verified: True'),
     ('HeteropolymerIsingEquilibrium', 'Verified: False'),
     ('HomozipperIsingEquilibrium', 'Verified: True'),
     ('ParallelTwoStateChevron', 'Verified: False'),
     ('ParallelTwoStateUnfoldingChevron', 'Verified: False'),
     ('TemplateModel', 'Verified: False'),
     ('ThreeStateChevron', 'Verified: True'),
     ('ThreeStateDimericIEquilibrium', 'Verified: True'),
     ('ThreeStateEquilibrium', 'Verified: True'),
     ('ThreeStateFastPhaseChevron', 'Verified: True'),
     ('ThreeStateMonoIEquilibrium', 'Verified: True'),
     ('ThreeStateSequentialChevron', 'Verified: True'),
     ('TwoStateChevron', 'Verified: True'),
     ('TwoStateChevronMovingTransition', 'Verified: True'),
     ('TwoStateDimerEquilibrium', 'Verified: True'),
     ('TwoStateEquilibrium', 'Verified: True'),
     ('TwoStateEquilibriumSloping', 'Verified: True')]

---
#### Print individual equations

If we want to know any further information, we can print each equation using the *Print_equation*:

```python
# The function has two parts:
# part1 states the model you want to print: "******()"
# part2 prints equation: ".print_equation()"
# as an example I have printed a two-state kinetic folding model:
# print the info
models.TwoStateEquilibrium().info()

```

$$F = \frac{\exp( m(x-d_{50})) / RT} { 1+\exp(m(x-d_{50}))/RT}$$

     Two state equilibrium denaturation curve - No sloping baseline.

    	Folding Scheme:
    		N <-> D

    	Params:
    		F = Fraction unfolded
    		m = m-value
    		x = denaturant concentration (M)
    		d50 = denaturant midpoint (M)
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Clarke and Fersht. Engineered disulfide bonds as probes of
    		the folding pathway of barnase: Increasing the stability
    		of proteins against the rate of denaturation.
    		Biochemistry (1993) vol. 32 (16) pp. 4322-4329

```python
# To print all equations at once
# create a list

folding_models = [TwoStateEquilibrium,
                  TwoStateDimerEquilibrium,
                  ThreeStateEquilibrium,
                  ThreeStateMonoIEquilibrium,
                  ThreeStateDimericIEquilibrium,
                  HomozipperIsingEquilibrium,
                  HeteropolymerIsingEquilibrium]

#feed the list into print_equation()

for model in folding_models:
    model().info()

```

$$F = \frac{\exp( m(x-d_{50})) / RT} { 1+\exp(m(x-d_{50}))/RT}$$

     Two state equilibrium denaturation curve - No sloping baseline.

    	Folding Scheme:
    		N <-> D

    	Params:
    		F = Fraction unfolded
    		m = m-value
    		x = denaturant concentration (M)
    		d50 = denaturant midpoint (M)
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Clarke and Fersht. Engineered disulfide bonds as probes of
    		the folding pathway of barnase: Increasing the stability
    		of proteins against the rate of denaturation.
    		Biochemistry (1993) vol. 32 (16) pp. 4322-4329

$$\begin{equation} \
				\begin{aligned} \
				& \Upsilon_{obs} = \Upsilon_N \cdot (1-F_D)  + \Upsilon_D \cdot F_D \\ \
				\text{where} \\ \
				& \Upsilon_N = \alpha_N+\beta_N x \\ \
				& \Upsilon_D = \alpha_D+\beta_D x \\ \
				& F_D = \frac{\sqrt{((K_U^2 + (8 K_U Pt)) - K_U}}  {4 Pt} \\ \
				& K_U = \exp \frac{(RT \ln(Pt - m(d_{50} - x))} {RT}\
				\end{aligned}\
				\end{equation}$$

     Two State model for a dimer denaturation Equilibrium - No Intermediate.

    	Folding Scheme:
    		N2 <-> 2D

    	Params:
    		Y_obs = spectroscopic signal at a given concentration of urea
    		Y_N = spectroscopic signal for native monomeric subunits at a concentration of Pt
    		Y_D = spectroscopic signal for denatured monomeric subunits at a concentration of Pt
    		alpha_N = intercept of the native baseline at low denaturation concentrations
    		beta_N = slope/gradient of the native baseline at low denaturation concentrations
    		alpha_D = intercept of the denatured baseline at high denaturation concentrations
    		beta_D = slope/gradient of the denatured baseline at high denaturation concentrations
    		F_D = fraction of unfolded monomers
    		K_U = Equilibrium Constant for Unfolding of dimer.
    		Pt = total protein concentration. This variable needs to be set per denaturation curve.
    		m = m-value
    		x = denaturant concentration (M)
    		d50 = denaturant midpoint (M)
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Mallam and Jackson. Folding studies on a knotted protein.
    		Journal of Molecular Biology (2005) vol. 346 (5) pp. 1409-1421

$$\begin{equation} \
				\begin{aligned} \
				& \Upsilon_{obs} = \Upsilon_N F_N + \Upsilon_I F_I + \Upsilon_D F_D \ \\ \
				\text{where:}  \\ \
				& F_N = \frac{1} {1 + K_{NI} + K_{NI} K_{ID}}\\ \
				& F_I = \frac{K_{NI}} {1 + K_{NI} + K_{NI} K_{ID}}\\ \
				& F_D = \frac{K_{NI} K_{ID}} {1 + K_{NI} + K_{NI} K_{ID}}\\ \
				\text{and:}  \\ \
				& K_{NI} = \exp \frac{\Delta G_{NI}^{H_2O} + m_{NI} x} {RT}\\ \
				& K_{ID} = \exp \frac{\Delta G_{ID}^{H_2O} + m_{ID} x} {RT}\\ \
				\\ \
				\text{thus:} \\ \
				& \Upsilon_{obs} = \frac{ \Upsilon_N + \Upsilon_I \exp \frac {\Delta G_{NI}^{H_2O} + m_{NI} x} {RT} + \
				\Upsilon_D \exp \frac{\Delta G_{NI}^{H_2O} + m_{NI} x} {RT} \cdot \exp \frac{\Delta G_{ID}^{H_2O} + m_{ID} x} {RT}} {1 + \exp \
				\frac{\Delta G_{NI}^{H_2O} + m_{NI} x} {RT} + \exp \frac{\Delta G_{NI}^{H_2O} + m_{NI} x} {RT} \cdot  \
 			    \exp \frac{\Delta G_{ID}^{H_2O} + m_{ID} x} {RT}}\
				\end{aligned}\
				\end{equation}$$

     Three state equilbrium denaturation curve.

    	Folding Scheme:
    		N <-> I <-> D

    	Params:
    		Y_obs = The spectroscopic signal maximum as a function of denaturant concentration
    		Y_N = spectroscopic signals of the  native state
    		Y_D = spectroscopic signals of the denatured state
    		F_D = fraction denatured
    		F_N = fraction native
    		F_I = fraction intermediate
    		Kni = equilibrium contstant of unfolding native to intermediate state
    		Kid = equilibrium contstant of unfolding intermediate to denatured state
    		DGni = stability of native state relative to intermediate state
    		m_ni = m-value of native to intermediate transition
    		DGid = stability of intermediate state relative to denatured state
    		m_id = m-value of intermediate to denatured transition
    		x = denaturant concentration (M)
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Hecky J, Muller K.M. Structural Perturbation and Compensation by Directed
    		Evolution at Physiological Temperature Leads to Thermostabilization of
    		beta-Lactamase. (2005) Biochemistry 44. pp. 12640-12654

$$\begin{equation} \
				\begin{aligned} \
				& \Upsilon_{rel} = \Upsilon_N F_N + \Upsilon_I F_I + \Upsilon_D F_D \\ \
				\text{expanded:} \\ \
				& \Upsilon_{rel} = \Upsilon_N \cdot \frac{2PtF_I^2} {K_1} + \Upsilon_I F_I + \Upsilon_D * K_2F_I \\ \
				\\ \
				\text{where:} \\ \
				& F_I = \frac {- K_1 (1+K_2) + \sqrt{(K_1^2 (1+K_2)^2 + (8 Pt K_1))}} {4Pt} \\ \
				& K_1 = \exp \frac{-\Delta G_{H_20}^1 + m_1 x} {RT} \\ \
				& K_2 = \exp \frac{-\Delta G_{H_20}^2 + m_2 x} {RT}\
				\end{aligned}\
				\end{equation}$$

     Three State model for a dimer denaturation Equilibrium - Monomeric intermediate.

    	Folding Scheme:
    		N2 <-> 2I <-> 2D

    	Params:
    		Y_rel = spectroscopic signal at a given concentration of urea
    		Y_N = spectroscopic signal for native state
    		Y_D = spectroscopic signal for denatured state
    		Y_I = spectroscopic signal for intermediate state
    		F_D = fraction denatured monomers
    		F_N = fraction native dimers
    		F_I = fraction intermediate dimers
    		Pt = total protein concentration. This variable needs to be set per denaturation curve.
    		K1 = equilibrium constant of unfolding for native to intermediate state
    		K2 = equilibrium constant of unfolding for intermediate to denatured state
    		DG1 = stability of native state relative to intermediate state
    		m1 = m-value of native to intermediate transition
    		DG2 = stability of intermediate state relative to denatured state
    		m2 = m-value of intermediate to denatured transition
    		x = denaturant concentration (M)
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Mallam and Jackson. Folding studies on a knotted protein.
    		Journal of Molecular Biology (2005) vol. 346 (5) pp. 1409-1421

$$\begin{equation} \
				\begin{aligned} \
				& \Upsilon_{rel} = \Upsilon_N F_N + \Upsilon_I F_I + \Upsilon_D F_D \\ \
				\text{expanded:} \\ \
				& \Upsilon_{rel} = \Upsilon_N \cdot \frac{2PtF_D^2} {K_1 K_2} + \Upsilon_I \frac{2PtF_D^2} {K_2} + \Upsilon_D * (F_D) \\ \
				\\ \
				\text{where:} \\ \
				& F_D = \frac {- K_1 K_2 + \sqrt{((K_1 K_2)^2 + 8(1+K_1)(K_1 K_2)Pt)}} {4Pt (1 + K_1)} \\ \
				& K_1 = \exp \frac{-\Delta G_{H_20}^1 + m_1 x} {RT} \\ \
				& K_2 = \exp \frac{-\Delta G_{H_20}^2 + m_2 x} {RT}\
				\end{aligned}\
				\end{equation}$$

     Three State model for a dimer denaturation Equilibrium - Dimeric Intermediate.

    	Folding Scheme:
    		N2 <-> I2 <-> 2D

    	Params:
    		Y_rel = spectroscopic signal at a given concentration of urea
    		Y_N = spectroscopic signal for native state
    		Y_D = spectroscopic signal for denatured state
    		Y_I = spectroscopic signal for intermediate state
    		F_D = fraction denatured monomers
    		F_N = fraction native dimers
    		F_I = fraction intermediate dimers
    		Pt = total protein concentration. This variable needs to be set per denaturation curve.
    		K1 = equilibrium contstant of unfolding native to intermediate state
    		K2 = equilibrium contstant of unfolding intermediate to denatured state
    		DG1 = stability of native state relative to intermediate state
    		m1 = m-value of native to intermediate transition
    		DG2 = stability of intermediate state relative to denatured state
    		m2 = m-value of intermediate to denatured transition
    		x = denaturant concentration (M)
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Mallam and Jackson. Folding studies on a knotted protein.
    		Journal of Molecular Biology (2005) vol. 346 (5) pp. 1409-1421

$$\text{the partition function } (q) \text{ and thus fraction of folded protein } (f) \text{ of n arrayed repeats are given by:}\\ \
				\begin{equation} \\ \
				\begin{aligned} \
				& q = 1 + \frac{\kappa([\kappa \tau]^{n+1} - [n+1]\kappa \tau - n)} {(\kappa \tau + 1)^2}  \\ \
				\\ \
				& f = \frac{1} {n} \sum^{n}_{i=0}i\frac{(n-i+1)\kappa^i\tau^{i-1}} {q} \\ \
				\\ \
				\text{where:} \\ \
				& \kappa (x) = \exp\frac{-G_i} {RT} = \exp\frac{-G_{i,H_20} + m_i x} {RT}  \\ \
				\\ \
				& \tau (x) = \exp\frac{-G_{i,i+1}} {RT} \
				\end{aligned}\
				\end{equation}$$

     Homopolymer Zipper Ising model

    	Params:
    		q = partition function
    		f = fraction of folded protein
    		Kappa = equilibrium constant of folding for a given repeating unit
    		Tau = equilibrium constant of association between 2 repeating units
    		n = number of repeating units
    		x = denaturant concentration (M)
    		Gi = intrinsic stability (folding energy) of a repeating unit i
    		mi = denaturant sensitivity of the intrinsic stability of a repeating unit i
    		Gi,i+1 = interface interaction energy between 2 repeating units
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Reference:
    		Aksel and Barrick. Analysis of repeat-protein folding using
    		nearest-neighbor statistical mechanical models.
    		Methods in enzymology (2009) vol. 455 pp. 95-125

$$\text{the partition function } (q) \text{ and thus fraction of folded protein } (f) \text{ of n arrayed repeats are given by:}  \\ \
			\begin{equation} \\ \
			\begin{aligned} \\ \
			\kappa(x) &= \exp(-(\Delta G_{intrinsic} - m_{intrinsic}x) / RT) \\ \
			\tau(x) &= \exp(-\Delta G_{interface}) / RT) \\ \
			q(i) &=  \
			\begin{bmatrix} 0 & 1\end{bmatrix}  \
			\begin{bmatrix} \kappa_1\tau_{-1} & 1\\ \kappa & 1 \end{bmatrix} \
			\ldots  \
			\begin{bmatrix} \kappa_n\tau_{n-1} & 1\\ \kappa & 1 \end{bmatrix} \
			\begin{bmatrix} 1 \\ 1 \end{bmatrix} \\ \
			\theta &=  \frac{1}{nq(n)} \sum_{i=0}^{n}{q(i)} \
			\end{aligned} \
			\end{equation}$$

     Heteropolymer Ising model

    	Params:
    		q = partition function
    		f = fraction of folded protein
    		Kappa = equilibrium constant of folding for a given repeating unit
    		Tau = equilibrium constant of association between 2 repeating units
    		n = number of repeating units
    		x = denaturant concentration (M)
    		DG_intrinsic = intrinsic stability (folding energy) of a repeating unit i
    		m_intrinsic = denaturant sensitivity of the intrinsic stability of a repeating unit i
    		DG_interface = interface interaction energy between 2 repeating units
    		R = Universal Gas Constant (kcal.mol-1.K-1)
    		T = Temperature (Kelvin)

    	Notes:
    		Aksel and Barrick. Analysis of repeat-protein folding using
    		nearest-neighbor statistical mechanical models.
    		Methods in enzymology (2009) vol. 455 pp. 95-125

---
#### Variables in a given Equation
If we would like to get further details, for example the variables of an equation, we can get those easily by request the fit function arguments:

```python
# The function has two parts:
# part1 states the model you want to find about: "******()"
# part2 prints the variables: ".fit_func_args"
# as an example I have used the two-state kinetic folding model:

TwoStateChevron().fit_func_args
```

    ['kf', 'mf', 'ku', 'mu']

---

### End of this Notebook.

---

max_stars_repo_name: quantumjot/PyFolding-Notebooks

max_issues_repo_name: quantumjot/PyFolding-Notebooks

max_forks_repo_name: quantumjot/PyFolding-Notebooks


text:
# Medición de tiempo con Python

El propósito de estos apuntes es presentar varias formas de medir el tiempo que toma una determinada función. La función de prueba estará basada en cuatro bucles `for`.

```python
def func(a,b,c,d):
    s = 0
    for i in range(a):
        for j in range(b):
            for k in range(c):
                for l in range(d):
                    s += i*j*k*l
    return s
```

```python
import random
```

```python
a = random.randrange(1,10)
b = random.randrange(1,10)
c = random.randrange(1,10)
d = random.randrange(1,10)
```

## Módulo time

Este módulo proporciona varios métodos para medir tiempo, todos presentan un resolución de orden similar.

```python
import time
```

```python
tic = time.time()
func(a,b,c,d)
toc = time.time()
print('tiempo de func con time.time %f'%(toc-tic))
```

    tiempo de func con time.time 0.000108

```python
tic = time.perf_counter()
func(a,b,c,d)
toc = time.perf_counter()
print('tiempo de func con time.perf_counter %f'%(toc-tic))
```

    tiempo de func con time.perf_counter 0.000108

```python
tic = time.process_time()
func(a,b,c,d)
toc = time.process_time()
print('tiempo de func con time.process_time %f'%(toc-tic))
```

    tiempo de func con time.process_time 0.000106

Se observa que las mediciones realizadas tienen resultados parecidos. Para ver si hay muchas diferencias se van a realizar varias medidas, las cuales serán guardadas en una lista.

```python
timelist = []
perflist = []
proclist = []
veces = 100
```

```python
for i in range(veces):
    tic = time.time()
    func(a,b,c,d)
    toc = time.time()
    timelist.append(toc-tic)

    tic = time.perf_counter()
    func(a,b,c,d)
    toc = time.perf_counter()
    perflist.append(toc-tic)

    tic = time.process_time()
    func(a,b,c,d)
    toc = time.process_time()
    proclist.append(toc-tic)
```

Ahora se graficarán las mediciones

```python
import matplotlib.pyplot as plt
```

```python
x = list(range(veces))
plt.plot(x,timelist,label='time.time')
plt.plot(x,perflist,label='time.perf_counter')
plt.plot(x,proclist,label='time.process_time')
plt.legend()
```

En la gráfica se observa que las tres curvas presentan, en general, la misma tendencia, pero con algunas distorsiones. ¿Cómo estará el promedio de los métodos de medición?

```python
import statistics
```

```python
noms = ['time','perf_counter','process_time']
proms1 = [
    statistics.mean(timelist),
    statistics.mean(perflist),
    statistics.mean(proclist)
]
plt.bar(noms,proms1,color='r')
plt.tight_layout()
```

Se observa que todos los promedios están en el mismo orden de magnitud, pero que no son exactamente iguales a pesar de tratarse de la misma función, que la cantidad de pruebas es la misma y que los métodos de medición tienen una resolución similar. ¿Cuál podría ser la  causa de estas variaciones?

## Módulo timeit

El método timeit devuelve el tiempo acumulado luego de ejecutar la función indicada, la cantidad de veces en su opción `number`. Por esta razón, para obtener el tiempo promedio transcurrido se dividirá el resultado de `timeit` entre la cantidad de veces que se indique en `number`.

```python
import timeit
```

```python
setup = '''
from __main__ import func

from __main__ import a

from __main__ import b

from __main__ import c

from __main__ import d

'''
```

```python
stmt = '''
func(a,b,c,d)
'''
```

```python
number = 10000
timeitres = timeit.timeit(
    setup = setup,
    stmt = stmt,
    number = number
)
print('tiempo de func con timeit.timeit %f'%(timeitres/number))
```

    tiempo de func con timeit.timeit 0.000056

Para realizar varias repeticiones `timeit` proporciona el método `repeat` el cual devuelve una lista con la cantidad de veces indicada en la opción respectiva. Igual que en el caso anterior los resultados se tienen que dividir entre la cantidad de veces porque devuelve el tiempo acumulado.

```python
repeat = 100
timeitreps = timeit.repeat(
    setup = setup,
    stmt = stmt,
    number = number,
    repeat = repeat
)
repeatlist = [t/number for t in timeitreps]
```

Ahora se procederá a hacer una gráfica de las mediciones

```python
plt.plot(x,repeatlist,label='timeit.repeat')
plt.legend()
```

## Comparación de los módulos

Gráfica conjunta de las mediciones

```python
plt.plot(x,timelist,label='time.time')
plt.plot(x,perflist,label='time.perf_counter')
plt.plot(x,proclist,label='time.process_time')
plt.plot(x,repeatlist,label='timeit.repeat')
plt.legend()
```

Ahora se graficarán los tiempos promedio de todos los métodos.

```python
noms = ['time','perf_counter','process_time','repeat']
proms2 = [
    statistics.mean(timelist),
    statistics.mean(perflist),
    statistics.mean(proclist),
    statistics.mean(repeatlist)
]
plt.bar(noms,proms2,color='r')
plt.tight_layout()
```

Se observa que las mediciones realizadas con los métodos proporcionados por los dos módulos son irregulares, es decir, ambas mediciones presentan picos. Sumar estos picos que se ha introducido en nuestras mediciones y luego promediar implicaría que estas irregularidades están incluidas en el cómputo del promedio y esto desviaría el resultado final de lo real. ¿Cómo medir esta distorsión? ¿Cómo eliminar el ruido?

## Filtro mediana unidimensional

El filtro de mediana es un algoritmo empleado para eliminar ruido.

```python
def filtro_mediana(sign, ventana):
    signc = []
    tam = len(sign)
    offs = int((ventana-1)/2)
    for i in range(tam):
        inicio = i - offs if (i - offs > 0) else i
        fin = i + offs if (i + offs < tam) else tam
        signc.append(statistics.median(sign[inicio:fin]))
    return signc
```

```python
ventana = 7
timelistsr = filtro_mediana(timelist,ventana)
```

```python
# gráfico de las mediciones realizadas con time.time con y sin ruido
plt.plot(x,timelist, label='con ruido')
plt.plot(x,timelistsr, label='sin ruido')
# para mostrar la leyenda
plt.legend()
```

```python
perflistsr = filtro_mediana(perflist,ventana)
```

```python
# gráfico de las mediciones realizadas con time.perf_counter con y sin ruido
plt.plot(x,perflist, label='con ruido')
plt.plot(x,perflistsr, label='sin ruido')
# para mostrar la leyenda
plt.legend()
```

```python
proclistsr = filtro_mediana(proclist,ventana)
```

```python
# gráfico de las mediciones realizadas con time.process_time con y sin ruido
plt.plot(x,proclist, label='con ruido')
plt.plot(x,proclistsr, label='sin ruido')
# para mostrar la leyenda
plt.legend()
```

```python
repeatlistsr = filtro_mediana(repeatlist,ventana)
```

```python
# gráfico de las mediciones realizadas con time.process_time con y sin ruido
plt.plot(x,repeatlist, label='con ruido')
plt.plot(x,repeatlistsr, label='sin ruido')
# para mostrar la leyenda
plt.legend()
```

¿Cómo se ven los tiempos sin ruido?

```python
plt.plot(x,timelistsr,label='time.time')
plt.plot(x,perflistsr,label='time.perf_counter')
plt.plot(x,proclistsr,label='time.process_time')
plt.plot(x,repeatlistsr,label='timeit.repeat')
plt.legend()
```

¿Varían mucho los promedios luego de la eliminación del ruido?

```python
noms = ['time','perf_counter','process_time','repeat']
promssr = [
    statistics.mean(timelistsr),
    statistics.mean(perflistsr),
    statistics.mean(proclistsr),
    statistics.mean(repeatlistsr)
]
plt.bar(noms,promssr,color='r')
plt.tight_layout()
```

¿Qué impacto tuvo eliminar el ruido de las señales en el valor final de los promedios? ¿Fueron significativas las distorsiones introducidas en las mediciones? ¿Cuál es la diferencia entre el promedio con ruido y sin ruido? Se deja estas preguntas para su reflexión.

# Implementación de funciones

En esta sección de implementará una función que calcula la suma de números enteros desde el $1$ hasta el número anterior indicado en su argumento.

$$
\begin{align}
\sum_{k=1}^{n-1}k
\end{align}
$$

Esta será implementada de múltiples formas con el propósito de mostrar como la implementación se puede traducir en tiempos de ejecución mayores o menores.

## Implementación con fórmula

Es probable que algunos usuarios, sea por cultura general, por formación, o por alguna razón cualquiera,  ya sepan que esa sumatoria tiene un resultado notable.

$$
\begin{align}
\frac{(n-1)n}{2}
\end{align}
$$

Está será la primera implementación

```python
def sumar_enteros_formula(n=100_000):
    return int((n*(n-1))/2)
```

```python
timeit.timeit(sumar_enteros_formula,number=1)
```

    2.6859997888095677e-06

```python
timeit.timeit(sumar_enteros_formula,number=10)/10
```

    4.5020001380180476e-07

```python
timeit.timeit(sumar_enteros_formula,number=100)/100
```

    1.9192000308976275e-07

## Implementación con funciones internas

Imagine que usted desconocía que el resultado de la sumatoria era notable. Un camino alternativo, es emplear las funciones internas que proporciona el lenguaje.

```python
def sumar_enteros_idiomatic(n=100_000):
    return sum(range(n))
```

```python
timeit.timeit(sumar_enteros_idiomatic,number=1)
```

    0.001114943000175117

```python
timeit.timeit(sumar_enteros_idiomatic,number=10)/10
```

    0.0011519827000029182

```python
timeit.timeit(sumar_enteros_idiomatic,number=100)/100
```

    0.001145526640002572

## Implementación for

Ahora suponga que usted es un usuario poco familiarizado con Python, y no sabía de estas funciones internas, pero está relativamente habituado a la sintaxis, y decide usar un bucle `for`.

```python
def sumar_enteros_for(n=100_000):
    s = 0
    for i in range(n):
        s += i
    return s
```

```python
timeit.timeit(sumar_enteros_for,number=1)
```

    0.003450418999818794

```python
timeit.timeit(sumar_enteros_for,number=10)/10
```

    0.0034587901999657334

```python
timeit.timeit(sumar_enteros_for,number=100)/100
```

    0.003610340300001553

Un usuario aún más inexperto quizás interaría algo con la siguiente forma:

```python
def sumar_enteros_for_incremento(n=100_000):
    s = 0
    c = 0
    for i in range(n):
        s += c
        c += 1
    return s
```

```python
timeit.timeit(sumar_enteros_for_incremento,number=1)
```

    0.006406569000318996

```python
timeit.timeit(sumar_enteros_for_incremento,number=10)/10
```

    0.0063905099000294285

```python
timeit.timeit(sumar_enteros_for_incremento,number=100)/100
```

    0.00630345029000182

Bastó con incluir un contador distinto del índice y acumularlo para que el tiempo sea mayor a la implementación con anterior. ¿Qué pasaría si se incluye una operación de comparación porque desconoce del funcionamiento del bucle `for`?

```python
def sumar_enteros_for_incremento_comp(n=100_000):
    s = 0
    c = 0
    for i in range(n):
        if i < n:
            s += c
            c += 1
    return s
```

```python
timeit.timeit(sumar_enteros_for_incremento_comp,number=1)
```

    0.008187398999780271

```python
timeit.timeit(sumar_enteros_for_incremento_comp,number=10)/10
```

    0.008175903100027426

```python
timeit.timeit(sumar_enteros_for_incremento_comp,number=100)/100
```

    0.0076293956099971185

Al añadir la operación de comparación la función ahora es incluso más lenta, y eso que los añadidos son operaciones muy sencillas.

## Implementación while

Quizás usted no entendió el uso del bucle `for` en la documentación y el bucle `while` le resulta más parecido a otros lenguajess de programación que ya ha empleado.

```python
def sumar_enteros_while(n=100_000):
    i = 0
    s = 0
    while (i < n):
        s += i
        i += 1
    return s
```

```python
timeit.timeit(sumar_enteros_while,number=1)
```

    0.007319155999994109

```python
timeit.timeit(sumar_enteros_while,number=10)/10
```

    0.007463831899985962

```python
timeit.timeit(sumar_enteros_while,number=100)/100
```

    0.006493551399998978

Se observa que la primera implmentación de la función que emplea el bucle `for` es ligeramente más rápida que la función implementada con bucle `while`. Esto se debe a que la implementación con bucle `for` está más optimizada a nivel de ensamblador. Se debe tener en cuenta estos detalles para implementar funciones.

## Comparación de las funciones implementadas

Se definen funciones con lo aprendido para poder usarlo de forma modular

```python
def medir_timeit(setup,stmt,number,repeat):
    timeitreps = timeit.repeat(
        setup = setup,
        stmt = stmt,
        number = number,
        repeat = repeat
    )
    return [t/number for t in timeitreps]
```

```python
def medir_time(f,veces):
    times = []
    tic = 0
    toc = 0
    e = 0
    for i in range(veces):
        tic = time.time()
        f()
        toc = time.time()
        e = toc - tic
        times.append(e)
    return times
```

Mediciones con timeit

```python
setup = '''
from __main__ import sumar_enteros_formula
'''
stmt='''
sumar_enteros_formula()
'''
number = 100
repeat = 100
t0 = medir_timeit(setup,stmt,number,repeat)
```

```python
setup = '''
from __main__ import sumar_enteros_idiomatic
'''
stmt='''
sumar_enteros_idiomatic()
'''
t1 = medir_timeit(setup,stmt,number,repeat)
```

```python
setup = '''
from __main__ import sumar_enteros_for
'''
stmt='''
sumar_enteros_for()
'''
t2 = medir_timeit(setup,stmt,number,repeat)
```

```python
setup = '''
from __main__ import sumar_enteros_for_incremento
'''
stmt='''
sumar_enteros_for_incremento()
'''
t3 = medir_timeit(setup,stmt,number,repeat)
```

```python
setup = '''
from __main__ import sumar_enteros_for_incremento_comp
'''
stmt='''
sumar_enteros_for_incremento_comp()
'''
t4 = medir_timeit(setup,stmt,number,repeat)
```

```python
setup = '''
from __main__ import sumar_enteros_while
'''
stmt='''
sumar_enteros_while()
'''
t5 = medir_timeit(setup,stmt,number,repeat)
```

```python
x = list(range(repeat))
plt.plot(x, t0, label='formula')
plt.plot(x, t1, label='idiomatic')
plt.plot(x, t2, label='for')
plt.plot(x, t3, label='for-inc')
plt.plot(x, t4, label='for-inc-comp')
plt.plot(x, t5, label='while')
# para mostrar la leyenda
plt.legend()
```

Luego de pasar filtro mediano a las mediciones

```python
t0sr = filtro_mediana(t0, 7)
```

```python
t1sr = filtro_mediana(t1, 7)
```

```python
t2sr = filtro_mediana(t2, 7)
```

```python
t3sr = filtro_mediana(t3, 7)
```

```python
t4sr = filtro_mediana(t4, 7)
```

```python
t5sr = filtro_mediana(t5, 7)
```

```python
plt.plot(x, t0sr, label='formula')
plt.plot(x, t1sr, label='idiomatic')
plt.plot(x, t2sr, label='for')
plt.plot(x, t3sr, label='for-inc')
plt.plot(x, t4sr, label='for-inc-comp')
plt.plot(x, t5sr, label='while')
# para mostrar la leyenda
plt.legend()
```

Mediciones con time

```python
r0 = medir_time(sumar_enteros_formula,100)
```

```python
r1 = medir_time(sumar_enteros_idiomatic,100)
```

```python
r2 = medir_time(sumar_enteros_for,100)
```

```python
r3 = medir_time(sumar_enteros_for_incremento,100)
```

```python
r4 = medir_time(sumar_enteros_for_incremento_comp,100)
```

```python
r5 = medir_time(sumar_enteros_while,100)
```

```python
plt.plot(x, r0, label='formula')
plt.plot(x, r1, label='idiomatic')
plt.plot(x, r2, label='for')
plt.plot(x, r3, label='for-inc')
plt.plot(x, r4, label='for-inc-comp')
plt.plot(x, r5, label='while')
# para mostrar la leyenda
plt.legend()
```

Luego de pasar filtro mediano a las mediciones

```python
r0sr = filtro_mediana(r0, 7)
```

```python
r1sr = filtro_mediana(r1, 7)
```

```python
r2sr = filtro_mediana(r2, 7)
```

```python
r3sr = filtro_mediana(r3, 7)
```

```python
r4sr = filtro_mediana(r4, 7)
```

```python
r5sr = filtro_mediana(r5, 7)
```

```python
plt.plot(x, r0sr, label='formula')
plt.plot(x, r1sr, label='idiomatic')
plt.plot(x, r2sr, label='for')
plt.plot(x, r3sr, label='for-inc')
plt.plot(x, r4sr, label='for-inc-comp')
plt.plot(x, r5sr, label='while')
# para mostrar la leyenda
plt.legend()
```

Para la función idiomática ¿Cómo quedan las mediciones de `timeit.repeat` luego de eliminar el ruido?

```python
plt.plot(x, t1, label='timeit.repeat')
plt.plot(x, t1sr, label='timeit.repeat-sr')
# para mostrar la leyenda
plt.legend()
```

¿Y para `time.time`?

```python
plt.plot(x, r1, label='time.time')
plt.plot(x, r1sr, label='time.time-sr')
# para mostrar la leyenda
plt.legend()
```

Para la misma función ¿Cómo varían los promedios de los dos métodos?

```python
noms = ['timeit.repeat','time.time']
proms = [
    statistics.mean(t1),
    statistics.mean(r1)
]
plt.bar(noms,proms,color='r')
plt.tight_layout()
```

¿Después de la eliminación de ruido?

```python
noms = ['timeit.repeat','time.time']
proms = [
    statistics.mean(t1sr),
    statistics.mean(r1sr)
]
plt.bar(noms,proms,color='r')
plt.tight_layout()
```

¿Será importante eliminar el ruido? ¿Cómo variaría el tiempo si cambiaramos el argumento a pasos de 100mil? ¿La tendencia sería la misma en todas las implementaciones? Se dejan estos dudas para el lector.

## Característica del equipo

Los resultados fueron obtenidos al ejecutar el cuaderno en una máquina virtual con cuatro núcleos y 4GB de RAM. Por esta razón, usted podría obtener resultados muy diferentes al ejecutarlo en su equipo.

ext: ipynb

max_stars_repo_name: isaiBr/sobre-python

max_issues_repo_name: isaiBr/sobre-python

max_forks_repo_name: isaiBr/sobre-python


text:
# Start-to-Finish Example: `GiRaFFE_NRPy` 1D tests

### Author: Patrick Nelson

### Adapted from [Start-to-Finish Example: Head-On Black Hole Collision](../Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide.ipynb)

## This module implements a basic GRFFE code to evolve one-dimensional GRFFE waves.

### NRPy+ Source Code for this module:
* [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Exact_Wald.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Exact_Wald.py) [\[**tutorial**\]](Tutorial-GiRaFFEfood_NRPy_Exact_Wald.ipynb) Generates Exact Wald initial data
* [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Aligned_Rotator.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Aligned_Rotator.py) [\[**tutorial**\]](Tutorial-GiRaFFEfood_NRPy_Aligned_Rotator.ipynb) Generates Aligned Rotator initial data
* [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_1D_tests.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_1D_tests.py) [\[**tutorial**\]](Tutorial-GiRaFFEfood_NRPy_1D_tests.ipynb) Generates Alfvén Wave initial data.
* [GiRaFFE_NRPy/GiRaFFE_NRPy_staggered_Afield_flux.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_staggered_Afield_flux.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy_staggered-Afield_flux.ipynb) Generates the expressions to find the flux term of the induction equation.
* [GiRaFFE_NRPy/GiRaFFE_NRPy_staggered_A2B.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_staggered_A2B.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy_staggered-A2B.ipynb) Generates the driver to compute the magnetic field from the vector potential/
* [GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy-BCs.ipynb) Generates the code to apply boundary conditions to the vector potential, scalar potential, and three-velocity.
* [GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy-C2P_P2C.ipynb) Generates the conservative-to-primitive and primitive-to-conservative solvers.
* [GiRaFFE_NRPy/GiRaFFE_NRPy_Metric_Face_Values.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Metric_Face_Values.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy-Metric_Face_Values.ipynb) Generates code to interpolate metric gridfunctions to cell faces.
* [GiRaFFE_NRPy/GiRaFFE_NRPy_PPM.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_PPM.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy-PPM.ipynb) Genearates code to reconstruct primitive variables on cell faces.
* [GiRaFFE_NRPy/GiRaFFE_NRPy_staggered_Source_Terms.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_staggered_Source_Terms.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy_staggered-Source_Terms.ipynb) Generates the expressions to find the flux term of the Poynting flux evolution equation.
* [GiRaFFE_NRPy/Stilde_flux.py](../../edit/in_progress/GiRaFFE_NRPy/Stilde_flux.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy-Stilde_flux.ipynb) Generates the expressions to find the flux term of the Poynting flux evolution equation.
* [../GRFFE/equations.py](../../edit/GRFFE/equations.py) [\[**tutorial**\]](../Tutorial-GRFFE_Equations-Cartesian.ipynb) Generates code necessary to compute the source terms.
* [../GRHD/equations.py](../../edit/GRHD/equations.py) [\[**tutorial**\]](../Tutorial-GRHD_Equations-Cartesian.ipynb) Generates code necessary to compute the source terms.

Here we use NRPy+ to generate the C source code necessary to set up initial data for an Alfvén wave (see [the original GiRaFFE paper](https://arxiv.org/pdf/1704.00599.pdf)). Then we use it to generate the RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration based on the [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4).

# Table of Contents
$$\label{toc}$$

This notebook is organized as follows

1. [Step 1](#initializenrpy): Set core NRPy+ parameters for numerical grids
1. [Step 2](#grffe): Output C code for GRFFE evolution
    1. [Step 2.a](#mol): Output macros for Method of Lines timestepping
1. [Step 3](#gf_id): Import `GiRaFFEfood_NRPy` initial data modules
1. [Step 4](#cparams): Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`
1. [Step 5](#mainc): `GiRaFFE_NRPy_standalone.c`: The Main C Code

# Step 1: Set up core functions and parameters for solving  GRFFE equations \[Back to [top](#toc)\]
$$\label{setup}$$

```python
import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions
# First, we'll add the parent directory to the list of directories Python will check for modules.
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

# Step P1: Import needed NRPy+ core modules:
from outputC import outCfunction, lhrh # NRPy+: Core C code output module
import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import grid as gri               # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface

# Step P2: Create C code output directory:
Ccodesdir = os.path.join("GiRaFFE_staggered_standalone_Ccodes/")
# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
# !rm -r ScalarWaveCurvilinear_Playground_Ccodes
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)

# Step P3: Create executable output directory:
outdir = os.path.join(Ccodesdir,"output/")
cmd.mkdir(outdir)

# Step P5: Set timestepping algorithm (we adopt the Method of Lines)
REAL      = "double"    # Best to use double here.
default_CFL_FACTOR= 0.5 # (GETS OVERWRITTEN WHEN EXECUTED.) In pure axisymmetry (symmetry_axes = 2 below) 1.0 works fine. Otherwise 0.5 or lower.

# Step P6: Set the finite differencing order to 2.
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER",4)

thismodule = "Start_to_Finish-GiRaFFE_NRPy-1D_tests"
TINYDOUBLE = par.Cparameters("REAL", thismodule, "TINYDOUBLE", 1e-100)

import GiRaFFE_NRPy.GiRaFFE_NRPy_Main_Driver_staggered as md
# par.set_paramsvals_value("GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C::enforce_speed_limit_StildeD = False")
par.set_paramsvals_value("GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C::enforce_current_sheet_prescription = False")
```

# Step 2: Output C code for GRFFE evolution \[Back to [top](#toc)\]
$$\label{grffe}$$

We will first write the C codes needed for GRFFE evolution. We have already written a module to generate all these codes and call the functions in the appropriate order, so we will import that here. We will take the slightly unusual step of doing this before we generate the initial data functions because the main driver module will register all the gridfunctions we need. It will also generate functions that, in addition to their normal spot in the MoL timestepping, will need to be called during the initial data step to make sure all the variables are appropriately filled in.

All of this is handled with a single call to `GiRaFFE_NRPy_Main_Driver_generate_all()`, which will register gridfunctions, write all the C code kernels, and write the C code functions to call those.

```python
md.GiRaFFE_NRPy_Main_Driver_generate_all(Ccodesdir)
```

    Output C function calculate_StildeD0_source_term() to file GiRaFFE_staggered_standalone_Ccodes/RHSs/calculate_StildeD0_source_term.h
    Output C function calculate_StildeD1_source_term() to file GiRaFFE_staggered_standalone_Ccodes/RHSs/calculate_StildeD1_source_term.h
    Output C function calculate_StildeD2_source_term() to file GiRaFFE_staggered_standalone_Ccodes/RHSs/calculate_StildeD2_source_term.h
    Output C function calculate_Stilde_rhsD() to file GiRaFFE_staggered_standalone_Ccodes/RHSs/calculate_Stilde_rhsD.h
    Output C function GiRaFFE_NRPy_cons_to_prims() to file GiRaFFE_staggered_standalone_Ccodes/C2P/GiRaFFE_NRPy_cons_to_prims.h
    Output C function GiRaFFE_NRPy_prims_to_cons() to file GiRaFFE_staggered_standalone_Ccodes/C2P/GiRaFFE_NRPy_prims_to_cons.h
    Output C function Workaround_ADM_to_BSSN() to file GiRaFFE_staggered_standalone_Ccodes/Workaround_ADM_to_BSSN.h
    Output C function Workaround_BSSN_to_ADM() to file GiRaFFE_staggered_standalone_Ccodes/Workaround_BSSN_to_ADM.h

## Step 2.a: Output macros for Method of Lines timestepping \[Back to [top](#toc)\]
$$\label{mol}$$

Now, we generate the code to implement the method of lines using the fourth-order Runge-Kutta algorithm.

```python
RK_method = "RK4"

# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.
#       As described above the Table of Contents, this is a 3-step process:
#       3.A: Evaluate RHSs (RHS_string)
#       3.B: Apply boundary conditions (post_RHS_string, pt 1)
import MoLtimestepping.C_Code_Generation as MoL
from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict
RK_order  = Butcher_dict[RK_method][1]
cmd.mkdir(os.path.join(Ccodesdir,"MoLtimestepping/"))
MoL.MoL_C_Code_Generation(RK_method,
    RHS_string      = """
GiRaFFE_NRPy_RHSs(&params,auxevol_gfs,RK_INPUT_GFS,RK_OUTPUT_GFS);""",
    post_RHS_string = """
GiRaFFE_NRPy_post_step(&params,xx,auxevol_gfs,RK_OUTPUT_GFS,n+1);\n""",
    outdir = os.path.join(Ccodesdir,"MoLtimestepping/"))
```

# Step 3: Import `GiRaFFEfood_NRPy` initial data modules \[Back to [top](#toc)\]
$$\label{gf_id}$$

With the preliminaries out of the way, we will write the C functions to set up initial data. There are two categories of initial data that must be set: the spacetime metric variables, and the GRFFE plasma variables. We will set up the spacetime first.

```python
# There are several initial data routines we need to test. We'll control which one we use with a string option
initial_data = "FastWave" # Valid options: "ExactWald", "AlignedRotator", "AlfvenWave", "FastWave"
spacetime = "flat" # Valid options: "ShiftedKerrSchild", "flat"

if spacetime == "ShiftedKerrSchild":
    # Exact Wald is more complicated. We'll need the Shifted Kerr Schild metric in Cartesian coordinates.
    import BSSN.ShiftedKerrSchild as sks
    sks.ShiftedKerrSchild(True)
    import reference_metric as rfm
    par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
    rfm.reference_metric()
    # Use the Jacobian matrix to transform the vectors to Cartesian coordinates.
    drrefmetric__dx_0UDmatrix = sp.Matrix([[sp.diff(rfm.xxSph[0],rfm.xx[0]), sp.diff(rfm.xxSph[0],rfm.xx[1]), sp.diff(rfm.xxSph[0],rfm.xx[2])],
                                           [sp.diff(rfm.xxSph[1],rfm.xx[0]), sp.diff(rfm.xxSph[1],rfm.xx[1]), sp.diff(rfm.xxSph[1],rfm.xx[2])],
                                           [sp.diff(rfm.xxSph[2],rfm.xx[0]), sp.diff(rfm.xxSph[2],rfm.xx[1]), sp.diff(rfm.xxSph[2],rfm.xx[2])]])
    dx__drrefmetric_0UDmatrix = drrefmetric__dx_0UDmatrix.inv()
    gammaDD = ixp.zerorank2()
    for i in range(3):
        for j in range(3):
            for k in range(3):
                for l in range(3):
                    gammaDD[i][j] += drrefmetric__dx_0UDmatrix[(k,i)]*drrefmetric__dx_0UDmatrix[(l,j)]*sks.gammaSphDD[k][l].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
    betaU = ixp.zerorank1()
    for i in range(3):
        for j in range(3):
            betaU[i] += dx__drrefmetric_0UDmatrix[(i,j)]*sks.betaSphU[j].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])
    alpha = sks.alphaSph.subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])

    # Description and options for this initial data
    desc = "Generate a spinning black hole with Shifted Kerr Schild metric."
    loopopts_id ="AllPoints,Read_xxs"

elif spacetime == "flat":
    gammaDD = ixp.zerorank2(DIM=3)
    for i in range(3):
        for j in range(3):
            if i==j:
                gammaDD[i][j] = sp.sympify(1) # else: leave as zero
    betaU = ixp.zerorank1() # All should be 0
    alpha = sp.sympify(1)

    # Description and options for this initial data
    desc = "Generate a flat spacetime metric."
    loopopts_id ="AllPoints" # we don't need to read coordinates for flat spacetime.

# For testing: Also set inverse metric:
gammaUU, unused_gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)

name = "set_initial_spacetime_metric_data"
values_to_print = [
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD00"),rhs=gammaDD[0][0]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD01"),rhs=gammaDD[0][1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD02"),rhs=gammaDD[0][2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD11"),rhs=gammaDD[1][1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD12"),rhs=gammaDD[1][2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD22"),rhs=gammaDD[2][2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU0"),rhs=betaU[0]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU1"),rhs=betaU[1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU2"),rhs=betaU[2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","alpha"),rhs=alpha),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaUU00"),rhs=gammaUU[0][0]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaUU01"),rhs=gammaUU[0][1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaUU02"),rhs=gammaUU[0][2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaUU11"),rhs=gammaUU[1][1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaUU12"),rhs=gammaUU[1][2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaUU22"),rhs=gammaUU[2][2])
                  ]

outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs",
    body     = fin.FD_outputC("returnstring",values_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts = loopopts_id)

```

    Output C function set_initial_spacetime_metric_data() to file GiRaFFE_staggered_standalone_Ccodes/set_initial_spacetime_metric_data.h

Now, we will write out the initial data function for the GRFFE variables.

```python
if initial_data=="AlfvenWave":
    import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_1D_tests as gid
    gid.GiRaFFEfood_NRPy_1D_tests(stagger = True)
    desc = "Generate Alfven wave 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="FastWave":
    import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_1D_tests_fast_wave as gid
    gid.GiRaFFEfood_NRPy_1D_tests_fast_wave()
    desc = "Generate fast wave 1D initial test data for GiRaFFEfood_NRPy."
elif initial_data=="AlignedRotator":
    import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_Aligned_Rotator as gid
    gid.GiRaFFEfood_NRPy_Aligned_Rotator()
    desc = "Generate aligned rotator initial test data for GiRaFFEfood_NRPy."
elif initial_data=="ExactWald":
    import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_Exact_Wald as gid
    M,r0 = sp.symbols("M r0")
    gid.GiRaFFEfood_NRPy_Exact_Wald(gammaDD,M,r0)
    desc = "Generate exact Wald initial test data for GiRaFFEfood_NRPy."
else:
    print("Unsupported Initial Data string "+initial_data+"! Supported ID: AlfvenWave, FastWave, AlignedRotator, or ExactWald")

name = "initial_data"

values_to_print = [\
                   lhrh(lhs=gri.gfaccess("out_gfs","AD0"),rhs=gid.AD[0]),\
                   lhrh(lhs=gri.gfaccess("out_gfs","AD1"),rhs=gid.AD[1]),\
                   lhrh(lhs=gri.gfaccess("out_gfs","AD2"),rhs=gid.AD[2]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=gid.ValenciavU[0]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=gid.ValenciavU[1]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=gid.ValenciavU[2]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","BU0"),rhs=gid.BU[0]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","BU1"),rhs=gid.BU[1]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","BU2"),rhs=gid.BU[2]),\
                   lhrh(lhs=gri.gfaccess("out_gfs","psi6Phi"),rhs=sp.sympify(0))\
                  ]

outCfunction(
    outfile  = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *out_gfs",
    body     = fin.FD_outputC("returnstring",values_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints,Read_xxs")

```

    Output C function initial_data() to file GiRaFFE_staggered_standalone_Ccodes/initial_data.h

# Step 4: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \[Back to [top](#toc)\]
$$\label{cparams}$$

Based on declared NRPy+ Cparameters, first we generate `declare_Cparameters_struct.h`, `set_Cparameters_default.h`, and `set_Cparameters[-SIMD].h`.

Then we output `free_parameters.h`, which sets initial data parameters, as well as grid domain & reference metric parameters, applying `domain_size` and `sinh_width`/`SymTP_bScale` (if applicable) as set above

```python
# Step 3.e: Output C codes needed for declaring and setting Cparameters; also set free_parameters.h
# Step 3.e.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))

# Step 3.e.ii: Set free_parameters.h
with open(os.path.join(Ccodesdir,"free_parameters.h"),"w") as file:
    file.write("""// Override parameter defaults with values based on command line arguments and NGHOSTS.
params.Nxx0 = atoi(argv[1]);
params.Nxx1 = atoi(argv[2]);
params.Nxx2 = atoi(argv[3]);
params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;
// Step 0d: Set up space and time coordinates
// Step 0d.i: Declare \Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:
const REAL xxmin[3] = {-1.35,-0.09,-0.09};
const REAL xxmax[3] = { 1.65, 0.11, 0.11};

params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0+1);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1+1);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2+1);
printf("dxx0,dxx1,dxx2 = %.5e,%.5e,%.5e\\n",params.dxx0,params.dxx1,params.dxx2);
params.invdx0 = 1.0 / params.dxx0;
params.invdx1 = 1.0 / params.dxx1;
params.invdx2 = 1.0 / params.dxx2;

const int poison_grids = 0;
// Standard GRFFE parameters:
params.GAMMA_SPEED_LIMIT = 2000.0;
params.diss_strength = 0.1;
""")
if initial_data=="ExactWald":
    with open(os.path.join(out_dir,"free_parameters.h"),"a") as file:
        file.write("""params.r0 = 0.4;
params.a = 0.0;
""")

```

# Step 4: Set up boundary condition functions for chosen singular, curvilinear coordinate system \[Back to [top](#toc)\]
$$\label{bc_functs}$$

Next apply singular, curvilinear coordinate boundary conditions [as documented in the corresponding NRPy+ tutorial notebook](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb)

...But, for the moment, we're actually just using this because it writes the file `gridfunction_defines.h`.

```python
import CurviBoundaryConditions.CurviBoundaryConditions as cbcs
cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,"boundary_conditions/"),Cparamspath=os.path.join("../"),enable_copy_of_static_Ccodes=False)
```

    Wrote to file "GiRaFFE_staggered_standalone_Ccodes/boundary_conditions/parity_conditions_symbolic_dot_products.h"
    Evolved parity: ( AD0:1, AD1:2, AD2:3, StildeD0:1, StildeD1:2, StildeD2:3,
        psi6Phi:0 )

    AuxEvol parity: ( BU0:1, BU1:2, BU2:3, B_lU0:1, B_lU1:2, B_lU2:3, B_rU0:1,
        B_rU1:2, B_rU2:3, BstaggerU0:1, BstaggerU1:2, BstaggerU2:3,
        Bstagger_lU0:1, Bstagger_lU1:2, Bstagger_lU2:3, Bstagger_rU0:1,
        Bstagger_rU1:2, Bstagger_rU2:3, Stilde_flux_HLLED0:1,
        Stilde_flux_HLLED1:2, Stilde_flux_HLLED2:3, ValenciavU0:1,
        ValenciavU1:2, ValenciavU2:3, Valenciav_lU0:1, Valenciav_lU1:2,
        Valenciav_lU2:3, Valenciav_llU0:1, Valenciav_llU1:2, Valenciav_llU2:3,
        Valenciav_lrU0:1, Valenciav_lrU1:2, Valenciav_lrU2:3, Valenciav_rU0:1,
        Valenciav_rU1:2, Valenciav_rU2:3, Valenciav_rlU0:1, Valenciav_rlU1:2,
        Valenciav_rlU2:3, Valenciav_rrU0:1, Valenciav_rrU1:2, Valenciav_rrU2:3,
        alpha:0, alpha_face:0, betaU0:1, betaU1:2, betaU2:3, beta_faceU0:1,
        beta_faceU1:2, beta_faceU2:3, cmax_x:0, cmax_y:0, cmax_z:0, cmin_x:0,
        cmin_y:0, cmin_z:0, gammaDD00:4, gammaDD01:5, gammaDD02:6, gammaDD11:7,
        gammaDD12:8, gammaDD22:9, gammaUU00:4, gammaUU01:5, gammaUU02:6,
        gammaUU11:7, gammaUU12:8, gammaUU22:9, gamma_faceDD00:4,
        gamma_faceDD01:5, gamma_faceDD02:6, gamma_faceDD11:7, gamma_faceDD12:8,
        gamma_faceDD22:9, gamma_faceUU00:4, gamma_faceUU01:5, gamma_faceUU02:6,
        gamma_faceUU11:7, gamma_faceUU12:8, gamma_faceUU22:9, phi:0,
        phi_face:0, psi6_temp:0, psi6center:0 )
    Wrote to file "GiRaFFE_staggered_standalone_Ccodes/boundary_conditions/EigenCoord_Cart_to_xx.h"

# Step 5: `GiRaFFE_NRPy_standalone.c`: The Main C Code \[Back to [top](#toc)\]
$$\label{mainc}$$

```python
# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER),
#          and set the CFL_FACTOR (which can be overwritten at the command line)

with open(os.path.join(Ccodesdir,"GiRaFFE_NRPy_REAL__NGHOSTS__CFL_FACTOR.h"), "w") as file:
    file.write("""
// Part P0.a: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER
#define NGHOSTS """+str(3)+"""
#define NGHOSTS_A2B """+str(2)+"""
// Part P0.b: Set the numerical precision (REAL) to double, ensuring all floating point
//            numbers are stored to at least ~16 significant digits
#define REAL """+REAL+"""
// Part P0.c: Set the CFL Factor. Can be overwritten at command line.
REAL CFL_FACTOR = """+str(default_CFL_FACTOR)+";")
```

```python
%%writefile $Ccodesdir/GiRaFFE_NRPy_standalone.c
// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.
#include "GiRaFFE_NRPy_REAL__NGHOSTS__CFL_FACTOR.h"

#include "declare_Cparameters_struct.h"

const int NSKIP_1D_OUTPUT = 1;

// Step P1: Import needed header files
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "time.h"
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#ifndef M_PI
#define M_PI 3.141592653589793238462643383279502884L
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2 0.707106781186547524400844362104849039L
#endif

// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of
//           data in a 1D array. In this case, consecutive values of "i"
//           (all other indices held to a fixed value) are consecutive in memory, where
//           consecutive values of "j" (fixing all other indices) are separated by
//           Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of
//           "k" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )
#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )
#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \
  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)
#define LOOP_ALL_GFS_GPS(ii) _Pragma("omp parallel for") \
  for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)

// Step P3: Set gridfunction macros
#include "boundary_conditions/gridfunction_defines.h"

// Step P4: Include the RHS, BC, and primitive recovery functions
#include "GiRaFFE_NRPy_Main_Driver.h"

// Step P5: Include the initial data functions
#include "set_initial_spacetime_metric_data.h"
#include "initial_data.h"
// main() function:
// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates
// Step 1: Set up scalar wave initial data
// Step 2: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,
//         applying quadratic extrapolation outer boundary conditions.
// Step 3: Output relative error between numerical and exact solution.
// Step 4: Free all allocated memory
int main(int argc, const char *argv[]) {
    paramstruct params;
#include "set_Cparameters_default.h"
    // Step 0a: Read command-line input, error out if nonconformant
    if(argc != 4 || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < NGHOSTS) {
      printf("Error: Expected three command-line arguments: ./GiRaFFE_NRPy_standalone [Nx] [Ny] [Nz],\n");
      printf("where Nx is the number of grid points in the x direction, and so forth.\n");
      printf("Nx,Ny,Nz MUST BE larger than NGHOSTS (= %d)\n",NGHOSTS);
      exit(1);
    }
    // Step 0c: Set free parameters, overwriting Cparameters defaults
    //          by hand or with command-line input, as desired.
#include "free_parameters.h"
#include "set_Cparameters-nopointer.h"

    //          ... and then set up the numerical grid structure in time:
    const REAL t_final = 0.5;
    const REAL CFL_FACTOR = 0.5; // Set the CFL Factor

    // Step 0c: Allocate memory for gridfunctions
    const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
    // Step 0k: Allocate memory for gridfunctions
#include "MoLtimestepping/RK_Allocate_Memory.h"
    REAL *restrict auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);
    REAL *evol_gfs_exact = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot);
    REAL *auxevol_gfs_exact  = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);

    // For debugging, it can be useful to set everything to NaN initially.
    if(poison_grids) {
        for(int ii=0;ii<NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
            y_n_gfs[ii] = 1.0/0.0;
            y_nplus1_running_total_gfs[ii] = 1.0/0.0;
            //k_odd_gfs[ii] = 1.0/0.0;
            //k_even_gfs[ii] = 1.0/0.0;
            diagnostic_output_gfs[ii] = 1.0/0.0;
            evol_gfs_exact[ii] = 1.0/0.0;
        }
        for(int ii=0;ii<NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
            auxevol_gfs[ii] = 1.0/0.0;
            auxevol_gfs_exact[ii] = 1.0/0.0;
        }
    }

    // Step 0d: Set up coordinates: Set dx, and then dt based on dx_min and CFL condition
    // This is probably already defined above, but just in case...
    #ifndef MIN
    #define MIN(A, B) ( ((A) < (B)) ? (A) : (B) )
    #endif
    REAL dt = CFL_FACTOR * MIN(dxx0,MIN(dxx1,dxx2)); // CFL condition
    int Nt = (int)(t_final / dt + 0.5); // The number of points in time.
                                        //Add 0.5 to account for C rounding down integers.

    // Step 0e: Set up cell-centered Cartesian coordinate grids
    REAL *xx[3];
    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);
    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);
    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);
    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + (j-NGHOSTS+1)*dxx0;
    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS+1)*dxx1;
    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS+1)*dxx2;

    // Step 1: Set up initial data to be exact solution at time=0:
    REAL time = 0.0;
    set_initial_spacetime_metric_data(&params,xx,auxevol_gfs);
    initial_data(&params,xx,auxevol_gfs,y_n_gfs);

    /*
    // Code to perturb the initial data:
    for(int ii=0;ii<NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
        y_n_gfs[ii] += 1.0e-15;
        y_nplus1_running_total_gfs[ii] += 1.0e-15;
        //k_odd_gfs[ii] = 1.0/0.0;
        //k_even_gfs[ii] = 1.0/0.0;
        diagnostic_output_gfs[ii] += 1.0e-15;
        evol_gfs_exact[ii] += 1.0e-15;
    }
    for(int ii=0;ii<NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot;ii++) {
        auxevol_gfs[ii] += 1.0e-15;
        auxevol_gfs_exact[ii] += 1.0e-15;
    }
    */

    // Fill in  the remaining quantities
    GiRaFFE_compute_B_and_Bstagger_from_A(&params,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD00GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD01GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD02GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD11GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD12GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*GAMMADD22GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*PSI6_TEMPGF, /* Temporary storage,overwritten */
                                          y_n_gfs+Nxx_plus_2NGHOSTS_tot*AD0GF,
                                          y_n_gfs+Nxx_plus_2NGHOSTS_tot*AD1GF,
                                          y_n_gfs+Nxx_plus_2NGHOSTS_tot*AD2GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BU0GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BU1GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BU2GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BSTAGGERU0GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BSTAGGERU1GF,
                                          auxevol_gfs+Nxx_plus_2NGHOSTS_tot*BSTAGGERU2GF);
    //override_BU_with_old_GiRaFFE(&params,auxevol_gfs,0);
    GiRaFFE_NRPy_prims_to_cons(&params,auxevol_gfs,y_n_gfs);
    // Extra stack, useful for debugging:
    GiRaFFE_NRPy_cons_to_prims(&params,xx,auxevol_gfs,y_n_gfs);
    //GiRaFFE_NRPy_prims_to_cons(&params,auxevol_gfs,y_n_gfs);
    //GiRaFFE_NRPy_cons_to_prims(&params,xx,auxevol_gfs,y_n_gfs);
    //GiRaFFE_NRPy_prims_to_cons(&params,auxevol_gfs,y_n_gfs);
    //GiRaFFE_NRPy_cons_to_prims(&params,xx,auxevol_gfs,y_n_gfs);

    for(int n=0;n<=Nt;n++) { // Main loop to progress forward in time.
        // Step 1a: Set current time to correct value & compute exact solution
        time = ((REAL)n)*dt;

        /* Step 2: Validation: Output relative error between numerical and exact solution, */
        if((n)%NSKIP_1D_OUTPUT ==0) {
            // Step 2c: Output relative error between exact & numerical at center of grid.
            const int i0mid=Nxx_plus_2NGHOSTS0/2;
            const int i1mid=Nxx_plus_2NGHOSTS1/2;
            const int i2mid=Nxx_plus_2NGHOSTS2/2;
            char filename[100];
            sprintf(filename,"out%d-%08d.txt",Nxx0,n);
            FILE *out2D = fopen(filename, "w");
            for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
                const int idx = IDX3S(i0,i1mid,i2mid);
                fprintf(out2D,"%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",
                        xx[0][i0],
                        auxevol_gfs[IDX4ptS(BU0GF,idx)],auxevol_gfs[IDX4ptS(BU1GF,idx)],auxevol_gfs[IDX4ptS(BU2GF,idx)],
                        y_n_gfs[IDX4ptS(AD0GF,idx)],y_n_gfs[IDX4ptS(AD1GF,idx)],y_n_gfs[IDX4ptS(AD2GF,idx)],
                        y_n_gfs[IDX4ptS(STILDED0GF,idx)],y_n_gfs[IDX4ptS(STILDED1GF,idx)],y_n_gfs[IDX4ptS(STILDED2GF,idx)],
                        auxevol_gfs[IDX4ptS(VALENCIAVU0GF,idx)],auxevol_gfs[IDX4ptS(VALENCIAVU1GF,idx)],auxevol_gfs[IDX4ptS(VALENCIAVU2GF,idx)],
                        y_n_gfs[IDX4ptS(PSI6PHIGF,idx)],
                         xx[1][i1mid], xx[2][i2mid]);
            }
            fclose(out2D);
            // For convergence testing, we'll shift the grid x -> x-1 and output initial data again, giving the exact solution.
            LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,1,0,1) {
                xx[0][i0] += -time;
            }
            set_initial_spacetime_metric_data(&params,xx,auxevol_gfs_exact);
            initial_data(&params,xx,auxevol_gfs_exact,evol_gfs_exact);
            // Fill in  the remaining quantities
            //driver_A_to_B(&params,evol_gfs_exact,auxevol_gfs_exact);
            GiRaFFE_NRPy_prims_to_cons(&params,auxevol_gfs_exact,evol_gfs_exact);
            // And now, we'll set the grid back to rights.
            LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,1,0,1) {
                xx[0][i0] -= -time;
            }
            sprintf(filename,"out%d-%08d_exact.txt",Nxx0,n);
            FILE *out2D_exact = fopen(filename, "w");
            for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
                const int idx = IDX3S(i0,i1mid,i2mid);
                fprintf(out2D_exact,"%.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",
                        xx[0][i0],
                        auxevol_gfs_exact[IDX4ptS(BU0GF,idx)],auxevol_gfs_exact[IDX4ptS(BU1GF,idx)],auxevol_gfs_exact[IDX4ptS(BU2GF,idx)],
                        evol_gfs_exact[IDX4ptS(AD0GF,idx)],evol_gfs_exact[IDX4ptS(AD1GF,idx)],evol_gfs_exact[IDX4ptS(AD2GF,idx)],
                        evol_gfs_exact[IDX4ptS(STILDED0GF,idx)],evol_gfs_exact[IDX4ptS(STILDED1GF,idx)],evol_gfs_exact[IDX4ptS(STILDED2GF,idx)],
                        auxevol_gfs_exact[IDX4ptS(VALENCIAVU0GF,idx)],auxevol_gfs_exact[IDX4ptS(VALENCIAVU1GF,idx)],auxevol_gfs_exact[IDX4ptS(VALENCIAVU2GF,idx)],
                        evol_gfs_exact[IDX4ptS(PSI6PHIGF,idx)]);
            }
            fclose(out2D_exact);

        }

        // Step 3.b: Step forward one timestep (t -> t+dt) in time using
        //           chosen RK-like MoL timestepping algorithm
#include "MoLtimestepping/RK_MoL.h"
    } // End main loop to progress forward in time.

    // Step 4: Free all allocated memory
#include "MoLtimestepping/RK_Free_Memory.h"
    free(auxevol_gfs);
    free(auxevol_gfs_exact);
    free(evol_gfs_exact);
    for(int i=0;i<3;i++) free(xx[i]);
    return 0;
}
```

    Writing GiRaFFE_staggered_standalone_Ccodes//GiRaFFE_NRPy_standalone.c

```python
cmd.C_compile(os.path.join(Ccodesdir,"GiRaFFE_NRPy_standalone.c"),
              os.path.join(Ccodesdir,"output","GiRaFFE_NRPy_standalone"),compile_mode="safe")
# !gcc -g -O2 -fopenmp GiRaFFE_standalone_Ccodes/GiRaFFE_NRPy_standalone.c -o GiRaFFE_NRPy_standalone -lm

# Change to output directory
os.chdir(outdir)
# Clean up existing output files
cmd.delete_existing_files("out*.txt")
cmd.delete_existing_files("out*.png")
# cmd.Execute(os.path.join(Ccodesdir,"output","GiRaFFE_NRPy_standalone"), "640 16 16", os.path.join(outdir,"out640.txt"))
# !taskset -c 0-7 ./GiRaFFE_NRPy_standalone 119 7 7
cmd.Execute("GiRaFFE_NRPy_standalone", "149 9 9","out150.txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "1280, 32, 32","out119.txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "119 119 119","out119.txt")
# cmd.Execute("GiRaFFE_NRPy_standalone", "239 15 15","out239.txt")
# !OMP_NUM_THREADS=1 valgrind --track-origins=yes -v ./GiRaFFE_NRPy_standalone 1280 32 32
# Return to root directory
os.chdir(os.path.join("../../"))
```

    Compiling executable...
    (EXEC): Executing `gcc -O2 -g -fopenmp GiRaFFE_staggered_standalone_Ccodes/GiRaFFE_NRPy_standalone.c -o GiRaFFE_staggered_standalone_Ccodes/output/GiRaFFE_NRPy_standalone -lm`...
    (BENCH): Finished executing in 2.8143463134765625 seconds.
    Finished compilation.
    (EXEC): Executing `taskset -c 0,1,2,3,4,5 ./GiRaFFE_NRPy_standalone 149 9 9`...
    (BENCH): Finished executing in 6.624955415725708 seconds.

Now, we will load the data generated by the simulation and plot it in order to test for convergence.

```python
import numpy as np
import matplotlib.pyplot as plt

#          0   1     2     3     4     5     6     7          8          9          10            11            12             13          14
labels = ["x","BU0","BU1","BU2","AD0","AD1","AD2","StildeD0","StildeD1","StildeD2","ValenciavU0","ValenciavU1","ValenciavU2", "psi6Phi", "time"]
Data_numer = np.loadtxt(os.path.join(Ccodesdir,"output","out149-00000000.txt"))
Data_exact = np.loadtxt(os.path.join(Ccodesdir,"output","out149-00000000_exact.txt"))

column = 8
plt.figure()
# # plt.plot(Data_exact[2:-2,0],np.log2(np.absolute((Data_numer[2:-2,column]-Data_exact[2:-2,column])/\
# #                                     (Data_num_2[2:-2:2,column]-Data_exa_2[2:-2:2,column]))),'.')
plt.plot(Data_exact[:,0],Data_exact[:,column],'k-',label="Exact")
plt.scatter(Data_numer[:,0],Data_numer[:,column],s=4,label="Approximate")
# plt.xlim(-0.0,1.0)
# # plt.ylim(-1.0,5.0)
# # plt.ylim(-0.0005,0.0005)
plt.xlabel(labels[0])
plt.ylabel(labels[column])
plt.title('NRPy+ GiRaFFE')
plt.legend()
plt.show()

plt.plot(Data_exact[:,0],np.log10(np.abs((Data_exact[:,column]-Data_numer[:,column])/(Data_exact[:,column]))))
plt.xlabel('x')
plt.ylabel('Relative Error')
plt.title('NRPy+ GiRaFFE')
plt.show()
# print(Data_numer[:,-2])
```

```python
# ETK_data = np.loadtxt('/home/terrence/Projects/ETK/ET_2021_01/Cactus/exe/ABE-GiRaFFEfood_1D_FastWave-GFood/giraffe-grmhd_conservatives.x.asc')
# # ETK_data = np.loadtxt('/home/terrence/Projects/ETK/ET_2021_01/Cactus/exe/ABE-GiRaFFEfood_1D_FastWave-GFood/giraffe-em_az.x.asc')
# n = 0
# start = n*155
# stop = (n+1)*155
# x = ETK_data[start:stop, 9]
# y = ETK_data[start:stop, 13]

# plt.scatter(x, y, s=1, label='ETK')
# plt.scatter(Data_numer[:,0],Data_numer[:,column],s=1, label='NRPy')
# plt.ylabel(labels[column])
# plt.xlabel('x')
# plt.title('NRPy+ vs ETK')
# plt.legend()
# plt.show()

# E_rel_diff = np.log10(np.abs(0.5*(y-Data_numer[:,column])/(y+Data_numer[:,column])))
# plt.plot(x,E_rel_diff)
# plt.title("Relative Difference")
# plt.xlabel('x')
# plt.ylabel('E_rel')
# plt.show()
```

```python
# ETK_data = np.loadtxt('/home/terrence/Projects/ETK/ET_2021_01/Cactus/exe/ABE-GiRaFFEfood_1D_FastWave-GFood/giraffe-grmhd_conservatives.x.asc')
# # ETK_data = np.loadtxt('/home/terrence/Projects/ETK/ET_2021_01/Cactus/exe/ABE-GiRaFFEfood_1D_FastWave-GFood/giraffe-em_az.x.asc')
# n = 0
# start = n*155
# stop = (n+1)*155
# x = ETK_data[start:stop, 9]
# y = ETK_data[start:stop, 13]

# plt.scatter(x, y, s=1, label='ETK')
# plt.plot(Data_exact[:,0],Data_exact[:,column],'k-', label='NRPy-Exact')
# plt.ylabel(labels[column])
# plt.xlabel('x')
# plt.title('ETK vs NRPy-Exact')
# plt.legend()
# plt.show()

# E_rel_diff = np.log10(np.abs((y-Data_exact[:,column])/(Data_exact[:,column])))
# plt.plot(x,E_rel_diff)
# plt.title("Relative Error")
# plt.xlabel('x')
# plt.ylabel('E_rel')
# plt.show()
```

```python
# for i in range(81):
#     start = i*133
#     stop = (i+1)*133
#     x = ETK_data[start:stop, 9]
#     y = ETK_data[start:stop, 13]
#     plt.scatter(x,y, s=1)
#     plt.show()
```

```python
# plt.plot(ETK_data[-155:,9],np.log10(np.abs(ETK_data[-155:,9] - Data_numer[:,0])))
# plt.title('Difference between NRPy and ETK grids')
# plt.xlabel('x')
# plt.ylabel('Difference')
# plt.show()
```

This code will create an animation of the wave over time.

```python
# import matplotlib.pyplot as plt
from matplotlib.pyplot import savefig
from IPython.display import HTML
import matplotlib.image as mgimg

import glob
import sys
from matplotlib import animation

cmd.delete_existing_files("out149-00*.png")
globby = glob.glob(os.path.join(Ccodesdir,'output','out149-00*.txt'))
file_list = []
for x in sorted(globby):
    file_list.append(x)

number_of_files = int(len(file_list)/2)

for timestep in range(number_of_files):
    fig = plt.figure()
    numer_filename = file_list[2*timestep]
    exact_filename = file_list[2*timestep+1]
    Numer = np.loadtxt(numer_filename)
    Exact = np.loadtxt(exact_filename)

    plt.title("Fast Wave")
    plt.xlabel("x")
    plt.ylabel(labels[column])
#     plt.xlim(-0.5,0.5)
#     plt.ylim(1.0,1.7)

    plt.scatter(Numer[3:-3,0],Numer[3:-3,column],s=1,label="Numerical")
    plt.plot(Exact[3:-3,0],Exact[3:-3,column],'k-',label="Exact")
    plt.legend()
    savefig(numer_filename+".png",dpi=150)
    plt.close(fig)
    sys.stdout.write("%c[2K" % 27)
    sys.stdout.write("Processing file "+numer_filename+"\r")
    sys.stdout.flush()
```

```python
## VISUALIZATION ANIMATION, PART 2: Combine PNGs to generate movie ##
# https://stackoverflow.com/questions/14908576/how-to-remove-frame-from-matplotlib-pyplot-figure-vs-matplotlib-figure-frame
# https://stackoverflow.com/questions/23176161/animating-pngs-in-matplotlib-using-artistanimation
# !rm -f GiRaFFE_NRPy-1D_tests.mp4
cmd.delete_existing_files("GiRaFFE_NRPy-1D_tests.mp4")

fig = plt.figure(frameon=False)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')

myimages = []

for i in range(number_of_files):
    img = mgimg.imread(file_list[2*i]+".png")
    imgplot = plt.imshow(img)
    myimages.append([imgplot])

ani = animation.ArtistAnimation(fig, myimages, interval=100,  repeat_delay=1000)
plt.close()
ani.save('GiRaFFE_NRPy-1D_tests.mp4', fps=5,dpi=150)
```

```python
%%HTML

```

```python
import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish-GiRaFFE_NRPy-1D_tests-staggered_FastWave",location_of_template_file=os.path.join(".."))
```

    Created Tutorial-Start_to_Finish-GiRaFFE_NRPy-1D_tests-
        staggered_FastWave.tex, and compiled LaTeX file to PDF file Tutorial-
        Start_to_Finish-GiRaFFE_NRPy-1D_tests-staggered_FastWave.pdf

max_stars_repo_name: fedelopezar/nrpytutorial

max_issues_repo_name: fedelopezar/nrpytutorial

max_forks_repo_name: fedelopezar/nrpytutorial


text:
# 数値計算

単振り子の運動について考える．

まずは数値計算に必要なライブラリのインポートから．[NumPy](http://www.numpy.org/)は科学計算を行うための基本的なライブラリである．[SciPy](https://www.scipy.org/)も科学計算を行う上で用いられるライブラリであるが，NumPyを用いて数学や工学などへの応用性のあるクラスが多い．

[math](https://docs.python.jp/3/library/math.html) はいつでも使えるデフォルトの数学関数系のライブラリである．

```python
import numpy as np
from scipy.integrate import odeint
from math import sin
```

## `math.sin` と `numpy.sin` について
`numpy` ライブラリの `sin` 関数を使わず， `math` ライブラリから `sin` 関数をインポートしたのは，計算を高速化させるためである．`math`ライブラリで定義される正弦関数と`numpy`ライブラリで定義される正弦関数には違いがある．次の例で示すように，スカラー値の正弦の計算は `numpy.sin` よりも`math.sin` の方が早い．

```python
import time
import math
import numpy

''' Sine function of math library '''
start = time.time()
for i in range(1000000):
    math.sin(i)
elapsed_time = time.time() - start
print ('Sine calculation by math: {0:f} [sec]'.format(elapsed_time))

''' Sine function of numpy library '''
start = time.time()
for i in range(1000000):
    numpy.sin(i)
elapsed_time = time.time() - start
print ('Sine calculation by numpy: {0:f} [sec]'.format(elapsed_time))
```

    Sine calculation by math: 0.193746 [sec]
    Sine calculation by numpy: 1.112725 [sec]

## `scipy.integrate.odeint` を用いた常微分方程式の計算

まず，運動方程式で用いる定数 ($m$，$l$，$g$) を定義しておく．

```python
''' constants '''
m = 1 # mass of the pendulum [kg]
l = 1 # length of the pendulum [m]
g = 10 # Gravitational acceleration [m/s^2]
```

次に，シミュレーションを行う時間に関する定数を定義しておく． `t_fps` はグラフの描画やアニメーションの生成のために1秒を何分割するかを意味するパラメータである．

```python
''' time setting '''
t_end = 10 # simulation time [s]
t_fps = 50 # frame per second. This value means smoothness of produced graph and animation
t_step = 1/t_fps
t = np.arange(0, t_end, t_step)
```

$\theta$自身とその一階微分である$\dot{\theta}$を縦に並べた列ベクトルをつくり，これを$s$とする．（GiuHub上で行列が上手く表示されないのは仕方のないことなのかな？）
\begin{align}
  s =
  \begin{bmatrix}
    \theta \\ \dot{\theta}
  \end{bmatrix}
\end{align}

常微分方程式$\dot{s} = f(s)$を解くために，初期値 $s_\mathrm{init}$ を与える．

```python
''' initial value '''
theta_init = 0 # initial value of theta [rad]
dtheta_init = 1 # initial value of dot theta [rad/s]
s_init = np.r_[theta_init, dtheta_init]
np.r_[s_init, s_init]
```

    array([0, 1, 0, 1])

そして，$s$の時間の一階微分を式で表す．ここで，$\ddot{\theta}$は先ほど求めた運動方程式によって，$\ddot{\theta}$を含まない形で表すことができる．
\begin{align}
  \frac{d}{dt} s =
  \begin{bmatrix}
    \dot{\theta} \\ \ddot{\theta}
  \end{bmatrix} =
  \begin{bmatrix}
    \dot{\theta} \\ -\frac{g}{l}\sin\theta
  \end{bmatrix}
\end{align}
こうして得られた常微分方程式は，数値計算に適した形となっている．具体的には，$\theta = $ `s[0]` ，$\dot{\theta} = $ `s[1]` とすることで，

```
def odefunc(s, t):
    theta = s[0]
    dtheta = s[1]
    ddtheta = -g/l*sin(theta) # <- This is the equation of motion
    return np.r_[dtheta, ddtheta]
```
のようにして常微分方程式をPythonの関数として記述したうえで

```
s = odeint(odefunc, s_init, t)
```
を実行することで常微分方程式を数値的に解くことができる．ただし [`odeint`](https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html) は [SciPy](https://www.scipy.org/) ライブラリの `integrate` クラスで定義される，常微分方程式を解く関数である．

```python
def odefunc(s, t):
    theta = s[0]
    dtheta = s[1]
    ddtheta = -g/l*sin(theta) # <- Equation of motion
    return np.r_[dtheta, ddtheta]

s = odeint(odefunc, s_init, t)
print('ODE calculation finished.')
```

    ODE calculation finished.

## 注意点：`odeint` 中では **行ベクトル** で解を扱う

振り子のモデリングにおいては，解$s = \begin{bmatrix} \theta & \dot{\theta} \end{bmatrix}^\mathrm{T}$は **列ベクトル** として扱っている．しかし，`odeint` 関数の引数の `s_init` や `odefunc`　の返り値は **行ベクトル** でなくてはならない．

max_stars_repo_name: yfur/basic-mechanics-python

max_issues_repo_name: yfur/basic-mechanics-python

max_forks_repo_name: yfur/basic-mechanics-python