From 4d3769615814102e5a13a789507c5de8bd11dd06 Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Fri, 4 Jun 2021 22:58:27 -0400
Subject: [PATCH 1/9] Update heatmaps of overall data yield
---
Snakefile | 2 +
rules/reports.smk | 10 +-
..._yield_per_participant_per_time_segment.py | 122 ++-
src/visualization/plotly_color_utils.py | 835 ++++++++++++++++++
tools/config.schema.yaml | 5 +-
5 files changed, 905 insertions(+), 69 deletions(-)
create mode 100644 src/visualization/plotly_color_utils.py
diff --git a/Snakefile b/Snakefile
index a73f412e..cc1ed4ba 100644
--- a/Snakefile
+++ b/Snakefile
@@ -381,6 +381,8 @@ if config["HEATMAP_SENSOR_ROW_COUNT_PER_TIME_SEGMENT"]["PLOT"]:
files_to_compute.append("reports/data_exploration/heatmap_sensor_row_count_per_time_segment.html")
if config["HEATMAP_PHONE_DATA_YIELD_PER_PARTICIPANT_PER_TIME_SEGMENT"]["PLOT"]:
+ if not config["PHONE_DATA_YIELD"]["PROVIDERS"]["RAPIDS"]["COMPUTE"]:
+ raise ValueError("Error: [PHONE_DATA_YIELD][PROVIDERS][RAPIDS][COMPUTE] must be True in config.yaml to get heatmaps of overall data yield.")
files_to_compute.append("reports/data_exploration/heatmap_phone_data_yield_per_participant_per_time_segment.html")
if config["HEATMAP_FEATURE_CORRELATION_MATRIX"]["PLOT"]:
diff --git a/rules/reports.smk b/rules/reports.smk
index 3a4fa6bd..de58139c 100644
--- a/rules/reports.smk
+++ b/rules/reports.smk
@@ -49,11 +49,13 @@ rule merge_heatmap_sensor_row_count_per_time_segment:
rule heatmap_phone_data_yield_per_participant_per_time_segment:
input:
- phone_data_yield = expand("data/processed/features/{pid}/phone_data_yield.csv", pid=config["PIDS"]),
- participant_file = expand("data/external/participant_files/{pid}.yaml", pid=config["PIDS"]),
- time_segments_labels = expand("data/interim/time_segments/{pid}_time_segments_labels.csv", pid=config["PIDS"])
+ participant_files = expand("data/external/participant_files/{pid}.yaml", pid=config["PIDS"]),
+ time_segments_file = config["TIME_SEGMENTS"]["FILE"],
+ phone_data_yield = "data/processed/features/all_participants/all_sensor_features.csv",
params:
- time = config["HEATMAP_PHONE_DATA_YIELD_PER_PARTICIPANT_PER_TIME_SEGMENT"]["TIME"]
+ pids = config["PIDS"],
+ time = config["HEATMAP_PHONE_DATA_YIELD_PER_PARTICIPANT_PER_TIME_SEGMENT"]["TIME"],
+ time_segments_type = config["TIME_SEGMENTS"]["TYPE"]
output:
"reports/data_exploration/heatmap_phone_data_yield_per_participant_per_time_segment.html"
script:
diff --git a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
index 8f639a62..34897d6e 100644
--- a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
+++ b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
@@ -1,93 +1,87 @@
+from plotly_color_utils import sample_colorscale
import pandas as pd
-import numpy as np
-import plotly.graph_objects as go
+import plotly.express as px
import yaml
+def getPidAndLabel(participant_file_paths, pids):
+ pid2label, y_axis_labels = {}, []
+ for participant_file_path, pid in zip(participant_file_paths, pids):
+ with open(participant_file_path, "r", encoding="utf-8") as f:
+ participant_file = yaml.safe_load(f)
+ label = str(participant_file["PHONE"]["LABEL"])
-def getPhoneDataYieldHeatmap(data_for_plot, y_axis_labels, time_segment, type, time, html_file):
+ pid2label[pid] = label
+ y_axis_labels.append(pid + "." + label)
+ return pid2label, y_axis_labels
- fig = go.Figure(data=go.Heatmap(z=data_for_plot.values.tolist(),
- x=data_for_plot.columns.tolist(),
- y=y_axis_labels,
- hovertext=data_for_plot.values.tolist(),
- hovertemplate="Time since first segment: %{x}
Participant: %{y}
Ratiovalidyielded" + type + ": %{z}
Participant: %{y}
Ratiovalidyielded" + type + ": %{z}
y-axis shows participant information (format: pid.label).
x-axis shows the time since their first segment.
z-axis (color) shows valid yielded " + type + " ratio during a segment instance.")
- else:
- fig.update_layout(title="Heatmap of valid yielded " + type + " ratio for " + time_segment + " segments.
y-axis shows participant information (format: pid.label).
x-axis shows the time.
z-axis (color) shows valid yielded " + type + " ratio during a segment instance.")
+ # Number of minutes after the first start date time of local segments
+ phone_data_yield["local_segment_end_datetime"] = (phone_data_yield["local_segment_end_datetime"] - phone_data_yield["local_segment_start_datetime"].min()) + pd.Timestamp(2000,1,1)
+ phone_data_yield["local_segment_start_datetime"] = (phone_data_yield["local_segment_start_datetime"] - phone_data_yield["local_segment_start_datetime"].min()) + pd.Timestamp(2000,1,1)
- fig["layout"]["xaxis"].update(side="bottom")
- fig["layout"].update(xaxis_title="Time Since First Segment" if time == "RELATIVE_TIME" else "Time")
- fig["layout"].update(margin=dict(t=160))
-
- html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
+ for type in ["minutes", "hours"]:
+
+ column_name = "phone_data_yield_rapids_ratiovalidyielded" + type
+
+ fig = px.timeline(phone_data_yield,
+ x_start="local_segment_start_datetime",
+ x_end="local_segment_end_datetime",
+ y="y_axis_label",
+ color=column_name,
+ color_continuous_scale="Viridis",
+ opacity=0.7,
+ hover_data={'local_segment_start_datetime':False, 'local_segment_end_datetime':False, 'local_segment':True})
+
+ fig.update_layout(title="Heatmap of valid yielded " + type + " ratio for " + time_segment + " segments and " + time.lower().replace("_", " ") + ".
y-axis shows participant information (format: pid.label).
x-axis shows the time" + (" since their first segment" if time == "RELATIVE_TIME" else "") + ".
z-axis (color) shows valid yielded " + type + " ratio during a segment instance.",
+ xaxis=dict(side="bottom", title="Time Since First Segment" if time == "RELATIVE_TIME" else "Time"),
+ yaxis=dict(side="left", title="Participant information"),
+ margin=dict(t=160))
+
+ if time == "RELATIVE_TIME":
+ fig.update_layout(xaxis_tickformat="%y years %j days
%X")
+
+ html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
+
+ return html_file
+pid2label, y_axis_labels = getPidAndLabel(snakemake.input["participant_files"], snakemake.params["pids"])
+time_segments_type = snakemake.params["time_segments_type"] # FREQUENCY or PERIODIC or EVENT
+time = snakemake.params["time"] # ABSOLUTE_TIME or RELATIVE_TIME
+time_segments = pd.read_csv(snakemake.input["time_segments_file"])["label"].unique()
+phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], parse_dates=["local_segment_start_datetime", "local_segment_end_datetime"]).sort_values(by=["pid", "local_segment_start_datetime"])
+if time_segments_type == "FREQUENCY":
+ phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.split("\d+", expand=True, n=1)[0]
-time = snakemake.params["time"]
-y_axis_labels, phone_data_yield_minutes, phone_data_yield_hours = [], {}, {}
-for phone_data_yield_path, participant_file_path, time_segments_path in zip(snakemake.input["phone_data_yield"], snakemake.input["participant_file"], snakemake.input["time_segments_labels"]):
-
- # set pid.label as y_axis_label
- pid = phone_data_yield_path.split("/")[3]
- time_segments = pd.read_csv(time_segments_path, header=0)["label"]
-
- with open(participant_file_path, "r", encoding="utf-8") as f:
- participant_file = yaml.safe_load(f)
- label = participant_file["PHONE"]["LABEL"]
-
- y_axis_label = pid + "." + label
- y_axis_labels.append(y_axis_label)
-
-
- phone_data_yield = pd.read_csv(phone_data_yield_path, index_col=["local_segment_start_datetime"], parse_dates=["local_segment_start_datetime"])
- # make sure the phone_data_yield file contains "phone_data_yield_rapids_ratiovalidyieldedminutes" and "phone_data_yield_rapids_ratiovalidyieldedhours" columns
+html_file = open(snakemake.output[0], "w", encoding="utf-8")
+if phone_data_yield.empty:
+ html_file.write("There is no sensor data for the sensors in [PHONE_DATA_YIELD][SENSORS].")
+else:
+ # Make sure the phone_data_yield file contains both "phone_data_yield_rapids_ratiovalidyieldedminutes" and "phone_data_yield_rapids_ratiovalidyieldedhours" columns
if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
- if not phone_data_yield.empty:
+ phone_data_yield[["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]] = phone_data_yield[["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]].round(3).clip(upper=1)
+ phone_data_yield["y_axis_label"] = phone_data_yield["pid"].apply(lambda pid: pid + "." + str(pid2label[pid]))
+ if time_segments_type == "EVENT":
+ html_file = getPhoneDataYieldHeatmap(phone_data_yield, time, "event", html_file)
+ else: # FREQUENCY or PERIODIC
for time_segment in time_segments:
+
phone_data_yield_per_segment = phone_data_yield[phone_data_yield["local_segment_label"] == time_segment]
if not phone_data_yield_per_segment.empty:
- if time == "RELATIVE_TIME":
- # set number of minutes after the first start date time of local segments as x_axis_label
- phone_data_yield_per_segment.index = phone_data_yield_per_segment.index - phone_data_yield_per_segment.index.min()
- elif time == "ABSOLUTE_TIME":
- pass
- else:
- raise ValueError("[HEATMAP_PHONE_DATA_YIELD_PER_PARTICIPANT_PER_TIME_SEGMENT][TIME] can only be RELATIVE_TIME or ABSOLUTE_TIME")
+ html_file = getPhoneDataYieldHeatmap(phone_data_yield_per_segment, time, time_segment, html_file)
- phone_data_yield_minutes_per_segment = phone_data_yield_per_segment[["phone_data_yield_rapids_ratiovalidyieldedminutes"]].rename(columns={"phone_data_yield_rapids_ratiovalidyieldedminutes": y_axis_label})
- phone_data_yield_hours_per_segment = phone_data_yield_per_segment[["phone_data_yield_rapids_ratiovalidyieldedhours"]].rename(columns={"phone_data_yield_rapids_ratiovalidyieldedhours": y_axis_label})
-
- if time_segment not in phone_data_yield_minutes.keys():
- phone_data_yield_minutes[time_segment] = phone_data_yield_minutes_per_segment
- phone_data_yield_hours[time_segment] = phone_data_yield_hours_per_segment
- else:
- phone_data_yield_minutes[time_segment] = pd.concat([phone_data_yield_minutes[time_segment], phone_data_yield_minutes_per_segment], axis=1, sort=True)
- phone_data_yield_hours[time_segment] = pd.concat([phone_data_yield_hours[time_segment], phone_data_yield_hours_per_segment], axis=1, sort=True)
-
-
-html_file = open(snakemake.output[0], "a", encoding="utf-8")
-if len(phone_data_yield_minutes.keys()) == 0:
- html_file.write("There is no sensor data for the sensors in [PHONE_DATA_YIELD][SENSORS].")
-for time_segment in phone_data_yield_minutes.keys():
- minutes_data_for_plot = phone_data_yield_minutes[time_segment].transpose().reindex(pd.Index(y_axis_labels)).round(3)
- hours_data_for_plot = phone_data_yield_hours[time_segment].transpose().reindex(pd.Index(y_axis_labels)).round(3)
-
- getPhoneDataYieldHeatmap(minutes_data_for_plot, y_axis_labels, time_segment, "minutes", time, html_file)
- getPhoneDataYieldHeatmap(hours_data_for_plot, y_axis_labels, time_segment, "hours", time, html_file)
html_file.close()
diff --git a/src/visualization/plotly_color_utils.py b/src/visualization/plotly_color_utils.py
new file mode 100644
index 00000000..12ffb279
--- /dev/null
+++ b/src/visualization/plotly_color_utils.py
@@ -0,0 +1,835 @@
+"""
+Source: https://github.com/plotly/plotly.py/pull/3136
+=====
+colors
+=====
+Functions that manipulate colors and arrays of colors.
+-----
+There are three basic types of color types: rgb, hex and tuple:
+rgb - An rgb color is a string of the form 'rgb(a,b,c)' where a, b and c are
+integers between 0 and 255 inclusive.
+hex - A hex color is a string of the form '#xxxxxx' where each x is a
+character that belongs to the set [0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f]. This is
+just the set of characters used in the hexadecimal numeric system.
+tuple - A tuple color is a 3-tuple of the form (a,b,c) where a, b and c are
+floats between 0 and 1 inclusive.
+-----
+Colormaps and Colorscales:
+A colormap or a colorscale is a correspondence between values - Pythonic
+objects such as strings and floats - to colors.
+There are typically two main types of colormaps that exist: numerical and
+categorical colormaps.
+Numerical:
+----------
+Numerical colormaps are used when the coloring column being used takes a
+spectrum of values or numbers.
+A classic example from the Plotly library:
+```
+rainbow_colorscale = [
+ [0, 'rgb(150,0,90)'], [0.125, 'rgb(0,0,200)'],
+ [0.25, 'rgb(0,25,255)'], [0.375, 'rgb(0,152,255)'],
+ [0.5, 'rgb(44,255,150)'], [0.625, 'rgb(151,255,0)'],
+ [0.75, 'rgb(255,234,0)'], [0.875, 'rgb(255,111,0)'],
+ [1, 'rgb(255,0,0)']
+]
+```
+Notice that this colorscale is a list of lists with each inner list containing
+a number and a color. These left hand numbers in the nested lists go from 0 to
+1, and they are like pointers tell you when a number is mapped to a specific
+color.
+If you have a column of numbers `col_num` that you want to plot, and you know
+```
+min(col_num) = 0
+max(col_num) = 100
+```
+then if you pull out the number `12.5` in the list and want to figure out what
+color the corresponding chart element (bar, scatter plot, etc) is going to be,
+you'll figure out that proportionally 12.5 to 100 is the same as 0.125 to 1.
+So, the point will be mapped to 'rgb(0,0,200)'.
+All other colors between the pinned values in a colorscale are linearly
+interpolated.
+Categorical:
+------------
+Alternatively, a categorical colormap is used to assign a specific value in a
+color column to a specific color everytime it appears in the dataset.
+A column of strings in a panadas.dataframe that is chosen to serve as the
+color index would naturally use a categorical colormap. However, you can
+choose to use a categorical colormap with a column of numbers.
+Be careful! If you have a lot of unique numbers in your color column you will
+end up with a colormap that is massive and may slow down graphing performance.
+"""
+from __future__ import absolute_import
+
+import decimal
+from numbers import Number
+import six
+
+from _plotly_utils import exceptions
+
+
+# Built-in qualitative color sequences and sequential,
+# diverging and cyclical color scales.
+
+DEFAULT_PLOTLY_COLORS = [
+ "rgb(31, 119, 180)",
+ "rgb(255, 127, 14)",
+ "rgb(44, 160, 44)",
+ "rgb(214, 39, 40)",
+ "rgb(148, 103, 189)",
+ "rgb(140, 86, 75)",
+ "rgb(227, 119, 194)",
+ "rgb(127, 127, 127)",
+ "rgb(188, 189, 34)",
+ "rgb(23, 190, 207)",
+]
+
+PLOTLY_SCALES = {
+ "Greys": [[0, "rgb(0,0,0)"], [1, "rgb(255,255,255)"]],
+ "YlGnBu": [
+ [0, "rgb(8,29,88)"],
+ [0.125, "rgb(37,52,148)"],
+ [0.25, "rgb(34,94,168)"],
+ [0.375, "rgb(29,145,192)"],
+ [0.5, "rgb(65,182,196)"],
+ [0.625, "rgb(127,205,187)"],
+ [0.75, "rgb(199,233,180)"],
+ [0.875, "rgb(237,248,217)"],
+ [1, "rgb(255,255,217)"],
+ ],
+ "Greens": [
+ [0, "rgb(0,68,27)"],
+ [0.125, "rgb(0,109,44)"],
+ [0.25, "rgb(35,139,69)"],
+ [0.375, "rgb(65,171,93)"],
+ [0.5, "rgb(116,196,118)"],
+ [0.625, "rgb(161,217,155)"],
+ [0.75, "rgb(199,233,192)"],
+ [0.875, "rgb(229,245,224)"],
+ [1, "rgb(247,252,245)"],
+ ],
+ "YlOrRd": [
+ [0, "rgb(128,0,38)"],
+ [0.125, "rgb(189,0,38)"],
+ [0.25, "rgb(227,26,28)"],
+ [0.375, "rgb(252,78,42)"],
+ [0.5, "rgb(253,141,60)"],
+ [0.625, "rgb(254,178,76)"],
+ [0.75, "rgb(254,217,118)"],
+ [0.875, "rgb(255,237,160)"],
+ [1, "rgb(255,255,204)"],
+ ],
+ "Bluered": [[0, "rgb(0,0,255)"], [1, "rgb(255,0,0)"]],
+ # modified RdBu based on
+ # www.sandia.gov/~kmorel/documents/ColorMaps/ColorMapsExpanded.pdf
+ "RdBu": [
+ [0, "rgb(5,10,172)"],
+ [0.35, "rgb(106,137,247)"],
+ [0.5, "rgb(190,190,190)"],
+ [0.6, "rgb(220,170,132)"],
+ [0.7, "rgb(230,145,90)"],
+ [1, "rgb(178,10,28)"],
+ ],
+ # Scale for non-negative numeric values
+ "Reds": [
+ [0, "rgb(220,220,220)"],
+ [0.2, "rgb(245,195,157)"],
+ [0.4, "rgb(245,160,105)"],
+ [1, "rgb(178,10,28)"],
+ ],
+ # Scale for non-positive numeric values
+ "Blues": [
+ [0, "rgb(5,10,172)"],
+ [0.35, "rgb(40,60,190)"],
+ [0.5, "rgb(70,100,245)"],
+ [0.6, "rgb(90,120,245)"],
+ [0.7, "rgb(106,137,247)"],
+ [1, "rgb(220,220,220)"],
+ ],
+ "Picnic": [
+ [0, "rgb(0,0,255)"],
+ [0.1, "rgb(51,153,255)"],
+ [0.2, "rgb(102,204,255)"],
+ [0.3, "rgb(153,204,255)"],
+ [0.4, "rgb(204,204,255)"],
+ [0.5, "rgb(255,255,255)"],
+ [0.6, "rgb(255,204,255)"],
+ [0.7, "rgb(255,153,255)"],
+ [0.8, "rgb(255,102,204)"],
+ [0.9, "rgb(255,102,102)"],
+ [1, "rgb(255,0,0)"],
+ ],
+ "Rainbow": [
+ [0, "rgb(150,0,90)"],
+ [0.125, "rgb(0,0,200)"],
+ [0.25, "rgb(0,25,255)"],
+ [0.375, "rgb(0,152,255)"],
+ [0.5, "rgb(44,255,150)"],
+ [0.625, "rgb(151,255,0)"],
+ [0.75, "rgb(255,234,0)"],
+ [0.875, "rgb(255,111,0)"],
+ [1, "rgb(255,0,0)"],
+ ],
+ "Portland": [
+ [0, "rgb(12,51,131)"],
+ [0.25, "rgb(10,136,186)"],
+ [0.5, "rgb(242,211,56)"],
+ [0.75, "rgb(242,143,56)"],
+ [1, "rgb(217,30,30)"],
+ ],
+ "Jet": [
+ [0, "rgb(0,0,131)"],
+ [0.125, "rgb(0,60,170)"],
+ [0.375, "rgb(5,255,255)"],
+ [0.625, "rgb(255,255,0)"],
+ [0.875, "rgb(250,0,0)"],
+ [1, "rgb(128,0,0)"],
+ ],
+ "Hot": [
+ [0, "rgb(0,0,0)"],
+ [0.3, "rgb(230,0,0)"],
+ [0.6, "rgb(255,210,0)"],
+ [1, "rgb(255,255,255)"],
+ ],
+ "Blackbody": [
+ [0, "rgb(0,0,0)"],
+ [0.2, "rgb(230,0,0)"],
+ [0.4, "rgb(230,210,0)"],
+ [0.7, "rgb(255,255,255)"],
+ [1, "rgb(160,200,255)"],
+ ],
+ "Earth": [
+ [0, "rgb(0,0,130)"],
+ [0.1, "rgb(0,180,180)"],
+ [0.2, "rgb(40,210,40)"],
+ [0.4, "rgb(230,230,50)"],
+ [0.6, "rgb(120,70,20)"],
+ [1, "rgb(255,255,255)"],
+ ],
+ "Electric": [
+ [0, "rgb(0,0,0)"],
+ [0.15, "rgb(30,0,100)"],
+ [0.4, "rgb(120,0,100)"],
+ [0.6, "rgb(160,90,0)"],
+ [0.8, "rgb(230,200,0)"],
+ [1, "rgb(255,250,220)"],
+ ],
+ "Viridis": [
+ [0, "#440154"],
+ [0.06274509803921569, "#48186a"],
+ [0.12549019607843137, "#472d7b"],
+ [0.18823529411764706, "#424086"],
+ [0.25098039215686274, "#3b528b"],
+ [0.3137254901960784, "#33638d"],
+ [0.3764705882352941, "#2c728e"],
+ [0.4392156862745098, "#26828e"],
+ [0.5019607843137255, "#21918c"],
+ [0.5647058823529412, "#1fa088"],
+ [0.6274509803921569, "#28ae80"],
+ [0.6901960784313725, "#3fbc73"],
+ [0.7529411764705882, "#5ec962"],
+ [0.8156862745098039, "#84d44b"],
+ [0.8784313725490196, "#addc30"],
+ [0.9411764705882353, "#d8e219"],
+ [1, "#fde725"],
+ ],
+ "Cividis": [
+ [0.000000, "rgb(0,32,76)"],
+ [0.058824, "rgb(0,42,102)"],
+ [0.117647, "rgb(0,52,110)"],
+ [0.176471, "rgb(39,63,108)"],
+ [0.235294, "rgb(60,74,107)"],
+ [0.294118, "rgb(76,85,107)"],
+ [0.352941, "rgb(91,95,109)"],
+ [0.411765, "rgb(104,106,112)"],
+ [0.470588, "rgb(117,117,117)"],
+ [0.529412, "rgb(131,129,120)"],
+ [0.588235, "rgb(146,140,120)"],
+ [0.647059, "rgb(161,152,118)"],
+ [0.705882, "rgb(176,165,114)"],
+ [0.764706, "rgb(192,177,109)"],
+ [0.823529, "rgb(209,191,102)"],
+ [0.882353, "rgb(225,204,92)"],
+ [0.941176, "rgb(243,219,79)"],
+ [1.000000, "rgb(255,233,69)"],
+ ],
+}
+
+
+def color_parser(colors, function):
+ """
+ Takes color(s) and a function and applies the function on the color(s)
+ In particular, this function identifies whether the given color object
+ is an iterable or not and applies the given color-parsing function to
+ the color or iterable of colors. If given an iterable, it will only be
+ able to work with it if all items in the iterable are of the same type
+ - rgb string, hex string or tuple
+ """
+ if isinstance(colors, str):
+ return function(colors)
+
+ if isinstance(colors, tuple) and isinstance(colors[0], Number):
+ return function(colors)
+
+ if hasattr(colors, "__iter__"):
+ if isinstance(colors, tuple):
+ new_color_tuple = tuple(function(item) for item in colors)
+ return new_color_tuple
+
+ else:
+ new_color_list = [function(item) for item in colors]
+ return new_color_list
+
+
+def validate_colors(colors, colortype="tuple"):
+ """
+ Validates color(s) and returns a list of color(s) of a specified type
+ """
+ from numbers import Number
+
+ if colors is None:
+ colors = DEFAULT_PLOTLY_COLORS
+
+ if isinstance(colors, str):
+ if colors in PLOTLY_SCALES:
+ colors_list = colorscale_to_colors(PLOTLY_SCALES[colors])
+ # TODO: fix _gantt.py/_scatter.py so that they can accept the
+ # actual colorscale and not just a list of the first and last
+ # color in the plotly colorscale. In resolving this issue we
+ # will be removing the immediate line below
+ colors = [colors_list[0]] + [colors_list[-1]]
+ elif "rgb" in colors or "#" in colors:
+ colors = [colors]
+ else:
+ raise exceptions.PlotlyError(
+ "If your colors variable is a string, it must be a "
+ "Plotly scale, an rgb color or a hex color."
+ )
+
+ elif isinstance(colors, tuple):
+ if isinstance(colors[0], Number):
+ colors = [colors]
+ else:
+ colors = list(colors)
+
+ # convert color elements in list to tuple color
+ for j, each_color in enumerate(colors):
+ if "rgb" in each_color:
+ each_color = color_parser(each_color, unlabel_rgb)
+ for value in each_color:
+ if value > 255.0:
+ raise exceptions.PlotlyError(
+ "Whoops! The elements in your rgb colors "
+ "tuples cannot exceed 255.0."
+ )
+ each_color = color_parser(each_color, unconvert_from_RGB_255)
+ colors[j] = each_color
+
+ if "#" in each_color:
+ each_color = color_parser(each_color, hex_to_rgb)
+ each_color = color_parser(each_color, unconvert_from_RGB_255)
+
+ colors[j] = each_color
+
+ if isinstance(each_color, tuple):
+ for value in each_color:
+ if value > 1.0:
+ raise exceptions.PlotlyError(
+ "Whoops! The elements in your colors tuples "
+ "cannot exceed 1.0."
+ )
+ colors[j] = each_color
+
+ if colortype == "rgb" and not isinstance(colors, six.string_types):
+ for j, each_color in enumerate(colors):
+ rgb_color = color_parser(each_color, convert_to_RGB_255)
+ colors[j] = color_parser(rgb_color, label_rgb)
+
+ return colors
+
+
+def validate_colors_dict(colors, colortype="tuple"):
+ """
+ Validates dictionary of color(s)
+ """
+ # validate each color element in the dictionary
+ for key in colors:
+ if "rgb" in colors[key]:
+ colors[key] = color_parser(colors[key], unlabel_rgb)
+ for value in colors[key]:
+ if value > 255.0:
+ raise exceptions.PlotlyError(
+ "Whoops! The elements in your rgb colors "
+ "tuples cannot exceed 255.0."
+ )
+ colors[key] = color_parser(colors[key], unconvert_from_RGB_255)
+
+ if "#" in colors[key]:
+ colors[key] = color_parser(colors[key], hex_to_rgb)
+ colors[key] = color_parser(colors[key], unconvert_from_RGB_255)
+
+ if isinstance(colors[key], tuple):
+ for value in colors[key]:
+ if value > 1.0:
+ raise exceptions.PlotlyError(
+ "Whoops! The elements in your colors tuples "
+ "cannot exceed 1.0."
+ )
+
+ if colortype == "rgb":
+ for key in colors:
+ colors[key] = color_parser(colors[key], convert_to_RGB_255)
+ colors[key] = color_parser(colors[key], label_rgb)
+
+ return colors
+
+
+def convert_colors_to_same_type(
+ colors,
+ colortype="rgb",
+ scale=None,
+ return_default_colors=False,
+ num_of_defualt_colors=2,
+):
+ """
+ Converts color(s) to the specified color type
+ Takes a single color or an iterable of colors, as well as a list of scale
+ values, and outputs a 2-pair of the list of color(s) converted all to an
+ rgb or tuple color type, aswell as the scale as the second element. If
+ colors is a Plotly Scale name, then 'scale' will be forced to the scale
+ from the respective colorscale and the colors in that colorscale will also
+ be coverted to the selected colortype. If colors is None, then there is an
+ option to return portion of the DEFAULT_PLOTLY_COLORS
+ :param (str|tuple|list) colors: either a plotly scale name, an rgb or hex
+ color, a color tuple or a list/tuple of colors
+ :param (list) scale: see docs for validate_scale_values()
+ :rtype (tuple) (colors_list, scale) if scale is None in the function call,
+ then scale will remain None in the returned tuple
+ """
+ colors_list = []
+
+ if colors is None and return_default_colors is True:
+ colors_list = DEFAULT_PLOTLY_COLORS[0:num_of_defualt_colors]
+
+ if isinstance(colors, str):
+ if colors in PLOTLY_SCALES:
+ colors_list = colorscale_to_colors(PLOTLY_SCALES[colors])
+ if scale is None:
+ scale = colorscale_to_scale(PLOTLY_SCALES[colors])
+
+ elif "rgb" in colors or "#" in colors:
+ colors_list = [colors]
+
+ elif isinstance(colors, tuple):
+ if isinstance(colors[0], Number):
+ colors_list = [colors]
+ else:
+ colors_list = list(colors)
+
+ elif isinstance(colors, list):
+ colors_list = colors
+
+ # validate scale
+ if scale is not None:
+ validate_scale_values(scale)
+
+ if len(colors_list) != len(scale):
+ raise exceptions.PlotlyError(
+ "Make sure that the length of your scale matches the length "
+ "of your list of colors which is {}.".format(len(colors_list))
+ )
+
+ # convert all colors to rgb
+ for j, each_color in enumerate(colors_list):
+ if "#" in each_color:
+ each_color = color_parser(each_color, hex_to_rgb)
+ each_color = color_parser(each_color, label_rgb)
+ colors_list[j] = each_color
+
+ elif isinstance(each_color, tuple):
+ each_color = color_parser(each_color, convert_to_RGB_255)
+ each_color = color_parser(each_color, label_rgb)
+ colors_list[j] = each_color
+
+ if colortype == "rgb":
+ return (colors_list, scale)
+ elif colortype == "tuple":
+ for j, each_color in enumerate(colors_list):
+ each_color = color_parser(each_color, unlabel_rgb)
+ each_color = color_parser(each_color, unconvert_from_RGB_255)
+ colors_list[j] = each_color
+ return (colors_list, scale)
+ else:
+ raise exceptions.PlotlyError(
+ "You must select either rgb or tuple for your colortype variable."
+ )
+
+
+def convert_dict_colors_to_same_type(colors_dict, colortype="rgb"):
+ """
+ Converts a colors in a dictionary of colors to the specified color type
+ :param (dict) colors_dict: a dictionary whose values are single colors
+ """
+ for key in colors_dict:
+ if "#" in colors_dict[key]:
+ colors_dict[key] = color_parser(colors_dict[key], hex_to_rgb)
+ colors_dict[key] = color_parser(colors_dict[key], label_rgb)
+
+ elif isinstance(colors_dict[key], tuple):
+ colors_dict[key] = color_parser(colors_dict[key], convert_to_RGB_255)
+ colors_dict[key] = color_parser(colors_dict[key], label_rgb)
+
+ if colortype == "rgb":
+ return colors_dict
+ elif colortype == "tuple":
+ for key in colors_dict:
+ colors_dict[key] = color_parser(colors_dict[key], unlabel_rgb)
+ colors_dict[key] = color_parser(colors_dict[key], unconvert_from_RGB_255)
+ return colors_dict
+ else:
+ raise exceptions.PlotlyError(
+ "You must select either rgb or tuple for your colortype variable."
+ )
+
+
+def validate_scale_values(scale):
+ """
+ Validates scale values from a colorscale
+ :param (list) scale: a strictly increasing list of floats that begins
+ with 0 and ends with 1. Its usage derives from a colorscale which is
+ a list of two-lists (a list with two elements) of the form
+ [value, color] which are used to determine how interpolation weighting
+ works between the colors in the colorscale. Therefore scale is just
+ the extraction of these values from the two-lists in order
+ """
+ if len(scale) < 2:
+ raise exceptions.PlotlyError(
+ "You must input a list of scale values that has at least two values."
+ )
+
+ if (scale[0] != 0) or (scale[-1] != 1):
+ raise exceptions.PlotlyError(
+ "The first and last number in your scale must be 0.0 and 1.0 "
+ "respectively."
+ )
+
+ if not all(x < y for x, y in zip(scale, scale[1:])):
+ raise exceptions.PlotlyError(
+ "'scale' must be a list that contains a strictly increasing "
+ "sequence of numbers."
+ )
+
+
+def validate_colorscale(colorscale):
+ """Validate the structure, scale values and colors of colorscale."""
+ if not isinstance(colorscale, list):
+ # TODO Write tests for these exceptions
+ raise exceptions.PlotlyError("A valid colorscale must be a list.")
+ if not all(isinstance(innerlist, list) for innerlist in colorscale):
+ raise exceptions.PlotlyError("A valid colorscale must be a list of lists.")
+ colorscale_colors = colorscale_to_colors(colorscale)
+ scale_values = colorscale_to_scale(colorscale)
+
+ validate_scale_values(scale_values)
+ validate_colors(colorscale_colors)
+
+
+def make_colorscale(colors, scale=None):
+ """
+ Makes a colorscale from a list of colors and a scale
+ Takes a list of colors and scales and constructs a colorscale based
+ on the colors in sequential order. If 'scale' is left empty, a linear-
+ interpolated colorscale will be generated. If 'scale' is a specificed
+ list, it must be the same legnth as colors and must contain all floats
+ For documentation regarding to the form of the output, see
+ https://plot.ly/python/reference/#mesh3d-colorscale
+ :param (list) colors: a list of single colors
+ """
+ colorscale = []
+
+ # validate minimum colors length of 2
+ if len(colors) < 2:
+ raise exceptions.PlotlyError(
+ "You must input a list of colors that has at least two colors."
+ )
+
+ if scale is None:
+ scale_incr = 1.0 / (len(colors) - 1)
+ return [[i * scale_incr, color] for i, color in enumerate(colors)]
+
+ else:
+ if len(colors) != len(scale):
+ raise exceptions.PlotlyError(
+ "The length of colors and scale must be the same."
+ )
+
+ validate_scale_values(scale)
+
+ colorscale = [list(tup) for tup in zip(scale, colors)]
+ return colorscale
+
+
+def find_intermediate_color(lowcolor, highcolor, intermed, colortype="tuple"):
+ """
+ Returns the color at a given distance between two colors
+ This function takes two color tuples, where each element is between 0
+ and 1, along with a value 0 < intermed < 1 and returns a color that is
+ intermed-percent from lowcolor to highcolor. If colortype is set to 'rgb',
+ the function will automatically convert the rgb type to a tuple, find the
+ intermediate color and return it as an rgb color.
+ """
+ if colortype == "rgb":
+ # convert to tuple color, eg. (1, 0.45, 0.7)
+ lowcolor = unlabel_rgb(lowcolor)
+ highcolor = unlabel_rgb(highcolor)
+
+ diff_0 = float(highcolor[0] - lowcolor[0])
+ diff_1 = float(highcolor[1] - lowcolor[1])
+ diff_2 = float(highcolor[2] - lowcolor[2])
+
+ inter_med_tuple = (
+ lowcolor[0] + intermed * diff_0,
+ lowcolor[1] + intermed * diff_1,
+ lowcolor[2] + intermed * diff_2,
+ )
+
+ if colortype == "rgb":
+ # back to an rgb string, e.g. rgb(30, 20, 10)
+ inter_med_rgb = label_rgb(inter_med_tuple)
+ return inter_med_rgb
+
+ return inter_med_tuple
+
+
+def unconvert_from_RGB_255(colors):
+ """
+ Return a tuple where each element gets divided by 255
+ Takes a (list of) color tuple(s) where each element is between 0 and
+ 255. Returns the same tuples where each tuple element is normalized to
+ a value between 0 and 1
+ """
+ return (colors[0] / (255.0), colors[1] / (255.0), colors[2] / (255.0))
+
+
+def convert_to_RGB_255(colors):
+ """
+ Multiplies each element of a triplet by 255
+ Each coordinate of the color tuple is rounded to the nearest float and
+ then is turned into an integer. If a number is of the form x.5, then
+ if x is odd, the number rounds up to (x+1). Otherwise, it rounds down
+ to just x. This is the way rounding works in Python 3 and in current
+ statistical analysis to avoid rounding bias
+ :param (list) rgb_components: grabs the three R, G and B values to be
+ returned as computed in the function
+ """
+ rgb_components = []
+
+ for component in colors:
+ rounded_num = decimal.Decimal(str(component * 255.0)).quantize(
+ decimal.Decimal("1"), rounding=decimal.ROUND_HALF_EVEN
+ )
+ # convert rounded number to an integer from 'Decimal' form
+ rounded_num = int(rounded_num)
+ rgb_components.append(rounded_num)
+
+ return (rgb_components[0], rgb_components[1], rgb_components[2])
+
+
+def n_colors(lowcolor, highcolor, n_colors, colortype="tuple"):
+ """
+ Splits a low and high color into a list of n_colors colors in it
+ Accepts two color tuples and returns a list of n_colors colors
+ which form the intermediate colors between lowcolor and highcolor
+ from linearly interpolating through RGB space. If colortype is 'rgb'
+ the function will return a list of colors in the same form.
+ """
+ if colortype == "rgb":
+ # convert to tuple
+ lowcolor = unlabel_rgb(lowcolor)
+ highcolor = unlabel_rgb(highcolor)
+
+ diff_0 = float(highcolor[0] - lowcolor[0])
+ incr_0 = diff_0 / (n_colors - 1)
+ diff_1 = float(highcolor[1] - lowcolor[1])
+ incr_1 = diff_1 / (n_colors - 1)
+ diff_2 = float(highcolor[2] - lowcolor[2])
+ incr_2 = diff_2 / (n_colors - 1)
+ list_of_colors = []
+
+ for index in range(n_colors):
+ new_tuple = (
+ lowcolor[0] + (index * incr_0),
+ lowcolor[1] + (index * incr_1),
+ lowcolor[2] + (index * incr_2),
+ )
+ list_of_colors.append(new_tuple)
+
+ if colortype == "rgb":
+ # back to an rgb string
+ list_of_colors = color_parser(list_of_colors, label_rgb)
+
+ return list_of_colors
+
+
+def label_rgb(colors):
+ """
+ Takes tuple (a, b, c) and returns an rgb color 'rgb(a, b, c)'
+ """
+ return "rgb(%s, %s, %s)" % (colors[0], colors[1], colors[2])
+
+
+def unlabel_rgb(colors):
+ """
+ Takes rgb color(s) 'rgb(a, b, c)' and returns tuple(s) (a, b, c)
+ This function takes either an 'rgb(a, b, c)' color or a list of
+ such colors and returns the color tuples in tuple(s) (a, b, c)
+ """
+ str_vals = ""
+ for index in range(len(colors)):
+ try:
+ float(colors[index])
+ str_vals = str_vals + colors[index]
+ except ValueError:
+ if colors[index] == "," or colors[index] == ".":
+ str_vals = str_vals + colors[index]
+
+ str_vals = str_vals + ","
+ numbers = []
+ str_num = ""
+ for char in str_vals:
+ if char != ",":
+ str_num = str_num + char
+ else:
+ numbers.append(float(str_num))
+ str_num = ""
+ return (numbers[0], numbers[1], numbers[2])
+
+
+def hex_to_rgb(value):
+ """
+ Calculates rgb values from a hex color code.
+ :param (string) value: Hex color string
+ :rtype (tuple) (r_value, g_value, b_value): tuple of rgb values
+ """
+ value = value.lstrip("#")
+ hex_total_length = len(value)
+ rgb_section_length = hex_total_length // 3
+ return tuple(
+ int(value[i : i + rgb_section_length], 16)
+ for i in range(0, hex_total_length, rgb_section_length)
+ )
+
+
+def colorscale_to_colors(colorscale):
+ """
+ Extracts the colors from colorscale as a list
+ """
+ color_list = []
+ for item in colorscale:
+ color_list.append(item[1])
+ return color_list
+
+
+def colorscale_to_scale(colorscale):
+ """
+ Extracts the interpolation scale values from colorscale as a list
+ """
+ scale_list = []
+ for item in colorscale:
+ scale_list.append(item[0])
+ return scale_list
+
+
+def convert_colorscale_to_rgb(colorscale):
+ """
+ Converts the colors in a colorscale to rgb colors
+ A colorscale is an array of arrays, each with a numeric value as the
+ first item and a color as the second. This function specifically is
+ converting a colorscale with tuple colors (each coordinate between 0
+ and 1) into a colorscale with the colors transformed into rgb colors
+ """
+ for color in colorscale:
+ color[1] = convert_to_RGB_255(color[1])
+
+ for color in colorscale:
+ color[1] = label_rgb(color[1])
+ return colorscale
+
+
+def named_colorscales():
+ """
+ Returns lowercased names of built-in continuous colorscales.
+ """
+ from _plotly_utils.basevalidators import ColorscaleValidator
+
+ return [c for c in ColorscaleValidator("", "").named_colorscales]
+
+
+def get_colorscale(name):
+ """
+ Returns the colorscale for a given name. See `named_colorscales` for the
+ built-in colorscales.
+ """
+ from _plotly_utils.basevalidators import ColorscaleValidator
+
+ if not isinstance(name, str):
+ raise exceptions.PlotlyError("Name argument have to be a string.")
+
+ name = name.lower()
+ if name[-2:] == "_r":
+ should_reverse = True
+ name = name[:-2]
+ else:
+ should_reverse = False
+
+ if name in ColorscaleValidator("", "").named_colorscales:
+ colorscale = ColorscaleValidator("", "").named_colorscales[name]
+ else:
+ raise exceptions.PlotlyError(f"Colorscale {name} is not a built-in scale.")
+
+ if should_reverse:
+ colorscale = colorscale[::-1]
+ return make_colorscale(colorscale)
+
+
+def sample_colorscale(colorscale, samplepoints, low=0.0, high=1.0, colortype="rgb"):
+ """
+ Samples a colorscale at specific points.
+ Interpolates between colors in a colorscale to find the specific colors
+ corresponding to the specified sample values. The colorscale can be specified
+ as a list of `[scale, color]` pairs, as a list of colors, or as a named
+ plotly colorscale. The samplepoints can be specefied as an iterable of specific
+ points in the range [0.0, 1.0], or as an integer number of points which will
+ be spaced equally between the low value (default 0.0) and the high value
+ (default 1.0). The output is a list of colors, formatted according to the
+ specified colortype.
+ """
+ from bisect import bisect_left
+
+ try:
+ validate_colorscale(colorscale)
+ except exceptions.PlotlyError:
+ if isinstance(colorscale, str):
+ colorscale = get_colorscale(colorscale)
+ else:
+ colorscale = make_colorscale(colorscale)
+
+ scale = colorscale_to_scale(colorscale)
+ validate_scale_values(scale)
+ colors = colorscale_to_colors(colorscale)
+ colors = validate_colors(colors, colortype="tuple")
+
+ if isinstance(samplepoints, int):
+ samplepoints = [
+ low + idx / (samplepoints - 1) * (high - low) for idx in range(samplepoints)
+ ]
+ elif isinstance(samplepoints, float):
+ samplepoints = [samplepoints]
+
+ sampled_colors = []
+ for point in samplepoints:
+ high = bisect_left(scale, point)
+ low = high - 1
+ interpolant = (point - scale[low]) / (scale[high] - scale[low])
+ sampled_color = find_intermediate_color(colors[low], colors[high], interpolant)
+ sampled_colors.append(sampled_color)
+ return validate_colors(sampled_colors, colortype=colortype)
diff --git a/tools/config.schema.yaml b/tools/config.schema.yaml
index 01e7ec02..d1c30054 100644
--- a/tools/config.schema.yaml
+++ b/tools/config.schema.yaml
@@ -1187,10 +1187,13 @@ properties:
HEATMAP_PHONE_DATA_YIELD_PER_PARTICIPANT_PER_TIME_SEGMENT:
type: object
- required: [PLOT]
+ required: [PLOT, TIME]
properties:
PLOT:
type: boolean
+ TIME:
+ type: string
+ enum: [ABSOLUTE_TIME, RELATIVE_TIME]
HEATMAP_SENSORS_PER_MINUTE_PER_TIME_SEGMENT:
type: object
From f436f1f530e7d3f564c3b1face60e07f24e44412 Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Wed, 16 Jun 2021 17:10:08 -0400
Subject: [PATCH 2/9] Update heatmap of correlation matrix
---
rules/reports.smk | 1 +
src/visualization/heatmap_feature_correlation_matrix.py | 9 ++++++++-
2 files changed, 9 insertions(+), 1 deletion(-)
diff --git a/rules/reports.smk b/rules/reports.smk
index de58139c..7935888d 100644
--- a/rules/reports.smk
+++ b/rules/reports.smk
@@ -65,6 +65,7 @@ rule heatmap_feature_correlation_matrix:
input:
all_sensor_features = "data/processed/features/all_participants/all_sensor_features.csv" # before data cleaning
params:
+ time_segments_type = config["TIME_SEGMENTS"]["TYPE"],
min_rows_ratio = config["HEATMAP_FEATURE_CORRELATION_MATRIX"]["MIN_ROWS_RATIO"],
corr_threshold = config["HEATMAP_FEATURE_CORRELATION_MATRIX"]["CORR_THRESHOLD"],
corr_method = config["HEATMAP_FEATURE_CORRELATION_MATRIX"]["CORR_METHOD"]
diff --git a/src/visualization/heatmap_feature_correlation_matrix.py b/src/visualization/heatmap_feature_correlation_matrix.py
index 8fb01934..5deff5c5 100644
--- a/src/visualization/heatmap_feature_correlation_matrix.py
+++ b/src/visualization/heatmap_feature_correlation_matrix.py
@@ -17,12 +17,19 @@ def getCorrMatrixHeatmap(corr_matrix, time_segment, html_file):
html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
-
+time_segments_type = snakemake.params["time_segments_type"]
min_rows_ratio = snakemake.params["min_rows_ratio"]
corr_threshold = snakemake.params["corr_threshold"]
corr_method = snakemake.params["corr_method"]
features = pd.read_csv(snakemake.input["all_sensor_features"])
+
+
+if time_segments_type == "FREQUENCY":
+ features["local_segment_label"] = features["local_segment_label"].str.split("\d+", expand=True, n=1)[0]
+if time_segments_type == "EVENT":
+ features["local_segment_label"] = "event"
+
time_segments = set(features["local_segment_label"])
html_file = open(snakemake.output[0], "a", encoding="utf-8")
From e98a8ff7ca77f060963dfa528f18db0f861ab004 Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Thu, 17 Jun 2021 12:27:32 -0400
Subject: [PATCH 3/9] Update histogram of phone data yield
---
rules/reports.smk | 2 ++
src/visualization/histogram_phone_data_yield.py | 6 ++++++
2 files changed, 8 insertions(+)
diff --git a/rules/reports.smk b/rules/reports.smk
index 7935888d..caed8993 100644
--- a/rules/reports.smk
+++ b/rules/reports.smk
@@ -1,6 +1,8 @@
rule histogram_phone_data_yield:
input:
"data/processed/features/all_participants/all_sensor_features.csv"
+ params:
+ time_segments_type = config["TIME_SEGMENTS"]["TYPE"]
output:
"reports/data_exploration/histogram_phone_data_yield.html"
script:
diff --git a/src/visualization/histogram_phone_data_yield.py b/src/visualization/histogram_phone_data_yield.py
index cd15ec8d..a8b90d77 100644
--- a/src/visualization/histogram_phone_data_yield.py
+++ b/src/visualization/histogram_phone_data_yield.py
@@ -2,8 +2,14 @@ import pandas as pd
import plotly.express as px
+time_segments_type = snakemake.params["time_segments_type"]
phone_data_yield = pd.read_csv(snakemake.input[0])
+if time_segments_type == "FREQUENCY":
+ phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.split("\d+", expand=True, n=1)[0]
+if time_segments_type == "EVENT":
+ phone_data_yield["local_segment_label"] = "event"
+
# make sure the input file contains "phone_data_yield_rapids_ratiovalidyieldedminutes" and "phone_data_yield_rapids_ratiovalidyieldedhours" columns
if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
From bc06477d8902d4f2577d6a85253b64a01c4425ae Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Fri, 25 Jun 2021 13:14:07 -0400
Subject: [PATCH 4/9] Update heatmap of sensor row count
---
rules/reports.smk | 4 +-
...atmap_sensor_row_count_per_time_segment.py | 116 +++++++++---------
2 files changed, 62 insertions(+), 58 deletions(-)
diff --git a/rules/reports.smk b/rules/reports.smk
index caed8993..a2c24468 100644
--- a/rules/reports.smk
+++ b/rules/reports.smk
@@ -35,7 +35,9 @@ rule heatmap_sensor_row_count_per_time_segment:
participant_file = "data/external/participant_files/{pid}.yaml",
time_segments_labels = "data/interim/time_segments/{pid}_time_segments_labels.csv"
params:
- pid = "{pid}"
+ pid = "{pid}",
+ sensor_names = config["HEATMAP_SENSOR_ROW_COUNT_PER_TIME_SEGMENT"]["SENSORS"],
+ time_segments_type = config["TIME_SEGMENTS"]["TYPE"]
output:
"reports/interim/{pid}/heatmap_sensor_row_count_per_time_segment.html"
script:
diff --git a/src/visualization/heatmap_sensor_row_count_per_time_segment.py b/src/visualization/heatmap_sensor_row_count_per_time_segment.py
index 6b62e6e1..df89a5dc 100644
--- a/src/visualization/heatmap_sensor_row_count_per_time_segment.py
+++ b/src/visualization/heatmap_sensor_row_count_per_time_segment.py
@@ -1,89 +1,91 @@
import pandas as pd
import numpy as np
-import plotly.graph_objects as go
+import plotly.express as px
from importlib import util
from pathlib import Path
import yaml
-
-def getRowCountHeatmap(data_for_plot, scaled_data_for_plot, pid, time_segment, html_file):
-
- fig = go.Figure(data=go.Heatmap(z=scaled_data_for_plot.values.tolist(),
- x=data_for_plot.columns,
- y=data_for_plot.index,
- hovertext=data_for_plot.values.tolist(),
- hovertemplate="Segment start: %{x}
Sensor: %{y}
Row count: %{hovertext}
y-axis shows the included sensors.
x-axis shows the start (date and time) of a time segment.
z-axis (color) shows row count per sensor per segment instance.")
- fig["layout"].update(margin=dict(t=160))
-
- html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
-
-
-
-
# import filter_data_by_segment from src/features/utils/utils.py
spec = util.spec_from_file_location("util", str(Path(snakemake.scriptdir).parent / "features" / "utils" / "utils.py"))
mod = util.module_from_spec(spec)
spec.loader.exec_module(mod)
filter_data_by_segment = getattr(mod, "filter_data_by_segment")
+def getRowCountHeatmap(data_for_plot, pid, time_segment, html_file):
+
+ fig = px.timeline(data_for_plot,
+ x_start="local_segment_start_datetime",
+ x_end="local_segment_end_datetime",
+ y="sensor",
+ color="scaled_value",
+ color_continuous_scale="Peach", #"Viridis",
+ opacity=0.7,
+ hover_data={"local_segment_start_datetime":False, "local_segment_end_datetime":False, "local_segment":True, "value":True, "scaled_value":True})
+
+ fig.update_layout(title="Heatmap of sensor row count for " + time_segment + " segments. Pid: " + pid +". Label: " + label + "
y-axis shows the included sensors.
x-axis shows time segments.
z-axis (color) shows row count per sensor per segment instance.",
+ xaxis=dict(side="bottom", title="Time Segments"),
+ yaxis=dict(side="left", title="Sensors"),
+ margin=dict(t=160))
+
+ html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
+
+ return html_file
-
-phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], index_col=["local_segment_start_datetime"], parse_dates=["local_segment_start_datetime"])
-# make sure the phone_data_yield file contains "phone_data_yield_rapids_ratiovalidyieldedminutes" and "phone_data_yield_rapids_ratiovalidyieldedhours" columns
-if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
- raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
-phone_data_yield = phone_data_yield[["local_segment_label", "phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]]
-
-time_segments = pd.read_csv(snakemake.input["time_segments_labels"], header=0)["label"]
-pid = snakemake.params["pid"]
-
with open(snakemake.input["participant_file"], "r", encoding="utf-8") as f:
participant_file = yaml.safe_load(f)
label = participant_file["PHONE"]["LABEL"]
-sensor_names = []
-sensors_row_count = dict(zip(time_segments, [pd.DataFrame()] * len(time_segments)))
+pid = snakemake.params["pid"]
+sensor_names = [sensor_name.lower() for sensor_name in snakemake.params["sensor_names"]]
+time_segments_type = snakemake.params["time_segments_type"]
+time_segments_labels = pd.read_csv(snakemake.input["time_segments_labels"], header=0)["label"]
-for sensor_path in snakemake.input["all_sensors"]:
+phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], index_col=["local_segment"], parse_dates=["local_segment_start_datetime", "local_segment_end_datetime"]) #index_col=["local_segment_start_datetime"],
+
+# make sure the phone_data_yield file contains "phone_data_yield_rapids_ratiovalidyieldedminutes" and "phone_data_yield_rapids_ratiovalidyieldedhours" columns
+if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
+ raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
+
+# extract row count
+sensors_row_count = pd.DataFrame()
+for sensor_path, sensor_name in zip(snakemake.input["all_sensors"], sensor_names):
sensor_data = pd.read_csv(sensor_path, usecols=["assigned_segments"])
- sensor_name = sensor_path.split("/")[-1].replace("_with_datetime.csv", "")
- sensor_names.append(sensor_name)
-
+
+ sensor_row_count = pd.DataFrame()
if not sensor_data.empty:
- for time_segment in time_segments:
+ for time_segment in time_segments_labels:
sensor_data_per_segment = filter_data_by_segment(sensor_data, time_segment)
if not sensor_data_per_segment.empty:
- # extract local start datetime of the segment from "local_segment" column
- sensor_data_per_segment["local_segment_start_datetime"] = pd.to_datetime(sensor_data_per_segment["local_segment"].apply(lambda x: x.split("#")[1].split(",")[0]))
- sensor_row_count = sensor_data_per_segment.groupby("local_segment_start_datetime")[["local_segment"]].count().rename(columns={"local_segment": sensor_name})
- sensors_row_count[time_segment] = pd.concat([sensors_row_count[time_segment], sensor_row_count], axis=1, sort=False)
+ sensor_row_count = pd.concat([sensor_row_count, sensor_data_per_segment.groupby(["local_segment"])[["local_segment"]].count().rename(columns={"local_segment": sensor_name})], axis=0, sort=False)
+ sensors_row_count = pd.concat([sensors_row_count, sensor_row_count], axis=1, sort=False)
+
+sensors_row_count.index.name = "local_segment"
+sensors_row_count.index = sensors_row_count.index.str.replace(r"_RR\d+SS", "")
+data_for_plot = phone_data_yield.rename(columns={"phone_data_yield_rapids_ratiovalidyieldedminutes": "ratiovalidyieldedminutes","phone_data_yield_rapids_ratiovalidyieldedhours": "ratiovalidyieldedhours"}).merge(sensors_row_count, how="left", left_index=True, right_index=True).reset_index()
+
+
+if time_segments_type == "FREQUENCY":
+ data_for_plot["local_segment_label"] = data_for_plot["local_segment_label"].str[:-4]
+elif time_segments_type == "EVENT":
+ data_for_plot["local_segment_label"] = "event"
-# add phone data yield features and plot heatmap
-html_file = open(snakemake.output[0], "a", encoding="utf-8")
sensor_names.extend(["ratiovalidyieldedminutes", "ratiovalidyieldedhours"])
-for time_segment in time_segments:
- if not phone_data_yield.empty:
- phone_data_yield_per_segment = phone_data_yield[phone_data_yield["local_segment_label"] == time_segment].rename(columns={"phone_data_yield_rapids_ratiovalidyieldedminutes": "ratiovalidyieldedminutes","phone_data_yield_rapids_ratiovalidyieldedhours": "ratiovalidyieldedhours"}).round(3)
- if not phone_data_yield_per_segment.empty:
- sensors_row_count[time_segment] = pd.concat([sensors_row_count[time_segment], phone_data_yield_per_segment], axis=1, sort=True)
-
- # consider all the sensors
- data_for_plot = sensors_row_count[time_segment].transpose().reindex(pd.Index(sensor_names))
-
- if data_for_plot.empty:
+html_file = open(snakemake.output[0], "a", encoding="utf-8")
+for time_segment in set(data_for_plot["local_segment_label"]):
+ if not data_for_plot.empty:
+ data_for_plot_per_segment = data_for_plot[data_for_plot["local_segment_label"] == time_segment]
+ if data_for_plot_per_segment.empty:
html_file.write("There are no records of selected sensors in database for " + time_segment + " segments. Pid: " + pid + ". Label: " + label + ".
")
else:
+ data_for_plot_per_segment = data_for_plot_per_segment.reindex(columns=["local_segment", "local_segment_start_datetime", "local_segment_end_datetime"] + sensor_names).set_index(["local_segment", "local_segment_start_datetime", "local_segment_end_datetime"])
# except for phone data yield sensor, scale each sensor (row) to the range of [0, 1]
- scaled_data_for_plot = data_for_plot.copy()
- scaled_data_for_plot.loc[sensor_names[:-2]] = scaled_data_for_plot.fillna(np.nan).loc[sensor_names[:-2]].apply(lambda x: (x - np.nanmin(x)) / (np.nanmax(x) - np.nanmin(x)) if np.nanmax(x) != np.nanmin(x) else (x / np.nanmin(x)), axis=1)
-
- getRowCountHeatmap(data_for_plot, scaled_data_for_plot, pid, time_segment, html_file)
+ scaled_data_for_plot_per_segment = data_for_plot_per_segment.copy()
+ scaled_data_for_plot_per_segment[sensor_names[:-2]] = scaled_data_for_plot_per_segment.fillna(np.nan)[sensor_names[:-2]].apply(lambda x: (x - np.nanmin(x)) / (np.nanmax(x) - np.nanmin(x)) if np.nanmax(x) != np.nanmin(x) else (x / np.nanmin(x)), axis=0)
+ data_for_plot_processed = pd.concat([data_for_plot_per_segment.stack(dropna=False).to_frame("value"), scaled_data_for_plot_per_segment.stack(dropna=False).to_frame("scaled_value")], axis=1).reset_index().rename(columns={"level_3": "sensor"})
+ data_for_plot_processed[["value", "scaled_value"]] = data_for_plot_processed[["value", "scaled_value"]].round(3).clip(upper=1)
+ getRowCountHeatmap(data_for_plot_processed, pid, time_segment, html_file)
html_file.close()
From cefcb0635b4959d588da6c54f0b454f7d04a471a Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Fri, 25 Jun 2021 16:48:55 -0400
Subject: [PATCH 5/9] Update heatmap of recorded phone sensors
---
rules/reports.smk | 3 +-
.../heatmap_feature_correlation_matrix.py | 2 +-
..._yield_per_participant_per_time_segment.py | 3 +-
...map_sensors_per_minute_per_time_segment.py | 123 +--
.../histogram_phone_data_yield.py | 2 +-
src/visualization/plotly_color_utils.py | 835 ------------------
6 files changed, 68 insertions(+), 900 deletions(-)
delete mode 100644 src/visualization/plotly_color_utils.py
diff --git a/rules/reports.smk b/rules/reports.smk
index a2c24468..bd1b448c 100644
--- a/rules/reports.smk
+++ b/rules/reports.smk
@@ -14,7 +14,8 @@ rule heatmap_sensors_per_minute_per_time_segment:
participant_file = "data/external/participant_files/{pid}.yaml",
time_segments_labels = "data/interim/time_segments/{pid}_time_segments_labels.csv"
params:
- pid = "{pid}"
+ pid = "{pid}",
+ time_segments_type = config["TIME_SEGMENTS"]["TYPE"]
output:
"reports/interim/{pid}/heatmap_sensors_per_minute_per_time_segment.html"
script:
diff --git a/src/visualization/heatmap_feature_correlation_matrix.py b/src/visualization/heatmap_feature_correlation_matrix.py
index 5deff5c5..1a42d447 100644
--- a/src/visualization/heatmap_feature_correlation_matrix.py
+++ b/src/visualization/heatmap_feature_correlation_matrix.py
@@ -26,7 +26,7 @@ features = pd.read_csv(snakemake.input["all_sensor_features"])
if time_segments_type == "FREQUENCY":
- features["local_segment_label"] = features["local_segment_label"].str.split("\d+", expand=True, n=1)[0]
+ features["local_segment_label"] = features["local_segment_label"].str.replace(r"[0-9]{4}", "")
if time_segments_type == "EVENT":
features["local_segment_label"] = "event"
diff --git a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
index 34897d6e..1e6f79fc 100644
--- a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
+++ b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
@@ -1,4 +1,3 @@
-from plotly_color_utils import sample_colorscale
import pandas as pd
import plotly.express as px
import yaml
@@ -59,7 +58,7 @@ time_segments = pd.read_csv(snakemake.input["time_segments_file"])["label"].uniq
phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], parse_dates=["local_segment_start_datetime", "local_segment_end_datetime"]).sort_values(by=["pid", "local_segment_start_datetime"])
if time_segments_type == "FREQUENCY":
- phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.split("\d+", expand=True, n=1)[0]
+ phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.replace(r"[0-9]{4}", "")
html_file = open(snakemake.output[0], "w", encoding="utf-8")
if phone_data_yield.empty:
diff --git a/src/visualization/heatmap_sensors_per_minute_per_time_segment.py b/src/visualization/heatmap_sensors_per_minute_per_time_segment.py
index dd524322..7c69a756 100644
--- a/src/visualization/heatmap_sensors_per_minute_per_time_segment.py
+++ b/src/visualization/heatmap_sensors_per_minute_per_time_segment.py
@@ -5,6 +5,11 @@ from importlib import util
from pathlib import Path
import yaml
+# import filter_data_by_segment from src/features/utils/utils.py
+spec = util.spec_from_file_location("util", str(Path(snakemake.scriptdir).parent / "features" / "utils" / "utils.py"))
+mod = util.module_from_spec(spec)
+spec.loader.exec_module(mod)
+filter_data_by_segment = getattr(mod, "filter_data_by_segment")
def colors2colorscale(colors):
colorscale = []
@@ -16,85 +21,83 @@ def colors2colorscale(colors):
colorscale.append([1, colors[i]])
return colorscale
-def getSensorsPerMinPerSegmentHeatmap(phone_data_yield, pid, time_segment, html_file):
-
- x_axis_labels = [pd.Timedelta(minutes=x) for x in phone_data_yield.columns]
+def getDataForPlot(phone_data_yield_per_segment):
+ # calculate the length (in minute) of per segment instance
+ phone_data_yield_per_segment["length"] = phone_data_yield_per_segment["timestamps_segment"].str.split(",").apply(lambda x: int((int(x[1])-int(x[0])) / (1000 * 60)))
+ # calculate the number of sensors logged at least one row of data per minute.
+ phone_data_yield_per_segment = phone_data_yield_per_segment.groupby(["local_segment", "length", "local_date", "local_hour", "local_minute"])[["sensor", "local_date_time"]].max().reset_index()
+ # extract local start datetime of the segment from "local_segment" column
+ phone_data_yield_per_segment["local_segment_start_datetimes"] = pd.to_datetime(phone_data_yield_per_segment["local_segment"].apply(lambda x: x.split("#")[1].split(",")[0]))
+ # calculate the number of minutes after local start datetime of the segment
+ phone_data_yield_per_segment["minutes_after_segment_start"] = ((phone_data_yield_per_segment["local_date_time"] - phone_data_yield_per_segment["local_segment_start_datetimes"]) / pd.Timedelta(minutes=1)).astype("int")
- fig = go.Figure(data=go.Heatmap(z=phone_data_yield.values.tolist(),
- x=x_axis_labels,
- y=phone_data_yield.index,
- zmin=0, zmax=16,
- colorscale=colors2colorscale(colors),
- colorbar=dict(thickness=25, tickvals=[1/2 + x for x in range(16)],ticktext=[x for x in range(16)])))
+ # impute missing rows with 0
+ columns_for_full_index = phone_data_yield_per_segment[["local_segment_start_datetimes", "length"]].drop_duplicates(keep="first")
+ columns_for_full_index = columns_for_full_index.apply(lambda row: [[row["local_segment_start_datetimes"], x] for x in range(row["length"] + 1)], axis=1)
+ full_index = []
+ for columns in columns_for_full_index:
+ full_index = full_index + columns
+ full_index = pd.MultiIndex.from_tuples(full_index, names=("local_segment_start_datetimes", "minutes_after_segment_start"))
+ phone_data_yield_per_segment = phone_data_yield_per_segment.set_index(["local_segment_start_datetimes", "minutes_after_segment_start"]).reindex(full_index).reset_index().fillna(0)
- fig.update_layout(title="Number of sensors with any data per minute for " + time_segment + " segments. Pid: "+pid+". Label: " + label + "
y-axis shows the start (date and time) of a time segment.
x-axis shows the time since the start of the time segment.
z-axis (color) shows how many sensors logged at least one row of data per minute.")
- fig["layout"].update(margin=dict(t=160))
-
- html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
-
-
-
-
-
-# import filter_data_by_segment from src/features/utils/utils.py
-spec = util.spec_from_file_location("util", str(Path(snakemake.scriptdir).parent / "features" / "utils" / "utils.py"))
-mod = util.module_from_spec(spec)
-spec.loader.exec_module(mod)
-filter_data_by_segment = getattr(mod, "filter_data_by_segment")
-
-
-
-
-
-
+ # transpose the dataframe per local start datetime of the segment and discard the useless index layer
+ phone_data_yield_per_segment = phone_data_yield_per_segment.groupby("local_segment_start_datetimes")[["minutes_after_segment_start", "sensor"]].apply(lambda x: x.set_index("minutes_after_segment_start").transpose())
+ phone_data_yield_per_segment.index = phone_data_yield_per_segment.index.get_level_values("local_segment_start_datetimes")
+ return phone_data_yield_per_segment
+def getSensorsPerMinPerSegmentHeatmap(phone_data_yield, pid, label, time_segment, html_file):
+ if phone_data_yield.empty:
+ html_file.write("There is no sensor data of " + time_segment + " segments for " + pid + " (pid) and " + label + " (label).
")
+ else:
+ phone_data_yield.sort_index(inplace=True)
+ x_axis_labels = [pd.Timedelta(minutes=x) for x in phone_data_yield.columns]
+
+ fig = go.Figure(data=go.Heatmap(z=phone_data_yield.values.tolist(),
+ x=x_axis_labels,
+ y=phone_data_yield.index,
+ zmin=0, zmax=16,
+ colorscale=colors2colorscale(colors),
+ colorbar=dict(thickness=25, tickvals=[1/2 + x for x in range(16)],ticktext=[x for x in range(16)])))
+
+ fig.update_layout(title="Number of sensors with any data per minute for " + time_segment + " segments. Pid: "+pid+". Label: " + label + "
y-axis shows the start (date and time) of a time segment.
x-axis shows the time since the start of the time segment.
z-axis (color) shows how many sensors logged at least one row of data per minute.")
+ fig["layout"].update(margin=dict(t=160))
+
+ html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
+ return
colors = ["red", "#3D0751", "#423176", "#414381", "#3F5688", "#42678B", "#42768C", "#45868B", "#4A968A", "#53A485", "#5FB57E", "#76C170", "#91CF63", "#B4DA55", "#D9E152", "#F8E755", "#DEE00F"]
pid = snakemake.params["pid"]
-time_segments_labels = pd.read_csv(snakemake.input["time_segments_labels"], header=0)
+time_segments_type = snakemake.params["time_segments_type"]
+time_segments_labels = pd.read_csv(snakemake.input["time_segments_labels"])
with open(snakemake.input["participant_file"], "r", encoding="utf-8") as f:
participant_file = yaml.safe_load(f)
label = participant_file["PHONE"]["LABEL"]
phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], parse_dates=["local_date_time"])
+if time_segments_type == "FREQUENCY":
+ phone_data_yield["assigned_segments"] = phone_data_yield["assigned_segments"].str.replace(r"[0-9]{4}#", "#")
+ time_segments_labels["label"] = time_segments_labels["label"].str.replace(r"[0-9]{4}", "")
+if time_segments_type == "PERIODIC":
+ phone_data_yield["assigned_segments"] = phone_data_yield["assigned_segments"].str.replace(r"_RR\d+SS#", "#")
+ time_segments_labels["label"] = time_segments_labels["label"].str.replace(r"_RR\d+SS", "")
html_file = open(snakemake.output[0], "a", encoding="utf-8")
if phone_data_yield.empty:
html_file.write("There is no sensor data for " + pid + " (pid) and " + label + " (label).")
else:
- for time_segment in time_segments_labels["label"]:
+ data_for_plot = pd.DataFrame()
+ for time_segment in set(time_segments_labels["label"]):
phone_data_yield_per_segment = filter_data_by_segment(phone_data_yield, time_segment)
-
- if phone_data_yield_per_segment.empty:
- html_file.write("There is no sensor data of " + time_segment + " segments for " + pid + " (pid) and " + label + " (label).
")
- else:
- # calculate the length (in minute) of per segment instance
- phone_data_yield_per_segment["length"] = phone_data_yield_per_segment["timestamps_segment"].str.split(",").apply(lambda x: int((int(x[1])-int(x[0])) / (1000 * 60)))
- # calculate the number of sensors logged at least one row of data per minute.
- phone_data_yield_per_segment = phone_data_yield_per_segment.groupby(["local_segment", "length", "local_date", "local_hour", "local_minute"])[["sensor", "local_date_time"]].max().reset_index()
- # extract local start datetime of the segment from "local_segment" column
- phone_data_yield_per_segment["local_segment_start_datetimes"] = pd.to_datetime(phone_data_yield_per_segment["local_segment"].apply(lambda x: x.split("#")[1].split(",")[0]))
- # calculate the number of minutes after local start datetime of the segment
- phone_data_yield_per_segment["minutes_after_segment_start"] = ((phone_data_yield_per_segment["local_date_time"] - phone_data_yield_per_segment["local_segment_start_datetimes"]) / pd.Timedelta(minutes=1)).astype("int")
-
- # impute missing rows with 0
- columns_for_full_index = phone_data_yield_per_segment[["local_segment_start_datetimes", "length"]].drop_duplicates(keep="first")
- columns_for_full_index = columns_for_full_index.apply(lambda row: [[row["local_segment_start_datetimes"], x] for x in range(row["length"] + 1)], axis=1)
- full_index = []
- for columns in columns_for_full_index:
- full_index = full_index + columns
- full_index = pd.MultiIndex.from_tuples(full_index, names=("local_segment_start_datetimes", "minutes_after_segment_start"))
- phone_data_yield_per_segment = phone_data_yield_per_segment.set_index(["local_segment_start_datetimes", "minutes_after_segment_start"]).reindex(full_index).reset_index().fillna(0)
-
- # transpose the dataframe per local start datetime of the segment and discard the useless index layer
- phone_data_yield_per_segment = phone_data_yield_per_segment.groupby("local_segment_start_datetimes")[["minutes_after_segment_start", "sensor"]].apply(lambda x: x.set_index("minutes_after_segment_start").transpose())
- phone_data_yield_per_segment.index = phone_data_yield_per_segment.index.get_level_values("local_segment_start_datetimes")
-
- # get heatmap
- getSensorsPerMinPerSegmentHeatmap(phone_data_yield_per_segment, pid, time_segment, html_file)
-
+ if not phone_data_yield_per_segment.empty:
+ data_for_plot_per_segment = getDataForPlot(phone_data_yield_per_segment)
+ if time_segments_type == "EVENT":
+ data_for_plot = pd.concat([data_for_plot, data_for_plot_per_segment], axis=0)
+ else:
+ getSensorsPerMinPerSegmentHeatmap(data_for_plot_per_segment, pid, label, time_segment, html_file)
+ if time_segments_type == "EVENT":
+ getSensorsPerMinPerSegmentHeatmap(data_for_plot, pid, label, "event", html_file)
html_file.close()
diff --git a/src/visualization/histogram_phone_data_yield.py b/src/visualization/histogram_phone_data_yield.py
index a8b90d77..998d4fd9 100644
--- a/src/visualization/histogram_phone_data_yield.py
+++ b/src/visualization/histogram_phone_data_yield.py
@@ -6,7 +6,7 @@ time_segments_type = snakemake.params["time_segments_type"]
phone_data_yield = pd.read_csv(snakemake.input[0])
if time_segments_type == "FREQUENCY":
- phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.split("\d+", expand=True, n=1)[0]
+ phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.replace(r"[0-9]{4}", "")
if time_segments_type == "EVENT":
phone_data_yield["local_segment_label"] = "event"
diff --git a/src/visualization/plotly_color_utils.py b/src/visualization/plotly_color_utils.py
deleted file mode 100644
index 12ffb279..00000000
--- a/src/visualization/plotly_color_utils.py
+++ /dev/null
@@ -1,835 +0,0 @@
-"""
-Source: https://github.com/plotly/plotly.py/pull/3136
-=====
-colors
-=====
-Functions that manipulate colors and arrays of colors.
------
-There are three basic types of color types: rgb, hex and tuple:
-rgb - An rgb color is a string of the form 'rgb(a,b,c)' where a, b and c are
-integers between 0 and 255 inclusive.
-hex - A hex color is a string of the form '#xxxxxx' where each x is a
-character that belongs to the set [0,1,2,3,4,5,6,7,8,9,a,b,c,d,e,f]. This is
-just the set of characters used in the hexadecimal numeric system.
-tuple - A tuple color is a 3-tuple of the form (a,b,c) where a, b and c are
-floats between 0 and 1 inclusive.
------
-Colormaps and Colorscales:
-A colormap or a colorscale is a correspondence between values - Pythonic
-objects such as strings and floats - to colors.
-There are typically two main types of colormaps that exist: numerical and
-categorical colormaps.
-Numerical:
-----------
-Numerical colormaps are used when the coloring column being used takes a
-spectrum of values or numbers.
-A classic example from the Plotly library:
-```
-rainbow_colorscale = [
- [0, 'rgb(150,0,90)'], [0.125, 'rgb(0,0,200)'],
- [0.25, 'rgb(0,25,255)'], [0.375, 'rgb(0,152,255)'],
- [0.5, 'rgb(44,255,150)'], [0.625, 'rgb(151,255,0)'],
- [0.75, 'rgb(255,234,0)'], [0.875, 'rgb(255,111,0)'],
- [1, 'rgb(255,0,0)']
-]
-```
-Notice that this colorscale is a list of lists with each inner list containing
-a number and a color. These left hand numbers in the nested lists go from 0 to
-1, and they are like pointers tell you when a number is mapped to a specific
-color.
-If you have a column of numbers `col_num` that you want to plot, and you know
-```
-min(col_num) = 0
-max(col_num) = 100
-```
-then if you pull out the number `12.5` in the list and want to figure out what
-color the corresponding chart element (bar, scatter plot, etc) is going to be,
-you'll figure out that proportionally 12.5 to 100 is the same as 0.125 to 1.
-So, the point will be mapped to 'rgb(0,0,200)'.
-All other colors between the pinned values in a colorscale are linearly
-interpolated.
-Categorical:
-------------
-Alternatively, a categorical colormap is used to assign a specific value in a
-color column to a specific color everytime it appears in the dataset.
-A column of strings in a panadas.dataframe that is chosen to serve as the
-color index would naturally use a categorical colormap. However, you can
-choose to use a categorical colormap with a column of numbers.
-Be careful! If you have a lot of unique numbers in your color column you will
-end up with a colormap that is massive and may slow down graphing performance.
-"""
-from __future__ import absolute_import
-
-import decimal
-from numbers import Number
-import six
-
-from _plotly_utils import exceptions
-
-
-# Built-in qualitative color sequences and sequential,
-# diverging and cyclical color scales.
-
-DEFAULT_PLOTLY_COLORS = [
- "rgb(31, 119, 180)",
- "rgb(255, 127, 14)",
- "rgb(44, 160, 44)",
- "rgb(214, 39, 40)",
- "rgb(148, 103, 189)",
- "rgb(140, 86, 75)",
- "rgb(227, 119, 194)",
- "rgb(127, 127, 127)",
- "rgb(188, 189, 34)",
- "rgb(23, 190, 207)",
-]
-
-PLOTLY_SCALES = {
- "Greys": [[0, "rgb(0,0,0)"], [1, "rgb(255,255,255)"]],
- "YlGnBu": [
- [0, "rgb(8,29,88)"],
- [0.125, "rgb(37,52,148)"],
- [0.25, "rgb(34,94,168)"],
- [0.375, "rgb(29,145,192)"],
- [0.5, "rgb(65,182,196)"],
- [0.625, "rgb(127,205,187)"],
- [0.75, "rgb(199,233,180)"],
- [0.875, "rgb(237,248,217)"],
- [1, "rgb(255,255,217)"],
- ],
- "Greens": [
- [0, "rgb(0,68,27)"],
- [0.125, "rgb(0,109,44)"],
- [0.25, "rgb(35,139,69)"],
- [0.375, "rgb(65,171,93)"],
- [0.5, "rgb(116,196,118)"],
- [0.625, "rgb(161,217,155)"],
- [0.75, "rgb(199,233,192)"],
- [0.875, "rgb(229,245,224)"],
- [1, "rgb(247,252,245)"],
- ],
- "YlOrRd": [
- [0, "rgb(128,0,38)"],
- [0.125, "rgb(189,0,38)"],
- [0.25, "rgb(227,26,28)"],
- [0.375, "rgb(252,78,42)"],
- [0.5, "rgb(253,141,60)"],
- [0.625, "rgb(254,178,76)"],
- [0.75, "rgb(254,217,118)"],
- [0.875, "rgb(255,237,160)"],
- [1, "rgb(255,255,204)"],
- ],
- "Bluered": [[0, "rgb(0,0,255)"], [1, "rgb(255,0,0)"]],
- # modified RdBu based on
- # www.sandia.gov/~kmorel/documents/ColorMaps/ColorMapsExpanded.pdf
- "RdBu": [
- [0, "rgb(5,10,172)"],
- [0.35, "rgb(106,137,247)"],
- [0.5, "rgb(190,190,190)"],
- [0.6, "rgb(220,170,132)"],
- [0.7, "rgb(230,145,90)"],
- [1, "rgb(178,10,28)"],
- ],
- # Scale for non-negative numeric values
- "Reds": [
- [0, "rgb(220,220,220)"],
- [0.2, "rgb(245,195,157)"],
- [0.4, "rgb(245,160,105)"],
- [1, "rgb(178,10,28)"],
- ],
- # Scale for non-positive numeric values
- "Blues": [
- [0, "rgb(5,10,172)"],
- [0.35, "rgb(40,60,190)"],
- [0.5, "rgb(70,100,245)"],
- [0.6, "rgb(90,120,245)"],
- [0.7, "rgb(106,137,247)"],
- [1, "rgb(220,220,220)"],
- ],
- "Picnic": [
- [0, "rgb(0,0,255)"],
- [0.1, "rgb(51,153,255)"],
- [0.2, "rgb(102,204,255)"],
- [0.3, "rgb(153,204,255)"],
- [0.4, "rgb(204,204,255)"],
- [0.5, "rgb(255,255,255)"],
- [0.6, "rgb(255,204,255)"],
- [0.7, "rgb(255,153,255)"],
- [0.8, "rgb(255,102,204)"],
- [0.9, "rgb(255,102,102)"],
- [1, "rgb(255,0,0)"],
- ],
- "Rainbow": [
- [0, "rgb(150,0,90)"],
- [0.125, "rgb(0,0,200)"],
- [0.25, "rgb(0,25,255)"],
- [0.375, "rgb(0,152,255)"],
- [0.5, "rgb(44,255,150)"],
- [0.625, "rgb(151,255,0)"],
- [0.75, "rgb(255,234,0)"],
- [0.875, "rgb(255,111,0)"],
- [1, "rgb(255,0,0)"],
- ],
- "Portland": [
- [0, "rgb(12,51,131)"],
- [0.25, "rgb(10,136,186)"],
- [0.5, "rgb(242,211,56)"],
- [0.75, "rgb(242,143,56)"],
- [1, "rgb(217,30,30)"],
- ],
- "Jet": [
- [0, "rgb(0,0,131)"],
- [0.125, "rgb(0,60,170)"],
- [0.375, "rgb(5,255,255)"],
- [0.625, "rgb(255,255,0)"],
- [0.875, "rgb(250,0,0)"],
- [1, "rgb(128,0,0)"],
- ],
- "Hot": [
- [0, "rgb(0,0,0)"],
- [0.3, "rgb(230,0,0)"],
- [0.6, "rgb(255,210,0)"],
- [1, "rgb(255,255,255)"],
- ],
- "Blackbody": [
- [0, "rgb(0,0,0)"],
- [0.2, "rgb(230,0,0)"],
- [0.4, "rgb(230,210,0)"],
- [0.7, "rgb(255,255,255)"],
- [1, "rgb(160,200,255)"],
- ],
- "Earth": [
- [0, "rgb(0,0,130)"],
- [0.1, "rgb(0,180,180)"],
- [0.2, "rgb(40,210,40)"],
- [0.4, "rgb(230,230,50)"],
- [0.6, "rgb(120,70,20)"],
- [1, "rgb(255,255,255)"],
- ],
- "Electric": [
- [0, "rgb(0,0,0)"],
- [0.15, "rgb(30,0,100)"],
- [0.4, "rgb(120,0,100)"],
- [0.6, "rgb(160,90,0)"],
- [0.8, "rgb(230,200,0)"],
- [1, "rgb(255,250,220)"],
- ],
- "Viridis": [
- [0, "#440154"],
- [0.06274509803921569, "#48186a"],
- [0.12549019607843137, "#472d7b"],
- [0.18823529411764706, "#424086"],
- [0.25098039215686274, "#3b528b"],
- [0.3137254901960784, "#33638d"],
- [0.3764705882352941, "#2c728e"],
- [0.4392156862745098, "#26828e"],
- [0.5019607843137255, "#21918c"],
- [0.5647058823529412, "#1fa088"],
- [0.6274509803921569, "#28ae80"],
- [0.6901960784313725, "#3fbc73"],
- [0.7529411764705882, "#5ec962"],
- [0.8156862745098039, "#84d44b"],
- [0.8784313725490196, "#addc30"],
- [0.9411764705882353, "#d8e219"],
- [1, "#fde725"],
- ],
- "Cividis": [
- [0.000000, "rgb(0,32,76)"],
- [0.058824, "rgb(0,42,102)"],
- [0.117647, "rgb(0,52,110)"],
- [0.176471, "rgb(39,63,108)"],
- [0.235294, "rgb(60,74,107)"],
- [0.294118, "rgb(76,85,107)"],
- [0.352941, "rgb(91,95,109)"],
- [0.411765, "rgb(104,106,112)"],
- [0.470588, "rgb(117,117,117)"],
- [0.529412, "rgb(131,129,120)"],
- [0.588235, "rgb(146,140,120)"],
- [0.647059, "rgb(161,152,118)"],
- [0.705882, "rgb(176,165,114)"],
- [0.764706, "rgb(192,177,109)"],
- [0.823529, "rgb(209,191,102)"],
- [0.882353, "rgb(225,204,92)"],
- [0.941176, "rgb(243,219,79)"],
- [1.000000, "rgb(255,233,69)"],
- ],
-}
-
-
-def color_parser(colors, function):
- """
- Takes color(s) and a function and applies the function on the color(s)
- In particular, this function identifies whether the given color object
- is an iterable or not and applies the given color-parsing function to
- the color or iterable of colors. If given an iterable, it will only be
- able to work with it if all items in the iterable are of the same type
- - rgb string, hex string or tuple
- """
- if isinstance(colors, str):
- return function(colors)
-
- if isinstance(colors, tuple) and isinstance(colors[0], Number):
- return function(colors)
-
- if hasattr(colors, "__iter__"):
- if isinstance(colors, tuple):
- new_color_tuple = tuple(function(item) for item in colors)
- return new_color_tuple
-
- else:
- new_color_list = [function(item) for item in colors]
- return new_color_list
-
-
-def validate_colors(colors, colortype="tuple"):
- """
- Validates color(s) and returns a list of color(s) of a specified type
- """
- from numbers import Number
-
- if colors is None:
- colors = DEFAULT_PLOTLY_COLORS
-
- if isinstance(colors, str):
- if colors in PLOTLY_SCALES:
- colors_list = colorscale_to_colors(PLOTLY_SCALES[colors])
- # TODO: fix _gantt.py/_scatter.py so that they can accept the
- # actual colorscale and not just a list of the first and last
- # color in the plotly colorscale. In resolving this issue we
- # will be removing the immediate line below
- colors = [colors_list[0]] + [colors_list[-1]]
- elif "rgb" in colors or "#" in colors:
- colors = [colors]
- else:
- raise exceptions.PlotlyError(
- "If your colors variable is a string, it must be a "
- "Plotly scale, an rgb color or a hex color."
- )
-
- elif isinstance(colors, tuple):
- if isinstance(colors[0], Number):
- colors = [colors]
- else:
- colors = list(colors)
-
- # convert color elements in list to tuple color
- for j, each_color in enumerate(colors):
- if "rgb" in each_color:
- each_color = color_parser(each_color, unlabel_rgb)
- for value in each_color:
- if value > 255.0:
- raise exceptions.PlotlyError(
- "Whoops! The elements in your rgb colors "
- "tuples cannot exceed 255.0."
- )
- each_color = color_parser(each_color, unconvert_from_RGB_255)
- colors[j] = each_color
-
- if "#" in each_color:
- each_color = color_parser(each_color, hex_to_rgb)
- each_color = color_parser(each_color, unconvert_from_RGB_255)
-
- colors[j] = each_color
-
- if isinstance(each_color, tuple):
- for value in each_color:
- if value > 1.0:
- raise exceptions.PlotlyError(
- "Whoops! The elements in your colors tuples "
- "cannot exceed 1.0."
- )
- colors[j] = each_color
-
- if colortype == "rgb" and not isinstance(colors, six.string_types):
- for j, each_color in enumerate(colors):
- rgb_color = color_parser(each_color, convert_to_RGB_255)
- colors[j] = color_parser(rgb_color, label_rgb)
-
- return colors
-
-
-def validate_colors_dict(colors, colortype="tuple"):
- """
- Validates dictionary of color(s)
- """
- # validate each color element in the dictionary
- for key in colors:
- if "rgb" in colors[key]:
- colors[key] = color_parser(colors[key], unlabel_rgb)
- for value in colors[key]:
- if value > 255.0:
- raise exceptions.PlotlyError(
- "Whoops! The elements in your rgb colors "
- "tuples cannot exceed 255.0."
- )
- colors[key] = color_parser(colors[key], unconvert_from_RGB_255)
-
- if "#" in colors[key]:
- colors[key] = color_parser(colors[key], hex_to_rgb)
- colors[key] = color_parser(colors[key], unconvert_from_RGB_255)
-
- if isinstance(colors[key], tuple):
- for value in colors[key]:
- if value > 1.0:
- raise exceptions.PlotlyError(
- "Whoops! The elements in your colors tuples "
- "cannot exceed 1.0."
- )
-
- if colortype == "rgb":
- for key in colors:
- colors[key] = color_parser(colors[key], convert_to_RGB_255)
- colors[key] = color_parser(colors[key], label_rgb)
-
- return colors
-
-
-def convert_colors_to_same_type(
- colors,
- colortype="rgb",
- scale=None,
- return_default_colors=False,
- num_of_defualt_colors=2,
-):
- """
- Converts color(s) to the specified color type
- Takes a single color or an iterable of colors, as well as a list of scale
- values, and outputs a 2-pair of the list of color(s) converted all to an
- rgb or tuple color type, aswell as the scale as the second element. If
- colors is a Plotly Scale name, then 'scale' will be forced to the scale
- from the respective colorscale and the colors in that colorscale will also
- be coverted to the selected colortype. If colors is None, then there is an
- option to return portion of the DEFAULT_PLOTLY_COLORS
- :param (str|tuple|list) colors: either a plotly scale name, an rgb or hex
- color, a color tuple or a list/tuple of colors
- :param (list) scale: see docs for validate_scale_values()
- :rtype (tuple) (colors_list, scale) if scale is None in the function call,
- then scale will remain None in the returned tuple
- """
- colors_list = []
-
- if colors is None and return_default_colors is True:
- colors_list = DEFAULT_PLOTLY_COLORS[0:num_of_defualt_colors]
-
- if isinstance(colors, str):
- if colors in PLOTLY_SCALES:
- colors_list = colorscale_to_colors(PLOTLY_SCALES[colors])
- if scale is None:
- scale = colorscale_to_scale(PLOTLY_SCALES[colors])
-
- elif "rgb" in colors or "#" in colors:
- colors_list = [colors]
-
- elif isinstance(colors, tuple):
- if isinstance(colors[0], Number):
- colors_list = [colors]
- else:
- colors_list = list(colors)
-
- elif isinstance(colors, list):
- colors_list = colors
-
- # validate scale
- if scale is not None:
- validate_scale_values(scale)
-
- if len(colors_list) != len(scale):
- raise exceptions.PlotlyError(
- "Make sure that the length of your scale matches the length "
- "of your list of colors which is {}.".format(len(colors_list))
- )
-
- # convert all colors to rgb
- for j, each_color in enumerate(colors_list):
- if "#" in each_color:
- each_color = color_parser(each_color, hex_to_rgb)
- each_color = color_parser(each_color, label_rgb)
- colors_list[j] = each_color
-
- elif isinstance(each_color, tuple):
- each_color = color_parser(each_color, convert_to_RGB_255)
- each_color = color_parser(each_color, label_rgb)
- colors_list[j] = each_color
-
- if colortype == "rgb":
- return (colors_list, scale)
- elif colortype == "tuple":
- for j, each_color in enumerate(colors_list):
- each_color = color_parser(each_color, unlabel_rgb)
- each_color = color_parser(each_color, unconvert_from_RGB_255)
- colors_list[j] = each_color
- return (colors_list, scale)
- else:
- raise exceptions.PlotlyError(
- "You must select either rgb or tuple for your colortype variable."
- )
-
-
-def convert_dict_colors_to_same_type(colors_dict, colortype="rgb"):
- """
- Converts a colors in a dictionary of colors to the specified color type
- :param (dict) colors_dict: a dictionary whose values are single colors
- """
- for key in colors_dict:
- if "#" in colors_dict[key]:
- colors_dict[key] = color_parser(colors_dict[key], hex_to_rgb)
- colors_dict[key] = color_parser(colors_dict[key], label_rgb)
-
- elif isinstance(colors_dict[key], tuple):
- colors_dict[key] = color_parser(colors_dict[key], convert_to_RGB_255)
- colors_dict[key] = color_parser(colors_dict[key], label_rgb)
-
- if colortype == "rgb":
- return colors_dict
- elif colortype == "tuple":
- for key in colors_dict:
- colors_dict[key] = color_parser(colors_dict[key], unlabel_rgb)
- colors_dict[key] = color_parser(colors_dict[key], unconvert_from_RGB_255)
- return colors_dict
- else:
- raise exceptions.PlotlyError(
- "You must select either rgb or tuple for your colortype variable."
- )
-
-
-def validate_scale_values(scale):
- """
- Validates scale values from a colorscale
- :param (list) scale: a strictly increasing list of floats that begins
- with 0 and ends with 1. Its usage derives from a colorscale which is
- a list of two-lists (a list with two elements) of the form
- [value, color] which are used to determine how interpolation weighting
- works between the colors in the colorscale. Therefore scale is just
- the extraction of these values from the two-lists in order
- """
- if len(scale) < 2:
- raise exceptions.PlotlyError(
- "You must input a list of scale values that has at least two values."
- )
-
- if (scale[0] != 0) or (scale[-1] != 1):
- raise exceptions.PlotlyError(
- "The first and last number in your scale must be 0.0 and 1.0 "
- "respectively."
- )
-
- if not all(x < y for x, y in zip(scale, scale[1:])):
- raise exceptions.PlotlyError(
- "'scale' must be a list that contains a strictly increasing "
- "sequence of numbers."
- )
-
-
-def validate_colorscale(colorscale):
- """Validate the structure, scale values and colors of colorscale."""
- if not isinstance(colorscale, list):
- # TODO Write tests for these exceptions
- raise exceptions.PlotlyError("A valid colorscale must be a list.")
- if not all(isinstance(innerlist, list) for innerlist in colorscale):
- raise exceptions.PlotlyError("A valid colorscale must be a list of lists.")
- colorscale_colors = colorscale_to_colors(colorscale)
- scale_values = colorscale_to_scale(colorscale)
-
- validate_scale_values(scale_values)
- validate_colors(colorscale_colors)
-
-
-def make_colorscale(colors, scale=None):
- """
- Makes a colorscale from a list of colors and a scale
- Takes a list of colors and scales and constructs a colorscale based
- on the colors in sequential order. If 'scale' is left empty, a linear-
- interpolated colorscale will be generated. If 'scale' is a specificed
- list, it must be the same legnth as colors and must contain all floats
- For documentation regarding to the form of the output, see
- https://plot.ly/python/reference/#mesh3d-colorscale
- :param (list) colors: a list of single colors
- """
- colorscale = []
-
- # validate minimum colors length of 2
- if len(colors) < 2:
- raise exceptions.PlotlyError(
- "You must input a list of colors that has at least two colors."
- )
-
- if scale is None:
- scale_incr = 1.0 / (len(colors) - 1)
- return [[i * scale_incr, color] for i, color in enumerate(colors)]
-
- else:
- if len(colors) != len(scale):
- raise exceptions.PlotlyError(
- "The length of colors and scale must be the same."
- )
-
- validate_scale_values(scale)
-
- colorscale = [list(tup) for tup in zip(scale, colors)]
- return colorscale
-
-
-def find_intermediate_color(lowcolor, highcolor, intermed, colortype="tuple"):
- """
- Returns the color at a given distance between two colors
- This function takes two color tuples, where each element is between 0
- and 1, along with a value 0 < intermed < 1 and returns a color that is
- intermed-percent from lowcolor to highcolor. If colortype is set to 'rgb',
- the function will automatically convert the rgb type to a tuple, find the
- intermediate color and return it as an rgb color.
- """
- if colortype == "rgb":
- # convert to tuple color, eg. (1, 0.45, 0.7)
- lowcolor = unlabel_rgb(lowcolor)
- highcolor = unlabel_rgb(highcolor)
-
- diff_0 = float(highcolor[0] - lowcolor[0])
- diff_1 = float(highcolor[1] - lowcolor[1])
- diff_2 = float(highcolor[2] - lowcolor[2])
-
- inter_med_tuple = (
- lowcolor[0] + intermed * diff_0,
- lowcolor[1] + intermed * diff_1,
- lowcolor[2] + intermed * diff_2,
- )
-
- if colortype == "rgb":
- # back to an rgb string, e.g. rgb(30, 20, 10)
- inter_med_rgb = label_rgb(inter_med_tuple)
- return inter_med_rgb
-
- return inter_med_tuple
-
-
-def unconvert_from_RGB_255(colors):
- """
- Return a tuple where each element gets divided by 255
- Takes a (list of) color tuple(s) where each element is between 0 and
- 255. Returns the same tuples where each tuple element is normalized to
- a value between 0 and 1
- """
- return (colors[0] / (255.0), colors[1] / (255.0), colors[2] / (255.0))
-
-
-def convert_to_RGB_255(colors):
- """
- Multiplies each element of a triplet by 255
- Each coordinate of the color tuple is rounded to the nearest float and
- then is turned into an integer. If a number is of the form x.5, then
- if x is odd, the number rounds up to (x+1). Otherwise, it rounds down
- to just x. This is the way rounding works in Python 3 and in current
- statistical analysis to avoid rounding bias
- :param (list) rgb_components: grabs the three R, G and B values to be
- returned as computed in the function
- """
- rgb_components = []
-
- for component in colors:
- rounded_num = decimal.Decimal(str(component * 255.0)).quantize(
- decimal.Decimal("1"), rounding=decimal.ROUND_HALF_EVEN
- )
- # convert rounded number to an integer from 'Decimal' form
- rounded_num = int(rounded_num)
- rgb_components.append(rounded_num)
-
- return (rgb_components[0], rgb_components[1], rgb_components[2])
-
-
-def n_colors(lowcolor, highcolor, n_colors, colortype="tuple"):
- """
- Splits a low and high color into a list of n_colors colors in it
- Accepts two color tuples and returns a list of n_colors colors
- which form the intermediate colors between lowcolor and highcolor
- from linearly interpolating through RGB space. If colortype is 'rgb'
- the function will return a list of colors in the same form.
- """
- if colortype == "rgb":
- # convert to tuple
- lowcolor = unlabel_rgb(lowcolor)
- highcolor = unlabel_rgb(highcolor)
-
- diff_0 = float(highcolor[0] - lowcolor[0])
- incr_0 = diff_0 / (n_colors - 1)
- diff_1 = float(highcolor[1] - lowcolor[1])
- incr_1 = diff_1 / (n_colors - 1)
- diff_2 = float(highcolor[2] - lowcolor[2])
- incr_2 = diff_2 / (n_colors - 1)
- list_of_colors = []
-
- for index in range(n_colors):
- new_tuple = (
- lowcolor[0] + (index * incr_0),
- lowcolor[1] + (index * incr_1),
- lowcolor[2] + (index * incr_2),
- )
- list_of_colors.append(new_tuple)
-
- if colortype == "rgb":
- # back to an rgb string
- list_of_colors = color_parser(list_of_colors, label_rgb)
-
- return list_of_colors
-
-
-def label_rgb(colors):
- """
- Takes tuple (a, b, c) and returns an rgb color 'rgb(a, b, c)'
- """
- return "rgb(%s, %s, %s)" % (colors[0], colors[1], colors[2])
-
-
-def unlabel_rgb(colors):
- """
- Takes rgb color(s) 'rgb(a, b, c)' and returns tuple(s) (a, b, c)
- This function takes either an 'rgb(a, b, c)' color or a list of
- such colors and returns the color tuples in tuple(s) (a, b, c)
- """
- str_vals = ""
- for index in range(len(colors)):
- try:
- float(colors[index])
- str_vals = str_vals + colors[index]
- except ValueError:
- if colors[index] == "," or colors[index] == ".":
- str_vals = str_vals + colors[index]
-
- str_vals = str_vals + ","
- numbers = []
- str_num = ""
- for char in str_vals:
- if char != ",":
- str_num = str_num + char
- else:
- numbers.append(float(str_num))
- str_num = ""
- return (numbers[0], numbers[1], numbers[2])
-
-
-def hex_to_rgb(value):
- """
- Calculates rgb values from a hex color code.
- :param (string) value: Hex color string
- :rtype (tuple) (r_value, g_value, b_value): tuple of rgb values
- """
- value = value.lstrip("#")
- hex_total_length = len(value)
- rgb_section_length = hex_total_length // 3
- return tuple(
- int(value[i : i + rgb_section_length], 16)
- for i in range(0, hex_total_length, rgb_section_length)
- )
-
-
-def colorscale_to_colors(colorscale):
- """
- Extracts the colors from colorscale as a list
- """
- color_list = []
- for item in colorscale:
- color_list.append(item[1])
- return color_list
-
-
-def colorscale_to_scale(colorscale):
- """
- Extracts the interpolation scale values from colorscale as a list
- """
- scale_list = []
- for item in colorscale:
- scale_list.append(item[0])
- return scale_list
-
-
-def convert_colorscale_to_rgb(colorscale):
- """
- Converts the colors in a colorscale to rgb colors
- A colorscale is an array of arrays, each with a numeric value as the
- first item and a color as the second. This function specifically is
- converting a colorscale with tuple colors (each coordinate between 0
- and 1) into a colorscale with the colors transformed into rgb colors
- """
- for color in colorscale:
- color[1] = convert_to_RGB_255(color[1])
-
- for color in colorscale:
- color[1] = label_rgb(color[1])
- return colorscale
-
-
-def named_colorscales():
- """
- Returns lowercased names of built-in continuous colorscales.
- """
- from _plotly_utils.basevalidators import ColorscaleValidator
-
- return [c for c in ColorscaleValidator("", "").named_colorscales]
-
-
-def get_colorscale(name):
- """
- Returns the colorscale for a given name. See `named_colorscales` for the
- built-in colorscales.
- """
- from _plotly_utils.basevalidators import ColorscaleValidator
-
- if not isinstance(name, str):
- raise exceptions.PlotlyError("Name argument have to be a string.")
-
- name = name.lower()
- if name[-2:] == "_r":
- should_reverse = True
- name = name[:-2]
- else:
- should_reverse = False
-
- if name in ColorscaleValidator("", "").named_colorscales:
- colorscale = ColorscaleValidator("", "").named_colorscales[name]
- else:
- raise exceptions.PlotlyError(f"Colorscale {name} is not a built-in scale.")
-
- if should_reverse:
- colorscale = colorscale[::-1]
- return make_colorscale(colorscale)
-
-
-def sample_colorscale(colorscale, samplepoints, low=0.0, high=1.0, colortype="rgb"):
- """
- Samples a colorscale at specific points.
- Interpolates between colors in a colorscale to find the specific colors
- corresponding to the specified sample values. The colorscale can be specified
- as a list of `[scale, color]` pairs, as a list of colors, or as a named
- plotly colorscale. The samplepoints can be specefied as an iterable of specific
- points in the range [0.0, 1.0], or as an integer number of points which will
- be spaced equally between the low value (default 0.0) and the high value
- (default 1.0). The output is a list of colors, formatted according to the
- specified colortype.
- """
- from bisect import bisect_left
-
- try:
- validate_colorscale(colorscale)
- except exceptions.PlotlyError:
- if isinstance(colorscale, str):
- colorscale = get_colorscale(colorscale)
- else:
- colorscale = make_colorscale(colorscale)
-
- scale = colorscale_to_scale(colorscale)
- validate_scale_values(scale)
- colors = colorscale_to_colors(colorscale)
- colors = validate_colors(colors, colortype="tuple")
-
- if isinstance(samplepoints, int):
- samplepoints = [
- low + idx / (samplepoints - 1) * (high - low) for idx in range(samplepoints)
- ]
- elif isinstance(samplepoints, float):
- samplepoints = [samplepoints]
-
- sampled_colors = []
- for point in samplepoints:
- high = bisect_left(scale, point)
- low = high - 1
- interpolant = (point - scale[low]) / (scale[high] - scale[low])
- sampled_color = find_intermediate_color(colors[low], colors[high], interpolant)
- sampled_colors.append(sampled_color)
- return validate_colors(sampled_colors, colortype=colortype)
From 1c57320ab3e8f79cb240dc5c08ec21be8e6082f8 Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Fri, 25 Jun 2021 19:06:46 -0400
Subject: [PATCH 6/9] Update segment labels and fix the bug when we do not have
any labels for event segments
---
.../phone_locations/barnett/daily_features.R | 3 +-
src/features/utils/utils.R | 3 ++
src/features/utils/utils.py | 26 +++++++------
.../heatmap_feature_correlation_matrix.py | 2 +-
..._yield_per_participant_per_time_segment.py | 2 +-
...atmap_sensor_row_count_per_time_segment.py | 38 ++++++++++---------
...map_sensors_per_minute_per_time_segment.py | 4 +-
.../histogram_phone_data_yield.py | 2 +-
8 files changed, 44 insertions(+), 36 deletions(-)
diff --git a/src/features/phone_locations/barnett/daily_features.R b/src/features/phone_locations/barnett/daily_features.R
index 4a9009c5..86e87718 100644
--- a/src/features/phone_locations/barnett/daily_features.R
+++ b/src/features/phone_locations/barnett/daily_features.R
@@ -26,9 +26,8 @@ barnett_daily_features <- function(snakemake){
location <- location %>%
filter(accuracy < accuracy_limit) %>%
mutate(is_daily = str_detect(assigned_segments, paste0(".*#", datetime_start_regex, ",", datetime_end_regex, ".*")))
-
- if(nrow(location) == 0 || all(location$is_daily == FALSE) || (max(location$timestamp) - min(location$timestamp) < 86400000)){
+ if(nrow(segment_labels) == 0 || nrow(location) == 0 || all(location$is_daily == FALSE) || (max(location$timestamp) - min(location$timestamp) < 86400000)){
warning("Barnett's location features cannot be computed for data or time segments that do not span one or more entire days (00:00:00 to 23:59:59). Values below point to the problem:",
"\nLocation data rows within accuracy: ", nrow(location %>% filter(accuracy < accuracy_limit)),
"\nLocation data rows within a daily time segment: ", nrow(filter(location, is_daily)),
diff --git a/src/features/utils/utils.R b/src/features/utils/utils.R
index 02931503..9a3dabce 100644
--- a/src/features/utils/utils.R
+++ b/src/features/utils/utils.R
@@ -60,6 +60,9 @@ fetch_provider_features <- function(provider, provider_key, sensor_key, sensor_d
source(provider[["SRC_SCRIPT"]])
features_function <- match.fun(paste0(tolower(provider_key), "_features"))
time_segments <- time_segments_labels %>% pull(label)
+ if(length(time_segments) == 0){
+ time_segments <- c("")
+ }
for (time_segment in time_segments){
print(paste(rapids_log_tag,"Processing", sensor_key, provider_key, time_segment))
diff --git a/src/features/utils/utils.py b/src/features/utils/utils.py
index 7183a82d..9288d1b4 100644
--- a/src/features/utils/utils.py
+++ b/src/features/utils/utils.py
@@ -99,19 +99,21 @@ def fetch_provider_features(provider, provider_key, sensor_key, sensor_data_file
if provider["COMPUTE"] == True:
- feature_module = import_path(provider["SRC_SCRIPT"])
- feature_function = getattr(feature_module, provider_key.lower() + "_features")
-
- for time_segment in time_segments_labels["label"]:
- print("{} Processing {} {} {}".format(rapids_log_tag, sensor_key, provider_key, time_segment))
- features = feature_function(sensor_data_files, time_segment, provider, filter_data_by_segment=filter_data_by_segment, chunk_episodes=chunk_episodes)
- if not "local_segment" in features.columns:
- raise ValueError("The dataframe returned by the " + sensor_key + " provider '" + provider_key + "' is missing the 'local_segment' column added by the 'filter_data_by_segment()' function. Check the provider script is using such function and is not removing 'local_segment' by accident (" + provider["SRC_SCRIPT"] + ")\n The 'local_segment' column is used to index a provider's features (each row corresponds to a different time segment instance (e.g. 2020-01-01, 2020-01-02, 2020-01-03, etc.)")
- features.columns = ["{}{}".format("" if col.startswith("local_segment") else (sensor_key + "_"+ provider_key + "_"), col) for col in features.columns]
- sensor_features = pd.concat([sensor_features, features], axis=0, sort=False)
+ feature_module = import_path(provider["SRC_SCRIPT"])
+ feature_function = getattr(feature_module, provider_key.lower() + "_features")
+
+ if time_segments_labels["label"].empty:
+ time_segments_labels["label"] = [""]
+ for time_segment in time_segments_labels["label"]:
+ print("{} Processing {} {} {}".format(rapids_log_tag, sensor_key, provider_key, time_segment))
+ features = feature_function(sensor_data_files, time_segment, provider, filter_data_by_segment=filter_data_by_segment, chunk_episodes=chunk_episodes)
+ if not "local_segment" in features.columns:
+ raise ValueError("The dataframe returned by the " + sensor_key + " provider '" + provider_key + "' is missing the 'local_segment' column added by the 'filter_data_by_segment()' function. Check the provider script is using such function and is not removing 'local_segment' by accident (" + provider["SRC_SCRIPT"] + ")\n The 'local_segment' column is used to index a provider's features (each row corresponds to a different time segment instance (e.g. 2020-01-01, 2020-01-02, 2020-01-03, etc.)")
+ features.columns = ["{}{}".format("" if col.startswith("local_segment") else (sensor_key + "_"+ provider_key + "_"), col) for col in features.columns]
+ sensor_features = pd.concat([sensor_features, features], axis=0, sort=False)
else:
- for feature in provider["FEATURES"]:
- sensor_features[feature] = None
+ for feature in provider["FEATURES"]:
+ sensor_features[feature] = None
segment_colums = pd.DataFrame()
sensor_features['local_segment'] = sensor_features['local_segment'].str.replace(r'_RR\d+SS', '')
split_segemnt_columns = sensor_features["local_segment"].str.split(pat="(.*)#(.*),(.*)", expand=True)
diff --git a/src/visualization/heatmap_feature_correlation_matrix.py b/src/visualization/heatmap_feature_correlation_matrix.py
index 1a42d447..a3474409 100644
--- a/src/visualization/heatmap_feature_correlation_matrix.py
+++ b/src/visualization/heatmap_feature_correlation_matrix.py
@@ -26,7 +26,7 @@ features = pd.read_csv(snakemake.input["all_sensor_features"])
if time_segments_type == "FREQUENCY":
- features["local_segment_label"] = features["local_segment_label"].str.replace(r"[0-9]{4}", "")
+ features["local_segment_label"] = features["local_segment_label"].str[:-4]
if time_segments_type == "EVENT":
features["local_segment_label"] = "event"
diff --git a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
index 1e6f79fc..423b0051 100644
--- a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
+++ b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
@@ -58,7 +58,7 @@ time_segments = pd.read_csv(snakemake.input["time_segments_file"])["label"].uniq
phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], parse_dates=["local_segment_start_datetime", "local_segment_end_datetime"]).sort_values(by=["pid", "local_segment_start_datetime"])
if time_segments_type == "FREQUENCY":
- phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.replace(r"[0-9]{4}", "")
+ phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str[:-4]
html_file = open(snakemake.output[0], "w", encoding="utf-8")
if phone_data_yield.empty:
diff --git a/src/visualization/heatmap_sensor_row_count_per_time_segment.py b/src/visualization/heatmap_sensor_row_count_per_time_segment.py
index df89a5dc..cd34793f 100644
--- a/src/visualization/heatmap_sensor_row_count_per_time_segment.py
+++ b/src/visualization/heatmap_sensor_row_count_per_time_segment.py
@@ -11,6 +11,25 @@ mod = util.module_from_spec(spec)
spec.loader.exec_module(mod)
filter_data_by_segment = getattr(mod, "filter_data_by_segment")
+def getRowCount(sensor_paths, sensor_names, time_segments_labels):
+ sensors_row_count = pd.DataFrame()
+ for sensor_path, sensor_name in zip(sensor_paths, sensor_names):
+ sensor_data = pd.read_csv(sensor_path, usecols=["assigned_segments"])
+
+ sensor_row_count = pd.DataFrame()
+ if not sensor_data.empty:
+ for time_segment in time_segments_labels:
+ sensor_data_per_segment = filter_data_by_segment(sensor_data, time_segment)
+
+ if not sensor_data_per_segment.empty:
+ sensor_row_count = pd.concat([sensor_row_count, sensor_data_per_segment.groupby(["local_segment"])[["local_segment"]].count().rename(columns={"local_segment": sensor_name})], axis=0, sort=False)
+ sensors_row_count = pd.concat([sensors_row_count, sensor_row_count], axis=1, sort=False)
+
+ sensors_row_count.index.name = "local_segment"
+ sensors_row_count.index = sensors_row_count.index.str.replace(r"_RR\d+SS#", "#")
+
+ return sensors_row_count
+
def getRowCountHeatmap(data_for_plot, pid, time_segment, html_file):
fig = px.timeline(data_for_plot,
@@ -18,7 +37,7 @@ def getRowCountHeatmap(data_for_plot, pid, time_segment, html_file):
x_end="local_segment_end_datetime",
y="sensor",
color="scaled_value",
- color_continuous_scale="Peach", #"Viridis",
+ color_continuous_scale="Peach",
opacity=0.7,
hover_data={"local_segment_start_datetime":False, "local_segment_end_datetime":False, "local_segment":True, "value":True, "scaled_value":True})
@@ -48,22 +67,7 @@ phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], index_col=["
if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
-# extract row count
-sensors_row_count = pd.DataFrame()
-for sensor_path, sensor_name in zip(snakemake.input["all_sensors"], sensor_names):
- sensor_data = pd.read_csv(sensor_path, usecols=["assigned_segments"])
-
- sensor_row_count = pd.DataFrame()
- if not sensor_data.empty:
- for time_segment in time_segments_labels:
- sensor_data_per_segment = filter_data_by_segment(sensor_data, time_segment)
-
- if not sensor_data_per_segment.empty:
- sensor_row_count = pd.concat([sensor_row_count, sensor_data_per_segment.groupby(["local_segment"])[["local_segment"]].count().rename(columns={"local_segment": sensor_name})], axis=0, sort=False)
- sensors_row_count = pd.concat([sensors_row_count, sensor_row_count], axis=1, sort=False)
-
-sensors_row_count.index.name = "local_segment"
-sensors_row_count.index = sensors_row_count.index.str.replace(r"_RR\d+SS", "")
+sensors_row_count = getRowCount(snakemake.input["all_sensors"], sensor_names, time_segments_labels)
data_for_plot = phone_data_yield.rename(columns={"phone_data_yield_rapids_ratiovalidyieldedminutes": "ratiovalidyieldedminutes","phone_data_yield_rapids_ratiovalidyieldedhours": "ratiovalidyieldedhours"}).merge(sensors_row_count, how="left", left_index=True, right_index=True).reset_index()
diff --git a/src/visualization/heatmap_sensors_per_minute_per_time_segment.py b/src/visualization/heatmap_sensors_per_minute_per_time_segment.py
index 7c69a756..bc975da2 100644
--- a/src/visualization/heatmap_sensors_per_minute_per_time_segment.py
+++ b/src/visualization/heatmap_sensors_per_minute_per_time_segment.py
@@ -79,10 +79,10 @@ label = participant_file["PHONE"]["LABEL"]
phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], parse_dates=["local_date_time"])
if time_segments_type == "FREQUENCY":
phone_data_yield["assigned_segments"] = phone_data_yield["assigned_segments"].str.replace(r"[0-9]{4}#", "#")
- time_segments_labels["label"] = time_segments_labels["label"].str.replace(r"[0-9]{4}", "")
+ time_segments_labels["label"] = time_segments_labels["label"].str[:-4]
if time_segments_type == "PERIODIC":
phone_data_yield["assigned_segments"] = phone_data_yield["assigned_segments"].str.replace(r"_RR\d+SS#", "#")
- time_segments_labels["label"] = time_segments_labels["label"].str.replace(r"_RR\d+SS", "")
+ time_segments_labels["label"] = time_segments_labels["label"].str.replace(r"_RR\d+SS$", "")
html_file = open(snakemake.output[0], "a", encoding="utf-8")
if phone_data_yield.empty:
diff --git a/src/visualization/histogram_phone_data_yield.py b/src/visualization/histogram_phone_data_yield.py
index 998d4fd9..34e24795 100644
--- a/src/visualization/histogram_phone_data_yield.py
+++ b/src/visualization/histogram_phone_data_yield.py
@@ -6,7 +6,7 @@ time_segments_type = snakemake.params["time_segments_type"]
phone_data_yield = pd.read_csv(snakemake.input[0])
if time_segments_type == "FREQUENCY":
- phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str.replace(r"[0-9]{4}", "")
+ phone_data_yield["local_segment_label"] = phone_data_yield["local_segment_label"].str[:-4]
if time_segments_type == "EVENT":
phone_data_yield["local_segment_label"] = "event"
From bb3c6141358ce40bfddbb723fb09d0b434a4792a Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Mon, 28 Jun 2021 16:30:44 -0400
Subject: [PATCH 7/9] Update analysis workflow example
---
docs/workflow-examples/analysis.md | 2 +-
example_profile/Snakefile | 23 ++++++++++++++++++++++-
example_profile/example_config.yaml | 6 +++++-
3 files changed, 28 insertions(+), 3 deletions(-)
diff --git a/docs/workflow-examples/analysis.md b/docs/workflow-examples/analysis.md
index e11ae3e7..c5288e97 100644
--- a/docs/workflow-examples/analysis.md
+++ b/docs/workflow-examples/analysis.md
@@ -69,7 +69,7 @@ Note you will see a lot of warning messages, you can ignore them since they happ
??? info "6. Feature cleaning."
In this stage we perform four steps to clean our sensor feature file. First, we discard days with a data yield hour ratio less than or equal to 0.75, i.e. we include days with at least 18 hours of data. Second, we drop columns (features) with more than 30% of missing rows. Third, we drop columns with zero variance. Fourth, we drop rows (days) with more than 30% of missing columns (features). In this cleaning stage several parameters are created and exposed in `example_profile/example_config.yaml`.
- After this step, we kept 158 features over 11 days for the individual model of p01, 101 features over 12 days for the individual model of p02 and 106 features over 20 days for the population model. Note that the difference in the number of features between p01 and p02 is mostly due to iOS restrictions that stops researchers from collecting the same number of sensors than in Android phones.
+ After this step, we kept 161 features over 11 days for the individual model of p01, 101 features over 12 days for the individual model of p02 and 109 features over 20 days for the population model. Note that the difference in the number of features between p01 and p02 is mostly due to iOS restrictions that stops researchers from collecting the same number of sensors than in Android phones.
Feature cleaning for the individual models is done in the `clean_sensor_features_for_individual_participants` rule and for the population model in the `clean_sensor_features_for_all_participants` rule in `rules/models.smk`.
diff --git a/example_profile/Snakefile b/example_profile/Snakefile
index f969fdcb..65cea721 100644
--- a/example_profile/Snakefile
+++ b/example_profile/Snakefile
@@ -204,15 +204,28 @@ for provider in config["PHONE_LOCATIONS"]["PROVIDERS"].keys():
else:
raise ValueError("Error: Add PHONE_LOCATIONS (and as many PHONE_SENSORS as you have) to [PHONE_DATA_YIELD][SENSORS] in config.yaml. This is necessary to compute phone_yielded_timestamps (time when the smartphone was sensing data) which is used to resample fused location data (ALL_RESAMPLED and RESAMPLED_FUSED)")
+ if provider == "BARNETT":
+ files_to_compute.extend(expand("data/interim/{pid}/phone_locations_barnett_daily.csv", pid=config["PIDS"]))
+ if provider == "DORYAB":
+ files_to_compute.extend(expand("data/interim/{pid}/phone_locations_processed_with_datetime_with_doryab_columns.csv", pid=config["PIDS"]))
+
files_to_compute.extend(expand("data/raw/{pid}/phone_locations_raw.csv", pid=config["PIDS"]))
files_to_compute.extend(expand("data/interim/{pid}/phone_locations_processed.csv", pid=config["PIDS"]))
files_to_compute.extend(expand("data/interim/{pid}/phone_locations_processed_with_datetime.csv", pid=config["PIDS"]))
- files_to_compute.extend(expand("data/interim/{pid}/phone_locations_processed_with_datetime_with_home.csv", pid=config["PIDS"]))
files_to_compute.extend(expand("data/interim/{pid}/phone_locations_features/phone_locations_{language}_{provider_key}.csv", pid=config["PIDS"], language=get_script_language(config["PHONE_LOCATIONS"]["PROVIDERS"][provider]["SRC_SCRIPT"]), provider_key=provider.lower()))
files_to_compute.extend(expand("data/processed/features/{pid}/phone_locations.csv", pid=config["PIDS"]))
files_to_compute.extend(expand("data/processed/features/{pid}/all_sensor_features.csv", pid=config["PIDS"]))
files_to_compute.append("data/processed/features/all_participants/all_sensor_features.csv")
+for provider in config["FITBIT_CALORIES_INTRADAY"]["PROVIDERS"].keys():
+ if config["FITBIT_CALORIES_INTRADAY"]["PROVIDERS"][provider]["COMPUTE"]:
+ files_to_compute.extend(expand("data/raw/{pid}/fitbit_calories_intraday_raw.csv", pid=config["PIDS"]))
+ files_to_compute.extend(expand("data/raw/{pid}/fitbit_calories_intraday_with_datetime.csv", pid=config["PIDS"]))
+ files_to_compute.extend(expand("data/interim/{pid}/fitbit_calories_intraday_features/fitbit_calories_intraday_{language}_{provider_key}.csv", pid=config["PIDS"], language=get_script_language(config["FITBIT_CALORIES_INTRADAY"]["PROVIDERS"][provider]["SRC_SCRIPT"]), provider_key=provider.lower()))
+ files_to_compute.extend(expand("data/processed/features/{pid}/fitbit_calories_intraday.csv", pid=config["PIDS"]))
+ files_to_compute.extend(expand("data/processed/features/{pid}/all_sensor_features.csv", pid=config["PIDS"]))
+ files_to_compute.append("data/processed/features/all_participants/all_sensor_features.csv")
+
for provider in config["FITBIT_DATA_YIELD"]["PROVIDERS"].keys():
if config["FITBIT_DATA_YIELD"]["PROVIDERS"][provider]["COMPUTE"]:
files_to_compute.extend(expand("data/raw/{pid}/fitbit_heartrate_intraday_raw.csv", pid=config["PIDS"]))
@@ -271,6 +284,12 @@ for provider in config["FITBIT_STEPS_SUMMARY"]["PROVIDERS"].keys():
for provider in config["FITBIT_STEPS_INTRADAY"]["PROVIDERS"].keys():
if config["FITBIT_STEPS_INTRADAY"]["PROVIDERS"][provider]["COMPUTE"]:
+
+ if config["FITBIT_STEPS_INTRADAY"]["EXCLUDE_SLEEP"]["TIME_BASED"]["EXCLUDE"] or config["FITBIT_STEPS_INTRADAY"]["EXCLUDE_SLEEP"]["FITBIT_BASED"]["EXCLUDE"]:
+ if config["FITBIT_STEPS_INTRADAY"]["EXCLUDE_SLEEP"]["FITBIT_BASED"]["EXCLUDE"]:
+ files_to_compute.extend(expand("data/raw/{pid}/fitbit_sleep_summary_raw.csv", pid=config["PIDS"]))
+ files_to_compute.extend(expand("data/interim/{pid}/fitbit_steps_intraday_with_datetime_exclude_sleep.csv", pid=config["PIDS"]))
+
files_to_compute.extend(expand("data/raw/{pid}/fitbit_steps_intraday_raw.csv", pid=config["PIDS"]))
files_to_compute.extend(expand("data/raw/{pid}/fitbit_steps_intraday_with_datetime.csv", pid=config["PIDS"]))
files_to_compute.extend(expand("data/interim/{pid}/fitbit_steps_intraday_features/fitbit_steps_intraday_{language}_{provider_key}.csv", pid=config["PIDS"], language=get_script_language(config["FITBIT_STEPS_INTRADAY"]["PROVIDERS"][provider]["SRC_SCRIPT"]), provider_key=provider.lower()))
@@ -357,6 +376,8 @@ if config["HEATMAP_SENSOR_ROW_COUNT_PER_TIME_SEGMENT"]["PLOT"]:
files_to_compute.append("reports/data_exploration/heatmap_sensor_row_count_per_time_segment.html")
if config["HEATMAP_PHONE_DATA_YIELD_PER_PARTICIPANT_PER_TIME_SEGMENT"]["PLOT"]:
+ if not config["PHONE_DATA_YIELD"]["PROVIDERS"]["RAPIDS"]["COMPUTE"]:
+ raise ValueError("Error: [PHONE_DATA_YIELD][PROVIDERS][RAPIDS][COMPUTE] must be True in config.yaml to get heatmaps of overall data yield.")
files_to_compute.append("reports/data_exploration/heatmap_phone_data_yield_per_participant_per_time_segment.html")
if config["HEATMAP_FEATURE_CORRELATION_MATRIX"]["PLOT"]:
diff --git a/example_profile/example_config.yaml b/example_profile/example_config.yaml
index be3deb29..c2f269c7 100644
--- a/example_profile/example_config.yaml
+++ b/example_profile/example_config.yaml
@@ -198,7 +198,11 @@ PHONE_DATA_YIELD:
# See https://www.rapids.science/latest/features/phone-keyboard/
PHONE_KEYBOARD:
CONTAINER: keyboard
- PROVIDERS: # None implemented yet but this sensor can be used in PHONE_DATA_YIELD
+ PROVIDERS:
+ RAPIDS:
+ COMPUTE: False
+ FEATURES: ["sessioncount","averageinterkeydelay","averagesessionlength","changeintextlengthlessthanminusone","changeintextlengthequaltominusone","changeintextlengthequaltoone","changeintextlengthmorethanone","maxtextlength","lastmessagelength","totalkeyboardtouches"]
+ SRC_SCRIPT: src/features/phone_keyboard/rapids/main.py
# See https://www.rapids.science/latest/features/phone-light/
PHONE_LIGHT:
From 97ef8a836839a25e45575534824ea14d3bf6974a Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Mon, 28 Jun 2021 17:31:35 -0400
Subject: [PATCH 8/9] Set color range and avoid SettingWithCopyWarning
---
src/visualization/heatmap_feature_correlation_matrix.py | 3 ++-
...tmap_phone_data_yield_per_participant_per_time_segment.py | 5 +++--
.../heatmap_sensor_row_count_per_time_segment.py | 5 +++--
src/visualization/histogram_phone_data_yield.py | 2 ++
4 files changed, 10 insertions(+), 5 deletions(-)
diff --git a/src/visualization/heatmap_feature_correlation_matrix.py b/src/visualization/heatmap_feature_correlation_matrix.py
index a3474409..2a75a513 100644
--- a/src/visualization/heatmap_feature_correlation_matrix.py
+++ b/src/visualization/heatmap_feature_correlation_matrix.py
@@ -10,7 +10,8 @@ def getCorrMatrixHeatmap(corr_matrix, time_segment, html_file):
fig = go.Figure(data=go.Heatmap(z=corr_matrix.values.tolist(),
x=feature_names,
y=feature_names,
- colorscale="Viridis"))
+ colorscale="Viridis",
+ zmin=-1, zmax=1))
fig.update_layout(title="Correlation matrix of features of " + time_segment + " segments.")
diff --git a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
index 423b0051..d3466e24 100644
--- a/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
+++ b/src/visualization/heatmap_phone_data_yield_per_participant_per_time_segment.py
@@ -33,6 +33,7 @@ def getPhoneDataYieldHeatmap(phone_data_yield, time, time_segment, html_file):
y="y_axis_label",
color=column_name,
color_continuous_scale="Viridis",
+ range_color=[0, 1],
opacity=0.7,
hover_data={'local_segment_start_datetime':False, 'local_segment_end_datetime':False, 'local_segment':True})
@@ -68,7 +69,7 @@ else:
if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
- phone_data_yield[["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]] = phone_data_yield[["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]].round(3).clip(upper=1)
+ phone_data_yield.loc[:, ["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]] = phone_data_yield.loc[:, ["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]].round(3).clip(upper=1)
phone_data_yield["y_axis_label"] = phone_data_yield["pid"].apply(lambda pid: pid + "." + str(pid2label[pid]))
if time_segments_type == "EVENT":
@@ -76,7 +77,7 @@ else:
else: # FREQUENCY or PERIODIC
for time_segment in time_segments:
- phone_data_yield_per_segment = phone_data_yield[phone_data_yield["local_segment_label"] == time_segment]
+ phone_data_yield_per_segment = phone_data_yield[phone_data_yield["local_segment_label"] == time_segment].copy()
if not phone_data_yield_per_segment.empty:
diff --git a/src/visualization/heatmap_sensor_row_count_per_time_segment.py b/src/visualization/heatmap_sensor_row_count_per_time_segment.py
index cd34793f..ecea23cf 100644
--- a/src/visualization/heatmap_sensor_row_count_per_time_segment.py
+++ b/src/visualization/heatmap_sensor_row_count_per_time_segment.py
@@ -38,6 +38,7 @@ def getRowCountHeatmap(data_for_plot, pid, time_segment, html_file):
y="sensor",
color="scaled_value",
color_continuous_scale="Peach",
+ range_color=[0, 1],
opacity=0.7,
hover_data={"local_segment_start_datetime":False, "local_segment_end_datetime":False, "local_segment":True, "value":True, "scaled_value":True})
@@ -66,6 +67,7 @@ phone_data_yield = pd.read_csv(snakemake.input["phone_data_yield"], index_col=["
# make sure the phone_data_yield file contains "phone_data_yield_rapids_ratiovalidyieldedminutes" and "phone_data_yield_rapids_ratiovalidyieldedhours" columns
if ("phone_data_yield_rapids_ratiovalidyieldedminutes" not in phone_data_yield.columns) or ("phone_data_yield_rapids_ratiovalidyieldedhours" not in phone_data_yield.columns):
raise ValueError("Please make sure [PHONE_DATA_YIELD][RAPIDS][COMPUTE] is True AND [PHONE_DATA_YIELD][RAPIDS][FEATURES] contains [ratiovalidyieldedminutes, ratiovalidyieldedhours].")
+phone_data_yield.loc[:, ["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]] = phone_data_yield.loc[:, ["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]].round(3).clip(upper=1)
sensors_row_count = getRowCount(snakemake.input["all_sensors"], sensor_names, time_segments_labels)
data_for_plot = phone_data_yield.rename(columns={"phone_data_yield_rapids_ratiovalidyieldedminutes": "ratiovalidyieldedminutes","phone_data_yield_rapids_ratiovalidyieldedhours": "ratiovalidyieldedhours"}).merge(sensors_row_count, how="left", left_index=True, right_index=True).reset_index()
@@ -88,8 +90,7 @@ for time_segment in set(data_for_plot["local_segment_label"]):
# except for phone data yield sensor, scale each sensor (row) to the range of [0, 1]
scaled_data_for_plot_per_segment = data_for_plot_per_segment.copy()
scaled_data_for_plot_per_segment[sensor_names[:-2]] = scaled_data_for_plot_per_segment.fillna(np.nan)[sensor_names[:-2]].apply(lambda x: (x - np.nanmin(x)) / (np.nanmax(x) - np.nanmin(x)) if np.nanmax(x) != np.nanmin(x) else (x / np.nanmin(x)), axis=0)
- data_for_plot_processed = pd.concat([data_for_plot_per_segment.stack(dropna=False).to_frame("value"), scaled_data_for_plot_per_segment.stack(dropna=False).to_frame("scaled_value")], axis=1).reset_index().rename(columns={"level_3": "sensor"})
- data_for_plot_processed[["value", "scaled_value"]] = data_for_plot_processed[["value", "scaled_value"]].round(3).clip(upper=1)
+ data_for_plot_processed = pd.concat([data_for_plot_per_segment.stack(dropna=False).to_frame("value"), scaled_data_for_plot_per_segment.stack(dropna=False).round(3).to_frame("scaled_value")], axis=1).reset_index().rename(columns={"level_3": "sensor"})
getRowCountHeatmap(data_for_plot_processed, pid, time_segment, html_file)
html_file.close()
diff --git a/src/visualization/histogram_phone_data_yield.py b/src/visualization/histogram_phone_data_yield.py
index 34e24795..3f765b7b 100644
--- a/src/visualization/histogram_phone_data_yield.py
+++ b/src/visualization/histogram_phone_data_yield.py
@@ -18,6 +18,8 @@ html_file = open(snakemake.output[0], "a", encoding="utf-8")
if phone_data_yield.empty:
html_file.write("There is no sensor data for the sensors in [PHONE_DATA_YIELD][SENSORS].")
else:
+ phone_data_yield.loc[:, ["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]] = phone_data_yield.loc[:, ["phone_data_yield_rapids_ratiovalidyieldedminutes", "phone_data_yield_rapids_ratiovalidyieldedhours"]].round(3).clip(upper=1)
+
# plot ratio valid yielded minutes histogram
fig_ratiovalidyieldedminutes = px.histogram(phone_data_yield, x="phone_data_yield_rapids_ratiovalidyieldedminutes", color="local_segment_label")
fig_ratiovalidyieldedminutes.update_layout(title="Histogram of valid yielded minutes ratio per time segment.")
From 2e12a061c7972a7b2ed352587c3729d9c84e1890 Mon Sep 17 00:00:00 2001
From: Meng Li <34143965+Meng6@users.noreply.github.com>
Date: Tue, 29 Jun 2021 00:28:48 -0400
Subject: [PATCH 9/9] Update docs of visualization module
---
...data-yield-participants-absolute-time.html | 20 ++--
...-data-yield-participants-absolute-time.png | Bin 184589 -> 129699 bytes
...data-yield-participants-relative-time.html | 20 ++--
...-data-yield-participants-relative-time.png | Bin 194497 -> 136067 bytes
docs/img/hm-feature-correlations.html | 102 ++----------------
docs/img/hm-feature-correlations.png | Bin 295093 -> 619587 bytes
docs/img/hm-sensor-rows.html | 22 ++--
docs/img/hm-sensor-rows.png | Bin 188334 -> 102313 bytes
.../data-quality-visualizations.md | 4 +-
9 files changed, 39 insertions(+), 129 deletions(-)
diff --git a/docs/img/hm-data-yield-participants-absolute-time.html b/docs/img/hm-data-yield-participants-absolute-time.html
index e902c07c..527b5e69 100644
--- a/docs/img/hm-data-yield-participants-absolute-time.html
+++ b/docs/img/hm-data-yield-participants-absolute-time.html
@@ -1,11 +1,11 @@
yUKeHZ>X}1?R&qu^_yc(Z$t-e5DZ3LZ zagV#K@#_+!f%$VZ;G}~KmUb3prlH1Et%x#Gus!ap45T!JSCP?p=JD9!PZz)VFm7-2@}id+b7_nskwy4e{eMskns z{1%*^P;8tf9&Z6Xx|6)QJMW{H3khA&CgKTJ@nQzMe%YgjN$$BE36m=H?TV~H#zO7H8B4nYodk(YKk&4%gVlV>CV@ 6SQ;B*(Sa{^c&$n1-=X>D<)m*{ptm@c_WZ~X K$yJ*KrZE|nRp?gcCA29Lsh5neTpNLXvvEm9lzP=XLG67pt_p_{ z4tm37Sr}zn-cQ;em*-sTRm!}C5=xYGeakaEK!X*w9_&M~_{gMbf`n0CE7Q*sqsGZ5 zOPBI`Px~f0&6hfWs7s~qIzJ7B=hL8#Bz`0HdW&Dl+TMQVD*-}in}3S*XSz(@P>7IA zPoiEU>xnPQtd7-5E=6pI%7V1g=qptvbRY8Crw(>MCPHQ@NQbfH@uC(a |po)9ofq_LDJ&VlFlG!nvPU}w7C*&vqHZpm~8T)49) z=Qd(d(d5PG&*#Xslf>m10vvLBjz`F97Wypj+G}t%4A>@ZnwGo}q-^R%Rh2i(!LpS@ z4J5mE&K&1F>-hnsn ;9#>;g3%tc60fR0IWM&PYhRMnDcI@zHO zD=bhp4=cqujlYR~BHVz5AH%lu)eNj#RIFAm!5(Axt%^) EeB=}%)KgDR{wIXb;~^!{VIXM?lIpW&c~46Z*|sn{v?24#(ti8KnCgQ!VfE3 zmcKBZd0ml*3Msy{-13Aaz~1T$fvd$iWK;#foa(a)E9DSKiWBee1|SBQweld4<&^uc zJ^-hhYF%uTg?$Yv?3%mk&kLXAeP((vSk-czacuXO;}79oKJb$b$VY4Cj=GbQWJ+Eu z$Q$sn+_G;=xa3?4>=xE9wkkG{7qR#$f(40pR@zr*01Ld2`y#tdUHAqSUB%P$j&|$2 z@3Kuwwyg9ZrPSZ9hWV9So$I%#pS>5j6z$GwHf|llJCvb+Yf-bP+djlpOKr5=I%>i! z6wR#q3&bdS6sqg>xS)N`Y4N~&;fpL(Y9&Cf>-;JLOMTl *b>3mtwTER#iiN6r(R~ z5cv5ZS^oJ#E+z1Lsnlw&?SOyUrNLrdgGY@xEOF9rn&BHSVL!5X->zb?@aZHdYnk30 zC0knrTXXg_CGcw>%41}QP8?7pJ1%ihTbKbt^YZ&$Husmt-DQi7O#$6@xOKL+TIXN2 z*5b#~PN)FX4W8JpwpcR|l%BZ12jtsMiM h?1Kmjrpj*oM*DD*f6`GKuDvBdM;s*~D z?3#Aok}Vc@`m>0ZZ!2=7jK#Fom&6 7C3kSnOXD1#lnse?BSD;@|usCIhRATWF)^a2? zJ85)xSQR%G+H0$6KtmIh?Y(OVQWPD)Yl&>UbBs&@k+l2->d) Xa-R~U7rX7c&0MW+5X$C1>3hrh7Jz>$uSoDLdYYqZd*uJWQ8 zofGK#9OG9eo(hC>@Lhk*vd!6xo#|bm^*0{qhSg0Hxj|<*8tIktoG5X=#}DtSfat1r ziQ>HdBvf?hoTH*=J)d185wthhrT+l$^L_}uce&o3LfB*dad}L)!km;8eiQZmjN#pf z6@HSpYA|;o0}+P0)nZSe3!kSWtC#~0pE?*xJN`6;@nH@!UEo?1z$>jJlWjvo-k6WY zDxtB1&wX2-KB(G2^zcBlZk^H!{pY;^*5^;4U{Av!uKnNa8%-`!oo5^qSup_{8)AZ% zvJ3KIa4FD@KQ|50^scLB{N5|7Z>0AgyI8C)JcaPXetG0~$Es^ucc!`T`Tl?x5t;nh zQ>6)Wg#}KGcKaMLsk@#>Z1XuC9=`QKc7C;(Yu>((&GVG8t@@sP6$_%IfMT|3$)84O zd2WTE!`F7cay_RlL_ODYgZq?G4reodRbnF2=C-r(sTN?ypu(1!I3W|I!{7b;{rr?H z=gZRb(&`Q?rS>7ao{v$TPa$i|Zkgnl=8J;^{k`2)d5N%V?b-+;gMlEs a|}W*{X0i44@ChSX)bI>_l~kY(lb6;#SuSfajc` zJCJlJ)Kt2|n~|u0_ ZQC49d7pkQ)d2XL|gDLysvTeY*S<#V&v5*0J->Q-s^t9Zvg%O~zl z1zW!pAevQoKJ$&SFA9A)ogv#bLoR8VTw#9%ua|;`>yL9AhO_tv|Em!HLvgOeaQcAc zDSn(QVLT~kl#%?>i`H5ob!&mabeuaMq_6u1>PW5{(_y25Y zEm|9Y9oLl*EwBR~L4EhroW$R -@v#%5u4*; zZRY%i#$Ka*lEir?Vg4R-kAtX`?cRlDB|9o9W*7 zY@NGkwoiExjqWI?X1zUhcKE+>5_pj%XWt=5v3*EQJXaR-$8>RN&}t=vtcaa?;9u%K zXA6Ht;=S&r124{R^|)?9_Q!piQ3#6HM?1^I5m1x^V %BUl|TzSi9<5o7}GPyZ$uP( zdUTy+mN?x+8lYs1qpo)!UtbRd9X69fj`)*Qq}|&d$Vn$KZUppXx*RnZP0Xie;zp<5 zLwErowVT3(J=Z3c5=#IFtVNh zq7Y_DWobFsz^Fl#GSYLZaU6kr^J2P&`*+*^Ahm)EawKfMzNq6HuoHZRzOJM9>G?TY zW1MEMd$FB4i|6Ev^8Cgyl9<0}x;;65GJj?89mnMzAG1m9Cg-2 dt-sw|>5M zX_HSsJx~gL>C|7Hp8|wD0!8r%E%)a!`62RqQsT5seX#fsmpZM3eXkgvv!Qohgg2hl zWNdc;h#Ev*CT^b?eYV| c rQ z0hllXWv`MnvKgbhiWce5aJk5_ WOV!dyI#69 uQxZ7Py@B$2I$Jc#Z+_u=h2 zI4beTt>-XUN~yS|uuZ>`;5dw@bp_yIx54q6?j@YowTu~(G22gLGD|EkWd3Dw6+O;y zT$Sz @`5=a7%fw&Y_eb>+ifBlownncE7 DxbdZ2XF}W}C;Ej!X_O>O^~ b10w%PQufqv{CXo-fnP{mnG2gbdHk_P@O8j?$L|%~w}8 zPi$vqbD%WWoy<(X=-*D=le1xaFOZR_g0YCvT6rdVQ`S<5r@0UPAeeFtxbTH(G$2 z4jvoRmuEzu |D`qg&li8Y zD5yk7_#p2r@t3pyUn>Fajs`T_0d6 YfUO7|@IXL3R-{x3Avc4P~I~4FU2m8h-Wj zPsgd)v(Q3hb^9J6GdrQ=ZGwn-`Qpx8;XpzO%i86Cz1_c^-+yk+h|=z4bm(n(5Ko;X zn>kL5s?hF(a7{gO{tqtr^-;?g_&LmpJGN8NedUmzS4Q=F&YS7Gmuda`7gj6WMBf$h zV7=EX7I&ImN<*LtgQ0(wt^MOUb$;3hGDTnWV$c0>G4@v@e&m|*wW{~OjlyXWu?wS! zj-{StwH)4U9g6(jai%p_Zg!0U{~Fc!$D3}RSK~pVMsCZT^8CGHzRTFNru~ATt{Gb3 zI=Lx^d>LY%Y#xHd-sy`9 ZbO}tU)`LvR#73EXOIwyPkKzt zPvq0_k}X2JFzXsxW4y`aGDf98UuSi{8y9991doYoI#Bmsn-o3C%JWg8P*3Ci)&2Qq z1FgB?WblZVr%~%u>wJxMGnqM>AxpcJ#L +$v>Re%Pll)zGBOU(5y@Ae)Y@K1OJbZ_nD6dJ_^PJ8r$Gl16*TmP$?ZvE(2Ma z5^%!_V=bD!BX9)udD9J0c3cVYoi7!2Ox3s~JY}qEsrd)zjeMW8-l9~*IoJ1NfY>k* zcD{gz2AsyDdm8Ut)MB(CMUqCekUWYSAR2dD_itPA6RjUzePU%k_Xxn&@zMIG%|BYQ zvmlE@1!j+EG@ub2^-;JDK8u~a9Jj&vHzWP{>AP=0N@;I`3XWCMlJ=i77lXQp^onu8 z0%h4M^t7sI>pe@s#PC!|Mzh=m=vO6iFHXdqv(VwV5GLW)b BGS(+ZYeB@qHrQ+Ec3k{5 zRh^NHI2XQxH{~Ua9Q8CZVAy#~MCJ6H*fIDusk<^Rovvx-zkIqKqo7}k~!u%^ #bzJje+&<-?(+ctQ$E`=a#c7iPn#!vp%Vd3;XKF%;zv| z`mK9A=S8O5Hd)?NAc1eUpF#Uu?bqImfG+OfUH*oUmy|5v5jPqWnBT05crE_U8HjZF z>R5hxfwuZ;ew>T;GKOnztC|G(Io6s$=d69c9+eXpxIKXKdC`JFs2zfMt})z|#T=_L! EIjNx`3pz8aBfedgcPk+40=rtk+H zmYM6(2Gco~7dk4X`ohsAD*WA@ju$S%k~7j?Ws53%R}J0o!=d${18mCJmj$!!&t6-! zrU|EROR=i%0$j 8#|6+j`EH*;DV0E=lMAar22A&B8Y1qYzW7%R~q*zW3yiG0e7 zsiGm0f79I2ll9?kF7k~h{%J1$QG3)^$*k)QB)aFKBi{3TD;RiM1G8Eo0Z0Go3_Z*a z*PjomJQ;_nlCLHex56wVcPqkH-s$5usg%Y7>+5D@+CX#APa%|lXxDaVm 1MbYJnse~W3?k^YBz~h;B6V%E*5}~)A$hyj=M(n}!V1jue{^Q{3o|4X^`%qX zJnu)y)%N|vh*^9zx(hINQ+kCfVz8%eVkFk=SptZ?@OzJVi<6qAuSR|mZ1NSR-I<0# zu71K3na}*LbfG-vxwSB0=}5E}#u<1;5#%k(t~`jjHk9W1q6@G&O3-EjCn7$6nOaR& z)NVG`aN*YkJlXrg@nZ}YiXXjG4+w0X7EW#j*lrD%rBr8ro2&d0FzU?>pp@X6ykU)q zn{8y7+xaeDOA_84wCTFULR@V;Gd(exhVb!S@R&`TDhY8y+{9NWue#q5(jB-zPkBr! zzz3iwoTsnK<~X$6sGiPO^7D8Q$vK*I(`8VNcyD%`Cqk#D7GWc(8$VTI{dpdLQ{4;9 z|Hs#Ag+s=p3zJ6OPfvzLT~V|R=@Hb&Ob7IP#_GbV{fEo7&L?#vFcyI7;8v>7t;^J} z{ktlXQZ-ua32bclULwGub#3_e-=QM5FI&;bIoGkir>d8rp|mdDFOL}cX>jWR`68QR zK3by-*b3a8ljdiIRd ) z7MBMG;KTW{ugX1B%308Sk#y44^68o(rvK|6R*dT1_1`>JM Kx@X3X!GAz1QfzrSCHuZm3PcU-z_+_yuoR2Rgt9SRuLXi56< z8gTbq!q=FjMsgyTDvJRnFoKBF ^Mu^0(|1(?x(2rg<^bmW zmQ*95EbFc8;PWZA?uWWo;`ou*6E(1Y30cqSn121mKcVHX4dHQ@&8nHDHKw?NES1Rl zaRzo67m+~vN|k!i@ad~x{?>*d*n5u(afmT#sFa)o3-3L6Q9Iizsg&xxXJmZEvw*NB z?RCuBScT~eej215fq7WjOn=@?a;f3o`AO8%zbe#H+qic=I(9{7dlcAjNJ;0YwMEV) z8BY~;OgmCER1MI>fCS>ZFGEB^=_4)o(q64;R9xfh{6EU;GM;q5o3G60!jGv~yd{!( zhMj;v|6{ZDT)Df-HT|L%j&bP1X?^c)@eX13oq#>+GS}8(p8L$i95NE`0@_nvVyDf@ zl?c6A?!9&0z3aUEVwedCKzg61FQ@062`IA9dskuaYy12?(lSxZ0+0jFWW+egy^p4z zdF}sWcs2sjs5 G`^p(r*9mPBl}5}`vnceMzAa%fdDSxJ>(kIENogn45b&?-J5X^ zt`{`9q4+%F^UXSfNLF=ZBj(VL+(_ITsE7o2;Irf3&xGAOj_%fy+vgu(lAPK&O^4)A z(veprw2HGOUZ)NHGBLhDyyRd!4BIpf~z5T5Pq04{ByU|&l3Q2ovA>)m$>XJh>ybd5sz9e_k5 zm+jY_vuX6PMnB5Xi@pq45F1NJN E-=sBfp{|2}a?H$hb|#A32t(L$Ev-@QJsta!(L13+-L(|Aw#zZs6Tkm(FjN zZSS@W9_+p6v)*?u2S)FRs7qJYB^LT`f41_~Dpy@;_Jn#M=6cxh*Ui(-?TU7zBi=Um z9vEmO &0E zO;;P9cyRUGRcdknjbre@&Gt5WV$gCE1#Mr`6mH#da$V J0H50=n*=Zk#83XESSz z#n$;4r~mqmXXlTbNgnf`_5YByc2@G8C_JB-;>4mFq@AKhVsOH%838xW6CT(CjuZK{ z{;Xq}uSmai#~8tadfb$p1LU#h_=?iB_-lq?6uFZz6+tS|9$$|2W?iWA)^R3df?gsr z-6f}PR`|-wH~QT`*EsUL-$w0~Trcx{WNmQ?2V%C>d7 zTIA4@V#44}O&ooY9=_k!QYAFv%YUsHR9`e1-O;rq9k+Y3sYCbZ9~!H!`_nDb5M09p z60OSp2zVeRCZ;|8ME;rtS=JlAlQiUt-E6`B!-;{IH-OVxy>6O~6s30m8$;UKTH|f} zzC98!s JHXyFYRTS#R7z2?5wL1@P}=Mpc( F(4D^# zlhJq>?t0-l;zja>r>TH!`%7IYwb`M@HUelH90kC D z#%Xq#sY)F_@+S1na6sKFn&m}TF=p-FB6yFCnepCsKWI06*#N2V*DB%nJizV$#Ypq< zIO<=beonA%hDT4?spU^Hpv%wwZyvUIKQK`|G8DeB5$z}dHnC=x$imw~9so<0%cddQ z )xfWTe8z;5?jKgY8%K&tIsv^seZn75-muz5Fj%?OMF5KNsL7OW{orj@ zx|8Dmn6_n4_E^QVdX2P5vO `h6%LXdR$7!_$6WZt9r2@ CnH~kk> z$N$46;^4uuKSlPTQ!ngdj8%_Y^JNb7WgV1Jb5`#PEIHBDh5vZeI`B!_g@u-U@)h8u zn=kFoFjhH)L#p<}oK?w@^w)qL=JLm;^&h3J;hP8Pc-LblCpEhBP8TRW7wco$t#{vL z1c$zA`>gv^s;bB3O@(?1do$Y!8sq;rvijR}LS$TefF;u8ge6N!Z8M`TLS4Q+@@Qw^ z;E2U5pj?i|5+ 4t{CmFM1%!Wj9TmNhN}H1tmIYDcdk?L$}L;v za$AH&E?_7Hl$@th)xq|N*?r+^06O8bu~xWpr)gEFuewEj)86a#k%x^PG+19uLeAgi zjas-qyMA_6KJ!G`N+{# 1s+RK7J+BiaUvp{Xj8D!6^IWc7YW0nVi@5#F zb;xA(SYonLYGeXT<%%8u5Qfb=vDrQnTMq>=%7{;_Yr19DS`#T>ZiklAwQn@Mm}!U! z?EF{o $(40u8lkWaR|?y|Y~x6O!lcTzOGFd%`3pr~tC-%&x^&tu*}is>?HDpK5c|EK z7%uX^qZzT|oAv>qp3<1X>J6V-!wBYJ7X&Ap nQPZM&m)UD>KNFXOJ@WSY+W9?XIU}Z3%3K6@(x|?xjzL zJUugssa=;h8QJJ^!MLSkkw&A_5_Q;5;>($7$;^L-fOb(<{tN*n;ibk!eoV^+b3)4e ztnrnnXH(BrbK)4!UKOWlk%$2fXWq0jQ4q$^iorqlwvxjypyn&Opkd?2uDZzLDnb2h z@<{0|0rGpRGwo-@&IIs13o;Td2f}BXN6jxhZ?pI7d~W =?lNb-p3dTJTGD`OujAWRXS#VcK-<4?bG6tfUcIj)^@>4NkY~~6VYlnF aYtRJnZ|K+WIozU|(h!F><^8QiCIR7C7s*UkOfsV@EJuYc0B z%J#cuH=g!}oq7$I4FmJCm^_ z0KOsP-r;kk)}{CQcf0kyet({!wd3-Or;jf0-q>;f{(<}V4j;K}Wx93iku)n)Q`2*2 z@Bfmbrccdi-rCj`%UH1`&d-nLV8_U`SuQbJOA!UBhmT 5@T9lmeGGY$I~ zPi0R0-j8EG_$1E{uVJR0Ry#>hn0V~vvUtE0R6s1t_8mhcAtb|SaOK4V!6db>HyZO+ z`Xx(&Yksimd|TdgDd0dGP5zE%lICgG+A)N9f%1rsa8&t#*L+=Q1Es#sj>gW+YkMj? z$@Phi_yeIFH7#WKuw4bIHsb`QuM^;7uXX&mrQeocOm^>NkchY3@8y>{wq>fGPlsCM za3fw_I{#4H*1@2`n~+#PMPVlj#@QFwR)pH(^}UdFQ`AVvsHk%_WwK!LN&N1Wulj!W zkyRee>e9u$9~mHrlS>hLVHEVbR|9|A&@DYkVd3nsQD~RkTl+Bn;?uFtF_VhRhL@;O zZDW!`EnLyvl@vd ?wW56HHr@!Dk0X6 z(>O1}u8AnDhYIAbv|;#ior$mFHw=;^C`pEy<|A+(o1*L>RIyD4%9~4-4pZM(-Z5Ri z@|?xUoYi}ejfB9%Ad0K9bexwCM&<^5%=NT48@rHeH@<_sL7AAyj06RDzGQW(C4+co zwSp`a-l|C@x6GPb6yJ-q--3Mzh4*P@7_ok=Ow^ o*EsvzU#2?-C9^~u|@5ark=_$PC`)5^m%g^)h W51BhPjTB z^rWA7cX!2a**b1@EMmCrAmusep5+*#T=@;$C;5=E9w*RNDrKlD#uYN?(FD@aZ)vU? zxRB!OI>x==;+8{srF`o4?|QR~*rp+tISfpjm5p2>NebiP^ZIUc$MefRFdR)jdyL@t z&Icz4;rK3;iKi&hDc0;=Z22W<73^4a8hBfnub^6a80BC#c QKFye+m_?NCYMiB8*_fji;KUzRw+q$HRr?hUxDVB{L{oYZaK`XD}RCDWy)uJ81c z7O*la#v)=N{1#$U9UeD={&X&1c@8u)be>=&xIFyBPb7TQWze-#6Z=lN1`$x36Yt_H z+8S8P3@lnpqimn@H0<-55vi}g&~3O5shY5nDnI70&(^hr2>De6h}<0o8MALlGWOX7 zv0K8S-r%~x8tzN!Gp-D&6wb(vBj|+-;J8`tdAa*c2*{w5k>%Oe>q{Dm5GaQ ?9SqAo#8zKA33qZZU+9zQbMBNLjJ^u=Y~9@OG+f4A=twjiB4+m;&e$)@p^;W zboyr@p6-G`95M5bRK6$!hb8N{B`f0X$t`iVxrGemSS3;s6SVgA53zy!22jQr?VGjH z9lQMn495}&UBR^Y4@jda;llaEVJ$nb6dLWpGyAmUwb2THY%VpCPSR}lu2cO4&V;1% zn$Dwa87#j8nUH*Rmw-Fx$;UcT`mgDQC4ZP}{8V`jNo>cZqmgk=+`xslQea_qG6o-f zbpHt)^}*$X6R3H0O~#7Z3`;~#;OVd?Hw lHK>kf0p+4*>I{~P1=pKh8f*pLCtqA; zhUeNz?HQ}4>$(y!E_Vk(e(6pTT^7^M9}k{&pK)-!uA9ZRMPetbk&^wEpl5K4 OK~u2HqrR!knz&^hBL$~Xhs zEylFt4h1?wqVj#fxy2{%p0eNs<3t$)XZ$CZtDWLXwq;ziXhzH#4eV12Dw|$5DLHsH z5hecF?<70Mz(xFeky{)XT)=f8d{*e-8bY5h$bgG|zI@v^?lWfMO6f=l$ovdXJHQNJ zn!I6OnDh!;F|L56U!%~%6oblK^?G(U%PrU73ORLToxAPvbgmC+yOp?QY&%$Y0EeH$ z(FY88f1VH_`pU8aOYz!%uAQ{v>Zi)3u%+|Ow)yA~aCRY&i&zW#W{MMzHHXMPM|Z8> zjA2Z#7?6kzG-=he6q?Vy&Lvh2!Go$ej5i@6K@p>B=6~%^g?~HQZfejUxw!2DA{HmP zmt`T=a^|qBuiUjI85!3UPd^4_ts{@$92~Eem>G<)yE&VV0sW$c@s>V zkNtY{ZamieA;~2mPW~p`WJo5Gg9B@SRbGw7spHDd_ege(XoKgD%~&-K7yHy&X^Ur) zif8vrIjB|lG`N?zj3TmgSjAwPw%=qq-NyFWaA9s}G%Zsv*0qR_bdB_Y1b%bnk>cd) zZ3rXR!$%{eBf>>(SGPOx%3PXnpu~ST*&Hu4)e(vk)J|#O6GR-(wNU~&LABkmcO_gJ zP6WB@isr``;ydLC(I?w5v%*CqD=3)!J|)HqvxFxyWw*hU#)E1qRbLGV-$myIIAB4{ z(JQ#i1aEm_b#MI%A@J$r#|(QAJ_TOye7wy_LFJ!R)g!rABPjfDpI7elb*fXz4>4m( zH3;_t8k(x4u~IiC2Wq0a!kFsl{R3JI8}VU#5-BWut44+=YVE1ExR4`d-^n7Z;L`=^ zwL$c(QHk4`L6h~v 4JZoO+ALFxI7SH zfj2#6!bo1_WsjYzbR@1}uMw}ObfRALuUNHJ72dD}`=PD|_=`G?R0Z<8Xe>yhRfZ+q zuZOQ*9Axd~nvQa?mj6OD<-z9i+uDepJidIsJP0k%ZH1ocmYg2Z9Zw#K%okjf6=k=E z-ctGtO0aw8u<1SGCGQ L7Jn=@ z9 }hsoY5Ot& *>b`3s6I;O>0zY9HLXzwlL)Jm^EiY zz?A-0Qp34xGLW-yq>E8EIqjlfb q%;0jT(0Wk>aMcrfMWxOCCfq0;kL!% znuU91XRssPRK+5kji2oiQXqyCj=DACJ6Y9)hQUk3yCN^V8&v&j!z{Sv*qmGKEBgCD zL`}zaXo+_9zCdjaVr)$t&y6K3dfFzASnRXEJ#svW&@-S#YrW7PI`+Qm)58<$_yjrF z)+`dW+F__9wOaWrc6@H7D)tzIv8rE%hGOjG({UoO(=nBHGr6eZR+Gsv_^$5;XJ9bh z6()t&$6g6Ugc5xzF(%mwWSiaSCz@c{qHQ+dPau<($Pyh|O+!CEdJHDH|H+Eay7urQ zV-DCvIGO%qS3A*N^(N-}q b<<2qyHl;$Lyxo~yU!k };osdL0}1Pw4TpjXXg8g^-lN1Sv9K@j?qoa}eWgN~GfTRAA!GG~ zq@vB($|O?I4X;f@2nK^h!nNf3e7)ommT)GSkD4=RhwbKw$j4Zj*L9s7I}?uO9+jG( zN&Jd#7>Rv{JrkA>Yc{rE2PpbB`=Ue(KOXnKo?dirM_R%Uh139fP!Y^TTy}UagU0ne z)}!pWR(Nnx39`H%{M`7$#lfwq&Y_`e2;-GfvQAIvO%lQ{e+#1qUJrZ+3$BJO!M`ql z^8*Vl7-fLt96rF90l>ad0+v>7 msN|2mfemyQ2JGu6 zOgxilkE@>DXI&{V<(FCh_?iHSw==$VtUG^c*r7{u)$hX~%t2i+kdUkNe5l<74AtRs zWS{s^qMpw>)Bxwm5`8?%=<83E;d^FEzwRhHyR%_sL1e=4!p2u*P_7rQOOIE%eS|^_ z!cFp20t&i&LRU-_q>3K)Nq_*z4f+=5xRG;+B?{UYFE#lc`PQeP)7AnR!ox41N^Fni zca-ExyY?W(UMEBqg%_mvNWKMU4XFKYtt;MA7!fCQj1w0I?GN=@eNcZ+CPQVE+qie4 zR|Qd2GPL5YfnnuYRr8X+vaUO(S?uy8w?YRA`3Gs9nbEDB7IzC82^{3vT+{Bh*^zsB z7EbQ9#GkYkM2HE;XcVc*BR}S6-C&yg8CzQhKq0JmC^3zVNy>)o<4$!cBs@R-fw4>c z%^0P32;JN#H;g6ytE4yyH40y_sL9!ZFxLWSpXP|68tyy^w%vy_QH)T}Mcn`geYrew z4fA#z6PXE*rmrUy; UA@*6L3k;6fqB`{GLxRmgH)wX}|d+ z!2-*UYr}p&&Fl{`TED8`c}BgHhQVE!@%! 6sTc&yoTRAKLR`mUBe@Cm~HIOSeHm>vBZn zGUp3%+1K &L8ro J?o0jE_rU zzVD%UAxSbU@T1#7=B@^UhLU;hh_LwH6PlRcYnOv+Vz%V(F2V)WvYy~5h6^#r<~Qh* zO)^7{`ZFRWa&uRKI947%h+M`#2-=BJ48GkHD#)wevh)PM(P>sD16^wt`dSK>|61L} zDny;jtT989{A*hxEI#!eY|GrML-7t5yDeR&d@)(n&N;FuXnC+i*T2T6KYY7>YMN(8 zcj8R!=>4WDTiKiP5;xc#%ZFrE7B7I_3`i mmXej1i3|u=eMqckZ2g{QUs+k^NPUqe?JRO3UOm!S z&3`>!`$0li?AnqS^qoySnO0nXwSR2ak6Ov937^PO>bo+0{Cv8jcj>K3ZAcy|Q`-M{ zBX<+T2xY87%@zz9kc_b)zo@e87plCInAx}=-dV5%=qxDX *25@c)jJ_0&_$j*+|6|ozAM{= zqt}owo<{C6v2*eaLGJ-^-iU2PC%zSj6qLfX>U{v@CBMToo;MaNxqp5=Ec9!kWsU3h zLbxuYUN8V#P67>U=~9*%45V{xTv@koCucb;cDB_8Nt~m3l@W2Iy08_>@$Dh-M%Z!l zY?&E9?iCj0)|tPm5pFsCvSYFJT);iD5p(_uuEm;wRp`depkR(C4$tlTHY-x9qR8)B z98*Tjf{$>ruqf+@k-$Te)2k1zRQZ+S4{a B~8vWeO32DUifbA_Cf z)#+I^w_iwFmlEPy%Y`SmD2~uXrGHU`ek)Z>)N_PXc^EZXM};>wk%`sBhp51Ay$DM1 z@=ZM3#7AIlp&w1JbL9X5^<=G$nS7b}y#Ze=wZSisd+P_IKPG0`fL(l@gVF=~Ydl8G z1?+TxZynx&=+1htHYTH2H-E+U7RpV{+M4;Dj|n2`&fmrmsvcY-0)P`#gOpfOzVYhp z*75VwuZDERvkMpVL)N~QKCgjYu*j1-o#~t9FqnWN`%Y?enKLYx!Sl0whLHHTL{~@S z?1kdJEOzis^ZGe-tdtcpd?E#E*qfj@;(J1Eiv0&-XWdYTyX0fid=xBZwWjH&jI%5Z zN(xsxos--{qzl{5zH^wY-mV`AwsQ{SWKG!1-?nRh?E<5B;T&(yCF12G>;l>XbSqVf zcJ#V8x(`Fuhp#(ATXLIji*@Y{ i4*rNcO17#>RaR=_6Bf2T`)q(vZ>jLDkvV(1u`OeUDf|ytTj# zo1Zn#>^EHq3h19xkYmU|yKkx}?Upc{FcP!~^X@))kC=_v$ilGYZTs!cM9Y5R^)kaR zvah@*TuJJ~POi{7y0Pf~id*pLj ne-|vwuD|sa&9c#l#q(~To2>F1wX6Aai~c}c@o}^1UJ!i530G{lG)MYYlhJW^ zUW~s!`yMSA8NeUm#6RsUtvjCKdk`_g#xn z*)Q(6{`~~nF`RILHJY1Nue4xo?!Hh)?LA9KuA6I`Y?`k)qa1T1B~r6tqh1hvrCD=M zUyo@)FpmN7!KQ-i%8BvQ$O>DqEhT?&^;We7cVH<1!+kk%&akye(_lPN*n&2K#q$=p z)m`$^8Wz7!X9P_5=#Bbd0@6$5pM?5UQ!4!n0J=yQE6K#FN(xPYx!J16;xXF;MA6`q z3W1WVD286veou#3a|xW=!0hIl_!l(g1df;{a}NOtx)rQlorP*2h(y;(Vxo!l#IhCv zzh^oVQvH^YJ3l>q)J7Sa9Nkx%3oFCEDD`M&-)o)h+x@xTm8qa7 X4MZbmivl=b`=d*90&WX#TgTr y#D7gIFhj6`M)BCAYjxOywD8 zz{Wi$95Df7)^fTVv_LpN8L`gJ)hs1iw{?H2mZ|McI5#LcTvJ9J-s} k&m94Tk~S8C|%QP`m;5hmDa!qQ)I*tcvux9`I?^dr0QZHRuZ@>5Ye(n``m zHOW8oV?smvy~LrVY&4%sH!;dWc;t9x^gcQC!OPTy^EuS53;b4MrCZxE Q7!d|jK2y({`2~#$zr>|G1r8CEdF;j&66X* z;EHMOr2qWS|09n1ANTQU?Ut=B$*!gn%74WR{zsYf|J~pk4;#Za6ky##F0t-5-N`SS zPrPq*(z2ATWju)h>mwt!jjqtc?CpNvy0HG63%I*5a|=7xP;p4=`aKBlkzGotuS%ma z 9>@RV z%LNcL&i^Y5;BW8y|GP{7|NU|w?0#{9F}s4hn2|p?!644s>);mGz@3{Frd`|mr$puW zQ|_%{MtOEFenGJ-itP9{BrIzZmQerTxihr~Wl=EcYlH~Ql&_&h;FdU*+=flye>94V zRg4^KuZD?4j fR3Wx?7-91HVm zb2eU^(VCdC*6=1VzKWIZ&>bJr-{)KRAjBzTFuZ16-==6V5Atqz+h7z=`1|PU;d*;} zIMe_-t}ei7Zp_=k=W2=9s+oqhKTz^bQ^av6e1XmwojRyfkiP#blgwDuYWmc+@ke0L zq|SO{wq@e>=bZK00Hdpp7N@Two7?_=2RYd*YP2z| Ri)*b2)d3`Rn%aph>2DAT^8AU^F7}
CsQu(p`+;?3r}Hx=U7 z%z-{oAZl%?e{MKD(V}@My?BbI1pQtYyjm!eAfmoMJ+wNl$;B|>R~>S&wrk{*fEAV@ zBidsWz7VaahF`tyh8K`^)INTR#qsN8T*W~~%s+w&+mqvu$Yg^2F(iYDB>B<>ngMo$ zglg_y8=N}ZEU$5rlEJHYy*nj7F-mb#9FO9TogCcvI#HbHa!O&Dj*f@(sQlpO5_S|t zV7w9b&LeNOSrZm(>OF|J?ulfbvx`?jLQiD`&nQ(9QeP)X8x3!J?rw?EQ>bSUo4(9S zOCX=}&upxYf)f1rA&_@>K0hiOZ1syr2V&jc!RibB$a# Na`Lyk~=5c`@)G33s2`%&>j=il)7 bVMEzJ-*(2T{(V|ux3%qe~U3MbP@%afZW9~;~F+$jRRA87@*^biCCyd zN0Jz8i*VPmfradv_v2d6*rp`bO+mEW<>eeLxJ)U(P4R5XlZvl;&hLCt$Dt0$niOcJ z|I)g4hFj&>c1$dFt*nixBe~fbfyjY)X$Ppw-M|`T1d&0B1l?g%b6MK;-nHM+0@2sk z3?9Qva6m~oK}>a)xx=gq(Animfw>C`PLvRSytb5MjfSjzBacAKzkg|I<0xf?j>8Dq z&B^dGFZ-ytx>&MxH|{gzWo9+E!H)gBoO_oxCczF@QK*KzqnI&lT~7eZI1Y-hYXn)x zXK`HauCgPCp1Dv)OV!=!O7eXW+}&cAxT0pyU7COh@ovG{+bt!KH~sm+3+*Ba zD#+;L8{ZzfCU}@s1iuSfAd!`-@2f+*ZROEE6JOpdDIsoUS@>!>o@AGdRZEBwKs6$< zaxT`O8CRlNF2UZaz{=B0Ofv7??bDPTut3v;ie!=3aLaeEMD57l?S5A5R|@x*P3>^g z^L;bz+Q#PB;y#1ajg0*9nN5|}ma^B@NI}MjUq(Z%t(LnUSXSM0M+E~gv;*?ji2`el zVn$zt12TH?lviSx`k9eiB~ofv-~XnDUAyl_wjeGwE-TlJN=X c`A(-2 z!(w0sZczhfrd1Jbx!-*8hq`yrN~^ExeZQC4yrPMd=5q+}8sqND8|+7jkF}g9SK_Pe zZniedLv1T2o$vMvy?4RKQ1`<5m|!woc?X3I5dls2!=nomBSjXNQ)OT7hsGyqVoRmJ zjfJ&gC$id*J5Sgzh9wHwcP*rOvQ2r^%%#rGX4k;A #OhJX2r90N3f`ar&C?-W|avRTt(5?rD~Dv^q|V~??`4HGOp3aSA8aG-v@E=ZO#5# z1ic~ljxqUUqa}Yl)p{=V7n90^o>F&$u `D;O aEEPtER|MUvH-xCd;~!QzI|S1<=l5c@{$WT}YX*705a7JIK5 z`v$qz45?m84246(Yl*O7bnxXkN_@WZp6RquTS>Qly9*62Y( z1U=ABtub9ZCw3UI1C3jnj67 7p#$#Krv{8wpXaE?Z3rZdq@8cux8q zncIDe26m^T$L_yQ-te@CKb!J%jyOdCwsV_WF2BuaK+gRWt%896z7=6mg<<)q(4Sv9 zyP%?+%=&Ou3MP65RKpVW>3_};ePBIP*zH$hP(Hh);kWf>1ufOY)U7uKju~Xsbm=}x z$ fro0p>HxIEg7k&>X?DwCL-;e7%8d8*Dkji0~0F(Nel#Cad F_}@rN>!%3`J~caKxK1r*gG?&FN-+t2u Q=i7Qjz0vN2i-0-2;5XHr zY_jvULrY#M=?}5)y_;$hTa0q~l+8DB*&~KO(Ru%)zuEuv-u-Nf6F)f*v)>J@C1+qM zbw(4roRDZL%DCYe^EkfPfZBD2%y0&U5hHw2eK&-jtM=dMms|#Q4CQ-tjaLO=lt=xB zFT_7Rp^o|8zVZGJAtBK;hpIRhnv!fLCvIE|1lV%8Cju9HJb$lqq8w>hXF*J2bDO%z z#g|WMiz||b!vf`Owc;b6(J*Pq)`Z}z7T`DxoSn1yRZJ6NYDg{CT33ya>T)RZqb$^h zSt``BQ=1>#f(!0l{nP>I0#LNwK_}ZB;Sxb p~&b0d)k8ksc-cCm{u^Ra#NAKVEN zHN~+p0q Q^(2 zEkPR_gRBV#3a%w^1rW|8p7#@g;-62X=XX;8WAqLKYV`#9APg~~-SOikG(X8kNH^E` z!ZmXmc9VS@kzd#$Y0L=&2uBZ2r2pys0F|^muNIq;PBX|N!jM)|>3SR8B8Kf+Ir0Qq z<(qebPI4EpR2x_5b?${J>W$ZLtpHHYn$~RKI-6(m(9Hqu_08D3ZEd0;xd~wM_8KrP z$yG5#uK-rvU{kIgWgo6=yhF)XhXp>@S6`$q54+o~4D}|QlAboY)%F2}ZmnK@B^JZ% zC&0q%M9gDMGR%5I0#JedT2rf|r-VxI *t0wPLdgS0^9a6?{;+cH(l`@3AseNp1m zO7??h`ru(>e}V_(+sphWO5g^Xp@sGK?)J$~Q+498$<21vqX*O%*e3tl(CDSE_y@Az z(=L^q4$J|cb@9>w$~W{9z;INuj-#FoKz}m?Cp}j(&Z84MWLe=MXz^!;l+lWHc0qgb z6K!L!^6LJBh!5L>9-Srqge3n5N@B|1=H(pap|`N=Dk3hZZIKyTq*mpYX#e8j-XBWy ztQKvU#Ni=&s}Dn)(VAS0`cAieO{j@~ZNx6YWpQ$$oJy5!K2m36bq3n>QU$4OOA6CO z_oa2K4oZ-N@&|IKxE-D96Xpl=I_d{)?K`-rz6<8FJ{OI}yP`SKsquBxcSo+CTE9~T zJ#MuQu-H3K1T=WP^JAqohvPwY0fLc)fTKF&x1``M!xD~6N$;@w!acK#yhsUvsV^#P z(MW7{*Omt}a^G@}KO70Ab~PXOa^9m+NzGTkLlDmBNYJ=SvlhH%gd+=0{7M_b?3=d{ z=xn`?2ZmimgU(NS7Dn<30g+A~7%!mAJPRKu=SX*SYf;ts{^rh;*GiU@X2Tk~N?oxP z3rL;>$k-p(d63^56FAPs)tS>z?yC*FO{5W|3v#JH!t1>%L+(E6hpxC*P~{dvdfGRY z#(|$`KTKl1>COITDXH=f`(dx+ztij=Ia2ye458`+@G+nRkKp=5^7>9`a}6wM8_upc z(lcxLRR^L!a7bB^Lom;pLh^Yd4C*>HZ47xSZKzA9a~I30QM=iWZF2M}f3<%Rqa`Ux z$JvfU1kGiDz&qh?p+eq!-@5{rY=+%rMT*+ef--=>TDnvAW;P}TJ$*mZowgTIwk@>F z$m1I>y^MIzQrQJSm?oIh{YM+Nus0s!g%v(6iUt!T(Yi>T387l5%7*s09SL^#Ah0Ox zdgZYXTn)>B9~?|VsHMqCwcE|Qs%r0>P; 5yZ1i8`@jILxWsmp!ThK8*82qFp)48jkGNICA%CHgE+;lwA8oxkpXC{!ZE9| z!l0O-2+5pq`Z)M{p=ev3C(@&M;+D42k32&qfAgYJ#RzUk9z9}`HiwtK(MYVljv-uu zI@0sUm{-WplOw@#I8G~>2P}9(Z^*2E3e&Z-1ZR6bnLqGn3A(PSFu?GNy=on&YhpuE z`KwgrpR16#?f8qim BNoqCsKROv zNiI589an#p5m`!dC}GWg;-qwOkb#8D{U@K6jF40mcB5ESNjsjFLS>|Did_-kp|sHo z9GCgNsif8GTX?8ZNfU7&s)_@Iv_UH-FLQ0m-Zz2O`xkR%Zc4coZ5CZ4i*FC5mo9`g zJUe@5@|r;aW!QD4Yaz~)Czo1SUoyxO5XcJ>=mwVOa4P~!v?vld0L}gwH8;TF-#6yB zp~T*2$lfTlAAtQ!#jlOUifhAKYrBD4QF_;(n6f3s%K22U&kmtAddtDjckGN00%+#m z)W;7L%z)|mUw={ L1($)f*hQi7(qDKAL_MtTeSQbnO8^)>|85Cm&9)=}zzcPgwjvu-2`=+*J+! zCP@$#)UJ1bpt^vI1~Zh7Bm#EPz3HaWw+-F8W&R%hJp94yi*xj+k8aEDDp!)CHo^F# zwOgd3_HOJ6-txDV)85Os-lpZTlQmL~4*a@U2fVLw7{HI81vG7r_ {gZ(lE4cSZjqf8g7@l-J6B#Hvsy4U+ zl5T22K({)#nB(`4AI~Wo;ykpv9+ m#5?O ze(~I)^Cj0SQ1yO)#o}cbal7ufEeEcF!##u`@4HZom!-$SliQNGw+E=QNHo!;N;Cdp zObxEO(2lC@6qiKMt9090u{ntz${c7+UnJf;^mlXUsKHjkxLSkt3+n9qU(|xFzf8k# z_hyJvW1;HztdCKRE6?uY4IKA`RC{DBLR(;>E5+F*-YyTQn2F7NxGyOHt5d2WBJ$I` zSz!8HRoJGKe;i|7oA;*sz7e`$rL4 Z9#H{VF|RIzBo6Z7{M>qkFr)c%&?0@KD{mg=izY z!Ad(?w5kls<48ClH7F@`Ayqp}5PGg~mxE7gsOKH_0@m`my*-*rKmA^*e>{o(wp)#B zy3=z5rP2o;rGeu^q_I8@Q=xO&cn42SKCxyQHg~5J$li>k_Tw(5M1DL<;EtT8m^kL^ zZcgqOky}+4mtpaGIE2*Nt*_XL{tiQ9$Obm%+mtSG*CAD1*4zYpAtcpv0IorXi{l3| z)Sm1M{N;4o!(ZKNsfoPhSp?C7U<{gHt*5F?ayQB_c65Z@i4cM$2D`9+33$txF;E!b zqgjtfSN1+K;WDCez)Qgu7L>mv$x!4=BEEn^I)3ecSf}C|_M6?8&RNGVc^|sAHN4uU z$5U2^SkFukM#A3l;ka1U_>l27nG-n`IGxOd{?+hmQ?9jlK5;*vRAkVKytT%UIt$$D zWi``#snnw8kIc$OOm@)DY}I7;s{ZQtPXxVnlv`iN6gFzc82WXv?2&*Pqz{aZ9lsIv zN_$=eGrH1DclkezeRouoS n6)u+C58RM){^CQ;&{EtlFPnAt(B@A^D~~ zJw?*L=D7A-1x=!Y&(7+p>$*k`*$JUNHN%7QYGh>>Ye)eX5+G90CS?bofQ@EwxOWk$ zPCdI_Cn2FTWT3o*0ysce1!vy2&hHV#1+#Syf)u?m14>W^jBR^S6BTt8zM?P@6r<}C zV3eA4t}x6k(+8JF*bmz_;!>Gl)Sr_5t%wFsKhD9?5fA_Y)+Kr;&T?TP#)>?9sCEb7 ze -bdxTP}hqgkb4&Qn>0Y%_N(lZljNo<_BNzEx48Jw8hHAA1#i?b{JG zu^k6400;Z6&J+H~;C^|6#7PZ27#;ju$(6gbrJX4fn^DK`i4Og*qG>lPc4;wRUtD zkdeQ^>zo{^0np=oj5iv=zc#OVB5FTzUQ7MEPH*MCpDw&L!r7bo;7Qe-c7WNmV`yc* z;pN0M1poM45;(gpbC{}4lALHdXYUPnQJll9w1Wc2@ar^Bfz37B01pjo0F1-HdyK=j z{Otd{0sF6C0}@K78^?)$$@pqG &}5BwiAT=E9FO=ku&%>xs-EiP)&_ <&lC=ofuz+v=#e7OPxw4j&O8eTE+u%L~Z+X5jyNQ z<})-R&47@q5-evOD*`EVrk&LqKU{1Ew>hc(!i@6^7J*+-TK(&ZQ$Kg1Q^sH(XqUuG z_HTwmp7J>Mot*atEy-|rDe4M1+e{B08cbhp2&XgZw8Fy}ooBR;`%#i}<5-=Glie%4 zEp40-v*d)>%B;IP7gy5p0KIJ_o-!IW7bejz3Cv5sf`p--d@>P7VW>6-lNb=TxrRlC zUZsf}_LSWHxByzY@ ~K4@qsmc1CGY&(@;FphZ58MglTci(v>IfWChknoUNKJd;r zXF-%*YOaaOHfe>oJXE!bLwp!gHpdKTAHPPx>lt}i^(5ox*2lf_gYfTCl7x8q1oa+l z!yCuFdBefnze&O0p8~GygqiRCA^naJf24 @?|+y1?+jhN75Cy5yH8M~GrNPJ4^Y?crTnxW+A>0MdOM zWLM$HgA#~$I%p9~UG_8>$BW-Tb_~)aSoPv&3jCAhF=a?o;cx7*k6MxJIkS1+(Re@H z4~STA@q8)=2#_KvymC`!PY1a@Bl5zXL92jijIBIt&Ez8X^u_tCvXBpVB%ek-LuF_A zt77*_YDr^QbsbXcnX!VsPB!#JrNr!5dm>r@b=0iEr-IX0Qv>)?(=upa`fqdpuTKCa zy^s!R$MCBQbWeQ85)K%qCf2I1KWie~A;rc?_fy2&tYPrrwN(Ge>46Y?{&bYZE4Yw9 z{n^9Ie_3g7Ds7M1Gkb1)hW}iOtt>Sw-7q6jAazv8Up}AU42E{oBtu(K0kQUlC>iQG zb7m;PSa#p1jNs2QO~+6Q3la^+_(N}WSkM2HgYrwhdQw*~O?srU51Da_Y-EOqjW|O1 zC+4KVx*?-5Uw&pqbcDBQG DG=i0S*1~Q&)AzLT!-qiT (thu7D`&G#X zOHKiH)(HWkb~X^9dWqZGxxCgl<4`EGZiz%lGv*ZOp}IhiU?0r{mWq<0B%yt<=Hf z?BXsnd`BI_xeo5Eu=)JRHb;#Xtr!)Kb*eczt5W&IQErIu%Co22L=8J)_`9V|7iLV& zB1@4JYnt!Zkb4ooQHnIT>8^5JA6!bHRTW^|s`Kzu;eUAykV&e*Cr)r%_XdxHCp6F$ z<;*fMiHji-4_;c2_bK|4ueR^e(X;A-fJ%8vunl&b$(DE@1xrSx)EFwPY-1@y*~E3) zRw?-^5Fk`NHfQ!9_&H=AxPX3u*2iAbaPUMIAX^>T=yDhzGXsrpw4NJ&W*IwAOY)JW zf_&)~MHano=5bmd K{33_e!K` zCM)rZriEw&^?72N)`qmEKFf8yL_+{n7#2HffT!0FY9 wL|duX B{h%d<`z=(L^jH^KM$P5ufx|r( zhFjqF6sjm6UET>#94yRd9F$3e>oNK^mfG04r+x79aC+q8Xhw zCjD8TT-Wqk^?KexMMU01^}r2Ocvym|Hl9Qwjj1p0*%|C&4FG!B8b81e7U}MHe+FQj zV6z3;d1If-bkrjm!xYhm@VNqYfk6`+5skZ|9#Z(crYjXC*0yZ8X1b+x!Ut$l0{ z^ZQRRWXYpX5hPK0RT6JWo6XD1vr7FJ_LfU`M>PdkdMx?Mr?WST-joPLvt9QZjNche zt}T0QHt0|XH$R{Ko`Z{(r+~wi|C~V%)RmiriaY{G4@lpI2r|X}V@seywGRDOM2g0v z*6D-}J0CO#{548Lg65u}WsPZ}1Jj+%+YPVvRs{ENF*DHvzd^rfx~HEw!^JBf=-k{P z_=#98Oiuen7UU-UC1V?zU6kc)*pJ#3MAuu-zbKoDXFoL(@`dVe*s0zVOPMBbn)?ZC zk9yyL1kC#llwKpXA~d}Azq0G!!5emh|AIeosflR|8WI=#amak)s dkItF_R9Biy<=4|n@ht@!C{J%Y`|H|^^C z#IQFEn+m$3=JKE2Dr!~GFz9}4#*`ouf+hK$*WIg!kB-9pJ>VXjTkHVIVF#;1154b# zH7~zOeoFR$Z4$_LK{Vpa^UKyJPAlwmlgCr&?}ll&0FYpwzGKN$Lsi~YLkrZw^5+>g zYBoK&4-%3N0nE<;V1D<^-m Xz{dh_NvDB?!7sf5O5`~_UvkC5`S|xlfB-kpp zEsfbF5tJG3M^hIVbxpi&DpEWHmXpFm*04TwGw7f@PGXR~Z%`V~F*^g)R*bI@JKw;q z7&&G93y%O5C(6h$H?+0F)j}FAgRY^&1GTabN7vnTKGD3?;8mJMm)+^DHGyIJ71``V zW!xoeH!un;*__OsvTH|r&3IVU=hNJ#pG3jm(!7O>Mt2Z$4ZJd58$!%ew~->~1v>b< zzoUkQ=8!3Hg|vKxkXRDa;&WT2t9DVU*(~uz!Rf7Hst`coe@RjJEJC^g)GuKnwsJX0 zc$k N&d9c}cWHp~fpxH%wngSx0Uw~Fnq_>h9 zCmH^X)2G092ojMzD59TGqO65HN^(g09^sUsfu>&3jXGPLR}jQ1vg31=`H?^@e3U3N zj~lLYj4gQaUh5F}I;RciDF^KBs{UsluC_V+n}<1nFZdOKl0{#77a)qhRf^je1#{7I zsw@b@2+b{*zesfAvD!N)I7`$_Td!voDo9Q{$KC+8b`o^tio!p3;LfMKYk8yU;1c^~ z6rT{`m(Y>*0-s^MVv!doPxEXUPxamSqLEnxSw|;G&I|5!7Bhd%SnY6%+~TS>KzI2c zQxp2ZllpcT%ovMPhms^dvinE_#kZZ{dkUI`SEl2Oa*v&N^~f@M){g@DnwDaRzZ0I7 zkX`$3Y#x}Z-NVUG|A82o^$_1{Xsfkk-Y3Ob@OO%INmS{Fp?$I~n8zK%tTUYJ?wRwY zTRFxjoAo(eF2va?BQ^`Y$;Cgm(swpch2NLkw~Lz&R}fZXrv@da=>SYJImepdBU}|T z@MzMxNpP4dIbT4M8gI9;aO%Ja{cz}ZD#eD;vFaSM5hi-z)Ia;X?qOp~Ce-h1pSNdh zfkS-a`=SZ&`5PKa(ZdteVdyn1bwLH>n{e|q`eEg52W&`iaNEn>(|3Qrx-HC_vvhdo zXL)Y!0y_NkQ Qq?0SvX$GfmZyTO{f5M)tQ$E?DRDs|~qJ zB(Q OBQ QV{p6dh?Ez_Q2>e zT3e3HFri_C9lqxW&uj&(8PGYa?tF5VTNTbP)pyq1+zt!xS>{U+{+cmBA}EI=VyP0X zQDh=pI!Ab#DAdt#SPkTrkXsG&vF7J>iOS!0aY0I2;ymyG(4kW`Qmw=vCO=0?z@Dx#o9W(tlu_1wT$L8h6H9L)c)iX3w7GU&!W>63;zpna36 z<3xV2%%XVZna>_)9@KmWI@v1nHo>HyoKD-a)w?^qrXC`1)B=l7%L{#`-cpHTn7qKw zi;g*IX|@`!7i4reIs{Z<-*!n nd? VQ`)PRe1fh>`qT>t2io^XiZn%|1l%jKL(W@U0H}# zeT>0WD3XT>(UOu9OrV+mukMLHXgo3;I`K57ow%8CN-i_D>AkAWeIX<*fEW@`gi@FM z!VXS%DXK1zV;Kc&c0oMh!IFixs)aBAV{jd)(}?N=2SKWY>r_E*=*m&z!%D$ciAR(R ziG?4@8X3{TNT#d$ $Mvmz@*fHwK>`s>vL0NDFZ#@@yqb(>^J^Fo z-E&@72L2+VGS^~v68+Au%dE&RG|3uh;9X3ry6Q8tj~c)A&HueAYTR-~1788m|1jWt zg!ruYr1E<=3 et g&l1RB5N)>pT68PzW+DYJ7WC6dzR!)WkbpI=>{Li01%xWGf z8?_5I!qSo+D&x@ I}MfW6z-CL|U5$xyRdN;f|KDmW* z265b8ea}Qbt1x~Onj&(&dk#5omS(dCZu!QEuOQBFW@^0d8@D BKm1^FxIh0&yF7GIqCo$X&BJd&J^GO&uJ=f{f|Su6U44nP{en|li?4Bl_aBbT z>GbcgQ$pFPggXlt(LP3aEld{r@b}#!IPkC9xkdeOJscj)P(*>e*CprL8&^h31LW5- zdoOgys&D`k3jn5C<;4@iKezoU Zm%q zdS)L?89Siia6cd~Ji(;e{1r1SNyY?b4O`9ibFmkkB<^Er8U{kqZX-8S%C^kY!!ggs zA&&*eY8S7hSw`)am6`^Wxug3QrT=PSC9IQd!wsQs$`QRDY6#r|^iAj)4rm%aiD? z_1RI`5>vKvF~QuPvR FVJ!kntt!8$L-rQKm%T*f~~O6!W3}m zB2vR2 YT@n5D`vi8JvO7zaqZ} z=&Pxkh|~|kvNTMeVapj6a*s(E7E>d)Gr$~MSaib!CDX1f9_WsioG$pyBHc~5->T@W z01t&!D8H4{^K^Rg7k(_AQv1n$4u6()VRpMO7(-f3Z$`pfu!A2fV}outM@G2eN+|q$ z*1)NhSMdv-9m`tz1vyjk0Yx_11RW!T=r5GkL`G+F3^NMJsr!MZ7sxtCPGA>T>Hw*U zXf;V0CSwU+&0lA3L9NsG^xYb@YJ6C7WQ^vRQ4R7I1)1Mt3i*jH5JO-y=rklh&-?^< z3*(M+2($-ERgW8iHs+my`2}E5m%z{#*b7Tj0}%z) qvco6q z##pH45XIBNDN0GzQM!hZPKsT2PbjjyoMN8l+D$|`i^!iUST5-(|8K)vV7^hStf46| zx+dYn|HCB4?B!*22htKgnz$G}RXA< *;ETYi3A8oK)g`jWOkm>B58+%j1sai^fr`?{FYF*)41n8u=^X6TU zHEqjyPtqR)L_c@1@P1;mK^S^n4ijjd%pAX^69@8o`60yA!z{bY#^%6Pjb4UzfwZjL z3GwqlPaMQK5C1?U2Ey$q#;j~#-5ok%Opy ^62t zK8C=0+bNL>DMT+EJ#T*?T_9K39-RDVSzQQZ%6|?59iBw(auTnKP)|loRvYQ6aw1uH zK0HwLEAau|&P5YXFmwnH1?Y4;{8=)_{+)=_^2%^w!;v4$5nIQNe%uE3p+O|#D%X9G z>;pjhkcPnbf!24zdB{s5nX_N{f8>v;87)_+lYcFc)?dncYP4s~&1MMgJ1~z_{)o}_ zFM@!Og)dzHY;M4SNBrYKQ^~&4#7j!h=Fw=Gr!ZioV!GQQ`n!pTNS3;h%N6mR9vip8 zh?LLdt6KAo3$}+HRzCz8XTN4y>TjibIlk%dIFXz~4Nw1vKKy?rtKbctM^*t tKrpMIBQSa5n5-7tJ zbkZ$tnMGU0!s`zHTMyX9L5#jL9#6R;# F@4HFrSt?lByWOrZZM>>x~M6(GgO!YET<#N8JGD(RnKvQx)I)3oJ5KKX^ZDW zc6}IU4TI&(7auaT 0{X1R-miC;-7?LFpG7!&SdZ0caLOYa&P2uf z_{;zMiiAJw)(cq=B!}OYB7H1 z;_JKWM;$!oig0B zB3%f%i539fLt&o-U(F5C0>dXrQvaz2OMkG-EpJ2C(OtMcqr(vDrS7n&|B5cTaok43 ztWYrLqfcz}a%-aAh~Q4iRLbVktpb(*xeH*RBr&1jV(S?#mh|}G9z^UGGW|P!Nx5w( zv(pWz>VG5YG4*J- U}Mm8%TX=2yWpYuu5vjLvtHe#-&Q-$ zY&AaMfX)6zNLpLgC2h$LGU!i&QqvBB1z4^z=iOm8P&%zbZuW%`6S?D^1CsVCp-7um z@RoHo&~&_;1phN@O8rZF`*4H~^{NI4DgPYOHvBC~lJ&N7RCm j) z1@zr2E%z9;>Y!k)44;mmLZgiRa>IeJ9qDQ|7;4pvN`D3x#EF1YOx}!ZYRrhner|jk z(1X+*Qo*M*T~Qa)q10*!%m@zCuFeXUNQ`7U)3=u%`I-^kF<$C7pFQ^M4y`!jLLoag zJa9BM&`Ex4YWxzIv^32S`?-xU!`=fWfh#*eJ^X&5dD*+N7dI`<_OM;E4WyWq)_Bo2 zhs2NAUgB4E#D{shmBiahx3;Km7qP$K0kBQKq{nk^w&~v?h;zDqW6yAYtBm2uuLxB3 zvs&F5$y=S14>P2o5 zT*8~@EUQH^6!H?yMz%RJU8isb5F%=qcS7;Atm1m=^AB*kuXy(#Q%x&oh@NIuZwxRy zmW`cQE)eJ(-Gfm3Q%9w<$@9=FB=>|Um%`wIj{J^_1d^e{a}dH>21~j1Eh)lBIUXIz z>_=0n%Yi2pAyi*zV&_i#xf9@C(n7k2f=mJAB>MbA=sKTea3MQiu1#?3QU5yYH4am4 zPJ39Qx|H@%SuO!`0B(TKcuMM&D1n&6@uEW#p9u~?2cgnGU8MMhz93({mWQ#y{MJ}- znKPQ&R2typHu$RD{TWv-*iO7-fDZ^kMcgm?x@V* gjLkX%4`IhN@su A zYCGn25&Q$L@gSer5m+vKDnJX(sLV8#Od^CANdsAq=LR*1p}?vu9@y&C|9eYLhw6^b zaG27YSVQ`mDv5k3JyK89x^_V9A@&E!v&{Sk{5Z@XKCV|}bwXC*KH?-R^w?ukEG~$= zz2?d-&s)HtM+bm7UTYIW&bz+ke=P z!O7+=Z)wRNo| ;e8*?@!wD2layzPL`xQ{LFx<%{pW zb_mu=#IO&@?#E^vwaL0m(q~T?T)n$eAj9sHTbN-LWMlo6ZC7#rv)FLD<+}8d=Y*=~ zTbukJ>OajX*(biHQA~0U;V|?mRt1ERaBX_1D%6A2v=Yzx%{J+Z&eZ!?IUh~ry4#t9 zq#z`sAt8l>ZOg)dgD%eZkNng!K>c=w-+voC3 _&5N>1QK0;NPB@)f}$K zhf*^<4X;{{4I<0|<(A57W@Nc)b`41*{f%7luI^=e7x|5j4Rak~+HLCUqG@jgh$%ja za)KAo>`MMK+I`w#Ckp>Pu@Hy# ~EByPd=kWF5<5Q`X%NuiB+Pqm=AWbexeV*Z# zZFB{vCDl;xfKgi*K~imBa)zSp6L)iwI5WW+Cx_+{@JLjcrSfykt(pQcDBFo}gywXG zmMjw)TwwNJhZ|Rv_8~INHN-Aco-~>Oq ZYZ1b(1n1A-&CnH~sn#m?!h9hh#3~3%nbunEkPJ6HCtB zy;Xq!moqOFCbFi@Kp!(e&B2wwZ*9g`5*5{?InDmXP%5h}=Y-UP$$e~}k)}N~+2D2n z_O#02Ef#g-!e^-$xei;LViPFxm#eRxIL#%Y8LRSyBO<%n1#6H}8It6~V2kuSq+K<0 z&ftCcjD2ecR6q0=VoQN^PJ3dwk3(QSI#{R2dK|p3z>k%<#T0FT!e+Dc=~=EK#tV01 zd`!i$W2%)ahsP`Q3?)#8oz?}`ZXUA7{KX6coJBIH3ZyM=<+&S?Yau?#3)vwaVtLg? zA31LMA4C@Z1sD;RjplOm9*nqf5z<+txAFp6a=hX-DXzPqvLRK$-dIG6-ELFoKrfFQ zLc&?CbyvvjpqE=9_wUMV9-_xe=2s0RYM;Yl(dL-Jl6TgcB5{k^i`lNhzuTg+=YsFm zng ($<1etA^nx{EWIjinY1m8{&x~_ z!7o#@oAgE%^`2K1 {m+~ftbO3?h_;B q qqrf37M_gFodcrHb ziO7n?-sS-Fha#)Mi(Q&Hr`%tvu?8*zr}=I38_(~%K2fdT74(5RUOtzoHE0$W(d2sk z W+}&-uXMK LFVBnTjJ-%bxn(K}QV64<|=3 zt%B&Q677VKImfQuEjP9VdP1A|HTd`a%0i6g3f~+7lj ON zdeFDa^U=j2;D)PEMu3qGwk*TUI#BBQk(>L0o2Nd=S>4ajMTGq7zYuOV79jr&=!5LF z1_ULaj_2> Z5RJ#s)|d~*!hn2(1pqE);Tyjpnd&a z?y>M|ohyego9KLe+O|XD%!T)79$uIYGBR=HTXMai%X9GZ<%y?vkbDbIF6dnc2^#DC z!O%#J5S(nb*eqO#4|MaC+wp9s`G? f)qWBK#qES$y8pGOoPrMMRRiQ2jK z24629S|~?kGn;ec5OcL*-&Hy1KhYc=6FK44Zm!#d&U>)S73Hjs{{YbufqMNxLQZ3D zfx`o<)*PPr1?(zhmZUecg@?RLf5p^a7YV(Ye94gQoFq9rptsp{Hd^aG5M`5Y(;vQl zxK8^Jtu3FPtsbxHk%U+qC=5on{GQh{oT07RUBs9mDhITu ~`h*w^d4+s5=ttwm2FlKddPuA+dF);-$igHP zI4$4UN4+P+CSKm(Hq>=aM(NWqDbWm7A3{n(Dmx7 @9<`@ax8s~ z2n)?MI^pxWWg1u63EZYleW3~D$LTDfJ`U`7d; -IH$Z})B`{HAC90-Nv7;73 z42d9!`1cyGLgH-V%oCUu);8cXDz-DL{%$RHc`dj?o*{>g{v+TUVd8$}AcKT#I$z2n z1!|B>#w{=Girbh!X5mxc&NtQmv30UE)2LV`td8%vp(TYQvXX-C*M z(H|j35L+9t93LMiKtfjn6E)#Eh{Dm8Z-Y3O`}tNI Nehbe)sU`wG5cyPF0v9T_(* zESn=MF~ym<9OtM rUGL=EW{)8?y=`pjy1ZIReW|)p?%S1z?MlCSmi21; zLeWqxOmYP*4}S? ?(;S1xMb+g#1S)00ZsuC*n33#3ZlA{k{EUg7*M zU;}I31L`IV!N}$?uKW3X$E8q;7t8yHht=*rZ%mqEVYGh8#}jgrb89+gf*VeFJxwc! zkJX1A-Yq5Xe-C2+Zg;MC{AhqfZH1||UpLYvaQ&qP9ktlwCS`zfJ$wD-h-L4Ru1hzS ziSAJ|t6BWE$Z`@xP9?Gy$|ngiW(jSj^*kKwT>!cn!fEHcxKVFD>9;x>V& b(#M3{t0#9URon!rn2;K$p!G3XLwCGThP=A@_>E z&ShM+I+mZ^zxbQK$}js~0FQhMv{!`>{zCcx8pD783Es|Q)Op4eVw+T}%sLK(Vqs^z z*_GzNsLG{)s3jY#EY>xfxldpW=b%7{HXB1;^_{N@>7(RweuLHNzRh(WZt@f0KSY`7 zKs~=w%weh`w&>W#=Pipe_ z$)dIgX6QwS9*R2%EZ)vh2(E-)@8L#fTfKM@;%@xyCRQS+lXkhK@fuP6kM=lh6X8Mo zwY9|@99wqK)G%s7)GLaTky I9o3C}@X0*z8Fm_3mHfu-m__ z^Y^E@bRCn=$Q8_xTuiToC5134SOlda#XH%`Q^vJrstaDwt_(GT8VJe1eb-~}Z|A%1 z&DET@act8H0uH(Re2V)l!>wkiJeTrHwkTg3^U^21Bb_AWWclHsQAtzI&HJ5CM`ZI$ z_*MmZCFj_A@x!p>;hT>(f6WnD_ndt9ltX4yh5`PdwE3_^lS{r!@dI7*cTw>=+QK)r zR);GyG`rMFx#FoepjXB;o6SRhtHZ}k%LXpPMn$h{$*RM9O|vpG_PYW`(5VF7zQVRy z3)z@0c46T|t8G)SFV4@?PWtt_1<(B2@u|)mKsC|a26&5TQ0g^!wC2X+?}3|FuQI+J zCNR4YevEN}L~3YTEjOEjD<9y7aYG7T>o;r;K{h{+T()}b9&!Zg%|ELLLOs+BN8%U| zS?of3Ci_x(RxQMssht)I-a)=v?k_f OykpV1zQ9WO-c4^JC=?BDvUT zvO2EiwnzAB=6cbhnA&?Mw~+5eSFLTj>cdg+m%Z$UmupYfhkNt~Llw U2KE3+(yV18rRCBsivCj_rMJlLlxXH-2mt>H5p*K_l+@w=|m4|H09gte*? zdU;WE`qjyR=Ej9@i&ZpP?eSl`LzBNx98O-xXf92Em?0gmk0y{q)=<8UM{Nql(EjP) zk&D3?W497k!ub6>;3Wb+B~5$2SLL$Ty6PID dYf0PiF zl&-G{2ZFrxR^OK|zK}7REsOAlHJUM?-z!f`E_?+K@oR9rxrnzgyGCPtO}l&d`Ss6} zRel%*Vd%+1hl(qlY~uyh=xdeh=ELh~e7!2SZYCI7y%5XydB0X4qKLun%w6BF^+}dJ zBm>u!Z@?J9NQ~u<&MudtCg)<6Nv~W)(M*{yNjI*cLc?|aP#|l6B2C^#^LtM9=>G1G zUgiV8NyjD`jYiqhH8J0wXWIEb$(tuu=V>AhHq)EzVU#OgrOGo*&TF#beC3n?>v`{H zW9Hzsu5s(aaMtn)%XRek!xOHda;)B^glxMW(Q;RN8}Mj@@1$7H(-3^GM-(ZDY=tO8 zR5b?p_7g+L3-*Ckr=rce)U*t{p`AUhZ{jTIRUvacc|k*1t<96=17qX0_d3n0E7 K5G-%{z;v_V9is1r3t`N}vioG`oKekl7p;_UX`2sNyHGPgLW))Rdj)z7U zUJw4gCGFgE!uH;E3@#I~Ir{?f^@juCsyX@lI;=7&Bu+~D+3sDeegJp3KF$1=#H(>q z>M@(Y{MutL=$$x_cR+yAU7v5*9p3wFlofU@*BqL0Is+X?8jDzIgjvPQ0k>VA-y$*d zhtnM%VBt*?>5z9AadD&7Wi`<@)Lreh`qTwdaYemv?;>O&n?fr1Zc|=((!VPX!H&y1 z-jDL{3}%cKjT_2EyN@r5SUZYo&wl1(6;upr jiXrNZ_D 1n&ZTkCG&H>cz zmBzj5$}{ciWBOP+hiRM)tfEhrx?!>_+TI#&IuiPsQolava+2N^@T7+rxBeZy*aID| z7GwhcrxX|#@eauMxy8s=s=gEV5}Ez|T4L(;*0~)x-vN2C{Nrw=O AG@ettf;=Cw?UbT;d>rF1InVaXG@=%Yfuh7?(HCghxmQ<7N0 zv8p?cql0RwO%1oYguF2f+Kbu1$Qwey-o>3p*Mem#VTWqPkE4vDfh2y`xM|zeu4>iJ z?HSbfk>ic$S6#bhloboVdn~RObsUrGZs^?<)|g;pyOu*&GG)d;^=wX-U*8op_1$KO z>ZC#&GP8ys>+W=>Wf%7zCc`0Eqlv{1K9J8EOKj~Nt_g=z2SGWU7dXp q|82 zaDzqV1o5VUU*k~c=Xdv{<* 5=6UO>&22
2ZEmrcrtUIwzfl6(9%^gVLMiZZ06kMl4v$MJ^|BL0O8AnP1&GIJJUDw$d} zcHBmYuYL!`5v@17l|5>23rsuy>SL`mv&QbxX4|gTFNVOZUT$-5Kv+9R2fR|7-e`DE zBENs&)rTPyHc#6`df>Z)9)-ni6;?m1r-TpO3+xdtP00D|#R8W>UWov=njEUuPcr*D zCHT!D0ehKtJ^Y~juiugvDS71AVW{MLB-i7rrYE=3)(U_7rJC|}(&`$_eNiBwgRbs^ z7p2~zU=4~?CJNHl*NR1GX=ax8!2L&0`3io`I^dMNU7LBj>-)jB?2=l~)A$2!>4=^# z`DJasp@Ah2`}JV4Be={>vlIIQNs>}D>Sl$L@dFh~9H!k9g 4vpOMK2OlDuR2b;^s>v$2}2;CJ8qsQ;$XQefrO zU_<<{YX7uA7vK7MIS#M*X~v-N{iyD795L1G%(;o1q5m-{5kNVV2oM29Y%Y16)s7ib zhTg{_y`lPudf`-5ar!Z%Pw0@*1j&spW6B-?AV&qcoHn@Zlf3M2THP|0E?FwI>G46L z8RG7 6~b=bKW8# zW9HyXRyFxJV!Xu~Jd$9DxnbRLe;!xnObOlhK2ojFED`<8Z*~+;k9)bMBfSVK6ewLx zoE#q)_{BlVxt(v7v3r-#!{)HT{hMDlyLhGo-dOn#6)q6nHoB^&Ux?9}0>MUW(exT5 zCZG%1{ArKrM#g$a_tfMKWR#hL^1A4`IOn M+9qASP2)zsI zblrcx+f>d%`=R_iRESRPWVkkjvaiqN;JaM62QKwD*zGp9(wN6qFEXcF3bVg@4ZU0& zShO4cdR1Tv?{CV;_Ujo8>(xZzit)irk)nojZR+xlfSruq;0>5j^UKK67uIuNfyrxh zY}GU(uvzNBycYl?MG7IC=WVXB vMPW|k!)On4!ta&9;w6iqtb?3epJ z_BcJBY3Z!;-Pf;zDN!|%9?M5%Pl+ZgbuR>wu@k^;{=nP8mRJj2S?32BZq?&e%Z~km zT#EE}Haamr1_Cso?M<`TJ*9V-xsShlR?WB6?4;EX&!cs^cd%+_mZcCOl|zVw42Kl# zvUx}jbI9R*@?*}NDZNxjMrztOh *hbsCoAd@bK?>UQq5 zG=tL6Q<^<