<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Use examples for numerical calculations
This jupyter notebook can be found at
Step1:  Numerical evaluations
See here for detailed instructions on how to turn sympy expressions into code
Step3:  Substitution of equations and values into equations
The easiest way is to define a dictionary with all variables we want to substitute as keys. We start with the default variables and then add more. First, however, we will define a function to display the contents of a dictionary
Step4:  We can substitute a range of equations into each other by using the custom function subs_eq
Step5:  We can also use subs_eq to substitute equations into each other and a dictionary with values. We will first add an entry for T_a into the dictionary and then substitute
Step6:  Evaluation of equations for long lists of variable sets
Substitution of variables into equations takes a lot of time if they need to be evaluated for a large number of variables. We can use theano to speed this up
Step7:  We will now create two long lists of values representing T_g and n_g respectively and show how long it takes to compute ideal gas law values.
Step8:  Both approaches give identical results, but theano_function makes it a lot faster.
Numerical solution
Some equations cannot be solved analytically for a given variable, e.g. eq_Nu_forced_all cannot be solved analytically for Re if Nu is given, so we can use numerical solvers instead
Step9:  Now applying to a long list of Nu-values
Step10:  We will now again use a theano function to make it faster. First we import optimize from scipy and preapre the theano_function
Step11:  Using theano and scipy makes it 2 orders of magnitude faster and the results are different only by 10$^{-10}$%!
Note, however, that scipy gets slowed down for large arrays, so it is more efficient to re-run it repreatedly with subsections of the arra
Step12:  We will now test that we can process Nuvals bit by bit and re-create it consistently
Step13:  Now we will run fsolve for portions of Nuvals bit by bit
Step14:  It is strange that resvals1 and resvals2 are different at all, but anyway, it is clear that slicing the data in relatively small portions is important to keep scipy.optimize.fsolve time-efficient.
Generate code from sympy expressions and execute
Need to install gfortran system-wide first!
Step15:  Use of autowrap made the calculation 3 orders of magnitude faster than substitution of values into the original expression!
Another way is to use binary_function
<ASSISTANT_TASK:>
Python Code:
from IPython.display import display
from sympy import init_printing, latex
init_printing() 
from sympy.printing import StrPrinter
StrPrinter._print_Quantity = lambda self, expr: str(expr.abbrev)    # displays short units (m instead of meter)
%run -i 'test_equation_definitions.py'
Explanation: Use examples for numerical calculations
This jupyter notebook can be found at:
https://github.com/environmentalscience/essm/blob/master/docs/examples/examples_numerics.ipynb
Below, we will import variable and equation defintions that were previously exported from api_features.ipynb
by running the file test_equation_definitions.py:
End of explanation
for eq in Equation.__registry__.keys():
        print(eq.definition.name + ': ' + str(eq))
Explanation: Numerical evaluations
See here for detailed instructions on how to turn sympy expressions into code: https://docs.sympy.org/latest/modules/codegen.html
We will first list all equations defined in this worksheet:
End of explanation
def print_dict(vdict, list_vars=None):
    Print values and units of variables in vdict.
    if not list_vars:
        list_vars = vdict.keys()
    for var1 in list_vars:
        unit1 = var1.definition.unit
        if unit1 == 1:
            unit1 = ''
        if vdict[var1] is not None:
            print('{0}: {1} {2}'.format(var1.name, str(vdict[var1]), str(unit1)))
vdict = Variable.__defaults__.copy()
print_dict(vdict)
Explanation: Substitution of equations and values into equations
The easiest way is to define a dictionary with all variables we want to substitute as keys. We start with the default variables and then add more. First, however, we will define a function to display the contents of a dictionary:
End of explanation
from essm.variables.utils import subs_eq
subs_eq(eq_Le, [eq_alphaa, eq_Dva])
Explanation: We can substitute a range of equations into each other by using the custom function subs_eq:
End of explanation
vdict[T_a] = 300.
subs_eq(eq_Le, [eq_alphaa, eq_Dva], vdict)
Explanation: We can also use subs_eq to substitute equations into each other and a dictionary with values. We will first add an entry for T_a into the dictionary and then substitute:
End of explanation
#import theano
from sympy.printing.theanocode import theano_function
import numpy as np
Explanation: Evaluation of equations for long lists of variable sets
Substitution of variables into equations takes a lot of time if they need to be evaluated for a large number of variables. We can use theano to speed this up:
End of explanation
npoints = 10000
xmin = 290.
xmax = 310.
Tvals = np.arange(xmin, xmax, (xmax - xmin)/npoints)
xmin = 0.1
xmax = 0.5
nvals = np.arange(xmin, xmax, (xmax-xmin)/npoints)
%%time
# looping
expr = eq_ideal_gas_law.rhs.subs(Variable.__defaults__)
resvals0 = []
for i in range(len(Tvals)):
    resvals0.append(expr.subs({T_g: Tvals[i], n_g: nvals[i]}))
