Merge pull request #44 from wseis/main

evaluation

.gitignore

@@ -127,4 +127,5 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
-qmra.db
+qmra.db
+notebooks

tools/qmratool/helper_functions.py

@@ -1,6 +1,6 @@
 import plotly.express as px
 import plotly.graph_objs as go
-
+from django_pandas.io import read_frame
 
 def plot_comparison(df, df2):
     fig = px.box(
@@ -82,3 +82,120 @@ def plot_comparison(df, df2):
         marker_size=8, hovertemplate=None, hoverinfo="skip", line=dict(width=0)
     )
     return fig
+
+
+
+def simulate_risk(ra):
+    # Selecting inflow concentration based in source water type
+    df_inflow = read_frame(
+        Inflow.objects.filter(water_source=ra.source).values(
+            "min", "max", "pathogen__pathogen", "water_source__water_source_name"
+        )
+    )
+    df_inflow = df_inflow[
+        df_inflow["pathogen__pathogen"].isin(
+            ["Rotavirus", "Campylobacter jejuni", "Cryptosporidium parvum"]
+        )
+    ]
+    # Querying dose response parameters based on pathogen inflow
+    dr_models = read_frame(
+        DoseResponse.objects.filter(
+            pathogen__in=Pathogen.objects.filter(
+                pathogen__in=df_inflow["pathogen__pathogen"]
+            )
+        )
+    )
+
+    # Querying for Logremoval based on selected treatments
+    df_treatment = read_frame(
+        LogRemoval.objects.filter(treatment__in=ra.treatment.all()).values(
+            "min", "max", "treatment__name", "pathogen_group__pathogen_group"
+        )
+    )
+
+    # Summarizing treatment to treatment max and treatment min
+    # df_treatment_summary=df_treatment.groupby(["pathogen_group__pathogen_group"]).sum().reset_index()
+    df_treatment_summary = (
+        df_treatment.groupby(["pathogen_group__pathogen_group"])[["min", "max"]]
+        .apply(lambda x: x.sum())
+        .reset_index()
+    )
+
+    results = pd.DataFrame()
+
+    for index, row in df_inflow.iterrows():
+        d = df_inflow.loc[df_inflow["pathogen__pathogen"] == row["pathogen__pathogen"]]
+        dr = dr_models.loc[dr_models["pathogen"] == row["pathogen__pathogen"]]
+
+        if row["pathogen__pathogen"] == "Rotavirus":
+            selector = "Viruses"
+        elif row["pathogen__pathogen"] == "Cryptosporidium parvum":
+            selector = "Protozoa"
+        else:
+            selector = "Bacteria"
+        # result.append(selector)
+
+        df_treat = df_treatment_summary[
+            df_treatment_summary["pathogen_group__pathogen_group"] == selector
+        ]
+
+        risk_df = pd.DataFrame(
+            {
+                "inflow": np.random.normal(
+                    loc=(
+                        np.log10(float(d["min"]) + 10 ** (-8))
+                        + np.log10(float(d["max"]))
+                    )
+                    / 2,
+                    scale=(
+                        np.log10(float(d["max"]))
+                        - np.log10(float(d["min"]) + 10 ** (-8))
+                    )
+                    / 4,
+                    size=10000,
+                ),
+                "LRV": np.random.uniform(
+                    low=df_treat["min"], high=df_treat["min"], size=10000
+                ),
+                "LRVmax": np.random.uniform(
+                    low=df_treat["max"], high=df_treat["max"], size=10000
+                ),
+            }
+        )
+        risk_df["outflow"] = risk_df["inflow"] - risk_df["LRV"]
+        risk_df["outflow_min"] = risk_df["inflow"] - risk_df["LRVmax"]
+
+        risk_df["dose"] = (10 ** risk_df["outflow"]) * float(
+            ra.exposure.volume_per_event
+        )
+        risk_df["dose_min"] = (10 ** risk_df["outflow_min"]) * float(
+            ra.exposure.volume_per_event
+        )
+
+        if selector != "Protozoa":
+            risk_df["probs"] = 1 - (
+                1 + (risk_df["dose"]) * (2 ** (1 / float(dr.alpha)) - 1) / float(dr.n50)
+            ) ** -float(dr.alpha)
+            risk_df["probs_min"] = 1 - (
+                1
+                + (risk_df["dose_min"])
+                * (2 ** (1 / float(dr.alpha)) - 1)
+                / float(dr.n50)
+            ) ** -float(dr.alpha)
+
+        else:
+            risk_df["probs"] = 1 - np.exp(-float(dr.k) * (risk_df["dose"]))
+            risk_df["probs_min"] = 1 - np.exp(-float(dr.k) * (risk_df["dose_min"]))
+
+        results[row["pathogen__pathogen"] + "_MinimumLRV"] = [
+            annual_risk(int(ra.exposure.events_per_year), risk_df["probs"])
+            for _ in range(1000)
+        ]
+        results[row["pathogen__pathogen"] + "_MaximumLRV"] = [
+            annual_risk(int(ra.exposure.events_per_year), risk_df["probs_min"])
+            for _ in range(1000)
+        ]
+
+    results_long = pd.melt(results)
+    results_long["log probability"] = np.log10(results_long["value"])
+    return results_long

