/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

#ifdef HAVE_CONFIG_H
#include "config.h"  // restrict
#endif

#include <float.h>
#include <fenv.h>
#include <assert.h>

#include "compare.h"
#include "array.h"
#include "cimdOmp.h"

//#pragma STDC FENV_ACCESS ON

const char *
fpe_errstr(int fpeRaised)
{
  const char *errstr = nullptr;

  // clang-format off
  if      (fpeRaised & FE_DIVBYZERO) errstr = "division by zero";
  else if (fpeRaised & FE_INEXACT)   errstr = "inexact result";
  else if (fpeRaised & FE_INVALID)   errstr = "invalid result";
  else if (fpeRaised & FE_OVERFLOW)  errstr = "overflow";
  else if (fpeRaised & FE_UNDERFLOW) errstr = "underflow";
  // clang-format on

  return errstr;
}

MinMaxVal
varrayMinMaxMV(const size_t len, const double *array, const double missval, size_t &nvals)
{
  double vmin = DBL_MAX;
  double vmax = -DBL_MAX;

  nvals = 0;
  for (size_t i = 0; i < len; ++i)
    {
      if (!DBL_IS_EQUAL(array[i], missval))
        {
          if (array[i] < vmin) vmin = array[i];
          if (array[i] > vmax) vmax = array[i];
          nvals++;
        }
    }

  return MinMaxVal(vmin, vmax);
}

MinMaxVal
varrayMinMaxMV(const size_t len, const double *array, const double missval)
{
  size_t nvals;
  return varrayMinMaxMV(len, array, missval, nvals);
}

MinMaxVal
varrayMinMaxMV(const size_t len, const Varray<double> &array, const double missval)
{
  size_t nvals;
  return varrayMinMaxMV(len, array.data(), missval, nvals);
}

template <typename T>
void
varrayMinMaxSum(const size_t len, const Varray<T> &array, double &rmin, double &rmax, double &rsum)
{
  // rmin, rmax and rsum will be initialized in Info

  /* 
  gnu 9.2.0: correct result
  intel 18.0.5: wrong result
#ifdef HAVE_OPENMP4
#pragma omp simd reduction(+:rsum) reduction(min:rmin) reduction(max:rmax)
#endif
  */
  for (size_t i = 0; i < len; ++i)
    {
      const auto v = array[i];
      if (v < rmin) rmin = v;
      if (v > rmax) rmax = v;
      rsum += v;
    }
}

// Explicit instantiation
template void varrayMinMaxSum(const size_t len, const Varray<float> &array, double &rmin, double &rmax, double &rsum);
template void varrayMinMaxSum(const size_t len, const Varray<double> &array, double &rmin, double &rmax, double &rsum);

template <typename T>
size_t
varrayMinMaxSumMV(const size_t len, const Varray<T> &array, const T missval, double &rmin, double &rmax, double &rsum)
{
  // rmin, rmax and rsum will be initialized in Info

  size_t nvals = 0;
  for (size_t i = 0; i < len; ++i)
    {
      const auto v = array[i];
      if (!DBL_IS_EQUAL(v, missval))
        {
          if (v < rmin) rmin = v;
          if (v > rmax) rmax = v;
          rsum += v;
          nvals++;
        }
    }

  if (nvals == 0 && IS_EQUAL(rmin, DBL_MAX)) rmin = missval;
  if (nvals == 0 && IS_EQUAL(rmax, -DBL_MAX)) rmax = missval;

  return nvals;
}

// Explicit instantiation
template size_t varrayMinMaxSumMV(const size_t len, const Varray<float> &array, const float missval, double &rmin, double &rmax, double &rsum);
template size_t varrayMinMaxSumMV(const size_t len, const Varray<double> &array, const double missval, double &rmin, double &rmax, double &rsum);

void
varrayMinMaxMean(const size_t len, const Varray<double> &array, double &rmin, double &rmax, double &rmean)
{
  rmin = DBL_MAX;
  rmax = -DBL_MAX;
  double sum = 0.0;
  varrayMinMaxSum(len, array, rmin, rmax, sum);

  rmean = len ? sum / (double) len : 0;
}