%%time
# Using theano
f1 = theano_function([T_g, n_g], [eq_ideal_gas_law.rhs.subs(Variable.__defaults__)], dims={T_g:1, n_g:1})
resvals1 = f1(Tvals,nvals)
list(resvals0) == list(resvals1)
Explanation: We will now create two long lists of values representing T_g and n_g respectively and show how long it takes to compute ideal gas law values.
End of explanation
from sympy import nsolve
vdict = Variable.__defaults__.copy()
vdict[Pr] = 0.71
vdict[Re_c] = 3000.
vdict[Nu] = 1000.
expr = eq_Nu_forced_all.subs(vdict)
nsolve(expr, 1000.)
Explanation: Both approaches give identical results, but theano_function makes it a lot faster.
Numerical solution
Some equations cannot be solved analytically for a given variable, e.g. eq_Nu_forced_all cannot be solved analytically for Re if Nu is given, so we can use numerical solvers instead:
End of explanation
npoints = 100
xmin = 1000.
xmax = 1200.
Nuvals = np.arange(xmin, xmax, (xmax - xmin)/npoints)
%%time
# Solving for a range of Nu values
vdict = Variable.__defaults__.copy()
vdict[Pr] = 0.71
vdict[Re_c] = 3000.
resvals = []
for Nu1 in Nuvals:
    vdict[Nu] = Nu1
    resvals.append(nsolve(eq_Nu_forced_all.subs(vdict), 1000.))
Explanation: Now applying to a long list of Nu-values:
End of explanation
import scipy.optimize as sciopt
vdict = Variable.__defaults__.copy()
vdict[Pr] = 0.71
vdict[Re_c] = 3000.
expr = eq_Nu_forced_all.subs(vdict)
expr1 = expr.rhs - expr.lhs
fun_tf = theano_function([Re, Nu], [expr1], dims={Nu:1, Re:1})
x0vals = np.full(Nuvals.shape, fill_value=2000.) # array of same shape as Nuvals, with initial guess
%%time
# Solving for a range of Nu values
resvals1 = sciopt.fsolve(fun_tf, args=Nuvals, x0=x0vals)
np.mean(abs((resvals - resvals1)/resvals))
Explanation: We will now again use a theano function to make it faster. First we import optimize from scipy and preapre the theano_function:
End of explanation
npoints = 1000
xmin = 1000.
xmax = 1200.
Nuvals = np.arange(xmin, xmax, (xmax - xmin)/npoints)
x0vals = np.full(Nuvals.shape, fill_value=2000.)
%%time
# Solving for a range of Nu values
resvals1 = sciopt.fsolve(fun_tf, args=Nuvals, x0=x0vals)
Explanation: Using theano and scipy makes it 2 orders of magnitude faster and the results are different only by 10$^{-10}$%!
Note, however, that scipy gets slowed down for large arrays, so it is more efficient to re-run it repreatedly with subsections of the arra:
End of explanation
# Solving for a range of Nu values
imax = len(Nuvals)
i0 = 0
idiff = 100
i1 = i0
resvals2 = []
while i1 < imax - 1:
    i0 = i1    # note that resvals[0:2] + resvals[2:4] = resvals[0:4]
    i1 = min(i0+idiff, imax)
    resvals0 = Nuvals[i0:i1]
    resvals2 = np.append(resvals2,resvals0)
print(list(resvals2) == list(Nuvals))
Explanation: We will now test that we can process Nuvals bit by bit and re-create it consistently:
End of explanation
%%time
# Solving for a range of Nu values
imax = len(Nuvals)
i0 = 0
idiff = 100
i1 = i0
resvals2 = []
while i1 < imax - 1:
    i0 = i1    # note that resvals[0:2] + resvals[2:4] = resvals[0:4]
    i1 = min(i0+idiff, imax)
    resvals0 = sciopt.fsolve(fun_tf, args=Nuvals[i0:i1], x0=x0vals[i0:i1])
    resvals2 = np.append(resvals2,resvals0)
np.mean(abs((resvals1 - resvals2)/resvals1))
Explanation: Now we will run fsolve for portions of Nuvals bit by bit:
End of explanation
from sympy.utilities.autowrap import autowrap
from sympy import symbols
x, y, z = symbols('x y z')
expr = ((x - y + z)**(13)).expand()
autowrap_func = autowrap(expr)
%%time
autowrap_func(1, 4, 2)
%%time
expr.subs({x:1, y:4, z:2})
Explanation: It is strange that resvals1 and resvals2 are different at all, but anyway, it is clear that slicing the data in relatively small portions is important to keep scipy.optimize.fsolve time-efficient.
Generate code from sympy expressions and execute
Need to install gfortran system-wide first!
End of explanation
from sympy.utilities.autowrap import binary_function
f = binary_function('f', expr)
%%time
f(x,y,z).evalf(2, subs={x:1, y:4, z:2})
Explanation: Use of autowrap made the calculation 3 orders of magnitude faster than substitution of values into the original expression!
Another way is to use binary_function:
End of explanation
<END_TASK>