///|
pub fn Grid::window_summary(
  self : Grid,
  row : Int,
  col : Int,
  radius : Int,
) -> Summary {
  summary(self.window(row, col, Int::max(0, radius)))
}

///|
pub fn Grid::local_variance(self : Grid, radius : Int) -> Grid {
  let safe = Int::max(0, radius)
  let output : Array[Double] = []
  for row in 0.. Grid {
  self.local_variance(radius).map(fn(value) { value.sqrt() })
}

///|
pub fn Grid::local_z_score(self : Grid, radius : Int) -> Grid {
  let safe = Int::max(0, radius)
  let output : Array[Double] = []
  for row in 0.. Grid {
  let safe = Int::max(0, radius)
  let output : Array[Double] = []
  for row in 0.. Grid {
  let safe = Int::max(0, radius)
  let output : Array[Double] = []
  for row in 0.. Grid {
  let safe = Int::max(1, radius)
  let output = self.values.copy()
  for row in 0.. Grid {
  let output = self.values.copy()
  let points = self.cells()
  for i in 0.. Grid {
  self.map(fn(value) {
    if !is_finite(value) {
      value
    } else {
      value.clamp(min=range.minimum, max=range.maximum)
    }
  })
}

///|
pub fn Grid::rescale_values(self : Grid, target : NumericRange) -> Grid {
  let source = self.value_range()
  let width = source.maximum - source.minimum
  let target_width = target.maximum - target.minimum
  if width == 0.0 {
    self.map(fn(value) { if is_finite(value) { target.minimum } else { value } })
  } else {
    self.map(fn(value) {
      if is_finite(value) {
        target.minimum + (value - source.minimum) / width * target_width
      } else {
        value
      }
    })
  }
}

///|
pub fn Grid::rank_transform(self : Grid) -> Grid {
  let ranks = rank_values(self.values)
  { ..self, values: ranks }
}

///|
pub fn grid_local_extrema(grid : Grid, find_maximum : Bool) -> Array[GridCell] {
  let output : Array[GridCell] = []
  for row in 0.. {
          let neighbours = grid.window(row, col, 1)
          let mut extremum = true
          for neighbour in neighbours {
            if is_finite(neighbour) && is_finite(cell.value) {
              if find_maximum && neighbour > cell.value {
                extremum = false
              }
              if !find_maximum && neighbour < cell.value {
                extremum = false
              }
            }
          }
          if extremum && is_finite(cell.value) {
            output.push(cell)
          }
        }
        None => ()
      }
    }
  }
  output
}

///|
pub fn grid_total_variation(grid : Grid) -> Double {
  let mut total = 0.0
  for row in 0.. {
          if col + 1 < grid.cols {
            match grid.value_at(row, col + 1) {
              Some(right) =>
                if is_finite(value) && is_finite(right) {
                  total += (right - value).abs()
                }
              None => ()
            }
          }
          if row + 1 < grid.rows {
            match grid.value_at(row + 1, col) {
              Some(down) =>
                if is_finite(value) && is_finite(down) {
                  total += (down - value).abs()
                }
              None => ()
            }
          }
        }
        None => ()
      }
    }
  }
  total
}

///|
pub fn grid_surface_area_m2(grid : Grid) -> Double {
  let cell_width = local_degree_longitude_m(grid.origin_lat).abs() *
    grid.lon_step.abs()
  let cell_height = local_degree_latitude_m(grid.origin_lat).abs() *
    grid.lat_step.abs()
  let mut area = 0.0
  for row in 0..
          if is_finite(gradient.magnitude) {
            area += cell_width *
              cell_height *
              (1.0 + gradient.magnitude * gradient.magnitude).sqrt()
          }
        None => ()
      }
    }
  }
  area
}

///|
pub fn grid_ruggedness(grid : Grid, radius : Int) -> Grid {
  grid.local_standard_deviation(radius)
}

///|
pub fn grid_entropy(grid : Grid, bins : Int) -> Double {
  let histogram_bins = histogram(grid.values, bins)
  let total = count_finite(grid.values)
  if total == 0 {
    0.0
  } else {
    let mut entropy = 0.0
    for bin in histogram_bins {
      if bin.count > 0 {
        let probability = Double::from_int(bin.count) / Double::from_int(total)
        entropy -= probability * @math.ln(probability)
      }
    }
    entropy
  }
}