tools/qmratool/templates/qmratool/results.html

@@ -9,6 +9,15 @@
           <a class = "btn btn-primary mt-3" href="{%url 'export' ra.id%}">Export summary as csv</a>
         </div>
         {%endif %}
+        {% if risk%}
+        <div class="alert alert-danger m-4" role="alert">
+          At least one pathogen exeeds tolerable infection risk value of 1/10000 pppy (red line)
+        </div>
+        {%else %}
+        <div class="alert alert-primary m-4" role="alert">
+         All risk simulations stay well below tolerable limits of 1/10000 pppy (blue line)
+        </div>
+        {% endif %}
         <div id="hero-graph">
           {% autoescape off %}
           {{ risk_plot }}
@@ -22,7 +31,16 @@
            </p> 
         </div>
         {%endif%}
-        {%if ra.id %}   
+        {%if ra.id %}
+        {% if dalyrisk%}
+        <div class="alert alert-danger m-4" role="alert">
+          At least one pathogen exeeds tolerable DALY risk value of 1µDALY pppy (red line)
+        </div>
+        {%else %}
+        <div class="alert alert-primary m-4" role="alert">
+         All risk simulations stay well below tolerable limits of 1µDALY pppy (blue line)
+        </div>
+        {% endif %}   
         <div id="hero-graph">
           {% autoescape off %}
           {{ daly_plot }}

tools/qmratool/views.py

@@ -5,24 +5,30 @@
 from django.shortcuts import render
 from django.db import IntegrityError
 from django.contrib.auth.decorators import login_required
-
 from qmratool.helper_functions import plot_comparison
+from django.urls import reverse_lazy
+from django.shortcuts import render, redirect
+from .forms import *
+from formtools.wizard.views import SessionWizardView
 from .forms import RAForm, SourceWaterForm, TreatmentForm, ExposureForm
 from .forms import LogRemovalForm, InflowForm, ComparisonForm
-
 from .models import *
 from django.utils.encoding import force_str
-
 django.utils.encoding.force_text = force_str
-
 from django_pandas.io import read_frame
 from plotly.offline import plot
-
 import numpy as np
 import pandas as pd
 import plotly.express as px
 import plotly.graph_objs as go
 import markdown2 as md
+from django.contrib import messages
+from django.contrib.auth import update_session_auth_hash
+from django.contrib.auth.forms import PasswordChangeForm
+from django.shortcuts import render, redirect
+
+
+
 
 # Create your views here.
 
@@ -167,11 +173,6 @@ def new_assessment(request):
     return render(request, "qmratool/new_ra.html", {"form": form, "content": content})
 
 
-from django.urls import reverse_lazy
-from django.shortcuts import render, redirect
-from .forms import *
-
-from formtools.wizard.views import SessionWizardView
 
 TEMPLATES = {
     "0": "qmratool/step1.html",
@@ -431,17 +432,40 @@ def calculate_risk(request, ra_id):
     )
 
     results_long = simulate_risk(ra)
+    #results_long["value"] = round(results_long["value"],1)
 
     results_long["pathogen"] = results_long["variable"].str.split("_", expand=True)[0]
     results_long["stat"] = results_long["variable"].str.split("_", expand=True)[1]
 
+
+
+
+
     health = read_frame(Health.objects.all())
     results_long = pd.merge(results_long, health, on="pathogen")
     results_long["DALYs pppy"] = (
         results_long.value
         * results_long.infection_to_illness.astype(float)
         * results_long.dalys_per_case.astype(float)
     )
