#!/usr/bin/env python
# coding: utf-8

# This tutorial describes two fundamental user APIs:
# 
# * Operator.apply
# * Operator.arguments
# 
# We will use a trivial `Operator` that, at each time step, increments by 1 all points in the physical domain.

# In[1]:


from devito import Grid, TimeFunction, Eq, Operator

grid = Grid(shape=(4, 4))
u = TimeFunction(name='u', grid=grid, save=3)
op = Operator(Eq(u.forward, u + 1))


# To run `op`, we have to "`apply`" it.

# In[2]:


#NBVAL_IGNORE_OUTPUT
summary = op.apply()


# Under-the-hood, some code has been generated (`print(op)` to display the generated code), JIT-compiled, and executed. Since no additional arguments have been passed, `op` has used `u` as input. We can verify that the content of `u.data` is as expected

# In[3]:


u.dimensions, u.shape


# In[4]:


u.data


# In particular, we observe that:
# 
# * `u` has size 3 along the time dimension, since it was built with `save=3`. Therefore `op` could only execute 2 timesteps, namely time=0 and time=1; given `Eq(u.forward, u + 1)`, executing time=2 would cause out-of-bounds access errors. Devito figures this out automatically and sets appropriate minimum and maximum iteration points.
# * All 16 points in each timeslice of the 4x4 `Grid` have been computed.
# 
# To access all default arguments used by `op` *without* running the `Operator`, one can run

# In[5]:


#NBVAL_IGNORE_OUTPUT
op.arguments()


# `'u'` stores a pointer to the allocated data; `'timers'` stores a pointer to a data structure used for C-level performance profiling.
# 
# One may want to replace some of these default arguments. For example, we could increase the minimum iteration point along the spatial Dimensions `x` and `y`, and execute only the very first timestep:

# In[6]:


#NBVAL_IGNORE_OUTPUT
u.data[:] = 0.  # Explicit reset to initial value
summary = op.apply(x_m=2, y_m=2, time_M=0)


# We look again at the computed data to convince ourselves that everything went as intended to go

# In[7]:


u.data


# Given a generic `Dimension` `d`, the naming convention is such that:
# 
# * `d_m` is the minimum iteration point
# * `d_M` is the maximum iteration point
# 
# Hence, `op.apply(..., d_m=4, d_M=7, ...)` will run `op` in the compact interval `[4, 7]` along `d`. For historical reasons, `d=...` aliases to `d_M=...`; in many Devito examples it happens to see `op.apply(..., time=10, ...)` -- this is just equivalent to `op.apply(..., time_M=10, ...)`.
# 
# If we try to specify an invalid iteration extreme, Devito will raise an exception.

# In[8]:


from devito.exceptions import InvalidArgument
try:
    op.apply(time_M=2)
except InvalidArgument as e:
    print(e)


# The same `Operator` can be applied to a different `TimeFunction`. For example:

# In[9]:


#NBVAL_IGNORE_OUTPUT
u2 = TimeFunction(name='u', grid=grid, save=5)
summary = op.apply(u=u2)


# In[10]:


u2.data


# Note that this was the third call to `op.apply`, but code generation and JIT-compilation only occurred upon the very first call.
# 
# There is one relevant case in which the maximum iteration point along the time dimension must be specified -- whenever `save` is unset, as in such a case the `Operator` wouldn't know for how many iterations to run.

# In[11]:


v = TimeFunction(name='v', grid=grid)
op2 = Operator(Eq(v.forward, v + 1))
try:
    op2.apply()
except ValueError as e:
    print(e)


# In[12]:


#NBVAL_IGNORE_OUTPUT
summary = op2.apply(time_M=4)


# In[13]:


v.data


# The `summary` variable can be inspected to retrieve performance metrics.

# In[14]:


#NBVAL_IGNORE_OUTPUT
summary


# We observe that basically all entries except for the execution time are fixed at 0. This is because by default Devito avoids to compute performance metrics, to minimize the processing time before returning control to the user (in complex `Operators`, the processing time to retrieve, for instance, the operation count or the memory footprint could be significant). To compute all performance metrics, a user could either export the environment variable `DEVITO_PROFILING` to `'advanced'` or change the profiling level programmatically *before* the `Operator` is constructed

# In[15]:


#NBVAL_IGNORE_OUTPUT
from devito import configuration
configuration['profiling'] = 'advanced'

op = Operator(Eq(u.forward, u*u + 1.))
op.apply()


# A `PerformanceSummary` will contain as many entries as "sections" in the generated code. Currently, there is no way to automatically tie a compiler-generated section to the user-provided `Eq`s (in general, there can be more than one `Eq` in a section). The only option is to look at the generated code and search for bodies of code wrapped within C comments such as
# ```
# /* Begin section0 */
# <code>
# /* End section0 \*/"
# ```
# For example

# In[16]:


# Uncomment me and search for /* Begin section0 */ ... /* End section0 */
# print(op)


# In the `PerformanceSummary`, associated to `section0` is a `PerfEntry`, whose entries represent:
# 
# * time: The time, in seconds, to execute the section.
# * gflopss: Performance of the section in Gigaflops per second.
# * gpointss: Performance of the section in Gigapoints per second.
# * oi: Operational intensity.
# * ops: Floating point operations at each (innermost loop) iteration.
# * itershape: The shape of the iteration space.