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);
}
}