Dendrogram and fast line graph with constant colours

py/clustering.py

@@ -1,4 +1,5 @@
 from util import getQuery, pickleQuery
+import numpy as np
 import pandas as p
 import matplotlib.pyplot as plt
 import seaborn as sns
@@ -34,7 +35,8 @@ from scipy.cluster.hierarchy import dendrogram, linkage, cophenet, fcluster
 # 
 # plt.show()
 
-df = p.read_pickle('../data/jan19s.pkl')
+numclusts = 7
+df = p.read_pickle('../data/2017-20s.pkl')
 dforig = df
 
 print(df.info())
@@ -53,36 +55,62 @@ lobj = linkage(lmat, method = 'ward')
 print(lobj)
 print(cophenet(lobj, lmat))
 
-#plt.figure(figsize = (25, 10))
-#plt.title('ICP Clustering Dendrogram')
-#plt.xlabel('ICP ID/(Number of ICPs)')
-#plt.ylabel('distance')
-#dendrogram(
-#    lobj,
-#    labels = cmat.index.values,
-#    leaf_rotation=90,
-#    leaf_font_size=8,
-#    #show_leaf_counts = True,
-#    #truncate_mode = 'lastp',
-#    #p = 50,
-#    #show_contracted = True,
-#    color_threshold = 1.9
-#)
-#plt.show()
 
-clusts = fcluster(lobj, 6, criterion='maxclust')
+
+clabs = [x + 1 for x in range(numclusts)]
+cpal = dict(zip(clabs, sns.color_palette("colorblind", numclusts).as_hex()))
+
+clusts = fcluster(lobj, numclusts, criterion='maxclust')
 print(clusts)
 print(cmat.index.values)
-clustdf = p.DataFrame({'icp_id' : cmat.index.values, 'cluster' : [chr(x + ord('A') - 1) for x in clusts]})
+clustdf = p.DataFrame({'icp_id' : cmat.index.values, 'cluster' : clusts})
 print(clustdf)
 mdf = p.merge(clustdf, dforig, on = 'icp_id', how = 'left')
 print(mdf)
 print(mdf.info())
+qlow  = lambda x: x.quantile(0.250)
+qhigh = lambda x: x.quantile(0.750)
 print(mdf.cluster.describe())
+mdagg = mdf.groupby(['read_time', 'cluster']).agg({
+        'kwh_tot': ['median', 'mean', ('CI_low', qlow), ('CI_high', qhigh)]
+}, q = 0.025)
+mdagg.columns = ['_'.join(x) for x in mdagg.columns.ravel()]
+mdagg = mdagg.reset_index()
+print(mdagg)
+print(mdagg.info())
+print(mdagg.describe())
+# mdf.to_csv('~/windows/Documents/clusters-ward.csv')
 
-mdf.to_csv('~/windows/Documents/clusters-ward.csv')
+# Algorithm via 
+# <https://stackoverflow.com/questions/38153829/custom-cluster-colors-of-scipy-dendrogram-in-python-link-color-func>
+ldict = {icp_id:cpal[cluster] for icp_id, cluster in zip(clustdf.icp_id, clustdf.cluster)}
+link_cols = {}
+for i, i12 in enumerate(lobj[:,:2].astype(int)):
+  c1, c2 = (link_cols[x] if x > len(lobj) else ldict[clustdf.icp_id[x]]
+    for x in i12)
+  link_cols[i+1+len(lobj)] = c1 if c1 == c2 else '#000000'
 
-sns.set()
+plt.figure(figsize = (25, 10))
+plt.title('ICP Clustering Dendrogram')
+plt.xlabel('ICP ID/(Number of ICPs)')
+plt.ylabel('distance')
+dendrogram(
+    lobj,
+    labels = cmat.index.values,
+    leaf_rotation=90,
+    leaf_font_size=8,
+    # show_leaf_counts = True,
+    # truncate_mode = 'lastp',
+    # p = 50,
+    # show_contracted = True,
+    link_color_func = lambda x: link_cols[x],
+    color_threshold = None
+)
+plt.show()
 
-sns.lineplot(x = 'read_time', y = 'kwh_tot', hue = 'cluster', data = mdf)
+sns.set()
+ax = sns.lineplot(x = 'read_time', y = 'kwh_tot_mean', hue = 'cluster', data = mdagg, palette = cpal)
+for c in clabs:
+    fds = mdagg[mdagg.cluster == c]
+    ax.fill_between(fds.read_time.dt.to_pydatetime(), fds.kwh_tot_CI_low, fds.kwh_tot_CI_high, alpha = 0.1, color = cpal[c])
 plt.show()

py/downkwh.py

