From 757489abd6ef13a6c3549ee9854a476834b5165d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tobias=20M=C3=BCller?= <to.mueller@gmail.com>
Date: Wed, 27 Mar 2024 12:16:16 -0700
Subject: [PATCH] feat: additional statistical parameters (#109)

* added version tag to distinguish test and prod deploy updates

* fixed pdf display

* removed debug tag, removed outdated todo comments

* fixed header size

* fix pandas slice warning

* added calculation of mean average percentage error, fixed stats update calc

* added mape to pdf export
---
 frontend/rctool/functions/fit_linear_model.py |  8 +-
 .../templates/rctool/components/navbar.html   |  6 +-
 .../rctool/develop/develop_table_rating.html  |  6 ++
 .../rctool/develop/rctool_develop_js.html     | 81 ++++++++++++-------
 .../rctool/export/rctool_export_pdf.html      |  3 +
 frontend/rctool/views.py                      | 12 +--
 6 files changed, 76 insertions(+), 40 deletions(-)

diff --git a/frontend/rctool/functions/fit_linear_model.py b/frontend/rctool/functions/fit_linear_model.py
index ad13185b..90677d10 100644
--- a/frontend/rctool/functions/fit_linear_model.py
+++ b/frontend/rctool/functions/fit_linear_model.py
@@ -1,7 +1,7 @@
 import pandas as pd
 import numpy as np
 from lmfit import models, Parameters
-from sklearn.metrics import mean_squared_error
+from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
 import math
 
 
@@ -71,9 +71,11 @@ def fit_linear_model(df, offset, label, weighted=None, intersect_points=None, *a
     # calculate statistical parameters to analyze goodness of fit
     unw_mse = mean_squared_error(df_data["Q"], unw_best)
     unw_rmse = math.sqrt(unw_mse)
+    unw_mape = mean_absolute_percentage_error(df_data["Q"], unw_best) * 100.0
 
     wgt_mse = mean_squared_error(df_data["Q"], wgt_best)
     wgt_rmse = math.sqrt(wgt_mse)
+    wgt_mape = mean_absolute_percentage_error(df_data["Q"], wgt_best) * 100.0
 
     # Process and ship output
     mdl_param = {
@@ -84,6 +86,7 @@ def fit_linear_model(df, offset, label, weighted=None, intersect_points=None, *a
             "seg_bounds": unw_seg_nodes,
             "offset": offset,
             "rmse": unw_rmse,
+            "mape": unw_mape,
         },
         "wgt": {
             "label": label,
@@ -92,6 +95,7 @@ def fit_linear_model(df, offset, label, weighted=None, intersect_points=None, *a
             "seg_bounds": wgt_seg_nodes,
             "offset": offset,
             "rmse": wgt_rmse,
+            "mape": wgt_mape,
         },
     }
 
@@ -103,6 +107,7 @@ def fit_linear_model(df, offset, label, weighted=None, intersect_points=None, *a
             "seg_bounds": wgt_seg_nodes,
             "offset": offset,
             "rmse": wgt_rmse,
+            "mape": wgt_mape,
         }
         wgt_data = [
             [a, b, c] for a, b, c in zip(df_data["H"].tolist(), wgt_best, wgt_residual)
@@ -117,6 +122,7 @@ def fit_linear_model(df, offset, label, weighted=None, intersect_points=None, *a
             "seg_bounds": unw_seg_nodes,
             "offset": offset,
             "rmse": unw_rmse,
+            "mape": unw_mape,
         }
         unw_data = [
             [a, b, c] for a, b, c in zip(df_data["H"].tolist(), unw_best, unw_residual)
diff --git a/frontend/rctool/templates/rctool/components/navbar.html b/frontend/rctool/templates/rctool/components/navbar.html
index 1e9c776e..f68820d6 100644
--- a/frontend/rctool/templates/rctool/components/navbar.html
+++ b/frontend/rctool/templates/rctool/components/navbar.html
@@ -4,11 +4,11 @@
     <div class="header-title-div justify-items-start"></div>
       <a class="navbar-brand" href="{% url 'home' %}"><img src="{% static 'images/gov_logo.JPG' %}" class="logo-gov"></a>
       <div class="header-title-div">
-        <h1 class="header-title" 
-          style="color: white; font-size: 1.5em; font-weight: bold; font-family:sans-serif; margin-top: 0.5em; margin-bottom: 0.5em; margin-left: 1em; margin-right: 1em; text-align: center; text-shadow:1px 1px 10px #000, 1px 1px 10px #000;
+        <p class="header-title" 
+          style="color: white; font-size: 1.5em; font-weight: bold; font-family:sans-serif; margin-top: 0em; margin-bottom: 0em; margin-left: 1em; margin-right: 1em; text-align: center; text-shadow:1px 1px 10px #000, 1px 1px 10px #000;
           ">
             HydRA
-          </h1>
+        </p>
       </div>
     </div>
     <button class="navbar-toggler" type="button" data-toggle="collapse" data-target="#navbarNav" aria-controls="navbarNav" aria-expanded="false" aria-label="Toggle navigation">
diff --git a/frontend/rctool/templates/rctool/rctool/develop/develop_table_rating.html b/frontend/rctool/templates/rctool/rctool/develop/develop_table_rating.html
index 6640bd66..3c1d42b2 100644
--- a/frontend/rctool/templates/rctool/rctool/develop/develop_table_rating.html
+++ b/frontend/rctool/templates/rctool/rctool/develop/develop_table_rating.html
@@ -12,6 +12,7 @@
                 <th class="table-cell">start point<br>(stage in m)</th>
                 <th class="table-cell">end point<br>(stage in m)</th>
                 <th class="table-cell">RMSE</th>
+                <th class="table-cell">MAPE</th>
                 <th class="table-cell"><p id="comparison-header", style="color:#9ec1a3">compare</p></th>
             </tr>
             {% if rc.parameters %}
@@ -23,6 +24,7 @@
                     <td class="table-cell"><p id="table1-endpoint-Seg1H0">{{ rc.parameters.0.seg_bounds.0.0|floatformat:3 }}</p></td>
                     <td class="table-cell"><p id="table1-endpoint-Seg1H1">{{ rc.parameters.0.seg_bounds.1.0|floatformat:3 }}</p></td>
                     <td class="table-cell"><p id="table1-segment1-rmse">{{ rc.parameters.0.rmse|floatformat:3 }}</p></td>
+                    <td class="table-cell"><p id="table1-segment1-mape">{{ rc.parameters.0.mape|floatformat:1 }} %</p></td>
                     <td class="table-cell"><p></p></td>
                 </tr>
             {% endif %}
@@ -35,6 +37,7 @@
                     <td class="table-cell"><p id="table1-endpoint-Seg2H0">{{ rc.parameters.1.seg_bounds.0.0|floatformat:3 }}</p></td>
                     <td class="table-cell"><p id="table1-endpoint-Seg2H1">{{ rc.parameters.1.seg_bounds.1.0|floatformat:3 }}</p></td>
                     <td class="table-cell"><p id="table1-segment2-rmse">{{ rc.parameters.1.rmse|floatformat:3 }}</p></td>
+                    <td class="table-cell"><p id="table1-segment2-mape">{{ rc.parameters.1.mape|floatformat:1 }} %</p></td>
                     <td class="table-cell"><p></p></td>
                 </tr>
             {% endif %}
@@ -60,6 +63,9 @@
                 <td class="table-cell">
                     <p id="compare-curve-rmse" style="width: 45px;"></p>
                 </td>
+                <td class="table-cell">
+                    <p id="compare-curve-mape" style="width: 45px;"></p>
+                </td>
                 <td class="table-cell" style="display:block">
                     <div id="add-curve">
                         <button 
diff --git a/frontend/rctool/templates/rctool/rctool/develop/rctool_develop_js.html b/frontend/rctool/templates/rctool/rctool/develop/rctool_develop_js.html
index b247e6af..7fb0fdb3 100644
--- a/frontend/rctool/templates/rctool/rctool/develop/rctool_develop_js.html
+++ b/frontend/rctool/templates/rctool/rctool/develop/rctool_develop_js.html
@@ -143,6 +143,26 @@
         }
     }
 
+    function calculate_stats(q_field, q_model) {
+
+            // Calculate statistical parameters
+            var sumCumulativeDifferencePerc = 0;
+            var sumSquaredCumulativeDifference = 0;
+            q_model.forEach((item, index) => {
+                sumSquaredCumulativeDifference += (item - q_field[index]) ** 2
+                sumCumulativeDifferencePerc += Math.abs(item - q_field[index]) / q_field[index]
+            });
+
+
+            // Calculate root mean squared error (rmse)
+            rmse = Math.sqrt( (sumSquaredCumulativeDifference/ q_model.length) );
+
+            // Calculate mean absolute percentage error (mape)
+            mape = (sumCumulativeDifferencePerc / q_model.length) * 100;
+
+            return {'rmse': rmse, 'mape': mape}
+    }
+
     function adjustSeg(segmentBounds, offsets) {
         // INPUTS:
             // segmentBounds: dictionary where keys are segmentName and values are segmentData
@@ -177,25 +197,16 @@
 
             // 3.2 Calculate model data from new rc equation, collect and prepaire output
             var stage = fieldDataFiltered.map( item => (item[0]));
-            var oldQ = fieldDataFiltered.map( item => (item[1]));
-            var newQ = stage.map(s => C * (s - offset) ** slope);
+            var q_field = fieldDataFiltered.map( item => (item[1]));
+            var q_model = stage.map(s => C * (s - offset) ** slope);
             // Calculate statistical parameters; residuals
             var newResiduals = fieldDataFiltered.map(item => ( 100 * (C * (item[0] - offset) ** slope - item[1]) / item[1] ) );
 
             // Collect and prepaire output
             var newSegData = [];
-            for (let t = 0; t < newQ.length; t++) {
-                newSegData.push([stage[t], newQ[t], newResiduals[t]]);
-            };
-
-            // Calculate statistical parameters; RMSE
-            var sumSquaredCumulative = 0;
-            for (var v = 0; v < newQ.length; v++) {
-                var currSumSquared = (newQ[v] - oldQ[v]) ** 2;
-                sumSquaredCumulative = sumSquaredCumulative + currSumSquared;
+            for (let t = 0; t < q_model.length; t++) {
+                newSegData.push([stage[t], q_model[t], newResiduals[t]]);
             };
-            rmse = Math.sqrt( (sumSquaredCumulative/ newQ.length) );
-
 
             // sort list by stage (index)
             newSegData.sort(function(a,b) {
@@ -207,12 +218,15 @@
             // add boundary points to segment data (for plotting on rcChart)
             newSegData.unshift([segmentData[0][0], segmentData[0][1],  0]);
             newSegData.push([segmentData[1][0], segmentData[1][1], 0]);
+
+            stats = calculate_stats(q_field, q_model)
             // 3.4 Update rcData & rcParam variables
             rcData[idx+1].data = newSegData;
             rcParam[idx].const = C;
             rcParam[idx].exp = slope;
             rcParam[idx].offset = offset;
-            rcParam[idx].rmse = rmse;
+            rcParam[idx].rmse = stats['rmse'];
+            rcParam[idx].mape = stats['mape'];;
             rcParam[idx].seg_bounds = [[segmentData[0][0], segmentData[0][1]], [segmentData[1][0], segmentData[1][1]]];
             rcDict = {'data':rcData, 'parameters':rcParam};
 
@@ -230,6 +244,7 @@
             document.getElementById('table1-segment' + segNumStr + '-const').innerHTML = C.toFixed(2);
             document.getElementById('table1-segment' + segNumStr + '-exp').innerHTML = slope.toFixed(2);
             document.getElementById('table1-segment' + segNumStr + '-rmse').innerHTML = rmse.toFixed(3);
+            document.getElementById('table1-segment' + segNumStr + '-mape').innerHTML = mape.toFixed(1) + '%';
             document.getElementById("endpointSeg1H0").value = rcParam[0].seg_bounds[0][0].toFixed(3);
             document.getElementById("endpointSeg1Q0").value = rcParam[0].seg_bounds[0][1].toFixed(3);
             document.getElementById("endpointSeg1H1").value = rcParam[0].seg_bounds[1][0].toFixed(3);
@@ -282,11 +297,6 @@
 
         // check if exponent in realistic range
         checkExponent();
-
-        // update status
-        // displayStatus('<strong>status: </strong> manual adjustment applied...')
-        {% comment %} toggleAxisFormat(document.getElementById("toggle_axis_format").checked) {% endcomment %}
-
     };
 
     function updateChartFromGlobal(){
@@ -482,6 +492,19 @@
 
         // update residual chart
         for (var seg = 0; seg < rcParam.length; seg++) {
+            
+            // get field data for this segment
+            const fieldDataActive_seg = fieldDataActive.filter(item => ( item.y >= rcParam[seg].seg_bounds[0][0] && item.y <= rcParam[seg].seg_bounds[1][0]));
+            q_field_active_seg = fieldDataActive_seg.map(item => item.x);
+            h_field_active_seg = fieldDataActive_seg.map(item => item.y);
+
+            // calculate model data for this segment
+            q_model_active_seg = h_field_active_seg.map(s => calc_q_model(s, rcParam[seg].const, rcParam[seg].exp, rcParam[seg].offset));
+            stats_seg = calculate_stats(q_field_active_seg, q_model_active_seg)
+            // update HTML table with new stats
+            document.getElementById('table1-segment' + (seg + 1) + '-rmse').innerHTML = stats_seg['rmse'].toFixed(3);
+            document.getElementById('table1-segment' + (seg + 1) + '-mape').innerHTML = stats_seg['mape'].toFixed(1) + '%';
+
             var dataArray_residual = Object.keys(fielddatacsv.toggle_point).map((id, index) => {
                 return {
                     id: id,
@@ -966,14 +989,14 @@
         var fieldDataFiltered = rawFieldData.filter(item => ( item[0] >= compStartBounds && item[0] <= compEndBounds));
         // Calculate model data from new rc equation, collect and prepaire output
         var stage = fieldDataFiltered.map( item => (item[0]));
-        var fieldQ = fieldDataFiltered.map( item => (item[1]));
-        var compQ = stage.map(s => compConst * (s - compOffset) ** compExp);
+        var q_field_compare = fieldDataFiltered.map( item => (item[1]));
+        var q_model_compare = stage.map(s => compConst * (s - compOffset) ** compExp);
         // Calculate statistical parameters; residuals
         var compResiduals = fieldDataFiltered.map(item => ( 100 * (compConst * (item[0] - compOffset) ** compExp - item[1]) / item[1] ) );
         // Collect and prepaire output
         var compData = [];
-        for (let t = 0; t < compQ.length; t++) {
-            compData.push([stage[t], compQ[t], compResiduals[t]]);
+        for (let t = 0; t < q_model_compare.length; t++) {
+            compData.push([stage[t], q_model_compare[t], compResiduals[t]]);
         };
         // Add endpoints for RC chart (not for residual chart)
         lowerH = compStartBounds
@@ -985,12 +1008,8 @@
         compRcData.unshift([lowerH, lowerQ, 0])
         compRcData.push([upperH, upperQ, 0])
         // 3. Calculate statistical parameters; RMSE
-        var sumSquaredCumulative = 0;
-        for (var v = 0; v < compQ.length; v++) {
-            var currSumSquared = (compQ[v] - fieldQ[v]) ** 2;
-            sumSquaredCumulative = sumSquaredCumulative + currSumSquared;
-        };
-        compRmse = Math.sqrt((sumSquaredCumulative/ compQ.length)).toFixed(2);
+        stats_compare = calculate_stats(q_field_compare, q_model_compare)
+
         // 4. Prepaire output
         var compChartData = compRcData.map(item => ({ x: item[1], y: item[0] }));
         var compResChartData = compData.map(item => ({ x: item[2], y: item[0] }));
@@ -1001,7 +1020,8 @@
             addResidualPlot(compResChartData, compLabel, '#DDCC77')
         };
         // 6. Update table with new values and remove option to add new curve
-        document.getElementById('compare-curve-rmse').innerHTML = compRmse;
+        document.getElementById('compare-curve-rmse').innerHTML = stats_compare['rmse'].toFixed(3)
+        document.getElementById('compare-curve-mape').innerHTML = stats_compare['mape'].toFixed(1) + '%'
         document.getElementById('add-curve').style.display = 'none';
         document.getElementById('remove-curve').style.display = 'block';
         document.getElementById('comparison-header').style.color = '#000000';
@@ -1021,6 +1041,7 @@
         // document.getElementById('compareStartBounds').value = ''
         // document.getElementById('compareEndBounds').value = ''
         document.getElementById('compare-curve-rmse').innerHTML = '';
+        document.getElementById('compare-curve-mape').innerHTML = '';
         // 4. hide 'remove curve' button and display 'add curve button'
         document.getElementById('remove-curve').style.display = 'none';
         document.getElementById('add-curve').style.display = 'block';
diff --git a/frontend/rctool/templates/rctool/rctool/export/rctool_export_pdf.html b/frontend/rctool/templates/rctool/rctool/export/rctool_export_pdf.html
index b8a02b84..6fd4385c 100644
--- a/frontend/rctool/templates/rctool/rctool/export/rctool_export_pdf.html
+++ b/frontend/rctool/templates/rctool/rctool/export/rctool_export_pdf.html
@@ -143,12 +143,14 @@ <h6 class="pdf-section-heading-text">Model Parameters</h6>
                 <th class="table-header">equation</th>
                 <th class="table-header">stage range (m)</th>
                 <th class="table-header">RMSE</th>
+                <th class="table-header">MAPE</th>
             </tr>
             <tr class="table-row" style="height:25px;">
                 <td class="table-cell">{{ rc.parameters.0.label }}</td>
                 <td class="table-cell">Q = {{ rc.parameters.0.const }} ( H {{ offsets_val.0 }} )<sup>{{ rc.parameters.0.exp }}</sup></td>
                 <td class="table-cell">{{ rc.parameters.0.seg_bounds.0.0|floatformat:3 }} , {{ rc.parameters.0.seg_bounds.1.0|floatformat:3}}</td>
                 <td class="table-cell">{{ rc.parameters.0.rmse|floatformat:3 }}</td>
+                <td class="table-cell">{{ rc.parameters.0.mape|floatformat:1 }} %</td>
                 </td>
             </tr>
             {% if rc.parameters.1 %}
@@ -157,6 +159,7 @@ <h6 class="pdf-section-heading-text">Model Parameters</h6>
                 <td class="table-cell">Q = {{ rc.parameters.1.const }}( H {{ offsets_val.1 }})<sup>{{ rc.parameters.1.exp }}</sup></td>
                 <td class="table-cell">{{ rc.parameters.1.seg_bounds.0.0|floatformat:3 }} , {{ rc.parameters.1.seg_bounds.1.0|floatformat:3 }}</td>
                 <td class="table-cell">{{ rc.parameters.1.rmse|floatformat:3 }}</td>
+                <td class="table-cell">{{ rc.parameters.1.mape|floatformat:1 }} %</td>
                 </td>
             </tr>
             {% endif %}
diff --git a/frontend/rctool/views.py b/frontend/rctool/views.py
index 0d69162a..431b458e 100644
--- a/frontend/rctool/views.py
+++ b/frontend/rctool/views.py
@@ -891,8 +891,8 @@ def calc_error(const, exp, offset, stage, discharge):
         result = -100 * ((const * (stage - offset) ** exp) - discharge) / discharge
         return np.array(result.tolist()).flatten()
 
-    df_error1 = field_data_output_df[field_data_output_df["stage"] <= breakpoint]
-    df_error2 = field_data_output_df[field_data_output_df["stage"] > breakpoint]
+    df_error1 = field_data_output_df[field_data_output_df["stage"] <= breakpoint].copy()
+    df_error2 = field_data_output_df[field_data_output_df["stage"] > breakpoint].copy()
 
     err_result1 = calc_error(
         segment_parameters[0]["const"],
@@ -915,15 +915,15 @@ def calc_error(const, exp, offset, stage, discharge):
         err_result2 = np.zeros(len(df_error2))
 
     # add list as column to df
-    df_error1["Discharge Error (%)"] = err_result1
-    df_error2["Discharge Error (%)"] = err_result2
+    df_error1.loc[:, "Discharge Error (%)"] = err_result1
+    df_error2.loc[:, "Discharge Error (%)"] = err_result2
 
     # merge dfs
     field_data_output_df = pd.concat([df_error1, df_error2])
     # round discharge error to 2 decimals
-    field_data_output_df["Discharge Error (%)"] = field_data_output_df[
+    field_data_output_df.loc[:, "Discharge Error (%)"] = field_data_output_df[
         "Discharge Error (%)"
-    ].round(decimals=2)
+    ].round(2)
     return field_data_output_df