text:
# SHORT COMMUNICATION: FLOW AS DISTRIBUTED LINES WITHIN THE LANDSCAPE

### John J. Armitage

Dynamique des Fluides Géologiques, Institute de Physique du Globe de Paris, Paris, France

## Abstract

Landscape evolution models (LEMs) aim to capture an aggregation of the processes of erosion and deposition within the Earth's surface and predict the evolving topography. Over long time-scales, i.e. greater than one million years, the computational cost is such that numerical resolution is coarse and all small-scale properties of the transport of material cannot be captured. A key aspect therefore of such a long time-scale LEM is the algorithm chosen to route water down the surface. I explore the consequences of two end-member assumptions of how water flows over the surface of a LEM, either down a single flow direction (SFD) or down multiple flow directions (MFD), on model sediment flux and valley spacing. I find that by distributing flow along the edges of the mesh cells, node-to-node, the resolution dependence of the evolution of a LEM is significantly reduced. Furthermore, the flow paths of water predicted by this node-to-node MFD algorithm is significantly closer to that observed in nature. This reflects the observation that river channels are not necessarily fixed in space, and a distributive flow captures the sub-grid scale processes that create non-steady flow paths. Likewise, drainage divides are not fixed in time. By comparing results between the distributive transport-limited LEM and the stream power model "Divide And Capture", which was developed to capture the sub-grid migration of drainage divides, I find that in both cases the approximation for sub-grid scale processes leads to resolution independent valley spacing. I would therefore suggest that LEMs need to capture processes at a sub-grid scale to accurately model the earth's surface over long time-scales.

## Introduction