@@ -1,40 +1,47 @@
-from util import getQuery, pickleQuery
+from util import getQuery, pickleQuery, getkwh
 import pandas as p
 import matplotlib.pyplot as plt
 import seaborn as sns
 
-query = """
-SELECT comb.icp_id, comb.read_time, COALESCE(kwh_tot, 0) AS kwh_tot
-FROM
-(
-    SELECT read_time, icp_id
-    FROM
-    (
-        SELECT read_time 
-        FROM GENERATE_SERIES('2017-01-01 00:30:00'::timestamp, '2017-02-01 00:00:00'::timestamp, 
-            '30 minutes'::interval) read_time
-    ) AS tsdata CROSS JOIN
-    (
-        SELECT *
-        FROM
-        (
-            SELECT icp_id, COUNT(DISTINCT read_date) AS data_days 
-            FROM coup_prd.coupdatamaster
-            WHERE read_date >= to_date('01/01/2017','dd/mm/yyyy')
-                AND read_date <  to_date('01/01/2018','dd/mm/yyyy')
-                AND content_code = 'UN'
-                AND icp_id LIKE '%%19'
-            GROUP BY icp_id
-        ) AS cir 
-        WHERE data_days >= 360
-    ) AS qual_icp
-) AS comb
-LEFT JOIN
-(
-    SELECT *, read_date + CONCAT(period / 2, ':', period %% 2 * 30, ':00')::time AS read_time
-    FROM public.coup_tall_jan
-) AS tall_timestamp 
-ON comb.read_time = tall_timestamp.read_time AND comb.icp_id = tall_timestamp.icp_id;
-"""
+# query = """
+# SELECT comb.icp_id, comb.read_time, COALESCE(kwh_tot, 0) AS kwh_tot
+# FROM
+# (
+#     SELECT read_time, icp_id
+#     FROM
+#     (
+#         SELECT read_time 
+#         FROM GENERATE_SERIES('2017-01-01 00:30:00'::timestamp, '2017-02-01 00:00:00'::timestamp, 
+#             '30 minutes'::interval) read_time
+#     ) AS tsdata CROSS JOIN
+#     (
+#         SELECT *
+#         FROM
+#         (
+#             SELECT icp_id, COUNT(DISTINCT read_date) AS data_days 
+#             FROM coup_prd.coupdatamaster
+#             WHERE read_date >= to_date('01/01/2017','dd/mm/yyyy')
+#                 AND read_date <  to_date('01/01/2018','dd/mm/yyyy')
+#                 AND content_code = 'UN'
+#                 AND icp_id LIKE '%%19'
+#             GROUP BY icp_id
+#         ) AS cir 
+#         WHERE data_days >= 360
+#     ) AS qual_icp
+# ) AS comb
+# LEFT JOIN
+# (
+#     SELECT *, read_date + CONCAT(period / 2, ':', period %% 2 * 30, ':00')::time AS read_time
+#     FROM public.coup_tall_jan
+# ) AS tall_timestamp 
+# ON comb.read_time = tall_timestamp.read_time AND comb.icp_id = tall_timestamp.icp_id;
+# """
+# 
+# pickleQuery(query, "../data/jan19s.pkl")
 
-pickleQuery(query, "../data/jan19s.pkl")
+kwhdata = getkwh('2017-01-01', '2018-01-01', '2017-01-01 00:30:00', '2018-01-01 00:00:00', '%%20')
+
+print(kwhdata.info())
+
+print("Pickling")
+kwhdata.to_pickle("../data/2017-20s.pkl")

py/util.py

@@ -1,6 +1,7 @@
 import psycopg2 as pg
 from configparser import ConfigParser
 import pandas.io.sql as psql
+import datetime as dt
 
 
 def config(filename='database.ini', section='postgresql'):
@@ -39,6 +40,7 @@ def getQuery(query, qparams=[]):
         cur = conn.cursor()
 
         # Get table
+        print("Retrieving table")
         dataframe = psql.read_sql(query, conn, params=qparams)
 
         cur.close()
@@ -63,6 +65,70 @@ def pickleQuery(query, path, qparams=[]):
     print("Table pickled")
 
 
