/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      https://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  
  /*
   * This is not the original file distributed by the Apache Software Foundation
   * It has been modified by the Hipparchus project
   */
  package org.hipparchus.stat.inference;
  
  import java.util.Map;
  import java.util.TreeMap;
  import java.util.stream.LongStream;
  
  import org.hipparchus.distribution.continuous.NormalDistribution;
  import org.hipparchus.exception.LocalizedCoreFormats;
  import org.hipparchus.exception.MathIllegalArgumentException;
  import org.hipparchus.exception.MathIllegalStateException;
  import org.hipparchus.exception.NullArgumentException;
  import org.hipparchus.stat.LocalizedStatFormats;
  import org.hipparchus.stat.ranking.NaNStrategy;
  import org.hipparchus.stat.ranking.NaturalRanking;
  import org.hipparchus.stat.ranking.TiesStrategy;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.Precision;
  
  /**
   * An implementation of the Mann-Whitney U test.
   * <p>
   * The definitions and computing formulas used in this implementation follow
   * those in the article,
   * <a href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney U
   * Test</a>
   * <p>
   * In general, results correspond to (and have been tested against) the R
   * wilcox.test function, with {@code exact} meaning the same thing in both APIs
   * and {@code CORRECT} uniformly true in this implementation. For example,
   * wilcox.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE
   * correct = TRUE) will return the same p-value as mannWhitneyUTest(x, y,
   * false). The minimum of the W value returned by R for wilcox.test(x, y...) and
   * wilcox.test(y, x...) should equal mannWhitneyU(x, y...).
   */
  public class MannWhitneyUTest { // NOPMD - this is not a Junit test class, PMD false positive here
  
      /**
       * If the combined dataset contains no more values than this, test defaults to
       * exact test.
       */
      private static final int SMALL_SAMPLE_SIZE = 50;
  
      /** Ranking algorithm. */
      private final NaturalRanking naturalRanking;
  
      /** Normal distribution */
      private final NormalDistribution standardNormal;
  
      /**
       * Create a test instance using where NaN's are left in place and ties get
       * the average of applicable ranks.
       */
      public MannWhitneyUTest() {
          naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
                                              TiesStrategy.AVERAGE);
          standardNormal = new NormalDistribution(0, 1);
      }
  
      /**
       * Create a test instance using the given strategies for NaN's and ties.
       *
       * @param nanStrategy specifies the strategy that should be used for
       *        Double.NaN's
       * @param tiesStrategy specifies the strategy that should be used for ties
       */
      public MannWhitneyUTest(final NaNStrategy nanStrategy,
                              final TiesStrategy tiesStrategy) {
          naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
          standardNormal = new NormalDistribution(0, 1);
      }
  
      /**
       * Computes the
       * <a href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U">
       * Mann-Whitney U statistic</a> comparing means for two independent samples
       * possibly of different lengths.
       * <p>
       * This statistic can be used to perform a Mann-Whitney U test evaluating
      * the null hypothesis that the two independent samples have equal mean.
      * <p>
      * Let X<sub>i</sub> denote the i'th individual of the first sample and
      * Y<sub>j</sub> the j'th individual in the second sample. Note that the
      * samples can have different lengths.
      * <p>
      * <strong>Preconditions</strong>:
      * <ul>
      * <li>All observations in the two samples are independent.</li>
      * <li>The observations are at least ordinal (continuous are also
      * ordinal).</li>
      * </ul>
      *
      * @param x the first sample
      * @param y the second sample
      * @return Mann-Whitney U statistic (minimum of U<sup>x</sup> and
      *         U<sup>y</sup>)
      * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
      * @throws MathIllegalArgumentException if {@code x} or {@code y} are
      *         zero-length.
      */
     public double mannWhitneyU(final double[] x, final double[] y)
         throws MathIllegalArgumentException, NullArgumentException {
 
         ensureDataConformance(x, y);
 
         final double[] z = concatenateSamples(x, y);
         final double[] ranks = naturalRanking.rank(z);
 
         double sumRankX = 0;
 
         /*
          * The ranks for x is in the first x.length entries in ranks because x
          * is in the first x.length entries in z
          */
         for (int i = 0; i < x.length; ++i) {
             sumRankX += ranks[i];
         }
 
         /*
          * U1 = R1 - (n1 * (n1 + 1)) / 2 where R1 is sum of ranks for sample 1,
          * e.g. x, n1 is the number of observations in sample 1.
          */
         final double U1 = sumRankX - ((long) x.length * (x.length + 1)) / 2;
 
         /*
          * U1 + U2 = n1 * n2
          */
         final double U2 = (long) x.length * y.length - U1;
 
         return FastMath.min(U1, U2);
     }
 
     /**
      * Concatenate the samples into one array.
      *
      * @param x first sample
      * @param y second sample
      * @return concatenated array
      */
     private double[] concatenateSamples(final double[] x, final double[] y) {
         final double[] z = new double[x.length + y.length];
 
         System.arraycopy(x, 0, z, 0, x.length);
         System.arraycopy(y, 0, z, x.length, y.length);
 
         return z;
     }
 
     /**
      * Returns the asymptotic <i>observed significance level</i>, or
      * <a href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
      * p-value</a>, associated with a <a href=
      * "http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U">Mann-Whitney U
      * Test</a> comparing means for two independent samples.
      * <p>
      * Let X<sub>i</sub> denote the i'th individual of the first sample and
      * Y<sub>j</sub> the j'th individual in the second sample.
      * <p>
      * <strong>Preconditions</strong>:
      * <ul>
      * <li>All observations in the two samples are independent.</li>
      * <li>The observations are at least ordinal.</li>
      * </ul>
      * <p>
      * If there are no ties in the data and both samples are small (less than or
      * equal to 50 values in the combined dataset), an exact test is performed;
      * otherwise the test uses the normal approximation (with continuity
      * correction).
      * <p>
      * If the combined dataset contains ties, the variance used in the normal
      * approximation is bias-adjusted using the formula in the reference above.
      *
      * @param x the first sample
      * @param y the second sample
      * @return approximate 2-sized p-value
      * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
      * @throws MathIllegalArgumentException if {@code x} or {@code y} are
      *         zero-length
      */
     public double mannWhitneyUTest(final double[] x, final double[] y)
         throws MathIllegalArgumentException, NullArgumentException {
         ensureDataConformance(x, y);
 
         // If samples are both small and there are no ties, perform exact test
         if (x.length + y.length <= SMALL_SAMPLE_SIZE &&
             tiesMap(x, y).isEmpty()) {
             return mannWhitneyUTest(x, y, true);
         } else { // Normal approximation
             return mannWhitneyUTest(x, y, false);
         }
     }
 
     /**
      * Returns the asymptotic <i>observed significance level</i>, or
      * <a href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
      * p-value</a>, associated with a <a href=
      * "http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U">Mann-Whitney U
      * Test</a> comparing means for two independent samples.
      * <p>
      * Let X<sub>i</sub> denote the i'th individual of the first sample and
      * Y<sub>j</sub> the j'th individual in the second sample.
      * <p>
      * <strong>Preconditions</strong>:
      * <ul>
      * <li>All observations in the two samples are independent.</li>
      * <li>The observations are at least ordinal.</li>
      * </ul>
      * <p>
      * If {@code exact} is {@code true}, the p-value reported is exact, computed
      * using the exact distribution of the U statistic. The computation in this
      * case requires storage on the order of the product of the two sample
      * sizes, so this should not be used for large samples.
      * <p>
      * If {@code exact} is {@code false}, the normal approximation is used to
      * estimate the p-value.
      * <p>
      * If the combined dataset contains ties and {@code exact} is {@code true},
      * MathIllegalArgumentException is thrown. If {@code exact} is {@code false}
      * and the ties are present, the variance used to compute the approximate
      * p-value in the normal approximation is bias-adjusted using the formula in
      * the reference above.
      *
      * @param x the first sample
      * @param y the second sample
      * @param exact true means compute the p-value exactly, false means use the
      *        normal approximation
      * @return approximate 2-sided p-value
      * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
      * @throws MathIllegalArgumentException if {@code x} or {@code y} are
      *         zero-length or if {@code exact} is {@code true} and ties are
      *         present in the data
      */
     public double mannWhitneyUTest(final double[] x, final double[] y,
                                    final boolean exact)
         throws MathIllegalArgumentException, NullArgumentException {
         ensureDataConformance(x, y);
         final Map<Double, Integer> tiesMap = tiesMap(x, y);
         final double u = mannWhitneyU(x, y);
         if (exact) {
             if (!tiesMap.isEmpty()) {
                 throw new MathIllegalArgumentException(LocalizedStatFormats.TIES_ARE_NOT_ALLOWED);
             }
             return exactP(x.length, y.length, u);
         }
 
         return approximateP(u, x.length, y.length,
                             varU(x.length, y.length, tiesMap));
     }
 
     /**
      * Ensures that the provided arrays fulfills the assumptions.
      *
      * @param x first sample
      * @param y second sample
      * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
      * @throws MathIllegalArgumentException if {@code x} or {@code y} are
      *         zero-length.
      */
     private void ensureDataConformance(final double[] x, final double[] y)
         throws MathIllegalArgumentException, NullArgumentException {
 
         if (x == null || y == null) {
             throw new NullArgumentException();
         }
         if (x.length == 0 || y.length == 0) {
             throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
         }
     }
 
     /**
      * Estimates the 2-sided p-value associated with a Mann-Whitney U statistic
      * value using the normal approximation.
      * <p>
      * The variance passed in is assumed to be corrected for ties. Continuity
      * correction is applied to the normal approximation.
      *
      * @param u Mann-Whitney U statistic
      * @param n1 number of subjects in first sample
      * @param n2 number of subjects in second sample
      * @param varU variance of U (corrected for ties if these exist)
      * @return two-sided asymptotic p-value
      * @throws MathIllegalStateException if the p-value can not be computed due
      *         to a convergence error
      * @throws MathIllegalStateException if the maximum number of iterations is
      *         exceeded
      */
     private double approximateP(final double u, final int n1, final int n2,
                                 final double varU)
         throws MathIllegalStateException {
 
         final double mu = (long) n1 * n2 / 2.0;
 
         // If u == mu, return 1
         if (Precision.equals(mu, u)) {
             return 1;
         }
 
         // Force z <= 0 so we get tail probability. Also apply continuity
         // correction
         final double z = -Math.abs((u - mu) + 0.5) / FastMath.sqrt(varU);
 
         return 2 * standardNormal.cumulativeProbability(z);
     }
 
     /**
      * Calculates the (2-sided) p-value associated with a Mann-Whitney U
      * statistic.
      * <p>
      * To compute the p-value, the probability densities for each value of U up
      * to and including u are summed and the resulting tail probability is
      * multiplied by 2.
      * <p>
      * The result of this computation is only valid when the combined n + m
      * sample has no tied values.
      * <p>
      * This method should not be used for large values of n or m as it maintains
      * work arrays of size n*m.
      *
      * @param u Mann-Whitney U statistic value
      * @param n first sample size
      * @param m second sample size
      * @return two-sided exact p-value
      */
     private double exactP(final int n, final int m, final double u) {
         final double nm = m * n;
         if (u > nm) { // Quick exit if u is out of range
             return 1;
         }
         // Need to convert u to a mean deviation, so cumulative probability is
         // tail probability
         final double crit = u < nm / 2 ? u : nm / 2 - u;
 
         double cum = 0d;
         for (int ct = 0; ct <= crit; ct++) {
             cum += uDensity(n, m, ct);
         }
         return 2 * cum;
     }
 
     /**
      * Computes the probability density function for the Mann-Whitney U
      * statistic.
      * <p>
      * This method should not be used for large values of n or m as it maintains
      * work arrays of size n*m.
      *
      * @param n first sample size
      * @param m second sample size
      * @param u U-statistic value
      * @return the probability that a U statistic derived from random samples of
      *         size n and m (containing no ties) equals u
      */
     private double uDensity(final int n, final int m, double u) {
         if (u < 0 || u > m * n) {
             return 0;
         }
         final long[] freq = uFrequencies(n, m);
         return freq[(int) FastMath.round(u + 1)] /
                (double) LongStream.of(freq).sum();
     }
 
     /**
      * Computes frequency counts for values of the Mann-Whitney U statistc. If
      * freq[] is the returned array, freq[u + 1] counts the frequency of U = u
      * among all possible n-m orderings. Therefore, P(u = U) = freq[u + 1] / sum
      * where sum is the sum of the values in the returned array.
      * <p>
      * Implements the algorithm presented in "Algorithm AS 62: A Generator for
      * the Sampling Distribution of the Mann-Whitney U Statistic", L. C. Dinneen
      * and B. C. Blakesley Journal of the Royal Statistical Society. Series C
      * (Applied Statistics) Vol. 22, No. 2 (1973), pp. 269-273.
      *
      * @param n first sample size
      * @param m second sample size
      * @return array of U statistic value frequencies
      */
     private long[] uFrequencies(final int n, final int m) {
         final int max = FastMath.max(m, n);
         if (max > 100) {
             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
                                                    max, 100);
         }
         final int min = FastMath.min(m, n);
         final long[] out = new long[n * m + 2];
         final long[] work = new long[n * m + 2];
         for (int i = 1; i < out.length; i++) {
             out[i] = (i <= (max + 1)) ? 1 : 0;
         }
         work[1] = 0;
         int in = max;
         for (int i = 2; i <= min; i++) {
             work[i] = 0;
             in = in + max;
             int n1 = in + 2;
             long l = 1 + in / 2;
             int k = i;
             for (int j = 1; j <= l; j++) {
                 k++;
                 n1 = n1 - 1;
                 final long sum = out[j] + work[j];
                 out[j] = sum;
                 work[k] = sum - out[n1];
                 out[n1] = sum;
             }
         }
         return out;
     }
 
     /**
      * Computes the variance for a U-statistic associated with samples of
      * sizes{@code n} and {@code m} and ties described by {@code tiesMap}. If
      * {@code tiesMap} is non-empty, the multiplicity counts in its values set
      * are used to adjust the variance.
      *
      * @param n first sample size
      * @param m second sample size
      * @param tiesMap map of <value, multiplicity>
      * @return ties-adjusted variance
      */
     private double varU(final int n, final int m,
                         Map<Double, Integer> tiesMap) {
         final double nm = (long) n * m;
         if (tiesMap.isEmpty()) {
             return nm * (n + m + 1) / 12.0;
         }
         final long tSum = tiesMap.entrySet().stream()
             .mapToLong(e -> e.getValue() * e.getValue() * e.getValue() -
                             e.getValue())
             .sum();
         final double totalN = n + m;
         return (nm / 12) * (totalN + 1 - tSum / (totalN * (totalN - 1)));
 
     }
 
     /**
      * Creates a map whose keys are values occurring more than once in the
      * combined dataset formed from x and y. Map entry values are the number of
      * occurrences. The returned map is empty iff there are no ties in the data.
      *
      * @param x first dataset
      * @param y second dataset
      * @return map of <value, number of times it occurs> for values occurring
      *         more than once or an empty map if there are no ties (the returned
      *         map is <em>not</em> thread-safe, which is OK in the context of the callers)
      */
     private Map<Double, Integer> tiesMap(final double[] x, final double[] y) {
         final Map<Double, Integer> tiesMap = new TreeMap<>(); // NOPMD - no concurrent access in the callers context
         for (double value : x) {
             tiesMap.merge(value, 1, Integer::sum);
         }
         for (double v : y) {
             tiesMap.merge(v, 1, Integer::sum);
         }
         tiesMap.entrySet().removeIf(e -> e.getValue() == 1);
         return tiesMap;
     }
 }