+    summary = results_long.groupby(["pathogen", "stat"]).mean() 
+    risk = sum(summary["value"]>10**-4) > 0
+    dalyrisk= sum(summary["DALYs pppy"]>10**-6) > 0
+
+    if risk:
+        bgcolor = "rgba(245, 0, 0, 0.15)"
+        lcolor = "firebrick"
+    else:
+        bgcolor = None
+        lcolor = "#0003e2"
+
+    if dalyrisk:
+        dalybgcolor = "rgba(245, 0, 0, 0.15)"
+        dlcolor = "firebrick"
+    else:
+        dalybgcolor = None
+        dlcolor = "#0003e2"
+
 
     risk_colors = ["#78BEF9", "#8081F1", "#5F60B3"]
     risk_colors_extended = [
@@ -463,37 +487,40 @@ def calculate_risk(request, ra_id):
         log_y=True,
         title="Risk as probability of infection per year",
         color_discrete_sequence=risk_colors,
+        hover_data={'value': ':.2e'}
     )
 
     fig.update_layout(
         font_family="Helvetica Neue, Helvetica, Arial, sans-serif",
         font_color="black",
+        plot_bgcolor= bgcolor,
         title={"text": "Risk assessment as probability of infection per year"},
         xaxis_title="",
         yaxis_title="Probability of infection per year",
         # markersize= 12,
     )
-    fig.add_hline(y=0.0001, line_dash="dashdot", line=dict(color="#0003e2", width=3))
+    fig.add_hline(y=0.0001, line_dash="dashdot", line=dict(color=lcolor, width=3))
     fig.update_layout(
+       # margin=dict(l=50, r=50, t=100, b=100),  # Adjust the margins
         legend=dict(
             orientation="h",
             yanchor="top",
             # text= "Reference pathogen",
-            y=-0.1,
+            y=-0.3,
             xanchor="left",
             x=0,
         ),
-        annotations=[
-            go.Annotation(
-                y=-4,
-                x=1.2,
-                text="Tolerable risk level of 1/10000 infections pppy",
-                bgcolor="#0003e2",
-                bordercolor="white",
-                borderpad=5,
-                font=dict(color="white"),
-            )
-        ],
+      #  annotations=[
+      #      go.Annotation(
+      #          y=-4,
+      #          x=1.2,
+      #          text="Tolerable risk level of 1/10000 infections pppy",
+      #          bgcolor="#0003e2",
+      #          bordercolor="white",
+      #          borderpad=5,
+      #          font=dict(color="white"),
+      #      )
+      #  ],
     )
 
     # fig.update_annotations(y = 0.0001, x = 1,  text = "Tolerable risk level")
@@ -522,17 +549,18 @@ def calculate_risk(request, ra_id):
         },
         xaxis_title="",
         yaxis_title="DALYs pppy",
-        annotations=[
-            go.Annotation(
-                y=-6,
-                x=1.2,
-                text="Tolerable risk level of 1µDALY pppy",
-                bgcolor="#0003e2",
-                bordercolor="white",
-                borderpad=5,
-                font=dict(color="white"),
-            )
-        ]
+        plot_bgcolor= dalybgcolor,
+       # annotations=[
+       #     go.Annotation(
+       #         y=-6,
+       #         x=1.2,
+       #         text="Tolerable risk level of 1µDALY pppy",
+       #         bgcolor="#0003e2",
+       #         bordercolor="white",
+       #         borderpad=5,
+       #         font=dict(color="white"),
+       #     )
+       # ]
         # markersize= 12,
     )
 
@@ -541,12 +569,12 @@ def calculate_risk(request, ra_id):
             orientation="h",
             yanchor="top",
             # text= "Reference pathogen",
-            y=-0.1,
+            y=-0.3,
             xanchor="left",
             x=0,
         )
     )
-    fig.add_hline(y=0.000001, line_dash="dashdot", line=dict(color="#0003e2", width=3))
+    fig.add_hline(y=0.000001, line_dash="dashdot", line=dict(color=dlcolor, width=3))
     # annotation_text="tolerable risk level")
     # annotation_position = "top left",
     # annotation=dict(font_size=20, font_family="Times New Roman"))
@@ -641,6 +669,8 @@ def calculate_risk(request, ra_id):
             "daly_plot": daly_plot,
             "risk_plot": risk_plot,
             "ra": ra,
+            "risk": risk,
+            "dalyrisk": risk
         },
     )
 
@@ -701,10 +731,6 @@ def register(request):
         return render(request, "qmratool/register.html")
 
 
-from django.contrib import messages
-from django.contrib.auth import update_session_auth_hash
-from django.contrib.auth.forms import PasswordChangeForm
-from django.shortcuts import render, redirect
 
 
 @login_required(login_url="/login")