+def getkwh(datestart, dateend, timestart, timeend, subset):
+    query = """
+    SELECT comb.icp_id, comb.read_time, COALESCE(kwh_tot, 0) AS kwh_tot
+    FROM
+    (
+        SELECT read_time, icp_id
+        FROM
+        (
+            SELECT read_time 
+            FROM GENERATE_SERIES(%(tsstart)s::timestamp, %(tsend)s::timestamp, 
+                '30 minutes'::interval) read_time
+        ) AS tsdata CROSS JOIN
+        (
+            SELECT *
+            FROM
+            (
+                SELECT icp_id, COUNT(DISTINCT read_date) AS data_days 
+                FROM coup_prd.coupdatamaster
+                WHERE read_date >= to_date('01/01/2017','dd/mm/yyyy')
+                    AND read_date <  to_date('01/01/2018','dd/mm/yyyy')
+                    AND content_code = 'UN'
+                    AND icp_id LIKE %(subset)s
+                GROUP BY icp_id
+            ) AS cir 
+            WHERE data_days >= 360
+        ) AS qual_icp
+    ) AS comb
+    LEFT JOIN
+    (
+        SELECT *, read_date + CONCAT(period / 2, ':', period %% 2 * 30, ':00')::time AS read_time
+        FROM (
+            SELECT  a.icp_id
+                 , a.read_date
+                 , c.period
+                 , sum(c.read_kwh) as kwh_tot
+                 , sum(case when a.content_code = 'UN' then c.read_kwh else 0 end) as kwh_un
+                 , sum(case when a.content_code in ('CN','EG') then c.read_kwh else 0 end) as kwh_cn
+            FROM    coup_prd.coupdatamaster a,
+                unnest(a.read_array) WITH ORDINALITY c(read_kwh, period)
+            WHERE   a.read_date >= to_date(%(datestart)s,'yyyy-mm-dd')
+             and   a.read_date <  to_date(%(dateend)s,'yyyy-mm-dd')
+             and   a.content_code  ~ ('UN|CN|EG')
+             AND   a.icp_id LIKE %(subset)s
+            GROUP BY 1, 2, 3
+        ) AS coup_tall
+    ) AS tall_timestamp 
+    ON comb.read_time = tall_timestamp.read_time AND comb.icp_id = tall_timestamp.icp_id;
+    """
+    pdict = {
+        'datestart': datestart,
+        'dateend': dateend,
+        'tsstart': timestart,
+        'tsend': timeend,
+        'subset': subset
+    }
+    print("Getting data with parameters:")
+    print(pdict)
+    qdf = getQuery(query, pdict)
+    print("Done")
+    return(qdf)
+
+
+
+
 if __name__ == "__main__":
     dv = getQuery('SELECT version()').version[0]
     print('PostgreSQL database version:')

sql/queries.pgsql

@@ -298,3 +298,49 @@ WHERE   a.read_date >= to_date('01/01/2017','dd/mm/yyyy')
  and   a.content_code  ~ ('UN|CN|EG')
 GROUP BY 1, 2, 3
 ORDER BY 1, 2, 3;
+
+-- Really big query (jan)
+SELECT comb.icp_id, comb.read_time, COALESCE(kwh_tot, 0) AS kwh_tot
+FROM
+(
+    SELECT read_time, icp_id
+    FROM
+    (
+        SELECT read_time 
+        FROM GENERATE_SERIES('2017-01-01 00:30:00'::timestamp, '2017-02-01 00:00:00'::timestamp, 
+            '30 minutes'::interval) read_time
+    ) AS tsdata CROSS JOIN
+    (
+        SELECT *
+        FROM
+        (
+            SELECT icp_id, COUNT(DISTINCT read_date) AS data_days 
+            FROM coup_prd.coupdatamaster
+            WHERE read_date >= to_date('01/01/2017','dd/mm/yyyy')
+                AND read_date <  to_date('01/01/2018','dd/mm/yyyy')
+                AND content_code = 'UN'
+                AND icp_id LIKE '%%19'
+            GROUP BY icp_id
+        ) AS cir 
+        WHERE data_days >= 360
+    ) AS qual_icp
+) AS comb
+LEFT JOIN
+(
+    SELECT *, read_date + CONCAT(period / 2, ':', period %% 2 * 30, ':00')::time AS read_time
+    FROM (
+        SELECT  a.icp_id
+             , a.read_date
+             , c.period
+             , sum(c.read_kwh) as kwh_tot
+             , sum(case when a.content_code = 'UN' then c.read_kwh else 0 end) as kwh_un
+             , sum(case when a.content_code in ('CN','EG') then c.read_kwh else 0 end) as kwh_cn
+        FROM    coup_prd.coupdatamaster a,
+            unnest(a.read_array) WITH ORDINALITY c(read_kwh, period)
+        WHERE   a.read_date >= to_date('01/01/2017','dd/mm/yyyy')
+         and   a.read_date <  to_date('01/02/2017','dd/mm/yyyy')
+         and   a.content_code  ~ ('UN|CN|EG')
+        GROUP BY 1, 2, 3
+    ) AS coup_tall
+) AS tall_timestamp 
+ON comb.read_time = tall_timestamp.read_time AND comb.icp_id = tall_timestamp.icp_id;