defmodule Growth.Calc.ZScore do
  @moduledoc """
  Calculate z-score for a given measurement and the fitted values of Box-Cox by age/height:

  * **power** `t:l/0`
  * **median** `t:m/0`
  * **coefficientof variation** `t:s/0`

  This calculation is described in the [instructions provided by the World Health Organization](https://cdn.who.int/media/docs/default-source/child-growth/growth-reference-5-19-years/computation.pdf).

  ## Examples:

      iex> measures =
      ...>   [
      ...>     [30, -1.7862, 16.9392, 0.1107],
      ...>     [14, -1.3592, 20.4951, 0.12579],
      ...>     [19, -1.6318, 16.049, 0.10038]
      ...>   ]
      iex> Enum.map(measures, &apply(Growth.Calc.ZScore, :compute, &1))
      [3.35390255606726, -3.7985108865993493, 1.4698319520484722]

  """

  @typedoc """
  The **power** value in Box-Cox transformation.
  """
  @type l :: number()
  @typedoc """
  The **median** value in Box-Cox transformation.
  """
  @type m :: number()
  @typedoc """
  The **coefficient of variation** in Box-Cox transformation.
  """
  @type s :: number()

  @spec compute(number(), l(), m(), s()) :: number()
  @doc """
  Compute the adjusted z-score for a given measurement (`y`) and the Box-Cox values: power (`l`), median (`m`), and
  coefficient of variation (`s`).

  ## Examples:

      iex> Growth.Calc.ZScore.compute(30, -1.7862, 16.9392, 0.1107)
      3.35390255606726
      iex> Growth.Calc.ZScore.compute(14, -1.3592, 20.4951, 0.12579)
      -3.7985108865993493
      iex> Growth.Calc.ZScore.compute(19, -1.6318, 16.049, 0.10038)
      1.4698319520484722

  """
  def compute(y, l, m, s) do
    y
    |> raw(l, m, s)
    |> adjust(y, l, m, s)
  end

  @spec raw(number(), l(), m(), s()) :: number()
  @doc """
  Compute the raw z-score for a given measurement (`y`) and the Box-Cox values, power (`l`), median (`m`), and
  coefficient of variation (`s`).

  ## Examples

      iex> Growth.Calc.ZScore.raw(30, -1.7862, 16.9392, 0.1107)
      3.2353902101095855
      iex> Growth.Calc.ZScore.raw(14, -1.3592, 20.4951, 0.12579)
      -3.969714730727475
      iex> Growth.Calc.ZScore.compute(19, -1.6318, 16.049, 0.10038)
      1.4698319520484722

  """
  def raw(y, l, m, s) do
    (:math.pow(y / m, l) - 1) / (s * l)
  end

  @spec adjust(number(), number(), l(), m(), s()) :: number()
  @doc """
  Adjust the raw z-score considering that extreme values, beyond -3 and 3 standard deviation, must be handled differently.
  """
  def adjust(zscore, y, l, m, s) when zscore > 3 do
    [sd2, sd3, _, _] = cutoffs(l, m, s)
    sd_delta = sd3 - sd2
    3 + (y - sd3) / sd_delta
  end

  def adjust(zscore, y, l, m, s) when zscore < -3 do
    [_, _, sd2, sd3] = cutoffs(l, m, s)
    sd_delta = sd2 - sd3
    -3 + (y - sd3) / sd_delta
  end

  def adjust(zscore, _, _, _, _) do
    zscore
  end

  @spec cutoffs(number(), number(), number()) :: [number()]
  @doc """
  Calculate the standard deviations cutoffs (2, 3, -2, and -3) for a given set of fitted Box-Cox values: power (`l`),
  median (`m`), and coefficient of variation (`s`).

  These cutoffs are used to calculate the adjusted z-score.
  """
  def cutoffs(l, m, s) do
    Enum.map([2, 3, -2, -3], &measure_at_standard_deviation(&1, l, m, s))
  end

  @spec measure_at_standard_deviation(integer(), number(), number(), number()) :: number()
  @doc """
  Calculate the measure value at a given standard deviation (`sd`), considering the fitted Box-Cox values: power (`l`),
  median (`m`), and coefficient of variation (`s`).
  """
  def measure_at_standard_deviation(sd, l, m, s) do
    m * :math.pow(1 + l * s * sd, 1 / l)
  end
end