maths.util.ts

import * as d3 from 'd3-format';

import {XySequence} from '../data';
import {Imt} from '../gmm';
import {imtToPeriod} from '../gmm/imt.model';

/**
 * Normal complementary cumulative distribution function.
 *
 * @param μ mean
 * @param σ standard deviation
 * @param x variate
 */
const normalCcdf = (μ: number, σ: number, x: number) => {
  return (1.0 + erf((μ - x) / (σ * Math.sqrt(2)))) * 0.5;
};

/**
 * Error function approximation of Abramowitz and Stegun, formula 7.1.26 in
 * the <em>Handbook of Mathematical Functions with Formulas, Graphs, and
 * Mathematical Tables</em>. Although the approximation is only valid for
 * x ≥ 0, because erf(x) is an odd function,
 * erf(x) = −erf(−x) and negative values are supported.
 */
const erf = (x: number) => {
  return x < 0.0 ? -erfBase(-x) : erfBase(x);
};

/**
 * Interpolate.
 *
 * @param xySequence The Xy sequence
 * @param value The value to interpolate at.
 */
const interpolate = (xySequence: XySequence, value: number): number => {
  const xValues = xySequence.xs;
  const yValues = xySequence.ys;

  const index = closestIndex(yValues, value);

  const x0 = xValues[index - 1];
  const x1 = xValues[index];
  const y0 = yValues[index - 1];
  const y1 = yValues[index];
  const x = interpolateValue(x0, x1, y0, y1, value);
  return isNaN(x) ? x : round(x, 6);
};

/**
 * Round a number to specific format
 *
 * @param value Value to round
 * @param scale Format scale
 */
const round = (value: number, scale: number) => {
  const format = d3.format(`.${scale}f`);
  return Number(format(value));
};

/**
 * Calculate the response sepectrum for each IMT.
 *
 * @param hazards The hazard curves by imt
 * @param returnPeriod The return period (in years) to calculate at
 */
const responseSpectrum = (
  hazards: Map<Imt, XySequence>,
  returnPeriod: number
): XySequence => {
  const xs: number[] = [];
  const ys: number[] = [];

  for (const [imt, xySequence] of hazards) {
    xs.push(imtToPeriod(imt));
    ys.push(calculateResponseSpectrum(xySequence, returnPeriod));
  }

  return {
    xs,
    ys,
  };
};

const calculateResponseSpectrum = (
  xySequence: XySequence,
  returnPeriod: number
): number => {
  return interpolate(xySequence, 1 / returnPeriod);
};

/**
 * Return index of closest value.
 *
 * @param values Array of values
 * @param value Value to find
 */
const closestIndex = (values: number[], value: number): number => {
  const closestValue = values.reduce((prev, curr) =>
    Math.abs(curr - value) < Math.abs(prev - value) ? curr : prev
  );

  return values.findIndex(value => value === closestValue);
};

/**
 * Interpolate.
 *
 * @param x0 X value before value
 * @param x1 X value after value
 * @param y0 Y value before value
 * @param y1 Y value after value
 * @param value Value to interpolate at
 */
const interpolateValue = (
  x0: number,
  x1: number,
  y0: number,
  y1: number,
  value: number
): number => {
  return x0 + (Math.log10(value / y0) * (x1 - x0)) / Math.log10(y1 / y0);
};

const erfBase = (x: number) => {
  const P = 0.3275911;
  const A1 = 0.254829592;
  const A2 = -0.284496736;
  const A3 = 1.421413741;
  const A4 = -1.453152027;
  const A5 = 1.061405429;

  const t = 1 / (1 + P * x);
  const tsq = t * t;

  return (
    1 -
    (A1 * t + A2 * tsq + A3 * tsq * t + A4 * tsq * tsq + A5 * tsq * tsq * t) *
      Math.exp(-x * x)
  );
};

/**
 * Export functions
 */
export const Maths = {
  normalCcdf,
  erf,
  interpolate,
  round,
  responseSpectrum,
};