size_t
varrayMinMaxMeanMV(const size_t len, const Varray<double> &array, const double missval, double &rmin, double &rmax, double &rmean)
{
  rmin = DBL_MAX;
  rmax = -DBL_MAX;
  double sum = 0.0;
  size_t nvals = varrayMinMaxSumMV(len, array, missval, rmin, rmax, sum);

  rmean = nvals ? sum / (double) nvals : missval;

  return nvals;
}

void
arrayMinMaxMask(const size_t len, const double *array, const Varray<int> &mask, double &rmin, double &rmax)
{
  double zmin = DBL_MAX;
  double zmax = -DBL_MAX;

  if (!mask.empty())
    {
      for (size_t i = 0; i < len; ++i)
        {
          if (!mask[i])
            {
              if (array[i] < zmin) zmin = array[i];
              if (array[i] > zmax) zmax = array[i];
            }
        }
    }
  else
    {
      for (size_t i = 0; i < len; ++i)
        {
          if (array[i] < zmin) zmin = array[i];
          if (array[i] > zmax) zmax = array[i];
        }
    }

  rmin = zmin;
  rmax = zmax;
}

void
arrayAddArray(const size_t len, double *restrict array1, const double *restrict array2)
{
  //#ifdef  _OPENMP
  //#pragma omp parallel for default(none) shared(array1,array2)
  //#endif
  for (size_t i = 0; i < len; ++i) array1[i] += array2[i];
}

void
arrayAddArrayMV(const size_t len, double *restrict array1, const double *restrict array2, const double missval)
{
  if (std::isnan(missval))
    {
      for (size_t i = 0; i < len; i++)
        if (!DBL_IS_EQUAL(array2[i], missval))
          {
            if (!DBL_IS_EQUAL(array1[i], missval))
              array1[i] += array2[i];
            else
              array1[i] = array2[i];
          }
    }
  else
    {
      for (size_t i = 0; i < len; i++)
        if (IS_NOT_EQUAL(array2[i], missval))
          {
            if (IS_NOT_EQUAL(array1[i], missval))
              array1[i] += array2[i];
            else
              array1[i] = array2[i];
          }
    }
}

size_t
arrayNumMV(const size_t len, const double *restrict array, const double missval)
{
  size_t nmiss = 0;

  if (std::isnan(missval))
    {
      for (size_t i = 0; i < len; ++i)
        if (DBL_IS_EQUAL(array[i], missval)) nmiss++;
    }
  else
    {
      for (size_t i = 0; i < len; ++i)
        if (IS_EQUAL(array[i], missval)) nmiss++;
    }

  return nmiss;
}

template <typename T>
size_t
varrayNumMV(const size_t len, const Varray<T> &array, const T missval)
{
  size_t nmiss = 0;

  if (std::isnan(missval))
    {
      for (size_t i = 0; i < len; ++i)
        if (DBL_IS_EQUAL(array[i], missval)) nmiss++;
    }
  else
    {
      for (size_t i = 0; i < len; ++i)
        if (IS_EQUAL(array[i], missval)) nmiss++;
    }

  return nmiss;
}

// Explicit instantiation
template size_t varrayNumMV(const size_t len, const Varray<float> &array, const float missval);
template size_t varrayNumMV(const size_t len, const Varray<double> &array, const double missval);

MinMaxVal
varrayMinMax(const size_t len, const double *restrict array)
{
  double vmin = DBL_MAX;
  double vmax = -DBL_MAX;

#ifdef HAVE_OPENMP4
#pragma omp simd reduction(min:vmin) reduction(max:vmax)
#endif
  for (size_t i = 0; i < len; ++i)
    {
      if (array[i] < vmin) vmin = array[i];
      if (array[i] > vmax) vmax = array[i];
    }

  return MinMaxVal(vmin, vmax);
}

MinMaxVal
varrayMinMax(const size_t len, const Varray<double> &v)
{
  return varrayMinMax(len, v.data());
}

template <typename T>
MinMaxVal
varrayMinMax(const Varray<T> &v)
{
  double vmin = DBL_MAX;
  double vmax = -DBL_MAX;

  size_t len = v.size();
#ifdef HAVE_OPENMP4
#pragma omp simd reduction(min:vmin) reduction(max:vmax)
#endif
  for (size_t i = 0; i < len; ++i)
    {
      if (v[i] < vmin) vmin = v[i];
      if (v[i] > vmax) vmax = v[i];
    }

  return MinMaxVal(vmin, vmax);
}

