CDPHE_lowFlowStats.java [src/java/m/cfa/timeseries] Revision:   Date:
package m.cfa.timeseries;

import m.cfa.DoubleArray;
import m.cfa.DoubleMath;
import java.io.IOException;
import java.text.DateFormat;
import java.text.ParseException;
import java.util.ArrayList;
import java.util.Date;

/**
* Last Updated: 2-June-2017
* @author Tyler Wible
* @since 16-December-2014
*/
public class CDPHE_lowFlowStats {
    /**
     * Calculate the CDPHE "extreme value" design flow (see DFLOW user manual) 
     * based on an 'm'-day average low flow of each year fitted to a Log-Pearson
     * Type III distribution and interpolated for a return period of 'R' years
     * @param flowData  flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @param M  the number of days to average for annual low flow analysis (m-day average)
     * @param R  the desired return period of the m-day low flow (in years)
     * @param waterYearBegin  the month and day (MM-dd) of the start of the water year for this analysis
     * @return
     * @throws IOException
     * @throws ParseException 
     */
    public double CDPHE_ExtremeValue(String[][] flowData,
                                     int M,
                                     double R,
                                     String waterYearBegin) throws IOException, ParseException{
        if(flowData.length == 0){
            return -1;
        }
        int currentYear = Integer.parseInt(flowData[0][0].substring(0,4));
        int finalYear = Integer.parseInt(flowData[flowData.length - 1][0].substring(0,4));
//        //Start by pulling the first complete water year (determine if this is necessary
//        if(flowData[0][0].compareToIgnoreCase(String.valueOf(currentYear) + "-04-01") > 0){
//            //If the start date of data is after the start of the water year, reset the 
//            //begin date to next year so that only complete water years are analyzed
//            currentYear++;
//        }
        int nextYear = currentYear + 1;
        String startDate = String.valueOf(currentYear) + "-" + waterYearBegin;//String startDate = String.valueOf(currentYear) + "-04-01";
        String previousDay = DoubleArray.getDay(startDate, -1);
        String endDate = String.valueOf(nextYear) + previousDay.substring(4);//String endDate = String.valueOf(nextYear) + "-03-31";
        
        //Calculate Flow statistics for each water-year in time period
        ArrayList<Double> mDay_annualNonZero = new ArrayList<>();
        double mDay_annual = 0;
        boolean moreYears = finalYear > nextYear;
        while(moreYears){
            //Get current water year's data and calculate it's statistics
            String[][] partialData = DoubleArray.getPeriodData(flowData, startDate, endDate);
            if(partialData.length > 0){
                
                //Calculate m-day statistics
                Object[] resultArray = DoubleArray.getMdayData(partialData, M, "arithmetic");
                //ArrayList<String> average_Mday_date = (ArrayList<String>) resultArray[0];
                ArrayList<Double> average_Mday = (ArrayList<Double>) resultArray[1];

                //Calculate minimum flow and save it
                double mDay_min = DoubleMath.min(average_Mday);
                if(mDay_min > 0) mDay_annualNonZero.add(mDay_min);
                mDay_annual++;
            }
            
            //Check if this is the last year
            if(finalYear >= nextYear){
                currentYear++;
                nextYear++;
                startDate = String.valueOf(currentYear) + "-" + waterYearBegin;//startDate = String.valueOf(currentYear) + "-04-01";
                previousDay = DoubleArray.getDay(startDate, -1);
                endDate = String.valueOf(nextYear) + previousDay.substring(4);//endDate = String.valueOf(nextYear) + "-03-31";
            }else{
                moreYears = false;
            }
        }
        
        //Convert low flows to natural-log low-flows
        ArrayList<Double> lnMday_annualNonZero = new ArrayList<>();
        for(int i=0; i<mDay_annualNonZero.size(); i++){
            lnMday_annualNonZero.add( Math.log(mDay_annualNonZero.get(i)) );
        }
        
        //Determine properties of natural-log low-flows
        double U = DoubleMath.meanArithmetic(lnMday_annualNonZero);//Mean
        double S = DoubleMath.StandardDeviationSample(lnMday_annualNonZero);//Standard Deviation
        double G = DoubleMath.SkewnessSample(lnMday_annualNonZero);//Skewness
        
        //Determine fraction of m-day low flows that are zero
        double nonZeroFlows = Double.parseDouble(String.valueOf(mDay_annualNonZero.size()));
        double f0 = (mDay_annual - nonZeroFlows) / mDay_annual;
        
        //Determine cumulative probability of corresponding user-supplied return period
        double P = ((1./R) - f0) / (1 - f0);
        
        //Determine standard normal deviate from cumulative probability (Joiner and Rosenblatt, 1971)
        double Z = 4.91 * (Math.pow(P, 0.14) - Math.pow(1-P, 0.14));
        
        //Compute gamma deviate, K, using Wilson-Hilferty transformation (Loucks et al., 1981)
        double K = (2./G) * (Math.pow(1 + (G*Z/6) - (G*G/36), 3) - 1);
        
        //Compute design flow
        double designFlow = Math.exp(U + K*S);
        
        //Round design flow
        designFlow = DoubleMath.round(designFlow, 2);
        
        return designFlow;
    }
    /**
     * Calculate the CDPHE "biologically-based" design flow (see DFLOW user manual) 
     * which is an 'm'-day harmonic average low flow based on a certain excursion count (exceedance?), 
     * performs this calculation for the entire flow record as well as each month of the year
     * @param flowData  a String[][] of flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @param M  the number of days to average for annual low flow analysis (m-day average)
     * @param R  the desired return period of the m-day low flow (in years)
     * @param clusterLength  the length of time for the definition of a 'cluster' of excursion 
     * periods (default should be 120). All excursion periods within this amount of 
     * time of each other will be counted, subject to the clusterCountMax below.
     * @param clusterCountMax  the upper count limit on how many excursion periods 
     * will be counted within an excursion cluster (default should be 5)
     * @return
     * @throws IOException
     * @throws ParseException 
     */
    public double[] CDPHE_Biological(String[][] flowData,
                                     int M,
                                     int R,
                                     int clusterLength,
                                     int clusterCountMax) throws IOException, ParseException{
        //Calculate annual design low flow
        double DFLOW_annual = CDPHE_Biological_calcs(flowData, M, R, clusterLength, clusterCountMax);
        
        //Calculate monthly design low flows
        double DFLOW_jan = -9999, DFLOW_feb = -9999, DFLOW_mar = -9999, DFLOW_apr = -9999, DFLOW_may = -9999, DFLOW_jun = -9999;
        double DFLOW_jul = -9999, DFLOW_aug = -9999, DFLOW_sep = -9999, DFLOW_oct = -9999, DFLOW_nov = -9999, DFLOW_dec = -9999;
        if(M == 30 && R == 3){
            //If the design flow is 30E3, follow the special provisions set out by reg31
            double[] monthlyDFLOWs = calc30E3(flowData, DFLOW_annual);
            DFLOW_jan = monthlyDFLOWs[0];
            DFLOW_feb = monthlyDFLOWs[1];
            DFLOW_mar = monthlyDFLOWs[2];
            DFLOW_apr = monthlyDFLOWs[3];
            DFLOW_may = monthlyDFLOWs[4];
            DFLOW_jun = monthlyDFLOWs[5];
            DFLOW_jul = monthlyDFLOWs[6];
            DFLOW_aug = monthlyDFLOWs[7];
            DFLOW_sep = monthlyDFLOWs[8];
            DFLOW_oct = monthlyDFLOWs[9];
            DFLOW_nov = monthlyDFLOWs[10];
            DFLOW_dec = monthlyDFLOWs[11];
        }else{
            //Otherwise treat it as a normal biological flow
            DFLOW_jan = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "01-01", "01-31"), M, R, clusterLength, clusterCountMax);
            DFLOW_feb = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "02-01", "02-28"), M, R, clusterLength, clusterCountMax);
            DFLOW_mar = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "03-01", "03-31"), M, R, clusterLength, clusterCountMax);
            DFLOW_apr = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "04-01", "04-30"), M, R, clusterLength, clusterCountMax);
            DFLOW_may = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "05-01", "05-31"), M, R, clusterLength, clusterCountMax);
            DFLOW_jun = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "06-01", "06-30"), M, R, clusterLength, clusterCountMax);
            DFLOW_jul = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "07-01", "07-31"), M, R, clusterLength, clusterCountMax);
            DFLOW_aug = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "08-01", "08-31"), M, R, clusterLength, clusterCountMax);
            DFLOW_sep = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "09-01", "09-30"), M, R, clusterLength, clusterCountMax);
            DFLOW_oct = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "10-01", "10-31"), M, R, clusterLength, clusterCountMax);
            DFLOW_nov = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "11-01", "11-30"), M, R, clusterLength, clusterCountMax);
            DFLOW_dec = CDPHE_Biological_calcs(DoubleArray.getSeasonalData(flowData, "12-01", "12-31"), M, R, clusterLength, clusterCountMax);
        }
        
        //Round design flows
        double[] designFlows = {DFLOW_annual, DFLOW_jan, DFLOW_feb, DFLOW_mar, DFLOW_apr, DFLOW_may, DFLOW_jun, DFLOW_jul, DFLOW_aug, DFLOW_sep, DFLOW_oct, DFLOW_nov, DFLOW_dec};
        DoubleMath.roundColumn(designFlows, 2);
        
        return designFlows;
    }
    /**
     * Calculate the CDPHE "biologically-based" design flow (see DFLOW user manual) 
     * which is an 'm'-day harmonic average low flow based on a certain excursion count (exceedance?)
     * @param flowData  flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @param M  the number of days to average for annual low flow analysis (m-day average)
     * @param R  the desired return period of the m-day low flow (in years)
     * @param clusterLength  the length of time for the definition of a 'cluster' of excursion 
     * periods (default should be 120). All excursion periods within this amount of 
     * time of each other will be counted, subject to the clusterCountMax below.
     * @param clusterCountMax  the upper count limit on how many excursion periods 
     * will be counted within an excursion cluster (default should be 5)
     * @return
     * @throws IOException
     * @throws ParseException 
     */
    private double CDPHE_Biological_calcs(String[][] flowData,
                                          int M,
                                          int R,
                                          int clusterLength,
                                          int clusterCountMax) throws IOException, ParseException{
        //Calculate m-day statistics
        Object[] resultArray = DoubleArray.getMdayData(flowData, M, "harmonic");
        ArrayList<String> average_Mday_date = (ArrayList<String>) resultArray[0];
        ArrayList<Double> average_Mday = (ArrayList<Double>) resultArray[1];
        
        //Get extreme-value based design flow (as a trial design flow to start searching for the biologically based one, note this still uses an arithmetic mean)
        double trialDFLOW = CDPHE_ExtremeValue(flowData, M, R, "04-01");
        if(trialDFLOW == 0){
            //Catch zero-flow guess to avoid ending the convergence loop early and reporting a zero when it shouldn't be zero
            trialDFLOW = DoubleMath.meanHarmonic(average_Mday);
        }
        double trialDFLOW_excursions = countExcursions(average_Mday_date, average_Mday, M, trialDFLOW, clusterLength, clusterCountMax);
        
        //Get allowed number of excursions
        double A = flowData.length/(R * 365.25);
        
        //Loop until the Method of False Position (Carnahan et al., 1969) converges for the limits of the design flow
        double flow_lowerBound = 0, flow_upperBound = trialDFLOW, excursions_lowerBound = 0, excursions_upperBound = trialDFLOW_excursions;
        double designFlow = -1;
        boolean convergance = false;
        int ctr = 0;
        while(!convergance){
            
            //Check for convergance
            if( Math.abs(flow_upperBound - flow_lowerBound) <= (0.005*flow_lowerBound) ){
                designFlow = flow_upperBound;
                convergance = true;
            }else if( Math.abs(A - excursions_lowerBound) <= (0.005*A) ){
                designFlow = flow_lowerBound;
                convergance = true;
            }else if( Math.abs(A - excursions_upperBound) <= (0.005*A) ){
                designFlow = flow_upperBound;
                convergance = true;
            }else if( ctr > 200 ){
                designFlow = -9999;
                convergance = true;
                System.err.println("Did not converge to a solution after " + ctr + " iterations, current upper flow boundary: " + flow_upperBound+ ", current lower flow boundary: " + flow_lowerBound);
            }else{
                //If convergance is not met, interpolate a new trial design flow and count it's excursions
                trialDFLOW = flow_lowerBound + ( ((flow_upperBound - flow_lowerBound) * (A - excursions_lowerBound)) / (excursions_upperBound - excursions_lowerBound) );
                trialDFLOW_excursions = countExcursions(average_Mday_date, average_Mday, M, trialDFLOW, clusterLength, clusterCountMax);
                if(trialDFLOW_excursions <= A){
                    flow_lowerBound = trialDFLOW;
                    excursions_lowerBound = trialDFLOW_excursions;
                }else{
                    flow_upperBound = trialDFLOW;
                    excursions_upperBound = trialDFLOW_excursions;
                }
            }
            ctr++;
        }
        return designFlow;
    }
    /**
     * Counts excursions of m-day average flows below the trialDFLOW based on an 
     * excursion cluster length (excursionClusterLength) and a maximum number of 
     * excursions counted per excursion cluster (clusterCountMax)
     * @param averageMday_date  a list of begin and end dates of the m-day average flows, format: "yyyy-mm-dd to yyyy-mm-dd"
     * @param averageMday  a list of m-day averages corresponding to the dates in the averageMday_date array (aka must be same size)
     * @param m  the m-day averaging number
     * @param trialDFLOW  trial design flow, m-day averages below this will be counted as excursions
     * @param clusterLength  the size of a cluster (default is 120 days)
     * @param clusterCountMax  the maximum number of excursions per cluster to be counted (default is 5)
     * @return 
     */
    private double countExcursions(ArrayList<String> averageMday_date, ArrayList<Double> averageMday, int m, double trialDFLOW, int clusterLength, double clusterCountMax) throws ParseException{
        DateFormat desiredDateFormat = DoubleArray.getDateFormat("daily",  false);
        
        //Determine number of excursions periods
        ArrayList<String> periodList_dates = new ArrayList<>();
        ArrayList<Integer> periodList_lengths = new ArrayList<>();
        for(int i=0; i<averageMday.size(); i++){
            //Check if m-day average flow is below the design flow (aka the start of an excursion period)
            if(averageMday.get(i) < trialDFLOW){
                //Check if this excursion period is a new one or the extension of the previous one
                if(periodList_dates.size() > 0){
                    //Get the date of the last excursion period
                    String lastPeriodDates = periodList_dates.get(periodList_dates.size() - 1);
                    Date lastPeriod_endDate = desiredDateFormat.parse(lastPeriodDates.substring(lastPeriodDates.length() - 10));
                    Date newPeriod_startDate = desiredDateFormat.parse(averageMday_date.get(i).substring(0,10));
                    if(lastPeriod_endDate.after(newPeriod_startDate) || lastPeriod_endDate.equals(newPeriod_startDate)){
                        //If the lasted entered excursion period ends 'after' the new excursion period then this is just a continuation of it
                        //(aka excursion period 1 starts on day 1 and lasts until day 4, excursion period 2 starts on day 2 and lasts until day 5,
                        // therefore the whole thing is just 1 excursion period that lasts from day 1 to day 5)
                        String lastPeriod_startDate = lastPeriodDates.substring(0,10);
                        String newPeriod_endDate = averageMday_date.get(i).substring(averageMday_date.get(i).length() - 10);
                        
                        periodList_dates.set(periodList_dates.size() - 1, lastPeriod_startDate + " to " + newPeriod_endDate);
                        periodList_lengths.set(periodList_dates.size() - 1, getExcursionLength(lastPeriod_startDate, newPeriod_endDate));
                    }else{
                        //This is a new excursion period, so add it to the list
                        periodList_dates.add(averageMday_date.get(i));
                        periodList_lengths.add(m);
                    }
                    
                }else{
                    //If this is the first excursion period, add it to the list
                    periodList_dates.add(averageMday_date.get(i));
                    periodList_lengths.add(m);
                }
            }
        }
        
        //Group the excursions periods into excursion clusters
        ArrayList<String> clusterList_dates = new ArrayList<>();
        ArrayList<Double> clusterList_lengths = new ArrayList<>();
        double m_double = m;
        for(int i=0; i<periodList_dates.size(); i++){
            //Get the dates of the excursion period
            String period_startDate_String = periodList_dates.get(i).substring(0,10);
            Date period_startDate = desiredDateFormat.parse(period_startDate_String);
            
            //Check if this is the first excursion cluster or not
            if(clusterList_dates.size() > 0){
                //Check if this excursion period is within the cluster or not
                String clusterDates = clusterList_dates.get(clusterList_dates.size() - 1);
                Date cluster_endDate = desiredDateFormat.parse(clusterDates.substring(clusterDates.length() - 10));
                if(period_startDate.before(cluster_endDate) || period_startDate.equals(cluster_endDate)){
                    //If the lasted excursion period starts before the end of the current excursion cluster, add it's 
                    //length to the cluster length (up to the maximum cluster count limit)
                    double periodLength = periodList_lengths.get(i);
                    double clusterCount = clusterList_lengths.get(clusterList_lengths.size() - 1) + (periodLength / m_double);
                    if(clusterCount > clusterCountMax){
                        clusterList_lengths.set(clusterList_lengths.size() - 1, clusterCountMax);
                    }else{
                        clusterList_lengths.set(clusterList_lengths.size() - 1, clusterCount);
                    }
                }else{
                    //This is excursion period is in a new excursion cluster, so add it to the list
                    String clusterEndDate = DoubleArray.getDay(period_startDate_String, clusterLength);
                    double periodLength = periodList_lengths.get(i);
                    double clusterCount = periodLength / m_double;
                    if(clusterCount > clusterCountMax){
                        clusterList_lengths.add(clusterCountMax);
                    }else{
                        clusterList_lengths.add(clusterCount);
                    }
                    clusterList_dates.add(period_startDate_String + " to " + clusterEndDate);
                }
            }else{
                //If this is the first excursion cluster, determine it's end date
                String clusterEndDate = DoubleArray.getDay(period_startDate_String, clusterLength);
                double periodLength = periodList_lengths.get(i);
                double clusterCount = periodLength / m_double;
                if(clusterCount > clusterCountMax){
                    clusterList_lengths.add(clusterCountMax);
                }else{
                    clusterList_lengths.add(clusterCount);
                }
                clusterList_dates.add(period_startDate_String + " to " + clusterEndDate);
            }
        }
        
        //Sum the sizes of the excursion clusters
        double excursionCountTotal = DoubleMath.sum(clusterList_lengths);
        return excursionCountTotal;
    }
    /**
     * Counts how many days exist between the provided begin date and end date (startDate cannot equal endDate)
     * @param startDate  the begin date (format yyyy-MM-dd)
     * @param endDate  the end date (format yyyy-MM-dd)
     * @return 
     */
    private int getExcursionLength(String startDate, String endDate) throws ParseException{
        //Determine how many days pass before "nextDay" == "endDay"
        String nextDay = DoubleArray.getDay(startDate, 1);
        int excursionLength = 2;
        while(!nextDay.equals(endDate)){
            nextDay = DoubleArray.getDay(nextDay, 1);
            excursionLength++;
        }
        return excursionLength;
    }
    /**
     * Calculate the 30-day average low flow associated with a 3 year return period 
     * using the special provisions set out in Regulation 31.9.3 "Streams with Rapid Flow Changes"
     * @param flowData  flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @param DFLOW_annual  the biologically based 30E3 design flow for the period of record
     * @return  a list of monthly low flows in order from January to December
     * @throws IOException 
     */
    private double[] calc30E3(String[][] flowData, double DFLOW_annual) throws IOException, ParseException{
        //Calculate m-day statistics
        Object[] resultArray = DoubleArray.getMdayData(flowData, 30, "harmonic");
        ArrayList<String> mDayAverage_dates = (ArrayList<String>) resultArray[0];
        ArrayList<Double> mDayAverage = (ArrayList<Double>) resultArray[1];
        
        //Determine the year and 'month' of each m-day flow
        ArrayList<String> mDayAverage_month = new ArrayList<>();
        for(int i=0; i<mDayAverage.size(); i++){
            //Get the date for 14 days from the begin date, 
            //aka a total of 15 days of data and check if this is still in the same month
            String startDate = mDayAverage_dates.get(i).substring(0,10);
            String midDate = DoubleArray.getDay(startDate, 12);
            
            String beginMonth = startDate.substring(5,7);
            String midMonth = midDate.substring(5,7);
            if(beginMonth.equalsIgnoreCase(midMonth)){
                mDayAverage_month.add(convertMonth(Double.parseDouble(beginMonth)));
            }else{
                mDayAverage_month.add(convertMonth(Double.parseDouble(midMonth)));
            }
        }
        
        //Determine monthly low flows
        String[] months = {"January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"};
        double[] monthlyLowFlow = new double[12];
        for(int j=0; j<12; j++){
            //Get all flows for the current month
            ArrayList<Double> currentMonthFlows = new ArrayList<>();
            for(int i=0; i<mDayAverage.size(); i++){
                if(mDayAverage_month.get(i).equalsIgnoreCase(months[j])){
                    currentMonthFlows.add(mDayAverage.get(i));
                }
            }
            
            //Determine m-day minimum flow for this month
            double mDayAveage_MonthlyMin = DoubleMath.min(currentMonthFlows);
            
            //Check whether to return the annual 30E3 value or the monthly flow (Reg 31.9.3.e)
            if(currentMonthFlows.isEmpty() || DFLOW_annual > mDayAveage_MonthlyMin){
                mDayAveage_MonthlyMin = DFLOW_annual;
            }
            monthlyLowFlow[j] = mDayAveage_MonthlyMin;
        }
        
        return monthlyLowFlow;
    }
    private String convertMonth(double monthDouble){
        int monthInteger = (int) monthDouble;
        String monthString = "?";
        switch(monthInteger){
            case 1:  monthString = "January"; break;
            case 2:  monthString = "February"; break;
            case 3:  monthString = "March"; break;
            case 4:  monthString = "April"; break;
            case 5:  monthString = "May"; break;
            case 6:  monthString = "June"; break;
            case 7:  monthString = "July"; break;
            case 8:  monthString = "August"; break;
            case 9:  monthString = "September"; break;
            case 10: monthString = "October"; break;
            case 11: monthString = "November"; break;
            case 12: monthString = "December"; break;
            default: monthString = "Not a valid month integer"; break;
        }
        return monthString;
    }
    /**
     * Calculate the CDPHE "human health" design flow (see DFLOW user manual) 
     * which is the harmonic mean of the flows
     * @param flowData  flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @return
     * @throws IOException
     * @throws ParseException 
     */
    public double CDPHE_HumanHealth(String[][] flowData) throws IOException, ParseException{
        //Calculate design flows
        double[] flowOnlyData = convertSecondColumn(flowData);
        double designFlow = DoubleMath.meanHarmonic(flowOnlyData);
        
        //Round design flows
        designFlow = DoubleMath.round(designFlow, 2);
        
        return designFlow;
    }
    /**
     * Convert a string[][] array's second column to a double[] array, 
     * aka pull out the flow values from flowData
     * @param flowData  flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @return 
     */
    private double[] convertSecondColumn(String[][] flowData){
        //Pull data
        double[] flowOnlyData = new double[flowData.length];
        for(int i=0; i<flowData.length; i++){
            flowOnlyData[i] = Double.parseDouble(flowData[i][1]);
        }
        
        return flowOnlyData;
    }
    /**
     * 
     * @param flowData  a String[][] of flow data, column1 = dates (format yyyy-mm-dd), column2 = flow values
     * @return
     * @throws IOException
     * @throws ParseException 
     */
    public String CDPHE_REG31(String[][] flowData) throws IOException, ParseException{
        ArrayList<Integer> monthlyRowIndex = new ArrayList<>();
        monthlyRowIndex.add(1);
        monthlyRowIndex.add(2);
        monthlyRowIndex.add(3);
        monthlyRowIndex.add(4);
        monthlyRowIndex.add(5);
        monthlyRowIndex.add(6);
        monthlyRowIndex.add(7);
        monthlyRowIndex.add(8);
        monthlyRowIndex.add(9);
        monthlyRowIndex.add(10);
        monthlyRowIndex.add(11);
        monthlyRowIndex.add(12);
        
        //Calculate acute biologically based low flows (1-day 3-year)
        System.out.println("Calcualting Reg. 31 1E3 Acute Monthly Low Flows...");
        double[] monthlyBiological_1day = CDPHE_Biological(flowData, 1, 3, 120, 5);
        double monlyBiological_1day_min = DoubleMath.min(DoubleArray.getRows(monthlyBiological_1day, monthlyRowIndex));
        
        //Calculate chronic biologically based low flows (7-day 3-year)
        System.out.println("Calcualting Reg. 31 7E3 Chronic Monthly Low Flows...");
        double[] monthlyBiological_7day = CDPHE_Biological(flowData, 7, 3, 120, 5);
        double monlyBiological_7day_min = DoubleMath.min(DoubleArray.getRows(monthlyBiological_7day, monthlyRowIndex));
        
        //Calculate chronic biologically based low flows (30-day 3-year)
        System.out.println("Calcualting Reg. 31 30E3 Chronic Monthly Low Flows...");
        double[] monthlyBiological_30day = CDPHE_Biological(flowData, 30, 3, 120, 5);
        double monlyBiological_30day_min = DoubleMath.min(DoubleArray.getRows(monthlyBiological_30day, monthlyRowIndex));
        
        //Calculate annual median of data with a 5-year return period
        System.out.println("Calcualting Reg. 31 1E5 Annual Median of Daily Average Flows...");
        double median_R_yearFlow = annualMedianReturnPeriod(flowData, 5);
        median_R_yearFlow = DoubleMath.round(median_R_yearFlow, 2);
        
        //Format the results into a summary table
        String summaryTable = "Month\t1E3 Acute Monthly Low Flows [cfs]\t7E3 Chronic Monthly Low Flows [cfs]\t30E3 Chronic Monthly Low Flows [cfs]";
        String[] months = {"Entire Record", "January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December", "Lowest Monthly Flow"};
        for(int i=0; i<13; i++){
            summaryTable = summaryTable + "\r\n" + months[i] + "\t" + monthlyBiological_1day[i] + "\t" + monthlyBiological_7day[i] + "\t" + monthlyBiological_30day[i];
        }
        summaryTable = summaryTable + "\r\n" + months[13] + "\t" + monlyBiological_1day_min + "\t" + monlyBiological_7day_min + "\t" + monlyBiological_30day_min;
        summaryTable = summaryTable + "\r\n" + median_R_yearFlow;
        
        return summaryTable;
    }
    public double annualMedianReturnPeriod(String[][] flowData, int R){
        //Calculate the annual median of daily average flow values with a recurrance intervale of 1 in "R" years

        //Calculate median for each year in time period
        ArrayList<Double> medianFlows = new ArrayList<>();
        boolean moreYears = flowData.length > 0;
        String currentYear = flowData[0][0].substring(0,4);
        String finalYear = flowData[flowData.length - 1][0].substring(0,4);
        while(moreYears){
            //Get current year's data and calculate it's statistics
            String[][] partialData = DoubleArray.getYearsData(flowData, currentYear);
            if(partialData.length > 0){
                double[] flowOnlyData = convertSecondColumn(partialData);
                medianFlows.add(DoubleMath.median(flowOnlyData));
            }
            
            int nextYear = Integer.parseInt(currentYear) + 1;
            if(finalYear.compareToIgnoreCase(String.valueOf(nextYear)) >= 0){
                currentYear = String.valueOf(nextYear);
            }else{
                moreYears = false;
            }
        }
        
        //Calculate return period of data
        double median_RyearFlow = DoubleArray.calculateLowFlowReturnPeriod(medianFlows, R);
        return median_RyearFlow;
    }
}