///|
pub fn grid_histogram(grid : Grid, bins : Int) -> Array[HistogramBin] {
  histogram(grid.values, bins)
}

///|
pub fn grid_quantile_mask(grid : Grid, fraction : Double) -> Grid {
  let threshold = percentile(grid.values, fraction)
  grid.map(fn(value) {
    if is_finite(value) && value >= threshold {
      1.0
    } else {
      0.0
    }
  })
}

///|
pub fn grid_classify_quantiles(grid : Grid, classes : Int) -> Grid {
  let safe = Int::max(1, classes)
  let thresholds : Array[Double] = []
  for i in 1.. Grid {
  let output : Array[Double] = []
  for row in 0..
            if is_finite(value) && is_finite(other) {
              difference += (value - other).abs()
              count += 1
            }
          None => ()
        }
      }
      output.push(
        if count == 0 {
          0.0
        } else {
          difference / Double::from_int(count)
        },
      )
    }
  }
  { ..grid, values: output }
}

///|
pub fn grid_normalized_gradient(grid : Grid) -> Grid {
  let gradient = grid.gradient()
  gradient.rescale_values({ minimum: 0.0, maximum: 1.0 })
}

///|
pub fn grid_slope_radians(grid : Grid) -> Grid {
  grid.map_cells(fn(cell) {
    match grid.gradient_at(cell.row, cell.col) {
      Some(gradient) => @math.atan(gradient.magnitude)
      None => 0.0 / 0.0
    }
  })
}

///|
pub fn grid_aspect_radians(grid : Grid) -> Grid {
  grid.map_cells(fn(cell) {
    match grid.gradient_at(cell.row, cell.col) {
      Some(gradient) => deg_to_rad(gradient.azimuth_deg)
      None => 0.0 / 0.0
    }
  })
}

///|
pub fn grid_apply_polynomial(
  grid : Grid,
  coefficients : ArrayView[Double],
) -> Grid {
  grid.map(fn(value) { polynomial_evaluate(coefficients, value) })
}

///|
pub fn grid_difference_of_gaussians(
  grid : Grid,
  small_radius : Int,
  large_radius : Int,
) -> Grid {
  let small = grid.gaussian_blur(small_radius)
  let large = grid.gaussian_blur(large_radius)
  small.difference(large).unwrap_or(grid)
}

///|
pub fn grid_percentile_clip(
  grid : Grid,
  lower_fraction : Double,
  upper_fraction : Double,
) -> Grid {
  grid.clamp_values({
    minimum: percentile(grid.values, lower_fraction),
    maximum: percentile(grid.values, upper_fraction),
  })
}

///|
pub fn grid_valid_cells(grid : Grid) -> Array[GridCell] {
  grid.cells().filter(fn(cell) { is_finite(cell.value) })
}

///|
pub fn grid_centroid(grid : Grid) -> GeoPoint? {
  let cells = grid_valid_cells(grid)
  if cells.length() == 0 {
    None
  } else {
    let mut lat = 0.0
    let mut lon = 0.0
    let mut weight = 0.0
    for cell in cells {
      let cell_weight = cell.value.abs()
      lat += cell.point.lat * cell_weight
      lon += cell.point.lon * cell_weight
      weight += cell_weight
    }
    if weight == 0.0 {
      None
    } else {
      Some({ lat: lat / weight, lon: lon / weight, elevation_m: 0.0 })
    }
  }
}

///|
pub fn grid_profile_sum(grid : Grid, row : Int) -> Double {
  sum_finite(grid.row(row))
}

///|
pub fn grid_column_sum(grid : Grid, col : Int) -> Double {
  sum_finite(grid.column(col))
}

///|
pub fn grid_row_means(grid : Grid) -> Array[Double] {
  let output : Array[Double] = []
  for row in 0.. Array[Double] {
  let output : Array[Double] = []
  for col in 0.. Summary {
  match observed.difference(predicted) {
    Some(residual) => residual.summary()
    None => summary([])
  }
}

///|
pub fn grid_signal_to_noise(grid : Grid, noise : Grid) -> Grid? {
  grid.combine(noise, fn(signal, uncertainty) {
    if is_finite(signal) && is_finite(uncertainty) && uncertainty != 0.0 {
      signal / uncertainty
    } else {
      0.0 / 0.0
    }
  })
}