It is known that resolution impacts landscape evolution models (LEMs) ([Schoorl et al., 2000](#)). The resolution dependence of LEMs is caused by how run-off is routed down the model surface. It has been demonstrated that either distributing flow down all slopes (multiple flow direction, MFD), or simply allowing flow to descent down the steepest slope (single flow direction, SFD), gives different outcomes for landscape evolution models ([Schoorl et al., 2000](#), [Pelletier 2004](https://doi.org/10.1029/2004GL020802)). It has been noted that landscape potentially has a characteristic wavelength for the spacing of valleys ([Perron et al., 2008](https://doi.org/10.1029/2007JF000977)). Therefore, a landscape evolution model should be able to reproduce such regular topographic features independently of the model resolution. For a model of channelised flow it was however found that the routing of run-off led to a resolution dependence in the valley spacing, which could be overcome by the addition of a parameterised flow width that was less than the numerical grid spacing ([Perron et al., 2008](https://doi.org/10.1029/2007JF000977)).

There is a potential problem with parameterising the flow width to be fixed at a sub grid level. The response time of LEMs to a change in external forcing is strongly dependent of the surface run-off ([Armitage et al., 2018)](https://doi.org/10.5194/esurf-6-77-2018)). This means that the model response time becomes likewise dependent on the chosen flow width. Ideally the LEM would be independent of grid resolution without introducing a predefined length scale that impacts the model response.

Water is the primary agent of landscape erosion. There are multiple pathways within the hydrological cycle from evaporation, transpiration, and ground water flow, however for many landscapes the river network is the primary route through which water flows down slope. Mean river width varies from 5 km to a few meters ([Allen & Palvesky, 2018](https://doi.org/10.1126/science.aat0636)). The very wide rivers, greater than 1 km, are however outliers within this global data set, with the median of the distribution of mean river width being 124 m and the upper quartile at 432 m (Figure 1). In LEMs developed for understanding long-term landscape evolution the large time scales necessitate large spatial scales, where a single grid cell can be a kilometer wide or more ([Temme et al., 2017](https://doi.org/10.1002/esp.4162)). A spatial resolution of cells larger than a few meters becomes necessary when modelling at the scale of a continent (e.g. [Salles et al., 2017](https://doi.org/10.1002/2016GC006617)). This means that flow has a width at a subgrid level.

```python
%matplotlib inline
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

data = gpd.read_file('data/GRWL_summaryStats.shp')
fig = plt.figure()
ax = data['width_mean'].plot(kind='hist',bins=20)
ax.set_xlabel('Mean river width (m)')
plt.show()
```

**Figure 1**: *Distribution of mean river width taken from the Global River Widths from Landsat (GRWL) Database (Allen & Palvesky, 2018)*

If the width of the flow path for run-off is narrower than can be reasonably modelled, then can the flow paths be treated as lines, from model node-to-node (Figure 2), where water collects along these lines? To explore this idea and understand LEM sensitivity to resolution, I wish to explore how a simple LEM evolves under four scenarios (Figure 2): (1) simple single flow direction (SFD) from cell area to cell area, (2) a multiple flow direction (MFD) version of this cell-to-cell algorithm, (3) a node-to-node SFD, and (4) a node-to-node MFD.

## A landscape evolution model

In this study I will assume landscape evolution can be effectively simulated with the classic set of diffusive equations described in [Smith & Bretherton, 1972](https://doi.org/10.1029/WR008i006p01506):
\begin{equation}
\frac{\partial z}{\partial t} = \nabla \left[ \left(\kappa + c q_{w}^{n}\right) \nabla z \right] + U
\end{equation}
where $\kappa$ is a linear diffusion coefficient, $c$ is the fluvial diffusion coefficient, $q_{w}$ is the water flux, $n$ is the water flux exponent, and $U$ is uplift. This heuristic concentrative-diffusive equation is capable of generating realistic landscape morphology, with the slope-area relationships commonly observed ([Simpson & Schlunneger, 2003](https://doi.org/10.1029/2002JB002162), [Armitage et al., 2018](https://doi.org/10.5194/esurf-6-77-2018)). Strictly it assumes that there is always a layer of material to be transported by surface run-off, and as such it can be classed as a transport-limited model. It accounts for both erosion and deposition, and is therefore appropriate for modelling landscape evolution beyond mountain ranges and into the depositional setting (see models such as DIONISOS; [Granjeon & Joseph, 1999](#)). It differs from mixed erosion and deposition models such as [Kooi & Beaumont, 1994](#) and [Davy & Lague, 2009](https://doi.org/10.1029/2008JF001146) because those models split the divergence of the sediment flux into two terms, a rate of erosion and rate of deposition. Here instead I assume that the sediment flux is a function of water flux and slope.

Equation 1 is solved with a finite element scheme written using Python and the FEniCS libraries (I will call the code [fLEM](https://github.com/johnjarmitage/flem)). The equations are solved on a Delaunay mesh, where the mesh is made up of predominantly equilateral triangles with an opening angle of 60$^{\circ}$. Model boundary conditions are initially of fixed elevation on the sides normal to the x-axis and zero gradient on the sides in normal to the y-axis. The model aspect ratio is 4 to 1. Uplift is fixed at $U=10^{-4}\rm\,m\,yr^{-1}$, the linear diffusion coefficient is $\kappa = 1\rm\,m^{2}\,yr^{-1}$, the fluvial diffusion coefficient is $c=10^{-4}\rm\,(m^{2}\,yr^{-1})^{n-1}$, and the water flux exponent is $n=1.5$.

```python
# markdown is not rendering images in voila...
import matplotlib.image as mpimg
img=mpimg.imread('figures/MFDandSFD.png')
plt.figure(figsize=(16, 6))
plt.imshow(img)
plt.axis('off')
plt.show()
```

**Figure 2**: *Diagram of flow routing from cell-to-cell and node-to-node for either a single flow direction (SFD) and a multiple flow direction (MFD) algorithm weighted by the relative gradient.*

Water can be routed from cell-to-cell, where precipitation is collected over the area of each cell, sent downwards, and accumulates. In this cell-to-cell configuration the water flux has units of length squared per unit time and is given by:
\begin{equation}
q_{w}\mathrm{[cell]} = \frac{\alpha a}{l_{s}},
\end{equation}
where $\alpha$ is precipitation rate, $a$ is the cell area, and $l{s}$ is the length from cell center to cell center down the steepest slope (Figure 2a and b). This gives a water discharge per unit length, which has the advantage of not having to explicitly state the sub-grid width of the flow ([Simpson & Schlunneger, 2003](https://doi.org/10.1029/2002JB002162)). However, implicitly this implies that the flow is over the width of a cell. An alternative is to route water from node to node along cell edges and for it to accumulate. I assume that along the length of each cell edge water can be added to the flow line, assuming that the input is linearly related to the length of the flow line,
\begin{equation}
q_{w}\mathrm{[node]} = \alpha l,
\end{equation}
where $l$ is the length of the edge that joins the up-slope node to the down-slope node (Figure 2c and d). This means that the cell area is ignored and instead water enters the flow path uniformly along its length and accumulates down slope.

Equation 3 makes the assumption that water accumulates as a function of length. Water flux is observed to related to catchment area, $Q_{w} \propto A^{0.8}$ ([Syvitski & Milliman, 2007](https://doi.org/10.1086/509246)). The catchment length, $l$ is then related to area by, $l\propto A^{1/p}$, where $1.4<p<2.0$ ([Armtiage et al., 2018)](https://doi.org/10.5194/esurf-6-77-2018)). At the lower end of the range this gives $Q_{w} \propto l^{1.12}$, suggesting that accumulating water as a linear function of flow length is a reasonable simplification. A knock on effect of this assumption is that the magnitude of the water flux predicted for the node-to-node routing is less than the cell-to-cell, as in the latter water is accumulated over cell areas, which is naturally larger than the cells edges.

Both equations 2 and 3 do not attempt to capture the interaction between water flux and river width, rather these are two methods to approximate run-off within a coarse numerical grid. For both the cell-to-cell and node-to-node methods the flow can then be routed down a single flow direction (SFD) or routed down multiple flow directions (MFD) weighted by the relative gradient, as in for example [Schoorl et al., 2000](#). I run the numerical model with a uniform precipitation rate of $\alpha = 1\rm\,m\,yr^{-1}$.

Equation 1 is made dimensionless following ([Simpson & Schlunneger, 2003](https://doi.org/10.1029/2002JB002162)) using the linear diffusion time scale and the model length in the $x$-direction, $L$. This means that equation 1 can be re-written as,
\begin{equation}
\frac{\partial \tilde{z}}{\partial \tilde{t}} = \nabla \left[ \left(1 + D \tilde{q}_{w}^{n}\right) \tilde{\nabla} \tilde{z} \right] + U
\end{equation}
and,
\begin{equation}
\tilde{\nabla} \cdot \left(\frac{\tilde{\nabla}\tilde{z}}{|\tilde{\nabla}\tilde{z}|}\right) \tilde{q}_{w} = -1
\end{equation}
where $x = \tilde{x}L$, $y = \tilde{y}L$, $z = \tilde{z}L$, $t = \tilde{t}L^{2}/\kappa$, $q = \tilde{q}\alpha L$, and,
\begin{equation}
D =  \frac{c\alpha^{n}L^{n}}{\kappa}.
\end{equation}

## The effect of model resolution

```python
import os,sys
import pyvista as pv
from pyvista import set_plot_theme
set_plot_theme('document')
from flem import initialise, solve_flem
```

```python
# Domain dimensions
dem = 0
bounds = [0, 2e5, 0, 8e5]
res = 128
model_space, u_n, mesh, V, bc = initialise(dem, bounds, res)
```

```python
# Physical parameters
physical_space = [1e+0,1e-4,1.5,1,1e-4] # [kappa, c, nexp, alpha, U]

# 0 = MFD node-to-node; 1 = MFD cell-to-cell; 2 = SD node-to-node; 3 = SD cell-to-cell
flow = 3

# Time
dt = 1e4
num_steps = 501

# Output
out_time = 100

# Plot stuff
plot = 0

# Calculate valley to valley wavelength
statistics = 1
```

```python
name = 'SFDcelltocell'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
```

```python
#[sed_SDctc_128, time_SDctc_128, wav_SDctc_128] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('SFDcelltocell/u_solution_128_500000000.vtu')
flux = pv.read('SFDcelltocell/q_solution_128_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpSFDcelltocell = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
# 0 = MFD node-to-node; 1 = MFD cell-to-cell; 2 = SD node-to-node; 3 = SD cell-to-cell
flow = 2
name = 'SFDnodetonode'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_SDntn_128, time_SDntn_128, wav_SDntn_128] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('SFDnodetonode/u_solution_128_500000000.vtu')
flux = pv.read('SFDnodetonode/q_solution_128_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpSFDnodetonode = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
# 0 = MFD node-to-node; 1 = MFD cell-to-cell; 2 = SD node-to-node; 3 = SD cell-to-cell
flow = 1
name = 'MFDcelltocell'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_MFctc_128, time_MFctc_128, wav_MFctc_128] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('MFDcelltocell/u_solution_128_500000000.vtu')
flux = pv.read('MFDcelltocell/q_solution_128_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpMFDcelltocell = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
# 0 = MFD node-to-node; 1 = MFD cell-to-cell; 2 = SD node-to-node; 3 = SD cell-to-cell
flow = 0
name = 'MFDnodetonode'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_MFntn_128, time_MFntn_128, wav_MFntn_128] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('MFDnodetonode/u_solution_128_500000000.vtu')
flux = pv.read('MFDnodetonode/q_solution_128_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpMFDnodetonode = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
p = pv.Plotter(notebook=True, shape=(2, 2), border=False)
p.subplot(0, 0)
p.add_text("(a) SFD cell-to-cell, resolution 128 by 512", font_size=12)
p.add_mesh(warpSFDcelltocell, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.subplot(0, 1)
p.add_text("(b) SFD node-to-node", font_size=12)
p.add_mesh(warpSFDnodetonode, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.subplot(1, 0)
p.add_text("(c) MFD cell-to-cell", font_size=12)
p.add_mesh(warpMFDcelltocell, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.subplot(1, 1)
p.add_text("(d) MFD node-to-node", font_size=12)
p.add_mesh(warpMFDnodetonode, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.show(use_panel=False)
```

**Figure 3**: *Shaded relief with colours showing the magnitude of surface run-off from the single (SFD) and multiple flow direction (MFD) routing using the cell-to-cell and node-to-node aproach. Model resolution is 128 by 512 nodes.*

```python
from itkwidgets import view
view(geometries=warpMFDnodetonode)
```

    Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…

**Figure 4**: *Interactive plot of the multiple flow direction (MFD) node-to-node model for a resolution of 128 by 512 nodes.*

```python
res=512
model_space, u_n, mesh, V, bc = initialise(dem, bounds, res)
flow = 3
name = 'SFDcelltocell'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_SDctc_512, time_SDctc_512, wav_SDctc_512] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('SFDcelltocell/u_solution_512_500000000.vtu')
flux = pv.read('SFDcelltocell/q_solution_512_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpSFDcelltocell512 = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
res=512
flow = 2
name = 'SFDnodetonode'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_SDntn_512, time_SDntn_512, wav_SDntn_512] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('SFDnodetonode/u_solution_512_500000000.vtu')
flux = pv.read('SFDnodetonode/q_solution_512_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpSFDnodetonode512 = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
res=512
flow = 1
name = 'MFDcelltocell'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_MFDctc_512, time_MFDctc_512, wav_MFDctc_512] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('MFDcelltocell/u_solution_512_500000000.vtu')
flux = pv.read('MFDcelltocell/q_solution_512_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpMFDcelltocell512 = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
res=512
flow = 0
name = 'MFDnodetonode'
dirtest = './%s/' % name
directory = os.path.dirname(dirtest)
if not os.path.exists(directory) :
  os.makedirs(directory)
#[sed_MFntn_512, time_MFntn_512, wav_MFntn_512] = solve_flem(model_space,physical_space,flow,u_n, mesh, V, bc,dt,num_steps,out_time,plot,statistics,name)
```

```python
elevation = pv.read('MFDnodetonode/u_solution_512_500000000.vtu')
flux = pv.read('MFDnodetonode/q_solution_512_500000000.vtu')
flux['water flux (log10)'] = np.log10(flux['flx'])
flux['elevation'] = elevation['elv']
warpMFDnodetonode512 = flux.warp_by_scalar(scalars='elevation',factor=5e2)
```

```python
p = pv.Plotter(notebook=True, shape=(2, 2), border=False)
p.subplot(0, 0)
p.add_text("(a) SFD cell-to-cell, resolution 512 by 2048", font_size=12)
p.add_mesh(warpSFDcelltocell, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.subplot(0, 1)
p.add_text("(b) SFD node-to-node", font_size=12)
p.add_mesh(warpSFDnodetonode512, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.subplot(1, 0)
p.add_text("(c) MFD cell-to-cell", font_size=12)
p.add_mesh(warpMFDcelltocell, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.subplot(1, 1)
p.add_text("(d) MFD node-to-node", font_size=12)
p.add_mesh(warpMFDnodetonode512, scalars='water flux (log10)', cmap='viridis', lighting=True, show_scalar_bar=True)
p.show(use_panel=False)
```

**Figure 5**: *Shaded relief with colours showing the magnitude of surface run-off from the single (SFD) and multiple flow direction (MFD) using the cell-to-cell and node-to-node aproach. Model resolution is 512 by 2048 nodes.*

```python
view(geometries=warpMFDnodetonode512)
```

**Figure 6**: *Interactive shaded relief for the MFD node-to-node model wth a resolution of 512 by 2048 nodes.*

At a low model resolution, 128 by 512 cells, all four methods of flow routing give similar landscape morphology after 5 Myr of model evolution (Figures 3). However, elevations are significantly lower for the cell-to-cell flow routing model as the water flux term is lower for the node-to-node routing algorithm (e.g. compare Figure 3a and 3b). As the resolution is increased to 512 by 2048 cells, the landscape morphology starts to diverge. For the cell-to-cell SFD algorithm the landscape shows more small scale branching, as previously discussed by ([Braun & Sambridge, 1997](#)) (Figure 5a). For the SFD algorithm it can be seen that the high resolution model has multiple peaks along the ridges (Figure 5a). This roughness to the topography is removed if the flow is distributed down slope from cell to cell (MFD; Figure 5c).

For the node-to-node SFD algorithm, the increase in resolution has lead to significant branching of the valleys, which is clearly visible when the water flux is plotted (Figure 5b). For the node-to-node MFD algorithm, the morphology and distribution of water flux are similar for both the low and high resolution (Figure 5c and d). The two MFD models give a smoother topography, as by distributing flow local carving of the landscape is reduced.

```python
def pltsed(file,var,xlimit):
    sed = np.genfromtxt(file)
    NUMB = sed[:, 0]
    RESO = sed[:, 1]
    TIME = sed[:, 2]
    QS = sed[:, 3]

    j = 0
    for res in resolutions:
      times = TIME[np.where(RESO == res)]
      seds = QS[np.where(RESO == res)]
      nums = NUMB[np.where(RESO == res)]
      for num in number :
        plt.plot(times[np.where(nums == num)],seds[np.where(nums == num)],colors[j])
      j += 1

    Tr = 3*np.percentile(var['elevation'], 90)*ly/U*1e-6
    plt.plot((Tr, Tr), (0, 25), 'k--')
    plt.xlim((0, xlimit))
    plt.ylim((0, 25))
```

```python
# don't run this code block unless you are on a big computer.
# Here we can call the python script run_models.py to run multiple models at different resoltuions.
# As is it is set up for STD cell-to-cell. The code is just as in the cells above.

#!python run_models.py
```

```python
ly = bounds[3]
U = physical_space[4]
kappa = physical_space[0]

dtreal = 1e4*kappa/(ly*ly)

resolutions = [64, 128, 256, 512]
number = np.linspace(0, 9, 10)
colors = ['r', 'g', 'b', 'k']

plt.figure(figsize=(9, 6))
plt.subplot(221)
file = 'data/seds_SFDctc.txt'
var = warpSFDcelltocell
pltsed(file,var,2.5)

plt.subplot(222)
file = 'data/seds_SFDntn.txt'
var = warpSFDnodetonode
pltsed(file,var,7.5)

plt.subplot(223)
file = 'data/seds_MFDctc.txt'
var = warpMFDcelltocell
pltsed(file,var,2.5)
plt.ylabel('Sediment flux (m^{2}/yr)')
plt.xlabel('Time (My)')

plt.subplot(224)
file = 'data/seds_MFDntn.txt'
var = warpMFDnodetonode
pltsed(file,var,7.5)
plt.xlabel('Time (My)')

plt.show()
```

**Figure 7**: *Sediment flux that exits the model domain for (a) SFD cell-to-cell, (b) SFD node-to-node, (c) MFD cell-to-cell, (d) MFD node-to-node.*

To understand better how increasing resolution impacts the model evolution the total sediment flux eroded from the model domain is plotted against time, and the final valley spacing is calculated (Figure 7 and 8). To calculate the valley spacing I take horizontal swaths of the spatial distribution of water flux. For each swath profile a peak finding algorithm [Negri & Vestri](https://doi.org/10.5281/zenodo.887917) is used to find the distance from peak to peak in water flux. This distance is then averaged over the hundred swath profiles and over ten model runs to give the minimum, lower quartile, median, upper quartile, and maximum valley wavelength (Figure 8).

For the cell-to-cell SFD it can be seen that the evolution of the model is resolution dependent, as the wind-up time reduces as resolution is increased from 64 to 512 cells along the y-axis (Figure 7a). Furthermore, the mean valley spacing reduces with increasing resolution (Figure 8a). This behavior is not ideal, as it means that model behavior to perturbations in forcing might become resolution dependent. For the MFD cell-to-cell wind-up times remain resolution dependent, while the mean valley spacing is similar for the four different resolutions (Figure 7c and 8c).

```python
def pltwav(file,ylimit):
    wave  = np.genfromtxt(file)
    ress  = wave[:,0]
    width = wave[:,2]

    wav64  = np.zeros(200)
    wav128 = np.zeros(200)
    wav256 = np.zeros(200)
    wav512 = np.zeros(200)

    i = 0
    j = 0
    k = 0
    l = 0
    m = 0
    for res in ress:
      if res == 64:
        wav64[i] = width[m]
        i += 1
      if res == 128:
        wav128[j] = width[m]
        j += 1
      if res == 256:
        wav256[k] = width[m]
        k += 1
      if res == 512:
        wav512[l] = width[m]
        l += 1
      m += 1

    x  = [64,128,256,512]
    y  = np.zeros(len(x))
    y_ = np.zeros(len(x))

    y[0] = np.mean(wav64)
    y[1] = np.mean(wav128)
    y[2] = np.mean(wav256)
    y[3] = np.mean(wav512)

    y_[0] = np.std(wav64)
    y_[1] = np.std(wav128)
    y_[2] = np.std(wav256)
    y_[3] = np.std(wav512)

    ys = [wav64,wav128,wav256,wav512]
    plt.boxplot(ys)
    plt.ylim((0,ylimit))
    plt.xticks([1, 2, 3, 4], ['64x256', '128x512', '256x1028', '512x2056'])
```

```python
plt.figure(figsize=(9, 6))
plt.subplot(221)
file = 'data/wavs_SFDctc.txt'
pltwav(file,.5)

plt.subplot(222)
file = 'data/wavs_SFDntn.txt'
pltwav(file,.5)

plt.subplot(223)
file = 'data/wavs_MFDctc.txt'
pltwav(file,.5)
plt.xlabel('Resolution')
plt.ylabel('Normalised wavelength')

plt.subplot(224)
file = 'data/wavs_MFDntn.txt'
pltwav(file,.5)
plt.xlabel('Resolution')

plt.show()
```

**Figure 8**: *Box wisker plots of the dimensionless valley-to-valley wavelength for each mode resolution: (a) SFD cell-to-cell, (b) SFD node-to-node, (c) MFD cell-to-cell, (d) MFD node-to-node.*

The node-to-node SFD algorithm is no better than the cell-to-cell SFD. In this case wind up time is resolution dependent, and the valley spacing increases with increasing resolution (Figure 7b and 8b). For the node-to-node SFD, at a resolution of 256 cells or less along the y-axis there is an instability in the sediment flux output. This is due to the flow tipping between adjacent nodes due to small differences in relative elevation after each time iteration. This unstable behavior disappears for the higher resolution of 512 cells along the y-axis (Figure 7b).

It is only when node-to-node MFD is used that the LEM becomes significantly less resolution dependent (Figure 7d and 8d). For the node-to-node MFD the time evolution of sediment flux is similar for all resolutions, and the valley spacing is similar as resolution is increased. The steady state sediment flux is however not completely stable (Figure 7d). This is due to the migration of the flow across the valley floors created within the model topography. Even once a balance has been achieved between erosion and uplift, small lateral changes in elevation can be seen to create a negative to positive change in elevation of a few meters between time iterations, where the time step is 100 yrs. This is associated with an equivalent change in water flux.

The model that has the least resolution dependence is the node-to-node MFD. The difference between this model and the other three is that this version has the maximum possible flow directions available within my set up. By treating flow paths as lines within the numerical grid, from any node there are 6 paths, which is twice as many as in the cell-to-cell MFD. This means that there is greater distribution of the flow, and a reduced localising of flow paths within the node-to-node distributed model. For SFD increasing resolution however leads to multiple branches.

The grid cells in the models presented are large. At the highest resolution, 2048 by 512 cells, the width of each triangle is of the order of 200 m if I was modelling a landscape 100 km wide. The model is therefore some approximation of local processes that give rise to the large scale landscape. By distributing flow in multiple directions the model is in a sense approximating for the hydrological processes that operate on a sub-grid scale that give rise to the river network. The assumption of SFD is however too strong, and the sub-grid scale processes are ignored.

The transport-limited model that I explore has certain limitations. In particular the valleys floors are wide and not representative of V-shaped valleys that would be expected from fluvial incision into bedrock (e.g. Figure 5). In order to generate such valleys a detachment-limited model, such as the stream power law, would be more appropriate. However, many stream power law models also suffer resolution dependence, as they typically use a SFD to route water (e.g. [Braun & Willett, 2013](https://doi.org/10.1016/j.geomorph.2012.10.008)). [Pelletier, 2010](#) looked at using MFD routing for the stream power law, and found that there remained some spatial resolution dependence. The model of \cite{pelletier-2010} used a rectangular grid and removed resolution issues by using a predictor-corrector algorithm to adjust for resolution effects. However, for the transport-limited model used here, I find that with a triangular grid the MFD routing is resolution independent without addition corrections. This is likely related to the fact that the length of each cell face is equal, while for rectangular cells the diagonal flow direction is longer than the cell faces. The implication is that for LEMs a mesh that has cells with node-to-node spacing of equal length is preferable to a rectangular grid, however this hypothesis will require further exploration.

## Conclusions

From this exploration of the various flow routing algorithms, it can be seen that for the two measures, convergence in wind-up time and valley-to-valley wavelength, a multiple flow direction (MFD) node-to-node routing algorithm is preferable.

```python

```

max_stars_repo_name: johnjarmitage/flem-execarticle

max_issues_repo_name: johnjarmitage/flem-execarticle

max_forks_repo_name: johnjarmitage/flem-execarticle


text:
Concept
==============

$$
\newcommand{\sumN}{\sum_{n = 1}^N}
\newcommand{\sumn}{\sum_n}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\bbeta}{\boldsymbol{\beta}}
\newcommand{\btheta}{\boldsymbol{\theta}}
\newcommand{\bbetahat}{\boldsymbol{\hat{\beta}}}
\newcommand{\bthetahat}{\boldsymbol{\hat{\theta}}}
\newcommand{\dadb}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\bphi}{\boldsymbol{\phi}}
\newcommand{\bPhi}{\boldsymbol{\Phi}}
$$

## Model Structure

*Linear regression* is a relatively simple method that is extremely widely-used. It is also a great stepping stone for more sophisticated methods, making it a natural algorithm to study first.

In linear regression, the target variable $y$ is assumed to follow a linear function of one or more predictor variables, $x_1, \dots, x_D$, plus some random error. Specifically, we assume the model for the $n^\text{th}$ observation in our sample is of the form

$$
y_n = \beta_0 + \beta_1 x_{n1} + \dots + \beta_Dx_{nD} + \epsilon_n.
$$

Here $\beta_0$ is the intercept term, $\beta_1$ through $\beta_D$ are the coefficients on our feature variables, and $\epsilon$ is an error term that represents the difference between the true $y$ value and the linear function of the predictors. Note that the terms with an $n$ in the subscript differ between observations while the terms without (namely the $\beta\text{s}$) do not.

The math behind linear regression often becomes easier when we use vectors to represent our predictors and coefficients. Let's define $\bx_n$ and $\bbeta$ as follows:

$$
\begin{align}
\bx_n &= \begin{pmatrix} 1 & x_{n1} & \dots & x_{nD} \end{pmatrix}^\top \\
\bbeta &= \begin{pmatrix} \beta_0 & \beta_1 & \dots & \beta_D \end{pmatrix}^\top.
\end{align}
$$

Note that $\bx_n$ includes a leading 1, corresponding to the intercept term $\beta_0$. Using these definitions, we can equivalently express $y_n$ as

$$
y_n = \bbeta^\top \bx_n + \epsilon_n.
$$

Below is an example of a dataset designed for linear regression. The input variable is generated randomly and the target variable is generated as a linear combination of that input variable plus an error term.

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

# generate data
np.random.seed(123)
N = 20
beta0 = -4
beta1 = 2
x = np.random.randn(N)
e = np.random.randn(N)
y = beta0 + beta1*x + e
true_x = np.linspace(min(x), max(x), 100)
true_y = beta0 + beta1*true_x

# plot
fig, ax = plt.subplots()
sns.scatterplot(x, y, s = 40, label = 'Data')
sns.lineplot(true_x, true_y, color = 'red', label = 'True Model')
ax.set_xlabel('x', fontsize = 14)
ax.set_title(fr"$y = {beta0} + ${beta1}$x + \epsilon$", fontsize = 16)
ax.set_ylabel('y', fontsize=14, rotation=0, labelpad=10)
ax.legend(loc = 4)
sns.despine()
```

## Parameter Estimation

The previous section covers the entire structure we assume our data follows in linear regression. The machine learning task is then to estimate the parameters in $\bbeta$. These estimates are represented by $\hat{\beta}_0, \dots, \hat{\beta}_D$ or $\bbetahat$. The estimates give us *fitted values* for our target variable, represented by $\hat{y}_n$.

This task can be accomplished in two ways which, though slightly different conceptually, are identical mathematically. The first approach is through the lens of *minimizing loss*. A common practice in machine learning is to choose a loss function that defines how well a model with a given set of parameter estimates the observed data. The most common loss function for linear regression is squared error loss. This says the *loss* of our model is proportional to the sum of squared differences between the true $y_n$ values and the fitted values, $\hat{y}_n$. We then *fit* the model by finding the estimates $\bbetahat$ that minimize this loss function. This approach is covered in the subsection {doc}`s1/loss_minimization`.

The second approach is through the lens of *maximizing likelihood*. Another common practice in machine learning is to model the target as a random variable whose distribution depends on one or more parameters, and then find the parameters that maximize its likelihood. Under this approach, we will represent the target with $Y_n$ since we are treating it as a random variable. The most common model for $Y_n$ in linear regression is a Normal random variable with mean $E(Y_n) = \bbeta^\top \bx_n$. That is, we assume

$$
Y_n|\bx_n \sim \mathcal{N}(\bbeta^\top \bx_n, \sigma^2),
$$

and we find the values of $\bbetahat$ to maximize the likelihood. This approach is covered in subsection {doc}`s1/likelihood_maximization`.

Once we've estimated $\bbeta$, our model is *fit* and we can make predictions. The below graph is the same as the one above but includes our estimated line-of-best-fit, obtained by calculating $\hat{\beta}_0$ and $\hat{\beta}_1$.

```python
# generate data
np.random.seed(123)
N = 20
beta0 = -4
beta1 = 2
x = np.random.randn(N)
e = np.random.randn(N)
y = beta0 + beta1*x + e
true_x = np.linspace(min(x), max(x), 100)
true_y = beta0 + beta1*true_x

# estimate model
beta1_hat = sum((x - np.mean(x))*(y - np.mean(y)))/sum((x - np.mean(x))**2)
beta0_hat = np.mean(y) - beta1_hat*np.mean(x)
fit_y = beta0_hat + beta1_hat*true_x

# plot
fig, ax = plt.subplots()
sns.scatterplot(x, y, s = 40, label = 'Data')
sns.lineplot(true_x, true_y, color = 'red', label = 'True Model')
sns.lineplot(true_x, fit_y, color = 'purple', label = 'Estimated Model')
ax.set_xlabel('x', fontsize = 14)
ax.set_title(fr"Linear Regression for $y = {beta0} + ${beta1}$x + \epsilon$", fontsize = 16)
ax.set_ylabel('y', fontsize=14, rotation=0, labelpad=10)
ax.legend(loc = 4)
sns.despine()
```

## Extensions of Ordinary Linear Regression

There are many important extensions to linear regression which make the model more flexible. Those include {doc}`/content/c2/s1/regularized`—which balances the bias-variance tradeoff for high-dimensional regression models—{doc}`/content/c2/s1/bayesian`—which allows for prior distributions on the coefficients—and {doc}`/content/c2/s1/GLMs`—which introduce non-linearity to regression models. These extensions are discussed in the next chapter.

```{toctree}
:hidden:
:titlesonly:

s1/loss_minimization
s1/likelihood_maximization
```

ext: ipynb

max_stars_repo_name: curioushruti/mlbook

max_issues_repo_name: curioushruti/mlbook

max_forks_repo_name: curioushruti/mlbook


text:
```python
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from numpy import fft

```

# Computer Vision Lecture 1

## 1.1.0) 1D Fourier Transforms

### Basic periodic function: sine waves

$y(t) = A \sin(\omega t + \phi)$

$y(x) = A \sin(k t + \phi)$

where $A$ is the amplitude, $\omega$ is $\textit{time}$-frequency, $k$ is a $\textit{spatial}$-frequency and $\phi$ is the phase shift

```python
f_s = 50 # sample rate / "resolution"
f = 1 # the frequency of the signal
A = 1 # Amplitude
p = 0 # phase

x = np.arange(-f_s,(f_s+1)) # points for horizontal axis
# compute the sin wave
y_0= [A*np.sin(2*np.pi*f*(i/f_s) + p*np.pi/180) for i in np.arange(-f_s,(f_s+1))]
y_1 = [A*np.sin(2*2*np.pi*f*(i/f_s) + p*np.pi/180) for i in np.arange(-f_s,(f_s+1))]
y_2 = [A*np.sin(3*2*np.pi*f*(i/f_s) + p*np.pi/180) for i in np.arange(-f_s,(f_s+1))]

# plot the sine wave
plt.subplot(3,1,1)
plt.stem(x,y_0,'r',)
plt.plot(x,y_0)

plt.subplot(3,1,2)
plt.stem(x,y_1,'r',)
plt.plot(x,y_1)

plt.subplot(3,1,3)
plt.stem(x,y_2,'r',)
plt.plot(x,y_2)

plt.show()
```

### 1.1.1) Principle of Superposition:
#### When $n$ waves interact the resulting wave is the sum of the interacting waves.

```python
A3 = 1
A4 = 1
A5 = 1
p3 = 0
p4 = 0
p5 = 0

y_total = [A3*np.sin(2*np.pi*f*(i/f_s) + p3*np.pi/180) +
           A4*np.sin(2*2*np.pi*f*(i/f_s) + p4*np.pi/180) +
           A5*np.sin(3*2*np.pi*f*(i/f_s) + p5*np.pi/180)
           for i in np.arange(-f_s,(f_s+1))]

plt.stem(x,y_total,'r',)
plt.plot(x,y_total)
plt.title('summation of three sine waves above')
plt.show()
```

#### So, if we can add sine waves, then we can take them apart, ie a Fourier decomposition / transform
##### The goal: given a periodic function, find a series of sine waves that if added up would give that periodic function

```python
# parameter settings for the sine waves
A3 = 1.8
A4 = 2.4
A5 = .6
p3 = 0
p4 = 0
p5 = 0
f_2 = 3
f_s2 = 300
x2 = np.arange(-f_s2,(f_s2+1))

# addition of three sine waves
y_total = [A3*np.sin(2*np.pi*f_2*(i/f_s2) + p3*np.pi/180) +
           A4*np.sin(2*2*np.pi*f_2*(i/f_s2) + p4*np.pi/180) +
           A5*np.sin(3*2*np.pi*f_2*(i/f_s2) + p5*np.pi/180) + 4.2
           for i in np.arange(-f_s2,(f_s2+1))]

y_0 = [3.5 for i in np.arange(-f_s2,(f_s2+1))]

y_fund = [3.5+ 2*np.sin(6*np.pi*f*(i/f_s2) + 0*np.pi/180) for i in np.arange(-f_s2,(f_s2+1))]
y_2nd = [4+ 1*np.sin(12*np.pi*f*(i/f_s2) + 0*np.pi/180) for i in np.arange(-f_s2,(f_s2+1))]
y_3rd = [3+ .5*np.sin(18*np.pi*f*(i/f_s2) + 0*np.pi/180) for i in np.arange(-f_s2,(f_s2+1))]

plt.subplot(10,1,1)
plt.plot(x2,y_total)
plt.yticks(np.arange(0,10,2.5))
plt.title('sum of periodic functions')

plt.subplot(10,1,4)
plt.plot(x2, y_0)
plt.yticks(np.arange(0,10,2.5))
plt.title('Scalar value (measures avg value): $\lambda = \infty$')
plt.show()

plt.subplot(10,1,6)
plt.plot(x2, y_fund)
plt.yticks(np.arange(0,10,2.5))
plt.title('Fundamental (sine wave = 1 period (sample_rate/freq)): $\lambda$ = a')
plt.show()

plt.subplot(10,1,8)
plt.plot(x2, y_2nd)
plt.yticks(np.arange(0,10,2.5))
plt.title('sine wave = 2 periods (sample_rate/freq): $\lambda$ = a/2')
plt.show()

plt.subplot(10,1,10)
plt.plot(x2, y_3rd)
plt.yticks(np.arange(0,10,2.5))
plt.title('sine wave = 3 periods (sample_rate/freq): $\lambda$ = a/3')
plt.show()
```

#### So on and on....
##### It matters how many decomposition terms it takes to acheive the solution. The number of terms used is known as the resolution. How exact your Fourier decomposition depends on the number of terms used (resolution).

##### Given any continuous periodic function, it can be decomposed into a series of sine and cosine functions:

\begin{equation}
f(x) = \frac{A_0}{2} + \sum_{n=1}^{\infty} A_n cos(\frac{2\pi nx}{\lambda}) + \sum_{n=1}^{\infty} B_n sin(\frac{2\pi nx}{\lambda})
\end{equation}

where,

\begin{equation}
A_n = \frac{2}{\lambda} \int_{0}^{\lambda} f(x) cos(\frac{2\pi nx}{\lambda}) dx
\end{equation}

\begin{equation}
B_n = \frac{2}{\lambda} \int_{0}^{\lambda} f(x) sin(\frac{2\pi nx}{\lambda}) dx
\end{equation}

the amplitudes tell you how much of the sine-wave, $n$ is in the continuous periodic function, $f(x)$

## 1.2.0) 1D Reciprocal Space

##### calculate the FFT on a 1D array

```python
one_d_pixel_vals =  np.array([0.3, 0.4, 1.3, 3.4, 4.5, 4.2, 2.8, 2.4, 1.4, 0.1])
print(one_d_pixel_vals)

plt.plot(one_d_pixel_vals)
plt.yticks(np.arange(5))
plt.xticks(np.arange(9))
plt.grid(True)
plt.show()

```

##### Returns a list of 10 numbers (+2 for the "DC" component): $(A_0, P_0, A_1, P_1, A_2, P_2, A_3, P_3, A_4, P_4, A_5, P_5)$

$A_0$ = Amplitude of the average value, "DC" $\lambda = \infty$ component

$P_0$ = Flat wave, so phase doesn't say anything meaningful

$A_1$ = Amplitude of the fundamental, $\lambda = 1a$

$P_1$ = Phase of the fundamental frequency component

$A_2$ = Amplitude of the first harmonic $\lambda = 2a$

$P_2$ = Phase of the first harmonic frequency component

$A_3$ = Amplitude of the second harmonic $\lambda = 3a$

$P_3$ = Phase of the second harmonic frequency component

$A_4$ = Amplitude of the third harmonic $\lambda = 4a$

$P_4$ = Phase of the third harmonic frequency component

$A_5$ = Amplitude of the Nyquist frequency

$P_5$ = Phase of the Nyquist

where $a$ is the period. Nyquist frequency is the highest frequency component that is present in the image. An array of length 10 is 10 pixels. The fastest frequency possible is one in which the max is in pixel $i$ and the min is in pixel $i+1$. For images the Nyquist frequency ( $f_{Nyq} = pixel_{total}/2$ )

#### 1.2.1) Python np.fft.fft output

For numpy fft, the output for a given array, $a$; $a[0]$ is the mean. $a[1]$ = np.conjugate($a[-1]$)

where $a[k],a[-k]$ represent the $k_{th}$ harmonic

So the ouptput for a np.fft.fft($a$), where $a$ is 10 pixels contains amplitudes only:

$a[0]$ = mean

$a[1]$ =  $a[-1]$ = $a[9]$ = fundamental (at $\pm k$)

$a[2]$ =  $a[-2]$ = $a[8]$ = first harmonic (at $\pm k$)

$a[3]$ =  $a[-3]$ = $a[7]$ = second harmonic (at $\pm k$)

$a[4]$ =  $a[-4]$ = $a[6]$ = third harmonic (at $\pm k$)

$a[5]$ =  Nyquist Frequency (at $\pm k$)

```python
one_d_fft = np.fft.fft(one_d_pixel_vals, norm='ortho')
print("for a 1D array given by: ", one_d_pixel_vals)
print("The fft returns an array of length {} with values : ".format(len(one_d_fft)) )
print(one_d_fft)

```

    for a 1D array given by:  [ 0.3  0.4  1.3  3.4  4.5  4.2  2.8  2.4  1.4  0.1]
    The fft returns an array of length 10 with values :
    [ 6.57753753 +0.00000000e+00j -3.27588825 -6.42423464e-01j
      0.01065311 +6.25512087e-01j  0.22429031 -4.34214340e-01j
      0.25814049 -7.06019011e-02j -0.06324555 -1.12346671e-15j
      0.25814049 +7.06019011e-02j  0.22429031 +4.34214340e-01j
      0.01065311 -6.25512087e-01j -3.27588825 +6.42423464e-01j]

```python
# Plot FFT and explain reciprocal space here

amplitudes =[]
for i in range(round(one_d_fft.size/2)+1):
        amplitudes.append(np.real(one_d_fft[i]))
plt.bar(np.arange(len(amplitudes)), np.abs(amplitudes))
plt.title("Plot of the amplitudes of the fft for +k (reciprocal space)")
plt.show()
```

#### To find the frequencies of the fft, look at the code below

```python
x = np.array([0.3, 0.4, 1.3, 3.4, 4.5, 4.2, 2.8, 2.4, 1.4, 0.1])
w = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x))

for coef,freq in zip(w,freqs):
    if coef:
        print('{c:>6} * exp(2 pi i t * {f})'.format(c=round(coef,2),f=round(freq,2)))

```

    (20.8+0j) * exp(2 pi i t * 0.0)
    (-10.36-2.03j) * exp(2 pi i t * 0.1)
    (0.03+1.98j) * exp(2 pi i t * 0.2)
    (0.71-1.37j) * exp(2 pi i t * 0.3)
    (0.82-0.22j) * exp(2 pi i t * 0.4)
    (-0.2-0j) * exp(2 pi i t * -0.5)
    (0.82+0.22j) * exp(2 pi i t * -0.4)
    (0.71+1.37j) * exp(2 pi i t * -0.3)
    (0.03-1.98j) * exp(2 pi i t * -0.2)
    (-10.36+2.03j) * exp(2 pi i t * -0.1)

### HW question 1 - for common 512 by 512 image, what is the Nyquiest Frequency?

```python
# put answer here
```

   ### HW question 2 - What is the inverse Fourier Transform of the array "one_d_fft"?

```python
# put answer here
```

### HW question 3 - Find the fft amplitudes and frequencies for $A \cos(kx)$, $A(3+ \cos(3kx))$, your own function

```python
# put answer here
```

### HW question 4 - Plot the array "one_d_pixel_vals" and the plane waves of the FFT on the same plot

```python
# put answer here
```

### *HW question 5 - plot the sum of the three sine waves in section 1.1.1 for 100, 50, 25, and 12  pixels per period. Take the FFT for each different sampling rate / resolution. Plot the amplitudes in reciprocal space. What do you see?

```python
# put answer here
```


text:
# Spectrum

[Colour](http://en.wikipedia.org/wiki/Colour) is defined as the characteristic of visual perception that can be described by attributes of hue, brightness (or lightness) and colourfulness (or saturation or chroma).

When necessary, to avoid confusion between other meanings of the word, the term "perceived colour" may be used.

Perceived colour depends on the spectral distribution of the colour stimulus, on the size, shape, structure and surround of the stimulus area, on the state of adaptation of the observer's visual system, and on the observer's experience of the prevailing and similar situations of observation. [1]

[Light](http://en.wikipedia.org/wiki/Light) is the electromagnetic radiation that is considered from the point of view of its ability to excite the human visual system. [2]

The portion of the electromatic radiation frequencies perceived in the approximate wavelength range 360-780 nanometres (nm) is called the [visible spectrum](http://en.wikipedia.org/wiki/Visible_spectrum).

```python
import colour

from colour.plotting import *
```

```python
colour_style();
```

```python
# Plotting the visible spectrum.
plot_visible_spectrum();
```

The spectrum is defined as the display or specification of the monochromatic components of the radiation considered. [3]

At the core of [Colour](https://github.com/colour-science/colour/) is the `colour.colorimetry` sub-package, it defines the objects needed for spectral related computations and many others:

```python
from pprint import pprint

import colour.colorimetry as colorimetry

pprint(colorimetry.__all__)
```

    ['SpectralShape',
     'DEFAULT_SPECTRAL_SHAPE',
     'SpectralDistribution',
     'MultiSpectralDistributions',
     'sds_and_multi_sds_to_sds',
     'sd_blackbody',
     'blackbody_spectral_radiance',
     'planck_law',
     'LMS_ConeFundamentals',
     'RGB_ColourMatchingFunctions',
     'XYZ_ColourMatchingFunctions',
     'CMFS',
     'LMS_CMFS',
     'RGB_CMFS',
     'STANDARD_OBSERVERS_CMFS',
     'ILLUMINANTS',
     'D_ILLUMINANTS_S_SDS',
     'HUNTERLAB_ILLUMINANTS',
     'ILLUMINANTS_SDS',
     'LIGHT_SOURCES',
     'LIGHT_SOURCES_SDS',
     'LEFS',
     'PHOTOPIC_LEFS',
     'SCOTOPIC_LEFS',
     'sd_constant',
     'sd_zeros',
     'sd_ones',
     'SD_GAUSSIAN_METHODS',
     'sd_gaussian',
     'sd_gaussian_normal',
     'sd_gaussian_fwhm',
     'SD_SINGLE_LED_METHODS',
     'sd_single_led',
     'sd_single_led_Ohno2005',
     'SD_MULTI_LEDS_METHODS',
     'sd_multi_leds',
     'sd_multi_leds_Ohno2005',
     'SD_TO_XYZ_METHODS',
     'MULTI_SD_TO_XYZ_METHODS',
     'sd_to_XYZ',
     'multi_sds_to_XYZ',
     'ASTME308_PRACTISE_SHAPE',
     'lagrange_coefficients_ASTME2022',
     'tristimulus_weighting_factors_ASTME2022',
     'adjust_tristimulus_weighting_factors_ASTME308',
     'sd_to_XYZ_integration',
     'sd_to_XYZ_tristimulus_weighting_factors_ASTME308',
     'sd_to_XYZ_ASTME308',
     'multi_sds_to_XYZ_integration',
     'multi_sds_to_XYZ_ASTME308',
     'wavelength_to_XYZ',
     'BANDPASS_CORRECTION_METHODS',
     'bandpass_correction',
     'bandpass_correction_Stearns1988',
     'sd_CIE_standard_illuminant_A',
     'sd_CIE_illuminant_D_series',
     'daylight_locus_function',
     'sd_mesopic_luminous_efficiency_function',
     'mesopic_weighting_function',
     'LIGHTNESS_METHODS',
     'lightness',
     'lightness_Glasser1958',
     'lightness_Wyszecki1963',
     'lightness_CIE1976',
     'lightness_Fairchild2010',
     'lightness_Fairchild2011',
     'intermediate_lightness_function_CIE1976',
     'LUMINANCE_METHODS',
     'luminance',
     'luminance_Newhall1943',
     'luminance_ASTMD1535',
     'luminance_CIE1976',
     'luminance_Fairchild2010',
     'luminance_Fairchild2011',
     'intermediate_luminance_function_CIE1976',
     'dominant_wavelength',
     'complementary_wavelength',
     'excitation_purity',
     'colorimetric_purity',
     'luminous_flux',
     'luminous_efficiency',
     'luminous_efficacy',
     'RGB_10_degree_cmfs_to_LMS_10_degree_cmfs',
     'RGB_2_degree_cmfs_to_XYZ_2_degree_cmfs',
     'RGB_10_degree_cmfs_to_XYZ_10_degree_cmfs',
     'LMS_2_degree_cmfs_to_XYZ_2_degree_cmfs',
     'LMS_10_degree_cmfs_to_XYZ_10_degree_cmfs',
     'WHITENESS_METHODS',
     'whiteness',
     'whiteness_Berger1959',
     'whiteness_Taube1960',
     'whiteness_Stensby1968',
     'whiteness_ASTME313',
     'whiteness_Ganz1979',
     'whiteness_CIE2004',
     'YELLOWNESS_METHODS',
     'yellowness',
     'yellowness_ASTMD1925',
     'yellowness_ASTME313']

> Note: `colour.colorimetry` sub-package public API is directly available from `colour` namespace.

[Colour](https://github.com/colour-science/colour/) computations are based on a comprehensive dataset available in pretty much each sub-packages, for example `colour.colorimetry.dataset` defines the following data:

```python
import colour.colorimetry.datasets as datasets

pprint(datasets.__all__)
```

    ['CMFS',
     'LMS_CMFS',
     'RGB_CMFS',
     'STANDARD_OBSERVERS_CMFS',
     'ILLUMINANTS',
     'D_ILLUMINANTS_S_SDS',
     'HUNTERLAB_ILLUMINANTS',
     'ILLUMINANTS_SDS',
     'LIGHT_SOURCES',
     'LIGHT_SOURCES_SDS',
     'LEFS',
     'PHOTOPIC_LEFS',
     'SCOTOPIC_LEFS']

> Note: `colour.colorimetry.dataset` sub-package public API is directly available from `colour` namespace.

## Spectral Distribution

Whether it be a sample spectral distribution, colour matching functions or illuminants, spectral data is manipulated using an object built with the `colour.SpectralDistribution` class or based on it:

```python
import colour

# Defining a sample spectral distribution data.
sample_sd_data = {
    380: 0.048,
    385: 0.051,
    390: 0.055,
    395: 0.06,
    400: 0.065,
    405: 0.068,
    410: 0.068,
    415: 0.067,
    420: 0.064,
    425: 0.062,
    430: 0.059,
    435: 0.057,
    440: 0.055,
    445: 0.054,
    450: 0.053,
    455: 0.053,
    460: 0.052,
    465: 0.052,
    470: 0.052,
    475: 0.053,
    480: 0.054,
    485: 0.055,
    490: 0.057,
    495: 0.059,
    500: 0.061,
    505: 0.062,
    510: 0.065,
    515: 0.067,
    520: 0.070,
    525: 0.072,
    530: 0.074,
    535: 0.075,
    540: 0.076,
    545: 0.078,
    550: 0.079,
    555: 0.082,
    560: 0.087,
    565: 0.092,
    570: 0.100,
    575: 0.107,
    580: 0.115,
    585: 0.122,
    590: 0.129,
    595: 0.134,
    600: 0.138,
    605: 0.142,
    610: 0.146,
    615: 0.150,
    620: 0.154,
    625: 0.158,
    630: 0.163,
    635: 0.167,
    640: 0.173,
    645: 0.180,
    650: 0.188,
    655: 0.196,
    660: 0.204,
    665: 0.213,
    670: 0.222,
    675: 0.231,
    680: 0.242,
    685: 0.251,
    690: 0.261,
    695: 0.271,
    700: 0.282,
    705: 0.294,
    710: 0.305,
    715: 0.318,
    720: 0.334,
    725: 0.354,
    730: 0.372,
    735: 0.392,
    740: 0.409,
    745: 0.420,
    750: 0.436,
    755: 0.450,
    760: 0.462,
    765: 0.465,
    770: 0.448,
    775: 0.432,
    780: 0.421}

sd = colour.SpectralDistribution(sample_sd_data, name='Sample')
print(sd)
```

    [[  3.80000000e+02   4.80000000e-02]
     [  3.85000000e+02   5.10000000e-02]
     [  3.90000000e+02   5.50000000e-02]
     [  3.95000000e+02   6.00000000e-02]
     [  4.00000000e+02   6.50000000e-02]
     [  4.05000000e+02   6.80000000e-02]
     [  4.10000000e+02   6.80000000e-02]
     [  4.15000000e+02   6.70000000e-02]
     [  4.20000000e+02   6.40000000e-02]
     [  4.25000000e+02   6.20000000e-02]
     [  4.30000000e+02   5.90000000e-02]
     [  4.35000000e+02   5.70000000e-02]
     [  4.40000000e+02   5.50000000e-02]
     [  4.45000000e+02   5.40000000e-02]
     [  4.50000000e+02   5.30000000e-02]
     [  4.55000000e+02   5.30000000e-02]
     [  4.60000000e+02   5.20000000e-02]
     [  4.65000000e+02   5.20000000e-02]
     [  4.70000000e+02   5.20000000e-02]
     [  4.75000000e+02   5.30000000e-02]
     [  4.80000000e+02   5.40000000e-02]
     [  4.85000000e+02   5.50000000e-02]
     [  4.90000000e+02   5.70000000e-02]
     [  4.95000000e+02   5.90000000e-02]
     [  5.00000000e+02   6.10000000e-02]
     [  5.05000000e+02   6.20000000e-02]
     [  5.10000000e+02   6.50000000e-02]
     [  5.15000000e+02   6.70000000e-02]
     [  5.20000000e+02   7.00000000e-02]
     [  5.25000000e+02   7.20000000e-02]
     [  5.30000000e+02   7.40000000e-02]
     [  5.35000000e+02   7.50000000e-02]
     [  5.40000000e+02   7.60000000e-02]
     [  5.45000000e+02   7.80000000e-02]
     [  5.50000000e+02   7.90000000e-02]
     [  5.55000000e+02   8.20000000e-02]
     [  5.60000000e+02   8.70000000e-02]
     [  5.65000000e+02   9.20000000e-02]
     [  5.70000000e+02   1.00000000e-01]
     [  5.75000000e+02   1.07000000e-01]
     [  5.80000000e+02   1.15000000e-01]
     [  5.85000000e+02   1.22000000e-01]
     [  5.90000000e+02   1.29000000e-01]
     [  5.95000000e+02   1.34000000e-01]
     [  6.00000000e+02   1.38000000e-01]
     [  6.05000000e+02   1.42000000e-01]
     [  6.10000000e+02   1.46000000e-01]
     [  6.15000000e+02   1.50000000e-01]
     [  6.20000000e+02   1.54000000e-01]
     [  6.25000000e+02   1.58000000e-01]
     [  6.30000000e+02   1.63000000e-01]
     [  6.35000000e+02   1.67000000e-01]
     [  6.40000000e+02   1.73000000e-01]
     [  6.45000000e+02   1.80000000e-01]
     [  6.50000000e+02   1.88000000e-01]
     [  6.55000000e+02   1.96000000e-01]
     [  6.60000000e+02   2.04000000e-01]
     [  6.65000000e+02   2.13000000e-01]
     [  6.70000000e+02   2.22000000e-01]
     [  6.75000000e+02   2.31000000e-01]
     [  6.80000000e+02   2.42000000e-01]
     [  6.85000000e+02   2.51000000e-01]
     [  6.90000000e+02   2.61000000e-01]
     [  6.95000000e+02   2.71000000e-01]
     [  7.00000000e+02   2.82000000e-01]
     [  7.05000000e+02   2.94000000e-01]
     [  7.10000000e+02   3.05000000e-01]
     [  7.15000000e+02   3.18000000e-01]
     [  7.20000000e+02   3.34000000e-01]
     [  7.25000000e+02   3.54000000e-01]
     [  7.30000000e+02   3.72000000e-01]
     [  7.35000000e+02   3.92000000e-01]
     [  7.40000000e+02   4.09000000e-01]
     [  7.45000000e+02   4.20000000e-01]
     [  7.50000000e+02   4.36000000e-01]
     [  7.55000000e+02   4.50000000e-01]
     [  7.60000000e+02   4.62000000e-01]
     [  7.65000000e+02   4.65000000e-01]
     [  7.70000000e+02   4.48000000e-01]
     [  7.75000000e+02   4.32000000e-01]
     [  7.80000000e+02   4.21000000e-01]]

The sample spectral distribution can be easily plotted against the visible spectrum:

```python
# Plotting the sample spectral distribution.
plot_single_sd(sd);
```

With the sample spectral distribution defined, we can retrieve its shape:

```python
# Displaying the sample spectral distribution shape.
print(sd.shape)
```

    (380.0, 780.0, 5.0)

The shape returned is an instance of `colour.SpectralShape` class:

```python
repr(sd.shape)
```

    'SpectralShape(380.0, 780.0, 5.0)'

The `colour.SpectralShape` class is used throughout [Colour](https://github.com/colour-science/colour/) to define spectral dimensions and is instantiated as follows:

```python
# Using *colour.SpectralShape* with iteration.
shape = colour.SpectralShape(start=0, end=10, interval=1)
for wavelength in shape:
    print(wavelength)

# *colour.SpectralShape.range* method is providing the complete range of values.
shape = colour.SpectralShape(0, 10, 0.5)
shape.range()
```

    0.0
    1.0
    2.0
    3.0
    4.0
    5.0
    6.0
    7.0
    8.0
    9.0
    10.0

    array([  0. ,   0.5,   1. ,   1.5,   2. ,   2.5,   3. ,   3.5,   4. ,
             4.5,   5. ,   5.5,   6. ,   6.5,   7. ,   7.5,   8. ,   8.5,
             9. ,   9.5,  10. ])

[Colour](https://github.com/colour-science/colour/) defines three convenient objects to create constant spectral distributions:

*  `colour.sd_constant`
*  `colour.sd_zeros`
*  `colour.sd_ones`

```python
# Defining a constant spectral distribution.
constant_sd = colour.sd_constant(100)
print('"Constant Spectral Distribution"')
print(constant_sd.shape)
print(constant_sd[400])

# Defining a zeros filled spectral distribution.
print('\n"Zeros Filled Spectral Distribution"')
zeros_sd = colour.sd_zeros()
print(zeros_sd.shape)
print(zeros_sd[400])

# Defining a ones filled spectral distribution.
print('\n"Ones Filled Spectral Distribution"')
ones_sd = colour.sd_ones()
print(ones_sd.shape)
print(ones_sd[400])
```

    "Constant Spectral Distribution"
    (360.0, 780.0, 1.0)
    100.0

    "Zeros Filled Spectral Distribution"
    (360.0, 780.0, 1.0)
    0.0

    "Ones Filled Spectral Distribution"
    (360.0, 780.0, 1.0)
    1.0

By default the shape used by `colour.sd_constant`, `colour.sd_zeros` and `colour.sd_ones` is the one defined by `colour.DEFAULT_SPECTRAL_SHAPE` attribute using the *CIE 1931 2° Standard Observer* shape.

```python
print(repr(colour.DEFAULT_SPECTRAL_SHAPE))
```

    SpectralShape(360, 780, 1)

A custom shape can be passed to construct a constant spectral distribution with tailored dimensions:

```python
colour.sd_ones(colour.SpectralShape(400, 700, 5))[450]
```

    1.0

Often interpolation of the spectral distribution is needed, this is achieved with the `colour.SpectralDistribution.interpolate` method. Depending on the wavelengths uniformity, the default interpolation method will differ. Following *CIE 167:2005* recommendation: The method developed by Sprague (1880) should be used for interpolating functions having a uniformly spaced independent variable and a *Cubic Spline* method for non-uniformly spaced independent variable. [4]

We can check the uniformity of the sample spectral distribution:

```python
# Checking the sample spectral distribution uniformity.
print(sd.is_uniform())
```

    True

Since the sample spectral distribution is uniform the interpolation will be using the `colour.SpragueInterpolator` interpolator.

> Note: Interpolation happens in place and may alter your original data, use the `colour.SpectralDistribution.copy` method to produce a copy of your spectral distribution before interpolation.

```python
# Copying the sample spectral distribution.
sd_copy = sd.copy()

# Interpolating the copied sample spectral distribution.
sd_copy.interpolate(colour.SpectralShape(400, 770, 1))
sd_copy[401]
```

    0.065809599999999996

```python
# Comparing the interpolated spectral distribution with the original one.
plot_multi_sds([sd, sd_copy], bounding_box=[730,780, 0.1, 0.5]);
```

Extrapolation although dangerous can be used to help aligning two spectral distributions together. *CIE 015:2004 Colorimetry, 3rd Edition* recommends that unmeasured values may be set equal to the nearest measured value of the appropriate quantity in truncation: [5]

```python
# Extrapolating the copied sample spectral distribution.
sd_copy.extrapolate(colour.SpectralShape(340, 830))
sd_copy[340], sd_copy[830]
```

    (0.064999999999999947, 0.44800000000000018)

The underlying interpolator can be swapped for any of the *Colour* interpolators.

```python
pprint([
    export for export in colour.algebra.interpolation.__all__
    if 'Interpolator' in export
])
```

    ['KernelInterpolator',
     'NearestNeighbourInterpolator',
     'LinearInterpolator',
     'SpragueInterpolator',
     'CubicSplineInterpolator',
     'PchipInterpolator',
     'NullInterpolator']

```python
# Changing interpolator while trimming the copied spectral distribution.
sd_copy.interpolate(
    colour.SpectralShape(400, 700, 10), interpolator=colour.LinearInterpolator)
```

    SpectralDistribution([[  4.00000000e+02,   6.50000000e-02],
                          [  4.10000000e+02,   6.80000000e-02],
                          [  4.20000000e+02,   6.40000000e-02],
                          [  4.30000000e+02,   5.90000000e-02],
                          [  4.40000000e+02,   5.50000000e-02],
                          [  4.50000000e+02,   5.30000000e-02],
                          [  4.60000000e+02,   5.20000000e-02],
                          [  4.70000000e+02,   5.20000000e-02],
                          [  4.80000000e+02,   5.40000000e-02],
                          [  4.90000000e+02,   5.70000000e-02],
                          [  5.00000000e+02,   6.10000000e-02],
                          [  5.10000000e+02,   6.50000000e-02],
                          [  5.20000000e+02,   7.00000000e-02],
                          [  5.30000000e+02,   7.40000000e-02],
                          [  5.40000000e+02,   7.60000000e-02],
                          [  5.50000000e+02,   7.90000000e-02],
                          [  5.60000000e+02,   8.70000000e-02],
                          [  5.70000000e+02,   1.00000000e-01],
                          [  5.80000000e+02,   1.15000000e-01],
                          [  5.90000000e+02,   1.29000000e-01],
                          [  6.00000000e+02,   1.38000000e-01],
                          [  6.10000000e+02,   1.46000000e-01],
                          [  6.20000000e+02,   1.54000000e-01],
                          [  6.30000000e+02,   1.63000000e-01],
                          [  6.40000000e+02,   1.73000000e-01],
                          [  6.50000000e+02,   1.88000000e-01],
                          [  6.60000000e+02,   2.04000000e-01],
                          [  6.70000000e+02,   2.22000000e-01],
                          [  6.80000000e+02,   2.42000000e-01],
                          [  6.90000000e+02,   2.61000000e-01],
                          [  7.00000000e+02,   2.82000000e-01]],
                         interpolator=SpragueInterpolator,
                         interpolator_args={},
                         extrapolator=Extrapolator,
                         extrapolator_args={'method': 'Constant', 'left': None, 'right': None})

The extrapolation behaviour can be changed for `Linear` method instead of the `Constant` default method or even use arbitrary constant `left` and `right` values:

```python
# Extrapolating the copied sample spectral distribution with *Linear* method.
sd_copy.extrapolate(
    colour.SpectralShape(340, 830),
    extrapolator_args={'method': 'Linear',
                       'right': 0})
sd_copy[340], sd_copy[830]
```

    (0.046999999999999341, 0.0)

Aligning a spectral distribution is a convenient way to first interpolate the current data within its original bounds then if needed extrapolates any missing values to match the requested shape:

```python
# Aligning the cloned sample spectral distribution.
# We first trim the spectral distribution as above.
sd_copy.interpolate(colour.SpectralShape(400, 700))
sd_copy.align(colour.SpectralShape(340, 830, 5))
sd_copy[340], sd_copy[830]
```

    (0.064999999999999974, 0.28199999999999986)

The `colour.SpectralDistribution` class also supports various arithmetic operations like *addition*, *subtraction*, *multiplication*, *division* or *exponentiation* with *numeric* and *array_like* variables or other `colour.SpectralDistribution` class instances:

```python
sd = colour.SpectralDistribution({
    410: 0.25,
    420: 0.50,
    430: 0.75,
    440: 1.0,
    450: 0.75,
    460: 0.50,
    480: 0.25
})

print((sd.copy() + 1).values)
print((sd.copy() * 2).values)
print((sd * [0.35, 1.55, 0.75, 2.55, 0.95, 0.65, 0.15]).values)
print((sd * colour.sd_constant(2, sd.shape) * colour.sd_constant(3, sd.shape)).values)
```

    [ 1.25  1.5   1.75  2.    1.75  1.5   1.25]
    [ 0.5  1.   1.5  2.   1.5  1.   0.5]
    [ 0.0875  0.775   0.5625  2.55    0.7125  0.325   0.0375]
    [ 1.5  3.   4.5  6.   4.5  3.   nan  1.5]

The spectral distribution can be normalised with an arbitrary factor:

```python
print(sd.normalise().values)
print(sd.normalise(100).values)
```

    [ 0.25  0.5   0.75  1.    0.75  0.5   0.25]
    [  25.   50.   75.  100.   75.   50.   25.]

## Colour Matching Functions

In the late 1920's, Wright (1928) and Guild (1931) independently conducted a series of colour matching experiments to quantify the colour ability of an average human observer which laid the foundation for the specification of the [CIE XYZ colourspace](http://en.wikipedia.org/wiki/CIE_color_space#Definition_of_the_CIE_XYZ_color_space). The results obtained were summarized by the *Wright & Guild 1931 2° RGB CMFs* $\bar{r}(\lambda)$,$\bar{g}(\lambda)$,$\bar{b}(\lambda)$ colour matching functions: they represent the amounts of three monochromatic primary colours $\textbf{R}$,$\textbf{G}$,$\textbf{B}$ needed to match the test colour at a single wavelength of light.

> See Also: The [Colour Matching Functions](cmfs.ipynb) notebook for in-depth information about the colour matching functions.

```python
# Plotting *Wright & Guild 1931 2 Degree RGB CMFs* colour matching functions.
plot_single_cmfs('Wright & Guild 1931 2 Degree RGB CMFs');
```

With an RGB model of human vision based on *Wright & Guild 1931 2° RGB CMFs* $\bar{r}(\lambda)$,$\bar{g}(\lambda)$,$\bar{b}(\lambda)$ colour matching functions and for pragmatic reasons the *CIE* members developed a new colour space that would relate to the *CIE RGB* colourspace but for which all tristimulus values would be positive for real colours: *CIE XYZ* described with $\bar{x}(\lambda)$,$\bar{y}(\lambda)$,$\bar{z}(\lambda)$ colour matching functions.

```python
# Plotting *CIE XYZ 1931 2 Degree Standard Observer* colour matching functions.
plot_single_cmfs('CIE 1931 2 Degree Standard Observer');
```

In the 1960's it appeared that cones were present in a larger region of eye than the one initially covered by the experiments that lead to the *CIE 1931 2° Standard Observer* specification.

As a result, colour computations done with the *CIE 1931 2° Standard Observer* do not always correlate to the visual observation.

In 1964, the *CIE* defined an additional standard observer: the *CIE 1964 10° Standard Observer* derived from the work of Stiles and Burch (1959), and Speranskaya (1959). The *CIE 1964 10° Standard Observer* is believed to be a better representation of the human vision spectral response and recommended when dealing with a field of view of more than 4°.

For example and as per *CIE* recommendation, the *CIE 1964 10° Standard Observer* is commonly used with spectrophotometers for colour measurements whereas colorimeters generally use the *CIE 1931 2° Standard Observer* for quality control and other colour evaluation applications.

## CIE XYZ Tristimulus Values

The *CIE XYZ* tristimulus values specify a colour stimulus in terms of the visual system. Their values for colour of a surface with spectral reflectance $\beta(\lambda)$ under an illuminant of spectral $S(\lambda)$ are calculated using the following equations: [6]

$$
\begin{equation}
X=k\int_{\lambda}\beta(\lambda)S(\lambda)\bar{x}(\lambda)d\lambda\\
Y=k\int_{\lambda}\beta(\lambda)S(\lambda)\bar{y}(\lambda)d\lambda\\
Z=k\int_{\lambda}\beta(\lambda)S(\lambda)\bar{z}(\lambda)d\lambda
\end{equation}
$$
where
$$
\begin{equation}
k=\cfrac{100}{\int_{\lambda}S(\lambda)\bar{y}(\lambda)d\lambda}
\end{equation}
$$

However in virtually all practical computations of *CIE XYZ* tristimulus values, the integrals are replaced by summations:

$$
\begin{equation}
X=k\sum\limits_{\lambda=\lambda_a}^{\lambda_b}\beta(\lambda)S(\lambda)\bar{x}(\lambda)\Delta\lambda\\
Y=k\sum\limits_{\lambda=\lambda_a}^{\lambda_b}\beta(\lambda)S(\lambda)\bar{y}(\lambda)\Delta\lambda\\
Z=k\sum\limits_{\lambda=\lambda_a}^{\lambda_b}\beta(\lambda)S(\lambda)\bar{z}(\lambda)\Delta\lambda\\
\end{equation}
$$
where
$$
\begin{equation}
k=\cfrac{100}{\sum\limits_{\lambda=\lambda_a}^{\lambda_b}S(\lambda)\bar{y}(\lambda)\Delta\lambda}
\end{equation}
$$

Calculating the *CIE XYZ* tristimulus values of a colour stimulus is done using the `colour.sd_to_XYZ` definition which follows *ASTM E2022–11* and *ASTM E308–15* practises computation method:

```python
sd = colour.SpectralDistribution(sample_sd_data, name='Sample')
cmfs = colour.STANDARD_OBSERVERS_CMFS['CIE 1931 2 Degree Standard Observer']
illuminant = colour.ILLUMINANTS_SDS['A']

# Calculating the sample spectral distribution *CIE XYZ* tristimulus values.
colour.sd_to_XYZ(sd, cmfs, illuminant)
```

    array([ 14.78676625,  10.97815867,   1.99023459])

> Note: Output *CIE XYZ* colourspace matrix is in domain [0, 100].

*CIE XYZ* tristimulus values can be plotted into the *CIE 1931 Chromaticity Diagram*:

```python
import pylab

# Plotting the *CIE 1931 Chromaticity Diagram*.
# The argument *standalone=False* is passed so that the plot doesn't get displayed
# and can be used as a basis for other plots.
plot_chromaticity_diagram_CIE1931(standalone=False)

# Calculating the *xy* chromaticity coordinates.
# The output domain of *colour.spectral_to_XYZ* is [0, 100] and
# the input domain of *colour.XYZ_to_sRGB* is [0, 1].
# We need to take it in account and rescale the input *CIE XYZ* colourspace matrix.
x, y = colour.XYZ_to_xy(colour.sd_to_XYZ(sd, cmfs, illuminant) / 100)

# Plotting the *xy* chromaticity coordinates.
pylab.plot(x, y, 'o-', color='white')

# Annotating the plot.
pylab.annotate(sd.name,
               xy=(x, y),
               xytext=(-50, 30),
               textcoords='offset points',
               arrowprops=dict(arrowstyle='->', connectionstyle='arc3, rad=-0.2'))

# Displaying the plot.
render(standalone=True);
```

Retrieving the *CIE XYZ* tristimulus values of any wavelength from colour matching functions is done using the `colour.wavelength_to_XYZ` definition, if the value requested is not available, the colour matching functions will be interpolated following *CIE 167:2005* recommendation:

```python
colour.wavelength_to_XYZ(546.1, colour.STANDARD_OBSERVERS_CMFS['CIE 1931 2 Degree Standard Observer'])
```

    array([ 0.3755316 ,  0.98444552,  0.01220285])

## Bibliography

1. ^ CIE. (n.d.). 17-198 colour (perceived). Retrieved June 26, 2014, from http://eilv.cie.co.at/term/198
2. ^ CIE. (n.d.). 17-659 light. Retrieved June 26, 2014, from http://eilv.cie.co.at/term/659
3. ^ CIE. (n.d.). 17-1238 spectrum. Retrieved June 27, 2014, from http://eilv.cie.co.at/term/1238
4. ^ CIE TC 1-38. (2005). 9. INTERPOLATION. In CIE 167:2005 Recommended Practice for Tabulating Spectral Data for Use in Colour Computations (pp. 14–19). ISBN:978-3-901-90641-1
5. ^ CIE TC 1-48. (2004). CIE 015:2004 Colorimetry, 3rd Edition. CIE 015:2004 Colorimetry, 3rd Edition (pp. 1–82). ISBN:978-3-901-90633-6
6. ^ Wyszecki, G., & Stiles, W. S. (2000). Integration Replace by Summation. In *Color Science: Concepts and Methods, Quantitative Data and Formulae* (pp. 158–163). Wiley. ISBN:978-0471399186

max_stars_repo_name: colour-science/colour-notebooks

max_issues_repo_name: colour-science/colour-ipython

max_forks_repo_name: colour-science/colour-ipython


text:
# 1. Algoritmos de Ordenação

Nesta aula, estudaremos diferentes algoritmos de ordenação, os quais podem ser utilizados em diversas situações. Quando temos um conjunto de dados quaisquer, ou seja, não temos conhecimento prévio dos mesmos, podemos fazer uso de um conjunto de técnicas de ordenação cuja complexidade inferior está limitada a $\theta(n\log n)$. As seguintes técnicas serão abordadas:
* Insertion Sort
* Mergesort
* Quicksort
* Heapsort

No entanto, quando algum conhecimento dos dados a serem ordenados é dado, podemos fazer uso de algumas técnicas de ordenação linear, tais como:
* Counting Sort
* Radix Sort
* Bucket Sort
As próximas seções tratam de abordas as técnicas mencionadas.

## 1.1. Insertion Sort

O algoritmo de ordenação por Insertion sort é indicado em situações nas quais os dados chegam ao **longo do tempo**, dado que a sua premissa consiste em inserir cada elemento em sua posição correta no vetor de dados. A Figura 1 apresenta um exemplo de funcionamento do algoritmo Insertion sort, em que o elemento que está sendo comparado na iteração atual (variável $v$ no algoritmo) encontra-se destacado em amarelo.

Figura 1: Funcionamento da técnica Insertion Sort.

O algoritmo abaixo implementa a técnica Insertion sort, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(n)$
* Complexidade de tempo: $O(n^2)$
* Estável: Sim
* In-place: Sim

```python
import numpy
import math
from collections import deque
```

```python
def InsertionSort(A):
    n = len(A)
    for i in range(2,n):
        v = A[i]
        j = i-1

        while (A[j] > v) and (j >= 0):
            A[j+1] = A[j]
            j = j-1
        A[j+1] = v

```

Note que o laço mais **externo** executa $n-1$ vezes, ou seja, $\theta(n)$ vezes. A complexidade do algoritmo de Insertion sort é dominada pelo laço mais **interno**. Caso tenhamos um vetor ordenado de forma **crescente**, por exemplo, a condição `A[j+1] > v` nunca será satisfeita e, por conseguinte, o laço mais interno não será executado. Nesta situação hipotética, a complexidade seria então dada por $\theta(n)$, que seria então a situação de **melhor caso** para o Insertion sort.

Entretanto, no **pior caso**, ou seja, um vetor ordenado de forma **decrescente**, o laço mais interno será executado $i-1$ vezes para cada valor de $i$. Desta forma, teremos $i(i-1)\approx i^2$ execuções, resultando em uma complexidade final de $O(n^2)$.

Exemplo de funcionamento:

```python
A = numpy.random.randint(-20, 20, 10)
print('Vetor de entrada: '+ str(A))
InsertionSort(A)
print('Vetor ordenado: '+ str(A))
```

    Vetor de entrada: [-10   5 -12   6  10  -4  -6  -2  -8 -11]
    Vetor ordenado: [-12 -11 -10  -8  -6  -4  -2   5   6  10]

## 1.2. Mergesort

O algoritmo de ordenação por Mergesort é baseado na metodologia de **divisão e conquista**, em que o problema principal é dividido em subproblemas que são resolvidos e, posteriormente, agrupados com o intuito de gerar a solução final. A divisão é realizada em $\theta(1)$, ao passo que a conquista é realizada em $\theta(n)$. De maneira geral, a ténica realiza divisões no vetor pela metade, até que um único elemento é obtido. Em seguida, o processo de conquista é realizado, em que os elementos são agrupados de maneira ordenada. A Figura 2 apresenta esse processo, em que os elementos cinza denotam a primeira metade do vetor, e os elementos em azul correspondem ao final da etapa de divisão.

Figura 2: Funcionamento da técnica Mergesort.

O algoritmo abaixo implementa a técnica Mergesort, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(2n)$
* Complexidade de tempo: $\theta(n\log n)$
* Estável: Sim
* In-place: Não

A complexidade de tempo pode ser calculada por meio do **Teorema Mestre** da seguinte maneira:

\begin{equation}
\tag{1}
\label{e.recurrence_analysis_mergesort_part1}
T(n) =
    \begin{cases}
      1 & \text{se $n=1$}\\
      aT\left(\frac{n}{b}\right) + f(n)& \text{caso contrário.}
    \end{cases}
\end{equation}

Como pode ser observado, a função `MergeSort` é invocada duas vezes de maneira recursiva, em que cada chamada é realizada em metade do vetor. Desta forma, temos que $a=b=2$.

Entretanto, falta ainda analisarmos o termo **residual** $f(n)$ da Equação \ref{e.recurrence_analysis_mergesort_part1}, o qual é composto por três laços: (i) o primeiro deles copia a primeira metade do vetor de dados de entrada $A$ para o vetor **auxiliar** $B$ e possui complexidade de tempo $\theta\left(\frac{n}{n}\right)\in\theta(n)$; (ii) o segundo laço copia a segunda metade do vetor $A$ em $B$, porém de maneira invertida, e também possui complexidade de tempo $\theta\left(\frac{n}{n}\right)\in\theta(n)$; e (iii) o último laço copia os elementos de $B$ no vetor de dados de entrada $A$ de maneira ordenada, ou seja, nesta etapa é realizada a junção (*merge*) dos dados, e possui complexidade $\theta(n)$. Desta forma, o termo residual, ou seja, a conquista, possui complexidade $\theta(n)$.

Assim sendo, podemos reescrever a Equação \ref{e.recurrence_analysis_mergesort_part1} da seguinte maneira:

\begin{equation}
\tag{2}
\label{e.recurrence_analysis_mergesort_part2}
T(n) =
    \begin{cases}
      1 & \text{se $n=1$}\\
      2T\left(\frac{n}{2}\right) + n& \text{caso contrário.}
    \end{cases}
\end{equation}
O termo residual pode ser mapeado para a forma $cn^k$ assumindo $c=k=1$. Desta forma, de acordo com o Teorema Mestre, a complexidade final de tempo do Mergesort é dada por $\theta(n^k\log_bn)=\theta(n\log n)$.

```python
def MergeSort(A, e, d):
    if (d > e):

        # Divisão é realizada em theta(1) *****
        m = math.floor((e+d)/2)
        MergeSort(A, e, m)
        MergeSort(A, m+1, d)
        # Fim da divisão **********************

        # Conquista é realizada em theta(n) ***
        B = numpy.zeros(d+1);

        for k in range(0, m+1):
            B[k] = A[k]

        j = d
        for k in range(m+1, d+1):
            B[j] = A[k]
            j = j-1

        i = e
        j = d
        for k in range(e, d+1):
            if(B[i] < B[j]):
                A[k] = B[i]
                i = i+1
            else:
                A[k] = B[j]
                j = j-1
        # Fim da conquista ********************
```

Exemplo de funcionamento:

```python
A = numpy.random.randint(-20, 20, 10)
print('Vetor de entrada: '+ str(A))
MergeSort(A, 0, 9)
print('Vetor ordenado: '+ str(A))
```

    Vetor de entrada: [ 17  14  -3 -12 -10  15  11 -11   0  -5]
    Vetor ordenado: [-12 -11 -10  -5  -3   0  11  14  15  17]

 ## 1.3. Quicksort

O algoritmo de ordenação Quicksort também é baseado na metodologia de divisão e conquista, em que esta última é realizada em tempo $\theta(1)$, ao passo que o "gargalo" da técnica está justamente na etapa de divisão. O algoritmo faz uso de um elemento denominado **pivô**, e assume que os elementos à **esquerda** deste são menores ou iguais à ele; ao passo que os elementos à sua **direita** são maiores ou iguais ao pivô. Note que, não necessariamente, os elementos que estão à esquerda e à direita estarão ordenados entre si.

Desta forma, a complexidade de tempo do Quicksort depende, necessariamente, da **posição** do pivô durate a execução do algoritmo. Se conseguíssemos garantir que o pivô sempre dividisse o vetor de entrada ao meio, então a a complexidade de ordenação seria igual à do Mergesort, ou seja, $\theta(n\log n)$. Entretanto, essa garantia não é possível de ser assegurada. A Figura 3 apresenta o funcionamento da técnica Quicksort, em que os elementos à esquerda do pivô são representador pela cor cinza, os elementos azuis correspondem aos valores que são maiores ou iguais ao pivô e, finalmente, o pivô está representado em laranja.

Figura 3: Funcionamento da técnica Quicksort.

O algoritmo abaixo implementa a técnica Quicksort, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(n)$
* Complexidade de tempo: $O(n^2)$, porém sua complexidade no caso médio é de $O(n\log n)$
* Estável: Não
* In-place: Sim

```python
def Quicksort(A, e, d):
    if (d > e):
            index = Partition(A, e, d) #função que retorna o novo índice do pivô (elemento laranja)
            Quicksort(A, e, index-1) # executa o Quicksort na parte cinza
            Quicksort(A, index+1, d) # executa o Quicksort na parte laranja

def Partition(A, e, d):
    pivot = A[d]
    p_index = e

    for i in range(e,d):
        if (A[i] <= pivot):
            A[i], A[p_index] = A[p_index], A[i]
            p_index = p_index + 1

    A[p_index], A[d] = A[d], A[p_index]

    return p_index
```

Exemplo de funcionamento:

```python
A = numpy.random.randint(-20, 20, 6)
print('Vetor de entrada: '+ str(A))
Quicksort(A, 0, 5)
print('Vetor ordenado: '+ str(A))
```

    Vetor de entrada: [-18 -18   4  13 -18  -3]
    Vetor ordenado: [-18 -18 -18  -3   4  13]

## 1.4. Heapsort

O algoritmo de ordenação Heapsort é uma técnica baseada em comparações, cujo funcionamento é dividido em três partes: (i) construção de um heap (usualmente binário) a partir dos dados de entrada, (ii) remover a raiz e inserir na última posição do vetor, e (iii) mover último nó para a posição da raiz. Note que este último passo exige que o heap seja atualizado a cada remoção de tal forma a manter suas propriedades. Para uma ordenação **crescente**, devemos utilizar um **heap máximo**, ao passo que para uma ordenação **decrescente** deve-se empregar um **heap mínimo**.

A primeira parte, como citado anteriormente, consiste em transformar um vetor de entrada em um **heap binário**, o qual, por sua vez, utiliza uma estrutura de **árvore binária**. Isto pode ser feito, basicamente, por meio de uma associação de índices do próprio vetor, como apresenta a Figura 4.

Figura 4: Mapeando um vetor para uma árvore binária.

Neste caso, temos que a raiz é representada pelo elemento $A[0]$, cujos filhos são os elementos $A[1]$ e $A[2]$. Por conseguinte, os filhos do nó $A[1]$ são dados por $A[3]$ e $A[4]$, e o último nó representado por $A[5]$. Assim sendo, temos que:

* A[0]: possui como filhos $A[1]$ e $A[2]$, ou seja, $A[i+1]$ e $A[i+2]$ para $i=0$.
* A[1]: possui como filhos $A[3]$ e $A[4]$, ou seja, $A[i+2]$ e $A[i+3]$ para $i=1$.
* A[2]: possui como filho $A[5]$, ou seja, $A[i+3]$ para $i=2$.

Podemos, então, generalizar da seguinte forma: $A[i]$ possui como filhos os elementos $A[2i+1]$ e $A[2i+2]$.

Em seguida, precisamos verificar se o heap obedece à propriedade de um heap máximo, ou seja, o valor armazenado em um nó pai tem que ser maior ou igual ao(s) valor(es) armazenzado(s) em seu(s) nó(s) filho(s). A Figura 5 apresenta essa metodologia.

Figura 5: Mapeando uma árvore árvore binária para um heap máximo.

Em seguida, o nó raiz (maior valor) é removido e inserido na última posição do vetor, e o último nó do heap toma o seu lugar. Entretanto, após esta operação, precisamos verificar se o heap obedece à propriedade de heap máximo, ou seja, um nó filho tem sempre um valor menor ou igual ao seu pai (procedimento implementado pela função `Heapify` no algoritmo abaixo). Após atualizarmos a estrutura, o procedimento todo é repetido novamente até que todos os elementos tenham sido removidos do heap. A Figura 6 ilustra este procedimento.

Figura 6: Funcionamento da técnica Heapsort.

O algoritmo abaixo implementa a técnica Heapsort, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(n)$
* Complexidade de tempo: $O(n\log n)$
* Estável: Não
* In-place: Sim

```python
#Algoritmo baseado no código disponível em https://www.tutorialspoint.com/python-program-for-heap-sort

# função que verifica se um heap obedece às propriedades de heap máximo
def Heapify(A, n, i):

    maior = i # maior valor
    l = 2*i + 1 # filho à esquerda
    r = 2*i + 2 # filho à esquerda

    # verifica se filho à esquerda existe e quem possui o maior valor
    if (l < n) and (A[i] < A[l]):
        maior = l

    # verifica se filho à direita existe e quem possui o maior valor
    if (r < n) and (A[maior] < A[r]):
        maior = r

    # analisa raiz
    if (maior != i):
        A[i],A[maior] = A[maior],A[i] # faz a troca entre a raiz e o maior valor encontrado

        # executa todo o processo novamente
        Heapify(A, n, maior)

# Heapsort
def Heapsort(A):

    n = len(A)

    # cria o heap máximo
    for i in range(n-1, -1, -1):
        Heapify(A, n, i)

    # remoção da raiz e sua inserção na posição correta do vetor ordenado
    for i in range(n-1, 0, -1):

        A[i], A[0] = A[0], A[i]
        Heapify(A, i, 0)
```

Note que o primeiro laço da função `Heapsort`executa $n+1$ vezes, ou seja, $\theta(n)$ vezes, ao passo que o segundo laço executa, também, $\theta(n)$ vezes. Agora, precisamos calcular a complexidade da função `Heapify`, a qual objetiva manter a estrutura de acordo com as propriedades de um heap máximo, conforme mencionado anteriormente. Como a criação do heap objetiva a construção de uma árvore binária balanceada, sabemos que sua altura pode ser dada por $O(\log n)$. No pior caso, ou seja, quando um nó folha precisa ser movido até a raiz, temos que um número $O(\log n)$ de comparações deve ser efetuado até que a propriedade do heap máximo seja satisfeita. Desta forma, temos que a função `Heapify` possui complexidade $O(\log n)$, resultando em uma complexidade final de $O(n\log n)$ para o algoritmo do Heapsort.

Exemplo de funcionamento:

```python
A = numpy.random.randint(-20, 20, 6)
print('Vetor de entrada: '+ str(A))
Heapsort(A)
print('Vetor ordenado: '+ str(A))
```

    Vetor de entrada: [ -8  17   1 -20   9   5]
    Vetor ordenado: [-20  -8   1   5   9  17]

## 1.4. Counting Sort

Suponha que tenhamos um vetor $A$ com $n$ elementos naturais, isto é, $A[i]\geq 0$, $\forall i=1,2,\dots,n$, tal que seu maior elemento possui o valor $k$. Neste caso, podemos fazer uso da técnica de ordenação por Counting Sort a qual, basicamente, conta quantos elementos são menores do que $A[i]$. O mais interessante da técnica é que ela faz a ordenação de um conjunto de valores **sem realizar comparações**. A Figura 7 ilustra este procedimento.

Figura 7: Funcionamento da técnica Counting Sort.

O algoritmo abaixo implementa a técnica Counting Sort, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(2n+k)$
* Complexidade de tempo: $\theta(n)$
* Estável: Sim
* In-place: Sim

```python
def CountingSort(A, k):

    n = len(A)
    B = numpy.zeros(n+1, dtype=int)
    C = numpy.zeros(k+1, dtype=int)

    # C[A[i]] armazena o número de vezes que o elemento A[i] aparece no vetor a ser ordenado
    for i in range(0, n):
        C[A[i]] = C[A[i]] + 1

    # C[i] armazena
    for i in range(1, k+1):
        C[i] = C[i] + C[i-1]

    for i in range(n-1, -1, -1):
        B[C[A[i]]] = A[i]
        C[A[i]] = C[A[i]] - 1

    return B
```

O primeiro laço do algoritmo é executado $\theta(n)$ vezes, ao passo que o segundo laço é executado por $\theta(k)$ iterações. Finalmente, o terceiro laço também é executado $\theta(n)$. Desta forma, temos uma complexidade dada por

Exemplo de funcionamento:

```python
k = 10
A = numpy.random.randint(0, k, 20)
print('Vetor de entrada: '+ str(A))
B = CountingSort(A, k)
print('Vetor ordenado: '+ str(B))
```

    Vetor de entrada: [4 5 3 9 0 3 0 2 3 4 9 8 5 8 9 9 5 3 6 8]
    Vetor ordenado: [0 0 0 2 3 3 3 3 4 4 5 5 5 6 8 8 8 9 9 9 9]

## 1.5. Bucket Sort

O algoritmo de ordenação Bucket Sort supõe que os elementos do vetor a ser ordenado $A$ encontram-se distribuídos **uniformemente** no intervalo $[0,1)$. A ideia consiste em dividir os elementos em $k$ segmentos (buckets) de igual tamanho e distribuir os elementos nesses segmentos. Como partimos da premissa de que os elementos do vetor de entrada encontram-se distribuídos uniformemente, é esperado que os segmentos contenham o mesmo número de elementos. Posteriormente, os elementos de cada segmento são ordenados por um método qualquer e o vetor final ordenado é obtido **concatenando** os segmentos em **ordem crescente**.

Diversos algoritmos podem ser utilizados para a ordenação dos elementos de cada segmento, muito embora uma das abordagens mais empregadas é o **Insertion Sort**, dado que o seu tempo de ordenação é baseado no quão distante um elemento está de sua posição final no conjunto ordenado. Ao dividir o vetor em segmentos, na prática, estamos facilitando a tarefa do algoritmo de ordenação, pois teremos menos elementos a serem ordenados quando consideramos um segmento por vez. A Figura 8 ilustra o funcionamento do Bucket Sort.

Figura 8: Funcionamento da técnica Bucket Sort.

O algoritmo abaixo implementa a técnica Bucket Sort utilizando Insertion Sort para o procedimento de ordenação de cada segmento, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(n+k)$
* Complexidade de tempo:
* * Melhor caso: $O(n)$
* * Caso médio: $O(n)$ quando $k\approx n$
* * Pior caso: $O(n^2)$
* Estável: Sim
* In-place: Não

```python
def BucketSort(A):

    n = len(A)

    # criando um vetor de listas ligadas
    B = numpy.empty(n)
    B = B.tolist()

    for i in range(0, n):
        B[i] = deque()

    for i in range(0, n):
        B[math.floor(n*A[i])].append(A[i]) # insere elemento A[i] em seu segmento (bucket) específico

    z = 0
    for i in range(0, n):
        InsertionSort(B[i]) # ordena segmento i

        for j in range(0, len(B[i])): # gerando vetor ordenado final pela concatenação dos segmentos
            A[z] = B[i].popleft()
            z = z+1

```

O primeiro laço do algoritmo é executado $\theta(n)$ vezes, e consiste na criação dos apontadores para as listas ligadas. Já o segundo laço, o qual também é executado $\theta(n)$ vezes, insere cada elemento do vetor de entrada $A$ em sua posição (segmento ou bucket) correspondente, a qual é calculada como segue:

\begin{equation}
\label{e.bucketsort}
\tag{3}
B[\lfloor n* A[i]\rfloor] = A[i].
\end{equation}

Note que estamos seguindo a implementaçNo último laço, o qual executa $\theta(n)$ vezes, cada segmento é ordenado utilizando Insertion Sort para posterior geração do vetor final ordenado bastando, apenas, concatenar cada segmento ordenado. Falta, apenas, determinarmos a complexidade do Insertion Sort. No **melhor caso**, teríamos apenas um único elemento em cada segmento, bastando copiar (concatenar) os elementos do vetor auxiliar $B$ no vetor de saída $A$, o qual resultaria em $\theta(n)$ operações. Já o **pior caso** contemplaria a situação em que todos os elementos seriam mapeados para um mesmo segmento. Nesta situação, sabemos que Insertion Sort possui complexidade $O(n^2)$.

Exemplo de funcionamento:

```python
A = numpy.random.random(10)
print('Vetor de entrada: '+ str(A))
BucketSort(A)
print('Vetor ordenado: '+ str(A))
```

    Vetor de entrada: [0.89806448 0.81719929 0.78607081 0.30571191 0.91212116 0.30742524
     0.64817982 0.58592705 0.38073669 0.62501044]
    Vetor ordenado: [0.30571191 0.30742524 0.38073669 0.58592705 0.64817982 0.62501044
     0.78607081 0.89806448 0.81719929 0.91212116]

## 1.5. Radix Sort

Suponha que o conjunto de dados $A$ a ser ordenado é composto por números **inteiros** com até $d$ dígitos. A ideia, de maneira geral, é bastante simples: considerando cada um dos dígitos, da esquerda para a direita, ordene o conjunto de elementos utilizando alguma técnica de ordenação **estável**. Assim sendo, a complexidade final do algoritmo Radix Sort será dependente da técnica utilizada para a ordenação dos números de acordo com cada dígito. A estabilidade garante que os números não trocarão de posição após serem ordenados para cada dígito. A Figura 9 ilustra o funcionamento do Radix Sort.

Figura 9: Funcionamento da técnica Radix Sort.

O algoritmo abaixo implementa a técnica Radix Sort, cuja análise de complexidade é dada da seguinte maneira:
* Complexidade de espaço: $\theta(g(n))$, em que $g(n)$ denota a complexidade de espaço do algoritmo utilizado para a ordenação
* Complexidade de tempo: $O(df(n))$, em que $f(n)$ denota a complexidade de tempo do algoritmo utilizado para a ordenação
* Estável: Sim
* In-place: Sim/Não (depende do algoritmo utilizado para a ordenação)

No exemplo abaixo, utilizaremos uma heurística similar à técnica Counting Sort, em que os números são posicionados em seus respectivos segmentos (buckets) de acordo com o dígito em questão.

```python
def RadixSort(A):

    radix = 10
    aux = 1
    fim = False

    while not fim:

        fim = True
        B = [list() for _ in range(radix)] # cria os buckets
        for i in A:
            tmp = i/aux
            B[int(tmp % radix)].append(i)
            if(fim) and (tmp > 0):
                fim = False

        j = 0
        for b in range(radix):
            buck = B[b]

            for i in buck:
                A[j] = i
                j = j + 1

        aux = aux * radix
```

O laço principal (comando `while`) executa $\theta(d)$ vezes, em que $d$ corresponde ao número de dígitos do maior elemento. Já o primeiro laço interno executa $\theta(n)$ vezes, pois percorre todos os elementos do vetor $A$, e o segundo laço interno percorre cada um dos $10$ segmentos. Como estamos utilizando um conjunto de valores a serem ordenados que pertencem à um sistema decimal, ou seja, `radix = 10`, temos que esse laço executa $10$ vezes,ou seja, em tempo constante $\theta(1)$. Note que a palavra **radix** tem origem no latim e significa **raiz**, ou base. Assim, a técnica Radix Sort ordena os números de acordo com sua base decimal.

Desta forma, a **complexidade de tempo** do Radix Sort é dada por $\theta(dn)$, ao passo que sua **complexidade de espaço** é dada por $\theta(2n)$, dado que temos o vetor de antrada $A$ com $n$ elementos e uma estrutura de lista ligada temporária (`buckets`) que armazena os elementos. Segue, abaixo, um exemplo de funcionamento:

```python
d = 3 # número de dígitos do maior elemento a ser gerado aleatoriamente
A = numpy.random.randint(0, math.pow(10,d)-1, 20)
print('Vetor de entrada: '+ str(A))
RadixSort(A)
print('Vetor ordenado: '+ str(A))
```

    Vetor de entrada: [367 622  62 777 214 360 580 863 153 653 818 701 826 962 700 939 336 376
     667 996]
    Vetor ordenado: [ 62 153 214 336 360 367 376 580 622 653 667 700 701 777 818 826 863 939
     962 996]

```python

```

max_stars_repo_name: jppbsi/lectures

max_issues_repo_name: jppbsi/lectures

max_forks_repo_name: jppbsi/lectures