Utils.java [src/java/m/prms/model] Revision: default  Date:
/*
 * 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 m.prms.model;

import csip.utils.TextParser;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;
import java.io.PrintWriter;
import java.text.ParseException;
import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Date;
import java.util.HashMap;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
import java.util.TreeMap;
import oms3.Conversions;
import oms3.ObjectiveFunction;
import oms3.io.CSProperties;
import oms3.io.CSTable;
import oms3.io.DataIO;
import static oms3.io.DataIO.asArray;
import oms3.io.MemoryTable;
import org.apache.commons.lang3.StringUtils;

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

  static Map<String, Integer> t = new HashMap<>();


  static {
    t.put("ssstor_init_frac", 2);
    t.put("tstorm_mo", 1);
    t.put("va_open_exp", 2);
    t.put("basin_spring_frost", 1);
    t.put("nhm_seg", 1);
    t.put("jh_coef_hru", 2);
    t.put("pref_flow_den", 2);
    t.put("soil_moist_max", 2);
    t.put("snarea_curve", 2);
    t.put("covden_sum", 2);
    t.put("hru_percent_imperv", 2);
    t.put("soil2gw_max", 2);
    t.put("dprst_frac", 2);
    t.put("temp_units", 1);
    t.put("fall_frost", 1);
    t.put("epan_coef", 2);
    t.put("elev_units", 1);
    t.put("tmax_allrain_offset", 2);
    t.put("freeh2o_cap", 2);
    t.put("dday_slope", 2);
    t.put("print_type", 1);
    t.put("adjmix_rain", 2);
    t.put("rain_cbh_adj", 2);
    t.put("segment_type", 1);
    t.put("smidx_coef", 2);
    // t.put("soil_moist_init", 2);
    t.put("print_freq", 1);
    t.put("settle_const", 2);
    t.put("rad_trncf", 2);
    t.put("precip_units", 1);
    t.put("gwsink_coef", 2);
    t.put("cov_type", 1);
    t.put("hru_lat", 2);
    t.put("K_coef", 2);
    t.put("tmax_index", 2);
    t.put("den_init", 2);
    t.put("obsin_segment", 1);
    t.put("op_flow_thres", 2);
    t.put("transp_tmax", 2);
    t.put("hru_type", 1);
    t.put("hru_segment_nhm", 1);
    t.put("dprst_seep_rate_clos", 2);
    t.put("hru_x", 2);
    t.put("ssr2gw_exp", 2);
    t.put("smidx_exp", 2);
    t.put("dprst_frac_open", 2);
    t.put("poi_type", 1);
    t.put("radj_wppt", 2);
    t.put("snow_intcp", 2);
    t.put("hru_lon", 2);
    t.put("radadj_intcp", 2);
    // t.put("ssstor_init", 2);
    t.put("snowinfil_max", 2);
    t.put("sat_threshold", 2);
    t.put("tmax_allsnow", 2);
    t.put("cecn_coef", 2);
    t.put("dprst_et_coef", 2);
    t.put("hru_segment", 1);
    t.put("ppt_rad_adj", 2);
    t.put("dday_intcp", 2);
    t.put("radmax", 2);
    t.put("soil_rechr_init_frac", 2);
    t.put("basin_fall_frost", 1);
    t.put("slowcoef_lin", 2);
    t.put("den_max", 2);
    t.put("x_coef", 2);
    t.put("radj_sppt", 2);
    t.put("tosegment", 1);
    t.put("snowpack_init", 2);
    t.put("dprst_depth_avg", 2);
    t.put("sro_to_dprst_imperv", 2);
    t.put("gwstor_min", 2);
    t.put("hru_elev", 2);
    t.put("hru_deplcrv", 1);
    t.put("soil_moist_init_frac", 2);
    t.put("segment_flow_init", 2);
    t.put("outlet_sta", 1);
    t.put("fastcoef_sq", 2);
    t.put("carea_max", 2);
    t.put("gwflow_coef", 2);
    t.put("gwstor_init", 2);
    t.put("nhm_id", 1);
    t.put("hru_aspect", 2);
    t.put("tmax_cbh_adj", 2);
    t.put("nhm_deplcrv", 1);
    t.put("tosegment_nhm", 1);
    t.put("dprst_frac_init", 2);
    t.put("runoff_units", 1);
    t.put("tmin_cbh_adj", 2);
    t.put("poi_gage_id", 4);
    t.put("covden_win", 2);
    t.put("va_clos_exp", 2);
    t.put("jh_coef", 2);
    t.put("spring_frost", 1);
    t.put("imperv_stor_max", 2);
    t.put("wrain_intcp", 2);
    t.put("sro_to_dprst_perv", 2);
    t.put("melt_force", 1);
    t.put("dprst_flow_coef", 2);
    t.put("snarea_thresh", 2);
    t.put("potet_sublim", 2);
    t.put("ssr2gw_rate", 2);
    t.put("soil_rechr_max_frac", 2);
    t.put("slowcoef_sq", 2);
    t.put("emis_noppt", 2);
    t.put("srain_intcp", 2);
    t.put("radadj_slope", 2);
    t.put("hru_slope", 2);
    t.put("soil_type", 1);
    t.put("albset_sna", 2);
    t.put("albset_rna", 2);
    t.put("melt_look", 1);
    t.put("albset_rnm", 2);
    t.put("albset_snm", 2);
    t.put("fastcoef_lin", 2);
    t.put("snow_cbh_adj", 2);
    t.put("hru_y", 2);
    t.put("hru_area", 2);
    t.put("dprst_seep_rate_open", 2);
    t.put("transp_end", 1);
    t.put("transp_beg", 1);
    t.put("poi_gage_segment", 1);
    // t.put("soil_rechr_init", 2);
  }


  public static File deriveFileName(File file, String ext) {
    String name = file.getName();
    return new File(file.getParentFile(), name.substring(0, name.indexOf('.')) + ext);
  }


  public static File convertCSVToData(File csv, File dir) throws IOException, ParseException {
    String name = csv.getName();
    File parent = csv.getParentFile();
    if (dir != null) {
      parent = dir;
    }
    File result = new File(parent, name.substring(0, name.indexOf('.')) + ".data");
    CSTable t = DataIO.table(csv);

    PrintWriter out = new PrintWriter(result);
    StringBuilder header = new StringBuilder();
    // Patterson: column count is one short for some reason.
    for (int i = 0; i < t.getColumnCount() + 1; i++) {
      header.append(t.getColumnName(i));
      header.append(" ");
    }

    out.println(t.getName());

    String cols = createDataHeader(header.toString());
    out.print(cols);
    out.println("#########################################################################");

    String datefmt = t.getInfo().get(DataIO.DATE_FORMAT);

    SimpleDateFormat filefmt = new SimpleDateFormat(datefmt);
    SimpleDateFormat mmsfmt = new SimpleDateFormat("yyyy MM dd H m s");

    for (String[] row : t.rows()) {
      String time = mmsfmt.format(filefmt.parse(row[1]));
      out.print(time + " ");
      for (int i = 2; i < row.length; i++) {
        out.print(row[i]);
        out.print(' ');
      }
      out.println();
    }
    out.close();
    return result;
  }

  static final List<String> vars = Arrays.asList("precip", "runoff", "tmax", "tmin", "rhum");


  static String createDataHeader(String csvHeader) {
    Map<Integer, String> order = new TreeMap<>();
    // find the index in the string, sort it using the map keys
    for (String var : vars) {
      if (csvHeader.indexOf(var) > 0) {
        order.put(csvHeader.indexOf(var), var);
      }
    }
    StringBuilder b = new StringBuilder();
    for (Integer i : order.keySet()) {
      b.append(String.format("%s %s \n", order.get(i), StringUtils.countMatches(csvHeader, order.get(i))));
    }
    return b.toString();
  }


  public static File convertStatvarToCSV(File statvar) throws IOException {
    File csvFile = deriveFileName(statvar, ".csv");

    BufferedReader r = new BufferedReader(new FileReader(statvar));
    int vars = Integer.parseInt(r.readLine());

    Map<String, Integer> cols = new LinkedHashMap<>();

    for (int i = 0; i < vars; i++) {
      String v = r.readLine();
      String[] content = v.split("\\s+");
      if (content.length != 2) {
        throw new RuntimeException("invalid column length");
      }
      cols.put(content[0], new Integer(content[1]));
    }

    PrintWriter w = new PrintWriter(csvFile);
    w.println("@T, statvar conversion");
    w.println("created_at, " + new Date());
    w.println("converted_from, " + statvar);
    w.println("date_format, yyyy MM dd H m s");

    w.print("@H, date");
    for (Entry<String, Integer> entry : cols.entrySet()) {
      if (entry.getValue() == 1) {
        w.print(", " + entry.getKey());
      } else {
        for (int i = 0; i < entry.getValue(); i++) {
          w.print(", " + entry.getKey() + "[" + i + "]");
        }
      }
    }

    w.println();
    String line = null;
    while ((line = r.readLine()) != null) {
      String[] l = line.split("\\s+");
      // skip first row, then all dates
      w.print(',');
      for (int i = 1; i < 7; i++) {
        w.print(" " + l[i]);
      }
      for (int i = 7; i < l.length; i++) {
        w.print("," + l[i]);
      }
      w.println();
    }

    r.close();
    w.close();
    return csvFile;
  }


  /**
   * Converts a OMS CSV parameter file to a PRMS parameter file.
   *
   * @param csv
   * @return
   * @throws Exception
   */
  public static File convertToPRMSParam(File csv) throws Exception {
    String name = csv.getName();
    CSProperties p = readParameter(csv);
    fixParameters(p);
    File result = new File(csv.getParentFile(), name.substring(0, name.indexOf('.')) + ".params");
    writeParamFile(p, result);
    return result;
  }


  static void writeParamFile(CSProperties p, File result) throws FileNotFoundException, ParseException {
    PrintWriter w = new PrintWriter(result);
    w.println("Generated by csip-prms");
    w.println("Version: 1.7");

    // write dimensions
    List<String> keys1D = DataIO.keysByMeta(p, "role", "dimension");
    w.println("** Dimensions **");
    for (String k : keys1D) {
      w.println("####");
      w.println(k);
      w.println(p.get(k));
    }
    // write parameter
    w.println("** Parameters **");

    // scalars
    List<String> scalars = DataIO.keysByNotMeta(p, "bound");
    for (String s : scalars) {
      if (keys1D.contains(s)) {
        continue;
      }
      w.println("####");
      w.println(s);
      w.println("1");
      w.println("one");
      w.println("1");
      String o = p.get(s).toString();
      w.println(getParameterTypeKey(s));
      w.println(o);
    }

    // 1D parameter
    for (String s : p.keySet()) {
      if (DataIO.isBound(p, s, 1)) {
        String[] d = s.split("\\$");

        // Added check to make sure that the paramter exists
        if (t.get(d[1]) != null) {
          w.println("####");
          w.println(d[1]);
          w.println("1");
          w.println(d[0]);
          w.println(p.get(d[0]));

          String[] d1 = Conversions.convert(p.get(s), String[].class);

          w.println(getParameterTypeKey(d[1]));

          for (String elem : d1) {
            w.println(elem.trim());
          }
        }
      }
    }

    // 2D parameter
    for (String s : p.keySet()) {
      if (DataIO.isBound(p, s, 2)) {
        String m = p.getInfo(s).get("bound");
        w.println("####");
        w.println(s);
        w.println("2");
        String[] dims = m.split(",");

        final int rows = DataIO.getInt(p, dims[0].trim());
        final int cols = DataIO.getInt(p, dims[1].trim());

        w.println(dims[0].trim());
        w.println(dims[1].trim());
        w.println(rows * cols);

        String[][] d = Conversions.convert(p.get(s), String[][].class);
        w.println(getParameterTypeKey(s));
        for (String[] row : d) {
          for (String col : row) {
            w.println(col.trim());
          }
        }

        System.out.println("2dim : " + s);
      }
    }
    w.close();
  }


  static CSProperties readParameter(File file) throws IOException {
    CSProperties p = DataIO.properties();
    if (DataIO.propertyExists("Parameter", file)) {
      // original properties.                    
      p.putAll(DataIO.properties(new FileReader(file), "Parameter"));
    }
    // check for tables in the file.
    List<String> tables = DataIO.tables(file);
    if (!tables.isEmpty()) {
      for (String name : tables) {
        CSTable t = DataIO.table(file, name);
        // convert them to Properties.
        CSProperties prop;
        try {
          prop = DataIO.toProperties(t);
        } catch (ArrayIndexOutOfBoundsException ex) {
          throw new RuntimeException("Bad array index getting properties for " + name + ": " + ex);
        }
        //prop.().put("file", f);
        p.putAll(prop);
      }
    }
    return p;
  }


  static String getParameterTypeKey(String name) {
    Integer i = t.get(name);
    if (i == null) {
      throw new IllegalArgumentException("No type information for :" + name);
    }
    return i.toString();
  }


  static String getParameterType(String p) {
    try {
      Integer.parseInt(p);
      return "1";
    } catch (NumberFormatException E) {
      try {
        Double.parseDouble(p);
        return "2";
      } catch (NumberFormatException E1) {
        return "4";
      }
    }
  }


  /**
   * Look for parameters that need to be updated before they can be converted to
   * mms props.
   * @param p
   */
  static void fixParameters(CSProperties p) throws ParseException {
    // Convert to 1D by updating ndeplval to ndepl * ndeplval
    convertTo1D(p, "snarea_curve", "ndeplval");

    for (String s : p.keySet()) {
      if (s.contains("$spring_frost")) {
        // Check if this needs to be populated
        Map<String, String> info = p.getInfo(s);
        int len = Integer.parseInt(info.get("len"));
        int[] springFrost = new int[len];
        for (int i = 0; i < len; i++) {
          // Setting to day 120
          springFrost[i] = 120;
        }
        String sArray = Conversions.convert(springFrost, String.class);
        p.put(s, sArray);
      }
      if (s.equals("nssr")) {
        String nhru = (String) p.get("nhru");
        p.put(s, nhru);
      }
    }
  }


  static void convertTo1D(CSProperties p, String tableKey, String boundsKey) {
    CSTable as2Table = DataIO.as2DTable(p, tableKey);
    MemoryTable oneDTable = new MemoryTable();
    String m = p.getInfo(tableKey).get("bound");

    double[][] d = asArray(as2Table);
    List<Double> rowList = new ArrayList<>();
    for (int i = 0; i < d.length; i++) {
      for (int j = 0; j < d[0].length; j++) {
        rowList.add(d[i][j]);
      }
    }

    // Replace parameter
    Map<String, String> m1 = p.getInfo(tableKey);
    m1.put("bound", boundsKey);
    m1.put("len", Integer.toString(rowList.size()));
    p.setInfo(tableKey, m1);
    // p.put(s, Conversions.convert(rowList.toArray(new Double[rowList.size()]), String.class));
    // Need to get a string array
    Double[] dArray = rowList.toArray(new Double[rowList.size()]);
    double[] smallDArray = Conversions.convert(dArray, double[].class);
    String sArray = Conversions.convert(smallDArray, String.class);

    // Update parameter with new table length
    p.put(boundsKey, m1.get("len"));

    String fixedKey = boundsKey + "$" + tableKey;
    p.put(fixedKey, sArray);
    p.setInfo(fixedKey, m1);
    p.remove(tableKey);
  }

  /**
   *
   */
  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 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++;
        }
      }
      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])), 2);
          var += Math.pow(Math.abs(Math.log(obs[i]) - avg), 2);
        }
      }
      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;
    }
  }


  public static double[] extractCSVColumn(File file, String name, int skiplines) throws IOException {
    try (TextParser p = new TextParser(file).autoClose(false)) {
      String[] h = p.nextLine().tokens(",").asStringArray();
      int idx = -1;
      for (int i = 0; i < h.length; i++) {
        if (h[i].equalsIgnoreCase(name)) {
          idx = i;
          break;
        }
      }
      if (idx == -1)
        throw new IllegalArgumentException("Not found: " + name);

      p.nextLine(skiplines);
      List<Double> l = new ArrayList<>();
      while (p.nextLineSkipEmpty().notEOF()) {
        String[] row = p.tokens(",").asStringArray();
        l.add(Double.parseDouble(row[idx]));
      }
      return l.stream().mapToDouble(Double::doubleValue).toArray();
    }
  }


  public static double[] extractStatvarColumn(File file, String name) throws IOException {
    try (TextParser p = new TextParser(file).autoClose(false)) {
      int h = p.nextLine().asInteger();
      int idx = 6;
      for (int i = 0; i < h; i++) {
        String[] h1 = p.nextLine().tokens().asStringArray();
        idx += Integer.parseInt(h1[1]);
        if (h1[0].equals(name))
          break;
      }
      if (idx == 6)
        throw new IllegalArgumentException("Not found: " + name);

      p.toLineStaringWith("1");
      List<Double> l = new ArrayList<>();
      while (p.nextLineSkipEmpty().notEOF()) {
        String[] row = p.tokens().asStringArray();
        l.add(Double.parseDouble(row[idx]));
      }
      return l.stream().mapToDouble(Double::doubleValue).toArray();
    }
  }


  public static void main(String[] args) throws IOException, ParseException {
    String h1 = createDataHeader("date,runoff[0],precip[0],precip[1],precip[2],precip[3],precip[4],tmin[0],tmin[1],tmax[0],tmax[1]");
    System.out.println(h1);

    File f = convertStatvarToCSV(new File("/od/projects/csip-all/csip-prms/data/sagehen.statvar"));
    System.out.println(f);

//    File f = convertCSVToData(new File("/od/oms/oms_examples/oms3.prj.prms/data/data.csv"), null);
//    long start = System.currentTimeMillis();
//    File f = convertCSVToData(new File("/od/projects/nwcc2015.wsf/simulation/teton_kr.old/data/io/data-hru.csv"), new File("/tmp"));
//    long end = System.currentTimeMillis();
//    System.out.println("Time " + (end - start));
//    System.out.println(f);
  }

}