// Explicit instantiation
template MinMaxVal varrayMinMax(const Varray<float> &v);
template MinMaxVal varrayMinMax(const Varray<double> &v);

double
varrayMin(const size_t len, const Varray<double> &v)
{
  assert(v.size() >= len);

  double min = v[0];

  for (size_t i = 0; i < len; ++i)
    if (v[i] < min) min = v[i];

  return min;
}

double
varrayMax(const size_t len, const Varray<double> &v)
{
  assert(v.size() >= len);

  double max = v[0];

  for (size_t i = 0; i < len; ++i)
    if (v[i] > max) max = v[i];

  return max;
}

double
varrayRange(const size_t len, const Varray<double> &v)
{
  assert(v.size() >= len);

  double min = v[0];
  double max = v[0];

  for (size_t i = 0; i < len; ++i)
    {
      if (v[i] < min) min = v[i];
      if (v[i] > max) max = v[i];
    }

  return (max - min);
}

double
varrayMinMV(const size_t len, const Varray<double> &v, const double missval)
{
  assert(v.size() >= len);

  double min = DBL_MAX;

  for (size_t i = 0; i < len; ++i)
    if (!DBL_IS_EQUAL(v[i], missval))
      if (v[i] < min) min = v[i];

  if (IS_EQUAL(min, DBL_MAX)) min = missval;

  return min;
}

double
varrayMaxMV(const size_t len, const Varray<double> &v, const double missval)
{
  assert(v.size() >= len);

  double max = -DBL_MAX;

  for (size_t i = 0; i < len; ++i)
    if (!DBL_IS_EQUAL(v[i], missval))
      if (v[i] > max) max = v[i];

  if (IS_EQUAL(max, -DBL_MAX)) max = missval;

  return max;
}

double
varrayRangeMV(const size_t len, const Varray<double> &v, const double missval)
{
  assert(v.size() >= len);

  double min = DBL_MAX;
  double max = -DBL_MAX;

  for (size_t i = 0; i < len; ++i)
    {
      if (!DBL_IS_EQUAL(v[i], missval))
        {
          if (v[i] < min) min = v[i];
          if (v[i] > max) max = v[i];
        }
    }

  double range;
  if (IS_EQUAL(min, DBL_MAX) && IS_EQUAL(max, -DBL_MAX))
    range = missval;
  else
    range = max - min;

  return range;
}

double
varraySum(const size_t len, const Varray<double> &v)
{
  assert(v.size() >= len);

  double sum = 0;

  for (size_t i = 0; i < len; ++i) sum += v[i];

  return sum;
}

double
varraySumMV(const size_t len, const Varray<double> &v, const double missval)
{
  assert(v.size() >= len);

  double sum = 0;
  size_t nvals = 0;

  if (std::isnan(missval))
    {
      for (size_t i = 0; i < len; ++i)
        if (!DBL_IS_EQUAL(v[i], missval))
          {
            sum += v[i];
            nvals++;
          }
    }
  else
    {
      for (size_t i = 0; i < len; ++i)
        if (IS_NOT_EQUAL(v[i], missval))
          {
            sum += v[i];
            nvals++;
          }
    }

  if (!nvals) sum = missval;

  return sum;
}

double
varrayMean(const size_t len, const Varray<double> &v)
{
  assert(v.size() >= len);

  const double sum = varraySum(len, v);

  return sum / len;
}

double
varrayMeanMV(const size_t len, const Varray<double> &v, const double missval)
{
  assert(v.size() >= len);

  double sum = 0, sumw = 0;

  for (size_t i = 0; i < len; ++i)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        sum += v[i];
        sumw += 1;
      }

  double missval1 = missval, missval2 = missval;
  return DIVMN(sum, sumw);
}

