/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

package util;

import oms3.ObjectiveFunction;

/**
 *
 * @author od
 */
public class ObjFuncs {

  /**
   *
   */
  public static class PBIAS implements ObjectiveFunction {

    @Override
    public double calculate(double[] obs, double[] sim, double missing) {
      if (sim.length != obs.length) {
        throw new IllegalArgumentException("obs/sim length differ: " + obs.length + "!=" + sim.length);
      }
      double diffsum = 0;
      double obssum = 0;
      for (int i = 0; i < sim.length; i++) {
        if (obs[i] > missing) {
          diffsum += sim[i] - obs[i];
          obssum += obs[i];
        }
      }
      return (diffsum / obssum) * 100.0;
    }


    @Override
    public boolean positiveDirection() {
      return false;
    }
  }

  /**
   *
   */
  public static class NSLOG1P implements ObjectiveFunction {

    @Override
    public double calculate(double[] obs, double[] sim, double missing) {
      if (sim.length != obs.length) {
        throw new IllegalArgumentException("obs/sim length differ: " + obs.length + "!=" + sim.length);
      }

      /**
       * calculating logarithmic values of both data sets. Sets 0 if data is 0
       */
      int valid = 0;
      double avg = 0.0;
      for (int i = 0; i < sim.length; i++) {
        if (sim[i] >= 0.0 && obs[i] >= 0.0) {
          // summing up 
          avg += Math.log1p(obs[i]);
          valid++;
        }
      }

      if (valid < 2) {
        return Double.NEGATIVE_INFINITY;
      }

      // calculating mean 
      avg /= valid;

      // calculating mean pow deviations
      double rmse = 0.0;
      double e = 0.0;
      for (int i = 0; i < sim.length; i++) {
        if (sim[i] >= 0 && obs[i] >= 0) {
          double l1po = Math.log1p(obs[i]);
          rmse += Math.pow(Math.abs(l1po - Math.log1p(sim[i])), 2);
          e += Math.pow(Math.abs(l1po - avg), 2);
        }
      }
      double r = 1 - (rmse / e);
      return Double.isNaN(r) ? 0.0 : r;
    }


    @Override
    public boolean positiveDirection() {
      return true;
    }
  }

  /**
   *
   */
  public static class NSLOG2 implements ObjectiveFunction {

    @Override
    public double calculate(double obs[], double sim[], double missing) {
      if (obs.length != sim.length) {
        throw new IllegalArgumentException("obs/sim length differ: " + obs.length + "!=" + sim.length);
      }
      double pow = 2;
      double rsme = 0;
      double var = 0;
      double avg = 0;
      double count = 0;
      for (int i = 0; i < obs.length; i++) {
        if (obs[i] > 0 && obs[i] != missing) {
          avg += Math.log(obs[i]);
          count += 1;
        }
      }
      avg /= count;

      for (int i = 0; i < obs.length; i++) {
        if (obs[i] > 0 && sim[i] > 0 && obs[i] != missing) {
          rsme += Math.pow(Math.abs(Math.log(obs[i]) - Math.log(sim[i])), pow);
          var += Math.pow(Math.abs(Math.log(obs[i]) - avg), pow);
        }
      }
      double result = 1.0 - (rsme / var);
      if (Double.isNaN(result)) {
        result = 0;
      }
      return result;
    }


    @Override
    public boolean positiveDirection() {
      return true;
    }
  }

  /**
   * KGE 2012
   */
  public static class KGE implements ObjectiveFunction {

    @Override
    public double calculate(double obs[], double sim[], double missing) {
      if (obs.length != sim.length) {
        throw new IllegalArgumentException("obs/sim length differ: " + obs.length + "!=" + sim.length);
      }
      int contamedia = 0;
      double sommamediaoss = 0;
      double sommamediasim = 0;
      for (int i = 0; i < obs.length; i++) {
        if (obs[i] > missing) {
          contamedia++;
          sommamediaoss += obs[i];
          sommamediasim += sim[i];
        }
      }
      double mediaoss = sommamediaoss / contamedia;
      double mediasim = sommamediasim / contamedia;
      int count = 0;
      double numvaprev = 0;
      double coef1_den = 0;
      double numR = 0;
      double den1R = 0;
      double den2R = 0;
      for (int i = 0; i < obs.length; i++) {
        if (obs[i] > missing) {
          count++;
          coef1_den += (obs[i] - mediaoss) * (obs[i] - mediaoss);
          numR += (obs[i] - mediaoss) * (sim[i] - mediasim);
          den1R += (obs[i] - mediaoss) * (obs[i] - mediaoss);
          den2R += (sim[i] - mediasim) * (sim[i] - mediasim);
          numvaprev += (sim[i] - mediasim) * (sim[i] - mediasim);
        }
      }
      double sdosservati = Math.sqrt(coef1_den / (count - 1));
      double sdsimulati = Math.sqrt(numvaprev / (count - 1));
      double R = numR / (Math.sqrt(den1R) * Math.sqrt(den2R));
      //double alpha = sdsimulati / sdosservati; // 2009
      double beta = mediasim / mediaoss;
      double gamma = (mediaoss * sdsimulati) / (sdosservati * mediasim); // 2012
      //return 1 - Math.sqrt((R - 1) * (R - 1) + (alpha - 1) * (alpha - 1) + (beta - 1) * (beta - 1)); // 2009
      return 1 - Math.sqrt((R - 1) * (R - 1) + (gamma - 1) * (gamma - 1) + (beta - 1) * (beta - 1)); // 2012
    }


    @Override
    public boolean positiveDirection() {
      return true;
    }
  }

}