---
# System prepended metadata

title: A Frontend Specification for the Functional GT4Py Model (Deprecated)
tags: [shaping, frontend]

---

###### tags: `frontend` `shaping`

A Frontend Specification for the Functional GT4Py Model (Deprecated)
=======================================================
[TOC]

## Introduction

This document collects both the concepts exposed to GT4Py users to implement computations and the syntax used to represent the concepts.

GT4Py aims to to define an interface to specify computations between n-dimensional data arrays, called **Fields**, defined on top of a **Grid** (or *Mesh*), which is a discretized representation of a 3D physical space. Computations are defined locally in a subset of the grid points called the **Domain** of the computation.

The frontend component expose these concepts with a specific syntax to the users, who then define the data arrays, grids and computations that will be processed by the toolchain component to generate a valid implementation. The frontend is in charge of mapping these concepts to the IR concepts defined in the [Iterator IR spec](https://github.com/GridTools/concepts/wiki/Iterator-View) and thus, in general, both sets of concepts are very similar.


## Concepts

### Dimension
A *Dimension* is a just a *tag* or type for indices.

### Index
An integer number tagged with a [dimension](#Dimension), i.e. an element of an axis.

__Syntax proposal (WIP)__:
```
Index[Dim](int_idx)
```

### Axis
An *Axis* of a [dimension](#Dimension) *D* can be generally defined as an ordered set of [indices](#Index) of *D* (i.e. a spatially discretized version of a _dimension_). For this project, we will use a slightly more restricted definition of *Axis* as just a range (`(start, end, step)`) of index values.

### Position
A tuple of [indices](#Index) indicating a position on the grid / mesh.

__Syntax proposal (WIP)__:
```python=
Position[*Axis](*int_indices)
```

### Domain
Cartesian product of [axes](#Axis) (encoded as a tuple of *axes*).

__Syntax proposal (WIP)__:
```python=
Domain[*Axis]((i_start, i_stop), ...)
```

### Field
A mapping from [positions](#Position) to values of a certain datatype (e.g. int, float, tuple of int). GT4Py assumes that actual implementations of the *Field* concept use a contiguous memory buffer representation.

__Syntax proposal (WIP)__:
- Since `Field` is just a concept, the actual allocation interface could use the `Storage` class or other Array-like libraries with the `__gt_data_interface__`. TODO: revise the current implementation of `__gt_data_interface__` and adapt it to this proposal.
- If/When [PEP646](https://www.python.org/dev/peps/pep-0646/) gets approved the syntax for typing could be:
```python=
Axes = TypeVarTuple('Axes')
DTypeT = TypeVar('DTypeT')

class Field(Generic[DTypeT, *Axes]):
    ...

field: Field[DTypeT, *Axes]

# Alternative
class Space(Generic[*Axes]):
    ...
    
field: Field[Space[*Axes], DTypeT]
```
- Currently, typing could be done by something like:
```python=
Axis1T = TypeVar('Axis1T')
Axis2T = TypeVar('Axis2T')
Axis3T = TypeVar('Axis3T')
SpaceT = TypeVar('SpaceT')
...

DTypeT = TypeVar('DTypeT')

class Space1D(Generic[Axis1T]):
    ...
    
class Space2D(Generic[Axis1T, Axis2T]):
    ...

class Space3D(Generic[Axis1T, Axis2T, Axis3T]):
    ...

class Field(Generic[SpaceT, DTypeT]):
    ...
    
Field1D = Field[Space1D[Axis1T], DTypeT]
Field2D = Field[Space2D[Axis1T, Axis2T], DTypeT]
...

f: Field1D[VAxis, float]
```

### Connectivity

A *Connectivity* is a mapping from an [index](#Index) to a list of indices. The length of the list is either limited, e.g. for the connectivity from vertex to neighboring vertices, or fixed, e.g. from edge to vertex.

__Syntax proposal (WIP)__:

Each connectivity is its own type implementing the following protocol:
```python=
class Connectivity(Protocol[From_Axis, To_Axis]):
    max_neighbors: ClassVar[int]
    has_skip_values: ClassVar[bool]
   
    # TODO: Aliases
    LEFT = 0
    RIGHT = 1
    ...
```

### Local View

_(Alternative name: LocatedView? Accessor?)_

A local-view of a field is a representation of the value of the field at a specific position. Intuitively, it can be thought of as a tuple of a field together with a position, which can be converted to a value by index. It can be used as the actual value for any value operation and it additionally supports an extra operation to create another local-view placed on a different position, by mapping or shifting its current one. 

__Syntax proposal (WIP)__:

- Typing: `LocalView[PositionT, SpaceT, DType]`
- Arithmetic operations and tuple access: just as regular Python values
    ```python=
    new_value = view + 3.5 + tuple_view[2]
    ```
- Creating new views:
    + Mapping position using `Connectivity`: _Connectivity_ will map the view position (or the selected axis) to a tuple of indices that will be expanded in place.
        ```python=
        view: LocalView[...]
        V2E: Connectivity[...]
        
        # General syntax: View(Connectivity(NeighborId))
        new_view = view(V2E(1)) # new_view(position=V2E.map(view.position)[1])      
        # If there are named aliases in connectivity: View(Connectivity.NeighborAlias)
        new_view = view(V2E.SECOND)

        # Concatenation of mappings
        new_view = view(AB(0) | BC(2) | CD(1))  

        # Mapping only some components of current position (other components are left identical)
        new_view = view(AxisA >> AB(0))
        
        ```
    + Shifting one axis of the position using a literal offset:
        ```python=
        new_view = view(I + 3)  
        ```

### Operator

A pure functional _kernel_ computing a value (of any accepted data type, including tuples) from different [local-views](#Local-View) passed as arguments.

#### Local operator

An operator accepting point-wise [local-views](#Local-View).

__Syntax proposal (WIP)__:

- A Python function decorated with `@gt.local_operator(...)`
    - Keyword arguments: TODO
- Accepted data arguments are `LocalView` 
- If the current iteration position is needed, an extra keyword-only argument `position: Position` should be added to the signature.
    + Identification of the position argument?
        * [ ] by keyword name: `position`
        * [ ] by type hint: `Position` 
        * [ ] both: `position: Position` 
- Question: explicit or implicit `lift` and `deref` ??
- Accepted statements inside the body:
    - Assignments to scalar temporaries (using arithmetic and ternary operators)
    - _Shifting_/creation of new [views ](#Local-View)
    - Conditionals or ternary operators
        * [ ] Only ternary operators
            ```python=
            @local_operator
            def op(...):
                return ... if ... else ...
            ```
        * [ ] Every branch needs to return 
            ```python=
            @local_operator
            def op(...):
                if b:
                    return ...
                else:
                    return ...

            @local_operator
            def op(...):
                if b:
                    return ...
                return ...
            ```
        * [X] Symbols should always be defined in the code-path of their usage. In other words:
            - If a symbol is defined in both branches it may be used outside of the conditional
            - If a symbol is only defined in a single branch it may only be used inside that branch unless it was already defined previously (as this is equal to assigning the previous value to the again).
        ```python=
        @local_operator
        def op(...):
            a = 1
            if some_condition:
                a = 2
                b = 10
            else:
                tmp = 20
                b = tmp

            # equal to writing
            a = 1
            a, b = (2, 10) if some_condition else (a, 20)
        ```
        * [ ] Unrestricted `if`: Symbols may be used regardless of whether it is possible to statically prove they were assigned.
            ```python=
            @local_operator
            def op(...):
                if some_condition:
                    a = ...
                else:
                    b = ...

                if some_other_condition_equivalent_to_the_first_one:
                    return a
                else:
                    return b
            ```
        * [ ] ~~All symbols directly or indirectly involved in the return expression should always be defined regardless of the control-flow~~
            This option was withdrawn as it would allow symbols not being used in the return expression to have undefined values, which we would like to avoid to avoid errors further down the toolchain.
    - call builtins
    - return values
    - Restrictions for assignments and uses of variables:
       * Intuitively: Every time a variable is used, it has to be guaranteed (_statically_) that the variable is initialized at that point.
       * _More formally_: For each statement, there is a set of initialized variables before and after that statement. A simple forward flow analysis can be used to correctly deduce those sets for each statement. If-conditions (or more generally _path constraints_) are ignored (hence this is a static algorithm as it ignores runtime information like the conditions to if statements). Because our control flow graph is very simple (no loops, no goto, etc) this can be done in a single pass (no fixed-point iteration required).
       * Further reading:
           * Definite assignment analysis: https://en.wikipedia.org/wiki/Definite_assignment_analysis
           * Data-flow analysis: https://en.wikipedia.org/wiki/Data-flow_analysis
    - Weird behaviour with symbol binding & _Fields_:
        ```python=
        # Question: are the following two examples equivalent?
        
        @gt.local_operator
        def bindings1(input_a: Iterator, input_b: Iterator, condition: Iterator):
            x: Iterator = input_a
            evaluate_if[...]
            if condition:
                x: Iterator = input_b
            y: Iterator = x(I + 1)
            return x


        @gt.local_operator
        def bindings2(input_a: Iterator, input_b: Iterator, condition: Iterator):
            tmp_it: Iterator = x(I + 1)
            tmp_it2: Iterator = another_op[...](a,b)
            tmp_it: Value = x(I+1) + 2.
            #tmp_it3: Iterator = x if condition else y#
            

            x: Iterator = input_b(x+1) if condition else input_a(x+1)
            y: Iterator = x(I + 1)
            return y
        
        # Answer: No
        # Because in `bindings2` the condition `condition` will be lifted:
        y = input_b[I + 1] if condition[I + 1] else input_a[I + 1]
        # Whereas as in `bindings1` the condition `condition` won't be shifted:
        y = input_b[I + 1] if condition else input_a[I + 1]

        # Lessons:
        # * Conditions to if-statements are point-wise (have a local value)
        #   and have eager-evaluation.
        # * Symbol bindings are local.
        # * When lifting an expression it can change from eager-evaluation to
        #   lazy-evaluation.
        # * Implicit lifting can mean that eager-evaluation gets changed to
        #   lazy-evaluation.
        ```
    - Question regarding temporaries: Can they be accessed with off-sets?
        - Some simple examples in pseudo code for the structured cartesian use case:
            ```python=
            @local_operator
            def op(flat_field, column_field, full_field, scalar):
                a: value = flat_field * column_field + full_field * scalar
                #a: value = deref(flat_field) * deref(column_field) + full_field * scalar
                
                a:view = compute_blabla[...](flat_field, column_field, full_field, scalar)
                # a = lift(compute_something)(flat_field, column_field, full_field, scalar)
                
                # can we access temporaries with offset and if so what does it mean:
                # if ...:
                #     return a[i + 1]
                # elif ...:
                #     return a[k + 1]
                # elif ...:
                #     return a[i + 1, k + 1]
                # return a

##------------
            def op(flat: LocalView[I, J], full: LocalView[I, J, K]) -> value:
                return flat*full(K+1)
                
                # lift explicit todo EGP
                # all axis not part of the iteration domain but in a local view
                #  need to be explicitly defined (e.g. by shift)
                
            @local_operator(domain=Domain[Vertex])
            def some_other_unstructured_op(v: LocalView[Vertex], e: LocalView[Edge]):
                # domain: [Vertex]
                tmp = e(V2E[0]) * v
                tmp2 = e.field[UNDEFINED] + 1.
                return tmp + tmp2
                op_view = op[...](flat, full_at_k)
                return op_view(K+1) #+ flat(K+1)
                
            @local_operator(domain=Domain[Cell])
            def some_other_unstructured_op(c: LocalView[Cell], e: LocalView[Edge]):
                val_left = 2*e(C2E, 0)**2
                val_right = 2*e**2 for e_nb_vals in e(C2E)
                return 0.5 * (val_left + val_right)
                
            @local_operator(domain=Domain[I, J])
            def some_other_op(flat: LocalView[I, J], full: LocalView[I, J, K]):
                # domain: [I, J]
                field_full_at_k = field_full(K[3])
                op_view = op[...](flat, full_at_k)
                return op_view(K+1) #+ flat(K+1)

            @local_operator(domain=Domain[I, J])
            def some_other_op(normal: LocalView[I, J], staggered: LocalView[I_s, J]):
                # domain: [I, J]
                tmp1 = normal + 1
                tmp2 = staggered(I-1/2) + 1
                return normal + staggered(I-1/2)

            @program  
            def main(field_flat, field_full):
                some_other_op[IndexRange(I, 0, 10), IndexRange(J, 0, 10)](field_flat, field_full_at_k)
##------------
        

            @local_operator
            def op(a: LocalView[I, J], full_b: LocalView[I, J, K]):
                class LocalView:
                    field: Field
                    position: Position = UNDEFINED
                    
                    def __call__(self, offset_spec) -> LocalView:
                        return LocalView(self.field, shift(self.position, offset_spec))
                        
                    @property
                    def __value__(self):
                        return self.field[self.position]
                    
                a = a + full_b  # field_a.__value__ + full_field.__value__
                a_2 = a(I+1) * 2                
                
                return a_2
                
                
                another_view:view = compute_something(flat_field, column_field, full_field, scalar)
                # a = lift(compute_something)(flat_field, column_field, full_field, scalar)
                
                # can we access temporaries with offset and if so what does it mean:
                if ...:
                    return a[i + 1]
                elif ...:
                    return a[k + 1]
                elif ...:
                    return a[i + 1, k + 1]
                return a

```

__Examples__:
```python=
# TODO
def local_operator(...) -> Operator[]:
    ...

@gt.local_operator(...)
def operator_definition(a: LocalView, b: LocalView) -> Value:
    return a + b

@gt.local_operator(...)
def operator_definition(
    a: LocalView, b: LocalView, *, position: Position
) -> Value:
    return a + b if position[0] > 3 else 0
```

#### Column operator

An operator accepting and returning _ColumnViews_. Those are used for specialized processing of fields _columns_  (usually vertical columns). Currently the only way to directly create column operators is by using the `scan()` builtin, which applies a _scan pass_ function over a specific axis while carrying along state. Combined column operators can be also defined by grouping calls to column operators in a single function

__Syntax proposal (WIP)__:

`scan()` builtin:
- `scan(scan_pass_op, *, forward: bool = True, init=None)` builtin. Alternative: `scan(scan_pass_op, *, reverse: bool = False, init=None)` builtin
    - Signature of scan pass operator: `scan_pass_op(state: Any, *views: LocalView) -> Tuple[Any, Value]  # (state, output)` 
    - Accepted data arguments and statements inside the body as in [local operator](#Local-operator)

`column_operator(...)`:
- A Python function decorated with `@gt.column_operator(...)`
    - Keyword arguments: TODO
- Accepted data arguments are `LocalView` 


__Examples__:
```python=
@gt.local_operator
def tridiag_forward(
    state: Tuple[float, float], a: LocaLView[], b : ..., c, d
):
    if state is None:
        return (c / b, d / b)
    cp_km1, dp_km1 = state
    cp_k = c / (b - a * cp_km1)
    dp_k = (d - a * dp_km1) / (b - a * cp_km1)
    return (cp_k, dp_k)
    
@gt.local_operator
def tridiag_backward(x_kp1: float, cp: ..., dp):
    if x_kp1 is None:
        return dp
    return dp - cp * x_kp1
    
@gt.column_operator
def solve_tridiag(a, b, c, d):
    cp, dp = scan(tridiag_forward, True, None)(a, b, c, d)
    return scan(tridiag_backward, forward=False, init=None)(cp, dp)
        
@gt.program(sequential_axis=K)
def solve_tridiag(i_size, j_size, k_size, a, b, c, d, x):
    domain = IndexRange[IDim](0, i_size)*IndexRange[JDim](0, j_size)*IndexRange[KDim](0, k_size)
    solve_tridiag[domain](a, b, c, d, out=[x])
    
# --- Option 2: directly defined a scan pass using scan as a decorator ---
@gt.scan(forward=True, init=None)
@gt.local_operator
def tridiag_forward(
    state: Tuple[float, float], a: LocalView[], b : ..., c, d
):
    ...
   
@gt.column_operator
def solve_tridiag(a, b, c, d):
    cp, dp = tridiag_forward(a, b, c, d)
    return tridiag_backward(cp, dp)
```

### Program

A sequence of field operations (and entry point from the host language to GT4Py). Currently, the only allowed field operation is the application of local/column operators.

__Syntax proposal (WIP)__:
- Typing:
    ```python=
    P = ParamSpec("P")
    
    class Program(Protocol[P]):
        __call__: Callable[P, None]
        
        def build(**kwargs: Any) -> Program[P]:
            """Build a program using different compile/backend settings."""
            ...        
    ```
- A Python function decorated with `@gt.program(...)`.
    - If all the required build arguments are provided in the decorator, it builds a `Program` object directly, otherwise it would return a lazy program object that needs to be explicitly built using `.build(**kwargs)`.
    ```python
    # Lazy application
    @gt.program
    def prog(...):
        ...
    
    compiled_program = prog.build(backend=..., )    
    ```
    - Keyword arguments:
        + _connectivities_: list of `Connectivity` type definitions
        + _sequential_axis_: Axis  (for all sequential/column operators)
        + _backend_: str ?
        + _options_: Mapping[str, Any] ?
        + `gt.globals` vs `externals` run-time args??
    - Compile-time constants ?
        + [ ] externals dict-like passed as argument at compilation
            ```python=
            @gt.operator(...)
            def op(a, ...):
                from __externals__ import A
                return a*A

            @gt.program(externals={"A": 32.34})
            def program():
                ...
            ```
        + [ ] direct binding to python symbols (maybe only for _frozen_ objects?)
            ```python=
            # option 1
            A = 3.5
            # option 2 (safer)
            A = gt.global(3.5, name="A", dtype=float)
            
            @gt.operator(...)
            def op(a, ...):
                return a*A
            ```
        + [ ] frozen dataclass like
            ```python=
            class ModelConstants(gt.constant_container):
                b = 1
                
            @gt.operator(...)
            def op(a, ...):
                # Option 1:
                ModelConstants.b
            ```
        + [ ] partial
            ```python=
            # functools partial inspired
            # https://docs.python.org/3/library/functools.html#functools.partial
            @gt.operator(...)
            def op(A, a, ...):
                return a*A
            
            @gt.program
            def program(A, a, ...):
                op(A, a, ...)
                
            gt.partial(program, A=32.34)
            ```
        + [ ] both
    - Run-time values ?
        + [ ] global variables defined as `gt.global()`
        + [ ] direct binding to python symbols
        + [ ] run-time arguments
        ```python=
        @gt.operator(...)
        def op(dt: float, a, ...):
            return dt*a

        @gt.program
        def program(dt: float, ...):
            ...
            op[domain](dt, a)
        ```
        + [ ] both
- Accepted data arguments are `Field`s and scalars (run-time parameters)
    + Question: should be scalar tuples supported ?? How to enforce fixed length?
- Accepted statements inside the body:
    - Operator application (taken from Numba syntax):
    ```python=
    operator_name[domain_definition](*input_fields, out=[*output_fields])
    ```
    - Domain/var definitions: TODO

__Examples__:
```python=
def program(connectivities: List[Connectivity]) -> Program:
    ...
ASDF = gt.global(3)
@gt.program(connectivities=[Conn1, Conn2, A2B])
def my_program(domain_indication: Any, field_a: Field[], field_b: Field) -> None:
    # Operator application
    domain = ... # TODO
    operator[domain](field_a, field_b, out=[field_b])

# Option A
another_build = my_program.build(connectivities=OTHER_CONNS)

my_program((3,10), field_a, field_b, connectivities=[...])
```



## Alternative Positional Proposal

In this proposal positions and indices respectively are treated as first-class-citizens. The datastructures introduced here corresponding 1-to-1 to the mathematical concepts given at the end of the document.

### Field operator

*Type signature* `FieldOperator`
```python=
__init__(local_op: Union[LocalOperator, ColumnOperator], *, domain)
__call__(*args: Union[Field, Value], *, out: Union[Field, Tuple[Field, ...]]) -> Field
```

*Construction*
Field operators are constructed from local- or column-operators either by specification of a domain or by lifting them. (The `__init__` function was only given for completeness previously)

```python=
# only valid in @gt.program context
field_op: FieldOperator = op[domain]
# only valid in @gt.local_operator context
lazy_field_op: FieldOperator = op[...] # lift
```
*Call* The usual function call syntax is used to evaluate a field operator for a given set of fields and values.

At the `gt.program` level field operators need to be called with the additional `out` argument specifying field(s) into which the results should be materialized. As such a `gt.program` is stateful.
```python=
# context: gt.program
@gt.program(...)
def program(..., f_out):
    ...
    # 1. construct field operator from local operator `op`,
    # 2. evaluate in all points of the `domain` with given fields
    #  and materialize into `f_out`
    return op[domain](f_1, ..., f_n, out=f_out)
```



This syntax is inspired by the [numba kernel call syntax](https://numba.pydata.org/numba-doc/latest/cuda/kernels.html).

### Local operator

Local operators are defined by decoration of some (python) function `op` with `@gt.local_operator` decorator.

*Definition signature* (unstructured + cartesian)
```python=
op(*args: Union[Field, DT], *, position: Position) -> Value
```
- `args`: set of arguments being either a Field or a value.
- `position`: mandatory keyword argument

*Definition signature* (cartesian)  
For the cartesian case the local operator may be defined and called without a position argument.
```python=
op(*args: Union[Accessor, DT]) -> Value
op(*args: Union[Accessor, DT], *, position: Position) -> Value
```
- `args`: set of arguments being either an Accessor or a value.
- `position`: mandatory keyword argument

*Type signature* `LocalOperator`  
The decorator returns a `LocalOperator` object with the following signature. Note that all overloads of `__call__` may be used regardless of the signature of the decorated function.
```python=
__init__(op: Callable, domain: Optional[Type[DomainT]])

__call__(*args: Union[Field, DT], *, position: Position) -> Value

__call__(*args: Union[Accessor, DT], *, position: Position) -> Value
__call__(*args: Union[Accessor, DT]) -> Value
```

If a `LocalOperator` is called with accessors their focal point must agree (can be checked easily both at the IR level and during embedded execution). As such they - together with the type of the domain passed to the decorator - implicitly provide the position. This is an additional requirement not given by the semantic-model but introduced in order to guide the user to write comprehensible code that maintains the perspective of a *local operator*, i.e. a function defined on __one__ position.


#### Usage

__Definition__

```python=
@gt.local_operator
def sum_vnb(inp: Field[[V], DT], *, position: Position[V]):
    (v,) = position
    return sum(inp[v2] for v2 in neighbors(v, V2V))

@gt.local_operator(domain=Domain[I, J])
def sum_vnb(inp_sx: Accessor[[I+1/2, J], DT], inp_sy: Accessor[[I, J+1/2], DT]):
    # todo: define how the accessor relates to the domain
    # accessors can not be created they are (at least conceptually)
    # constructed by the decorator and their focal point is
    # the current position potentially shifted based on the index
    # convention (here from I, J to I+1/2, J or I, J+1/2)
    return inp_sx[I-1/2]+inp[I+1/2]+inp_sy[J-1/2]+inp_sy[J+1/2]
```

__Evaluation__ Local operators may either be called inside of another local operator (by direct call or after lifting)

```python=
@gt.local_operator
def sum_vnb_wrapper(inp: Field[[V], DT], *, position: Position[V]):
    (v,) = position
    
    # direct call
    return sum_vnb(inp, position=neighbor(v, (V2V, 0)))
    
    # lifting
    #  instructive version
    #   1. lift local operator to field operator
    sum_vnb_fop: FieldOperator = sum_vnb[...]
    #   2. field operator evaluation into a field
    sumvnb_field: Field = sum_vnb_fop(inp)
    #   3. access at a neighbor position
    return sumvnb_field[neighbor(v, (V2V, 0))]
    #  compact version
    return sum_vnb[...](inp)[neighbor(v, (V2V, 0))]
```
or at the program level by transforming them to a field operator.
```python=
@gt.program
def fencil(...):
    # 1. transformation into a field operator
    #  `sum_vnb[vertices]`
    # 2. materializing into `out` field
    sum_vnb[vertices](inp, out=out)
```

__Composition__

```python=
@gt.local_operator
def loc_op(inp: Field[[V], DT], *, position: Position[V]):
    # direct call + lift
    return loc_op1(local_op2[...](inp), position=position)
    # lift + lift + access
    return loc_op1[...](local_op2[...](inp))[position]
```

WIP:
```python=
@gt.local_operator
def compute_nb_vals_avg(nb_vals: Field[V, Field[[ConnT], DT]], *, position: Position[V]):
    (v,) = position
    
    return sum(nb_vals[v])/len(nb_vals[v])

@gt.local_operator
def field_getitem_with_nbs(self, nbs):
    return self[v] for v in nbs

@gt.local_operator
def loc_op(val: Field[[V], DT], *, position):    
    nbs1 = neighbors(v, V2V)
    nbs2 = neighbors(v, V2V, V2V)
    
    # creation of a sparse field
    nb_vals: Field[[V], Field[V2V, DT]] = field_getitem_with_nbs[...](val, nbsX)
    nb_vals = val[nbsX] # doesn't match numpy syntax
    nb_vals = val.remap(v, nbs)
    nb_vals = (lambda v: val[v] for v in nbsX)[...]() # inline local op
    
    
    return compute_nb_vals_avg(nb_vals, position=position)
```

<!--in which case the position has to be passed manually 

together with a domain at the program level where they are automatically lifted to a field operator such that the local operator is called for each position in the domain.-->

### Column operator

A Python function decorated with `@gt.column_operator(init: DT=init)`, where `init` is an initial value passed to `op`.

*Signature*
```python=
op(state: DT, *args: Union[Field, DT], *, position: Position) -> Field[[Dim], DT]
```
where `DIM` as the specified column dimension.
- `state`: value from the previous call to `op` else `init`
- `args`: set of arguments being either a Field or a value.
- `position`: mandatory keyword argument

```python=
@gt.column_operator
def column_sum(state: Union[DT, Tuple[DT]], column: Field[[K], DT]):
    *_, k = position
    
    return column[k]
```

### Accessing neighbors

Neighbors can be accessed using the built-in functions `neighbor` for accessing a single neighbor or `neighbors` for accessing a set of neighbors.

*Signature*
```python=
neighbor(idx: Index, conn: Type[ConnT], int) -> Index

neighbors(idx: Index, *shifts: Tuple[Type[ConnT], int]) -> List[Index]
```

with `ConnT = TypeVar("ConnT", bound=Connectivity)`

*Neighbor aliasing*  
Since indices are first-class-citizens there is no need to label or alias them.
```python=
v_left, v_right = vertices(e)
```

#### Field

`__init__(domain: Union[Domain[*Dim], None], image: np.array)`

`__getitem__(self, pos: Position)`

## Mathematical concepts

### Positional

This section contains the definition of all concepts in pure mathematical terms. 

[WIP] Not thoroughly connected to the rest of the document.

__Value__  
A *value* is either a number, a tuple of numbers, or a field.

__Dimension__  
A *dimension* is a just a tag or type for [indices](#index) ([wikipedia](https://en.wikipedia.org/wiki/Dimension)).

__Index__  
An index is an integer in $\mathbb{Z}$ tagged with a [dimension](#Dimension), i.e. an element of an axis.

__Neighbor Dimension__

Similar to a dimension a *neighbor dimension* is just a tag or type, but for [local indices](#local_index)

__Local index__  
A local index is an integer number tagged with a [neighbor dimension].

__Index range__  
An *index range* is the set of all integers between two integers $a, b \in \mathbb{Z}$, e.g $\{1, 2, 3\}$, $\{-1, 0, 1\}$ ([wikipedia](https://en.wikipedia.org/wiki/Interval_(mathematics)#Integer_intervals))

todo: index type

__Domain__  
A *domain* is a cartesian product of index ranges.

todo: index dims only once  
todo: enforce canonical representation  
todo: elements of domain are positions

<!--__Axis__  
The set of all possible [indices](#index) in a particular dimension.-->
  
__Position__  
A tuple of indices indicating a position on the grid / mesh.

__Field__  
A *field* is a function $f: D \to V, pos \mapsto v$, with $D$ a domain and $V$ a set of values.

__Connectivity__  
A *connectivity* is a function $f: pos \mapsto (nb\_pos_1, \ldots, nb\_pos_n)$, with $pos$ a position and $nb\_pos_i$ a set of positions neighboring $pos$.

__Field operator__  
A *field operator* is (higher-order) function $op: (f_1, \ldots, f_n) \mapsto f_{\text{out}}$, where $f_1, \ldots, f_n, f_{out}$ are a set fields. In other words a field operator is a function mapping fields to a field.

TODO: Apply local operator to all points

Note the difference in notation to the local operator originating from the higher-order function nature of the field operator. Instead of specifying the spaces of the functions arguments and value we are denoting the arguments and the value directly.

__Local operator__  
A *local operator* is a (higher-order) function $\hat{op}: (pos, f_1, \ldots, f_n) \mapsto v$, where $pos$ is a position, element of the domain $D$, $f_1, \ldots, f_n$ a set of fields, and $v \in V$ a value.

TODO: spaces

__Column operator__  
A *column operator* is a (higher-order) function $\hat{op_c}: (pos, f_1, \ldots, f_n) \mapsto f_c$, where $pos$ is a (todo: horizontal) position, $f_1, \ldots, f_n$ a set of fields, and $f_c$ a column, i.e. a field whose domain is one dimensional.

__Local view / iterator / accessor__   
A *local view* is a tuple $(pos, f)$ where $pos$ is a position and $f$ a field. It can be used as the actual value for any value operation and it additionally supports an extra operation to create another local-view placed on a different position, by mapping or shifting its current one.

<!--__Shift__  
$$\text{shift}((pos, f), offset) = ()$$

$$\text{shift}(f, dim, o)(pos) = f(pos_0, \ldots, pos_{dim}+o, \ldots, pos_n) \tag{Cartesian}$$ <-- imprecise --
$$\text{shift}(f, nb\_dim, l\_idx)(pos) = f((conn(nb\_dim)(pos))_{l\_idx}) \tag{Unstructured}$$-->


__Lift__  
The (higher-order) function $lift: \hat{op} \to op$ transforms a local operator or column operator $\hat{op}$ to a field operator $op$ by partial application of all arguments after the $pos$ argument of $\hat{op}$.

$$\underbrace{\text{lift}(\hat{op})}_{= op}(f_1, \ldots, f_n) = pos \mapsto\hat{op}(pos, f_1, \ldots, f_n)$$

Note that for $\hat{op}$ a column operator the resulting field operator $op$ returns a field of columns. To avoid this we can (optionally) define a different lift for column operators $\hat{op_c}$ as:
$$\underbrace{\text{lift}(\hat{op_c})}_{= op}(f_1, \ldots, f_n) = (pos, k) \mapsto \hat{op_c}(pos, f_1, \ldots, f_n)(k)$$

Let $\text{shift}$ be a function such that the following holds $\text{shift}(f, o)(pos) = f(pos+o)$ then:

$$\text{shift}(\text{lift}(\hat{op})(f_1, \ldots, f_n), o)(pos) = \hat{op}(pos, \text{shift}(f_1, o), \ldots, \text{shift}(f_n, o))$$

In other words shifting a lifted local operator is equal to shifting the arguments.

__Scan pass__  
A *scan pass* is a function $\hat{op_s}: (pos, state, f_1, \ldots, f_n) \mapsto v$, where $pos$ is a position, *state* is a value, $f_1, \ldots, f_n$ a set of fields, and $v \in V$ a value.

__Scan__  
The function $scan_C: (\hat{op_s}, dir, init) \mapsto \hat{op_c}$ transforms a scan pass $\hat{op_s}$, a direction $dir \in \{-1, 1\}$ (either forward $1$ or backward $-1$), and an initial value $init$ into a column operator $\hat{op_c}$, where $\hat{op_c}$ is defined as

$$\begin{aligned}&\hat{op_c}(pos, f_1, \ldots, f_n) \\ &= k \mapsto  \begin{cases}
    \hat{op_s}((pos, k_{0}), init, f_1, \ldots, f_n) & | \text{ if } k = k_{0} \\
    \hat{op_s}\Big((pos, k), \hat{op_c}\big((pos, k-dir), f_1, \ldots, f_n\big), f_1, \ldots, f_n\Big)
\end{cases}\end{aligned}$$

with $k_0$ being either the first, resp. last, position of the column $C$ for $dir = 1$, resp. $dir=-1$.

Alternatively we can make $scan$ transform into a field operator $op$ instead of a column operator:

$$\begin{aligned}&op(f_1, \ldots, f_1) = \\ &= (pos, k) \mapsto \begin{cases}
    \hat{op_s}((pos, k_{0}), init, f_1, \ldots, f_n) & | \text{ if } k = k_{0} \\
    \hat{op_s}\Big((pos, k), \hat{op_c}\big((pos, k-dir), f_1, \ldots, f_n\big), f_1, \ldots, f_n\Big)
\end{cases}\end{aligned}$$

__Map__

The (higher-order) function $map: (\hat{op}, D) \mapsto f_{out}$ takes a local operator $\hat{op}$ and a domain $D$ and returns a field $f_{out}: D \to V, pos \mapsto \hat{op}(pos, f_1, \ldots, f_n)$.

### Iterator view

__Iterator__
An *iterator* is a function $o \mapsto f(\text{shift}_{pos}(pos, o))$ from an offset $o$ to the value of some field $f$ at some position $pos$ shifted by $o$.

$\text{shift}_{it}(o) = it \mapsto (o \mapsto it(o))$

It can be used as the actual value for any value operation $f(pos)$ (`deref`) and it additionally supports an extra operation to create another iterator placed on a different position $\text{shift}_{it}$ (`shift`).

__Lift__
The (higher-order) function $lift: \hat{op} \to op$ transforms a stencil $\hat{op}$ returning a value to a stencil returning an iterator by shifting $\hat{op}$'s arguments.

$$
\text{lift}(\hat{op})(it_1, \ldots, it_n) = (it_1, \ldots, it_n) \mapsto (o \mapsto\hat{op}(\text{shift}_{it}(it_1)(o), \ldots, \text{shift}_{it}(it_n)(o)))
$$

As such $\text{lift}$ allows to evaluate a stencil returning values at an offset.

<!--__Relation to the iterator view__

By replacing each field $f_i: pos \mapsto f_i(pos)$ within the above representation (except for the map operation) by an iterator $it: o \mapsto \text{shift}(f_i, o)(pos_i)$ we arrive at the iterator view.-->

Notes:
- The lifted iterator of the iterator view needs to store different data than the operator passed from the apply_stencil call. As such iterators can not be encoded using a single (mathematical representation) easily.
- remove integer requirement for indices
- domain is not an isolated concept, but tied to a field
- the current axis concept does not match previous discussions with EGP (universe concept). Is it even useful as what we called universe must be tied to a field right?

Scratchpad:


```python=
var_offset=1 if some_cond else 2

var_offset=shift(K, 1) if some_cond else shift(K, 2)

var_offset(it)
```