double
varrayWeightedMean(const size_t len, const Varray<double> &v, const Varray<double> &w, const double missval)
{
  assert(v.size() >= len);
  assert(w.size() >= len);

  double sum = 0, sumw = 0;

  for (size_t i = 0; i < len; ++i)
    {
      sum += w[i] * v[i];
      sumw += w[i];
    }

  return IS_EQUAL(sumw, 0.) ? missval : sum / sumw;
}

double
varrayWeightedMeanMV(const size_t len, const Varray<double> &v, const Varray<double> &w, const double missval)
{
  assert(v.size() >= len);
  assert(w.size() >= len);

  const double missval1 = missval, missval2 = missval;
  double sum = 0, sumw = 0;

  for (size_t i = 0; i < len; ++i)
    if (!DBL_IS_EQUAL(v[i], missval1) && !DBL_IS_EQUAL(w[i], missval1))
      {
        sum += w[i] * v[i];
        sumw += w[i];
      }

  return DIVMN(sum, sumw);
}

double
varrayAvgMV(const size_t len, const Varray<double> &v, const double missval)
{
  assert(v.size() >= len);

  const double missval1 = missval, missval2 = missval;
  double sum = 0, sumw = 0;

  for (size_t i = 0; i < len; ++i)
    {
      sum = ADDMN(sum, v[i]);
      sumw += 1;
    }

  return DIVMN(sum, sumw);
}

double
varrayWeightedAvgMV(const size_t len, const Varray<double> &v, const Varray<double> &w, const double missval)
{
  assert(v.size() >= len);
  assert(w.size() >= len);

  const double missval1 = missval, missval2 = missval;
  double sum = 0, sumw = 0;

  for (size_t i = 0; i < len; ++i)
    if (!DBL_IS_EQUAL(w[i], missval))
      {
        sum = ADDMN(sum, MULMN(w[i], v[i]));
        sumw = ADDMN(sumw, w[i]);
      }

  return DIVMN(sum, sumw);
}

static void
varrayPrevarsum0(size_t len, const Varray<double> &v, double &rsum, double &rsumw)
{
  rsum = 0;
  for (size_t i = 0; i < len; i++)
    {
      rsum += v[i];
    }
  rsumw = len;
}

static void
varrayPrevarsum0MV(size_t len, const Varray<double> &v, double missval, double &rsum, double &rsumw)
{
  rsum = rsumw = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        rsum += v[i];
        rsumw += 1;
      }
}

static void
varrayPrevarsum(size_t len, const Varray<double> &v,
                double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = 0;
  for (size_t i = 0; i < len; i++)
    {
      rsum += v[i];
      rsumq += v[i] * v[i];
    }
  rsumw = len;
  rsumwq = len;
}

static void
varrayPrevarsumMV(size_t len, const Varray<double> &v, double missval,
                  double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = rsumw = rsumwq = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        rsum += v[i];
        rsumq += v[i] * v[i];
        rsumw += 1;
        rsumwq += 1;
      }
}

double
varrayVar(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayPrevarsumMV(len, v, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayPrevarsum(len, v, rsum, rsumw, rsumq, rsumwq);

  double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

double
varrayVar1(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayPrevarsumMV(len, v, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayPrevarsum(len, v, rsum, rsumw, rsumq, rsumwq);

  double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

static void
varrayWeightedPrevarsum(size_t len, const Varray<double> &v, const Varray<double> &w,
                        double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = rsumw = rsumwq = 0;
  for (size_t i = 0; i < len; i++)
    {
      rsum += w[i] * v[i];
      rsumq += w[i] * v[i] * v[i];
      rsumw += w[i];
      rsumwq += w[i] * w[i];
    }
}

static void
varrayWeightedPrevarsumMV(size_t len, const Varray<double> &v, const Varray<double> &w, double missval,
                          double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = rsumw = rsumwq = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval) && !DBL_IS_EQUAL(w[i], missval))
      {
        rsum += w[i] * v[i];
        rsumq += w[i] * v[i] * v[i];
        rsumw += w[i];
        rsumwq += w[i] * w[i];
      }
}

double
varrayWeightedVar(size_t len, const Varray<double> &v, const Varray<double> &w, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayWeightedPrevarsumMV(len, v, w, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayWeightedPrevarsum(len, v, w, rsum, rsumw, rsumq, rsumwq);

  double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

double
varrayWeightedVar1(size_t len, const Varray<double> &v, const Varray<double> &w, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayWeightedPrevarsumMV(len, v, w, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayWeightedPrevarsum(len, v, w, rsum, rsumw, rsumq, rsumwq);

  double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

static void
varrayPrekurtsum(size_t len, const Varray<double> &v, const double mean,
                 double &rsum3w, double &rsum4w, double &rsum2diff, double &rsum4diff)
{
  rsum2diff = rsum4diff = 0;
  for (size_t i = 0; i < len; i++)
    {
      const auto vdiff = v[i] - mean;
      rsum2diff += vdiff * vdiff;
      rsum4diff += vdiff * vdiff * vdiff * vdiff;
    }
  rsum3w = len;
  rsum4w = len;
}

static void
varrayPrekurtsumMV(size_t len, const Varray<double> &v, double missval, const double mean,
                   double &rsum3w, double &rsum4w, double &rsum2diff, double &rsum4diff)
{
  rsum3w = rsum4w = rsum2diff = rsum4diff = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        const auto vdiff = v[i] - mean;
        rsum2diff += vdiff * vdiff;
        rsum4diff += vdiff * vdiff * vdiff * vdiff;
        rsum3w += 1;
        rsum4w += 1;
      }
}

double
varrayKurt(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum3w;  // 3rd moment variables
  double rsum4w;  // 4th moment variables
  double rsum2diff, rsum4diff;
  double rsum, rsumw;
  if (nmiss)
    {
      varrayPrevarsum0MV(len, v, missval, rsum, rsumw);
      varrayPrekurtsumMV(len, v, missval, (rsum / rsumw), rsum3w, rsum4w, rsum2diff, rsum4diff);
    }
  else
    {
      varrayPrevarsum0(len, v, rsum, rsumw);
      varrayPrekurtsum(len, v, (rsum / rsumw), rsum3w, rsum4w, rsum2diff, rsum4diff);
    }

  if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum2diff, 0.0)) return missval;

  double rkurt = ((rsum4diff / rsum3w) / std::pow(rsum2diff / rsum3w, 2)) - 3.0;
  if (rkurt < 0 && rkurt > -1.e-5) rkurt = 0;

  return rkurt;
}

static void
varrayPreskewsum(size_t len, const Varray<double> &v, const double mean,
                 double &rsum3w, double &rsum4w, double &rsum3diff, double &rsum2diff)
{
  rsum3diff = rsum2diff = 0;
  for (size_t i = 0; i < len; i++)
    {
      const auto vdiff = v[i] - mean;
      rsum3diff += vdiff * vdiff * vdiff;
      rsum2diff += vdiff * vdiff;
    }
  rsum3w = len;
  rsum4w = len;
}

static void
varrayPreskewsumMV(size_t len, const Varray<double> &v, double missval, const double mean,
                   double &rsum3w, double &rsum4w, double &rsum3diff, double &rsum2diff)
{
  rsum3w = rsum4w = rsum3diff = rsum2diff = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        const auto vdiff = v[i] - mean;
        rsum3diff += vdiff * vdiff * vdiff;
        rsum2diff += vdiff * vdiff;
        rsum3w += 1;
        rsum4w += 1;
      }
}

double
varraySkew(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum3w;  // 3rd moment variables
  double rsum4w;  // 4th moment variables
  double rsum3diff, rsum2diff;
  double rsum, rsumw;
  if (nmiss)
    {
      varrayPrevarsum0MV(len, v, missval, rsum, rsumw);
      varrayPreskewsumMV(len, v, missval, (rsum / rsumw), rsum3w, rsum4w, rsum3diff, rsum2diff);
    }
  else
    {
      varrayPrevarsum0(len, v, rsum, rsumw);
      varrayPreskewsum(len, v, (rsum / rsumw), rsum3w, rsum4w, rsum3diff, rsum2diff);
    }

  if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum3w, 1.0) || IS_EQUAL(rsum2diff, 0.0)) return missval;

  double rskew = (rsum3diff / rsum3w) / std::pow((rsum2diff) / (rsum3w - 1.0), 1.5);
  if (rskew < 0 && rskew > -1.e-5) rskew = 0;

  return rskew;
}