distributed PCA

Author: Marigold

python-correlation-heatmap/README.md

@@ -6,13 +6,39 @@
 
 Calculate correlation heatmap, only works for **real variables**.
 
-You can run it on single node with `compute` or in a distributed way with
+You can run it on single node with `compute` and env variable `MODEL_PARAM_graph=correlation_heatmap`.
+
+1. `compute`
+
+Or in a distributed mode with
 
 1. `compute --mode intermediate`
 2. `compute --mode aggregate --job-ids 1 2 3`
 
 Intermediate mode calculates covariance matrix from a single node, while aggregate mode is used after intermediate to
-combine statistics from multiple jobs.
+combine statistics from multiple jobs and produce the final graph.
+
+
+# Python Distributed PCA
+
+Calculate PCA and return biplot visualization and screeplot. It only works for **real variables**.
+
+You can run it on single node with `compute` and env variable `MODEL_PARAM_graph=pca`.
+
+1. `compute`
+
+Or in a distributed mode with
+
+1. `compute --mode intermediate`
+2. `compute --mode aggregate --job-ids 1 2 3`
+
+<!--
+Proposal for distributed mode with PCA scores graph included. See https://trello.com/c/jfLav9K6/58-distributed-pca for
+discussion
+1. `compute --mode intermediate` (calculate covariance matrices on nodes)
+2. `compute --mode aggregate --job-ids 1 2 3` (calculate aggregated correlation matrix)
+3. `compute --mode intermediate --agg-job-id 4` (calculate covariance matrices and also sample scores for PCA)
+4. `compute --mode aggregate --job-ids 5 6 7 --graph pca` (produce plotly visualization) -->
 
 
 ## Build (for contributors)

python-correlation-heatmap/correlation_heatmap.py

@@ -107,7 +107,7 @@ def _compute_intermediate_result(inputs):
 
 
 @utils.catch_user_error
-def aggregate_stats(job_ids):
+def aggregate_stats(job_ids, graph_type=None):
     """Get all partial statistics from all nodes and aggregate them.
     :input job_ids: list of job_ids with intermediate results
     """
@@ -117,12 +117,14 @@ def aggregate_stats(job_ids):
 
     corr, columns = _aggregate_results(results)
 
-    graph_type = parameters.get_parameter('MODEL_PARAM_graph', str, 'correlation_heatmap')
+    graph_type = graph_type or parameters.get_parameter('graph', str, 'correlation_heatmap')
 
     if graph_type == 'correlation_heatmap':
         _save_corr_heatmap(corr, columns)
     elif graph_type == 'pca':
-        raise errors.UserError('MODEL_PARAM_graph=pca is not yet supported for distributed mode.')
+        # save PCA graphs, but leave out the one with PCA scores
+        logging.warning('Sample scores graph is not yet implemented for distributed PCA.')
+        _save_pca(corr, columns, X=None)
     else:
         raise errors.UserError('MODEL_PARAM_graph only supports values `correlation_heatmap` and `pca`')
 
@@ -168,7 +170,7 @@ def _save_corr_heatmap(corr, columns):
     logging.info("DONE")
 
 
-def _pca(corr, X):
+def _pca(corr, X=None):
     # calculate eigenvectors and eigenvalues
     eig_vals, eig_vecs = np.linalg.eig(corr)
 
@@ -187,43 +189,55 @@ def _pca(corr, X):
 
     # convert original data to scores
     # NOTE: since we are working with correlation matrix, original data must be standardized first!
-    X_std = (X - X.mean()) / X.std()
-    Y = X_std.dot(W)
+    if X is not None:
+        X_std = (X - X.mean()) / X.std()
+        Y = X_std.dot(W)
+    else:
+        Y = None
 
     return eig_vals, eig_vecs, Y
 
 
-def _figure(eig_vals, eig_vecs, Y, X):
+def _figure(eig_vals, eig_vecs, Y, columns):
+    show_scores = Y is not None
+
+    titles = ['Scree plot', 'Eigen-components', 'Variables scores']
+    titles.append('Samples scores' if show_scores else 'Samples scores <br>(not available in distributed mode)')
+
     # plotting
     fig = tools.make_subplots(
-        rows=2, cols=2, subplot_titles=('Samples scores', 'Variables scores', 'Scree plot', 'Eigen-components')
+        rows=2, cols=2, subplot_titles=titles
     )
 
-    for d in _biplot_samples(Y.values):
+    for d in _screeplot(eig_vals):
         fig.append_trace(d, 1, 1)
 
-    for d in _biplot_variables(eig_vecs, list(X.columns)):
+    for d in _eigencomponents(eig_vecs, columns):
         fig.append_trace(d, 1, 2)
 
-    for d in _screeplot(eig_vals):
+    for d in _biplot_variables(eig_vecs, columns):
         fig.append_trace(d, 2, 1)
 
-    for d in _eigencomponents(eig_vecs, list(X.columns)):
-        fig.append_trace(d, 2, 2)
+    # only show sample scores in single node mode
+    if show_scores:
+        for d in _biplot_samples(Y.values):
+            fig.append_trace(d, 2, 2)
 
     var_exp = _explained_variance(eig_vals)
 
-    fig['layout']['xaxis1'].update(title='PC1 ({:.1%})'.format(var_exp[0]))
-    fig['layout']['yaxis1'].update(title='PC2 ({:.1%})'.format(var_exp[1]))
-    fig['layout']['xaxis2'].update(title='PC1 ({:.1%})'.format(var_exp[0]), range=[-1.05, 1.05])
-    fig['layout']['yaxis2'].update(title='PC2 ({:.1%})'.format(var_exp[1]), range=[-1.05, 1.05])
-    fig['layout']['yaxis3'].update(title='Explained variance in percent')
+    fig['layout']['yaxis1'].update(title='Explained variance in percent')
+    fig['layout']['xaxis3'].update(title='PC1 ({:.1%})'.format(var_exp[0]), range=[-1.05, 1.05])
+    fig['layout']['yaxis3'].update(title='PC2 ({:.1%})'.format(var_exp[1]), range=[-1.05, 1.05])
+
+    if show_scores:
+        fig['layout']['xaxis4'].update(title='PC1 ({:.1%})'.format(var_exp[0]))
+        fig['layout']['yaxis4'].update(title='PC2 ({:.1%})'.format(var_exp[1]))
 
     # unit-circle for biplot
     circle = {
         'type': 'circle',
-        'xref': 'x2',
-        'yref': 'y2',
+        'xref': 'x3',
+        'yref': 'y3',
         'x0': -1,
         'y0': -1,
         'x1': 1,
@@ -238,10 +252,10 @@ def _figure(eig_vals, eig_vecs, Y, X):
     return fig
 
 
-def _save_pca(corr, columns, X):
+def _save_pca(corr, columns, X=None):
     """Generate PCA visualization in plotly format. Inspired by https://plot.ly/ipython-notebooks/principal-component-analysis/"""
     eig_vals, eig_vecs, Y = _pca(corr, X)
-    fig = _figure(eig_vals, eig_vecs, Y, X)
+    fig = _figure(eig_vals, eig_vecs, Y, columns)
 
     logging.info("Results:\n{}".format(fig))
     io_helper.save_results(json.dumps(fig), shapes.Shapes.PLOTLY)

python-correlation-heatmap/tests/Dockerfile

@@ -88,6 +88,15 @@ services:
       MODEL_PARAM_design: "factorial"
       MODEL_PARAM_graph: "correlation_heatmap"
 
+  correlation-heatmap-single:
+    extends: correlation-heatmap-base
+    container_name: "correlation-heatmap-single"
+    environment:
+      JOB_ID: '1'
+      PARAM_query: "SELECT lefthippocampus, minimentalstate, opticchiasm, subjectageyears FROM cde_features_a LIMIT 100"
+    links:
+      - "db:db"
+
   correlation-heatmap-a:
     extends: correlation-heatmap-base
     container_name: "correlation-heatmap-a"
@@ -114,7 +123,7 @@ services:
     links:
       - "db:db"
 
-  distributed-pca:
+  distributed-pca-base:
     image: "hbpmip/python-distributed-pca:latest"
     container_name: "distributed-pca"
     restart: "no"
@@ -136,9 +145,42 @@ services:
       PARAM_variables: "lefthippocampus"
       PARAM_covariables: "minimentalstate,opticchiasm,subjectageyears"
       PARAM_grouping: ""
-      PARAM_query: "SELECT lefthippocampus, minimentalstate, opticchiasm, subjectageyears FROM cde_features_a LIMIT 100"
       PARAM_meta: "{\"lefthippocampus\":{\"code\":\"lefthippocampus\",\"type\":\"real\",\"mean\":3.0,\"std\":0.35},\"minimentalstate\":{\"code\":\"minimentalstate\",\"type\":\"real\",\"mean\":24.0,\"std\":5.0},\"opticchiasm\":{\"code\":\"opticchiasm\",\"type\":\"real\",\"mean\":0.08,\"std\":0.009},\"subjectage\":{\"code\":\"subjectage\",\"type\":\"real\",\"mean\":71.0,\"std\":8.0},\"subjectageyears\":{\"description\":\"Subject age in years.\",\"methodology\":\"mip-cde\",\"label\":\"Age Years\",\"minValue\":0,\"code\":\"subjectageyears\",\"units\":\"years\",\"length\":3,\"maxValue\":130.0,\"type\":\"integer\"}}"
       MODEL_PARAM_design: "factorial"
       MODEL_PARAM_graph: "pca"
+
+  distributed-pca-single:
+    extends: distributed-pca-base
+    container_name: "distributed-pca-single"
+    environment:
+      JOB_ID: '1'
+      PARAM_query: "SELECT lefthippocampus, minimentalstate, opticchiasm, subjectageyears FROM cde_features_a LIMIT 100"
+    links:
+      - "db:db"
+
+  distributed-pca-a:
+    extends: distributed-pca-base
+    container_name: "distributed-pca-a"
+    environment:
+      JOB_ID: '1'
+      PARAM_query: "SELECT lefthippocampus, minimentalstate, opticchiasm, subjectageyears FROM cde_features_a LIMIT 100"
+    links:
+      - "db:db"
+
+  distributed-pca-b:
+    extends: distributed-pca-base
+    container_name: "distributed-pca-b"
+    environment:
+      JOB_ID: '2'
+      PARAM_query: "SELECT lefthippocampus, minimentalstate, opticchiasm, subjectageyears FROM cde_features_b LIMIT 100"
+    links:
+      - "db:db"
+
+  distributed-pca-agg:
+    extends: distributed-pca-base
+    container_name: "distributed-pca-agg"
+    environment:
+      JOB_ID: '3'
+      PARAM_query: "SELECT lefthippocampus, minimentalstate, opticchiasm, subjectageyears FROM cde_features_b LIMIT 100"
     links:
       - "db:db"

python-correlation-heatmap/tests/docker-compose.yml

@@ -46,10 +46,12 @@ echo "Initialise the databases..."
 $DOCKER_COMPOSE run sample_data_db_setup
 $DOCKER_COMPOSE run woken_db_setup
 
+## single-node mode
 echo
-echo "Run the distributed-pca algorithm..."
-$DOCKER_COMPOSE run distributed-pca compute
+echo "Run the correlation-heatmap algorithm on single node..."
+$DOCKER_COMPOSE run correlation-heatmap-single compute
 
+## distributed mode
 # echo
 # echo "Run the correlation-heatmap-a..."
 # $DOCKER_COMPOSE run correlation-heatmap-a compute --mode intermediate

python-correlation-heatmap/tests/test.sh

@@ -0,0 +1,65 @@
+#!/usr/bin/env bash
+
+set -o pipefail  # trace ERR through pipes
+set -o errtrace  # trace ERR through 'time command' and other functions
+set -o errexit   ## set -e : exit the script if any statement returns a non-true return value
+
+get_script_dir () {
+     SOURCE="${BASH_SOURCE[0]}"
+
+     while [ -h "$SOURCE" ]; do
+          DIR="$( cd -P "$( dirname "$SOURCE" )" && pwd )"
+          SOURCE="$( readlink "$SOURCE" )"
+          [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE"
+     done
+     cd -P "$( dirname "$SOURCE" )"
+     pwd
+}
+
+cd "$(get_script_dir)"
+
+if [[ $NO_SUDO || -n "$CIRCLECI" ]]; then
+  DOCKER_COMPOSE="docker-compose"
+elif groups $USER | grep &>/dev/null '\bdocker\b'; then
+  DOCKER_COMPOSE="docker-compose"
+else
+  DOCKER_COMPOSE="sudo docker-compose"
+fi
+
+function _cleanup() {
+  local error_code="$?"
+  echo "Stopping the containers..."
+  $DOCKER_COMPOSE stop | true
+  $DOCKER_COMPOSE down | true
+  $DOCKER_COMPOSE rm -f > /dev/null 2> /dev/null | true
+  exit $error_code
+}
+trap _cleanup EXIT INT TERM
+
+echo "Starting the databases..."
+$DOCKER_COMPOSE up -d --remove-orphans db
+$DOCKER_COMPOSE run wait_dbs
+$DOCKER_COMPOSE run create_dbs
+
+echo
+echo "Initialise the databases..."
+$DOCKER_COMPOSE run sample_data_db_setup
+$DOCKER_COMPOSE run woken_db_setup
+
+# single-node mode
+# echo
+# echo "Run the distributed-pca algorithm on single node..."
+# $DOCKER_COMPOSE run distributed-pca-single compute
+
+## distributed mode
+echo
+echo "Run the distributed-pca-a..."
+$DOCKER_COMPOSE run distributed-pca-a compute --mode intermediate
+echo "Run the distributed-pca-b..."
+$DOCKER_COMPOSE run distributed-pca-b compute --mode intermediate
+echo "Run the distributed-pca-agg..."
+$DOCKER_COMPOSE run distributed-pca-agg compute --mode aggregate --job-ids 1 2
+
+echo
+# Cleanup
+_cleanup

python-correlation-heatmap/tests/test_pca.sh

@@ -105,7 +105,7 @@ def intermediate_data_2():
 @mock.patch('correlation_heatmap.io_helper.fetch_data')
 @mock.patch('correlation_heatmap.io_helper.get_results')
 @mock.patch('correlation_heatmap.io_helper.save_results')
-def test_aggregate_stats(mock_save_results, mock_get_results, mock_fetch):
+def test_aggregate_stats_correlation_heatmap(mock_save_results, mock_get_results, mock_fetch):
 
     def mock_results(job_id):
         job_id = str(job_id)
@@ -116,7 +116,7 @@ def mock_results(job_id):
 
     mock_get_results.side_effect = mock_results
 
-    aggregate_stats([1, 2])
+    aggregate_stats([1, 2], graph_type='correlation_heatmap')
     results = json.loads(mock_save_results.call_args[0][0])
     assert round_dict(results) == [
         {
@@ -128,3 +128,22 @@ def mock_results(job_id):
             'zmax': 1
         }
     ]
+
+
+@mock.patch('correlation_heatmap.io_helper.fetch_data')
+@mock.patch('correlation_heatmap.io_helper.get_results')
+@mock.patch('correlation_heatmap.io_helper.save_results')
+def test_aggregate_stats_pca(mock_save_results, mock_get_results, mock_fetch):
+
+    def mock_results(job_id):
+        job_id = str(job_id)
+        if job_id == '1':
+            return mock.MagicMock(data=json.dumps(intermediate_data_1()))
+        elif job_id == '2':
+            return mock.MagicMock(data=json.dumps(intermediate_data_2()))
+
+    mock_get_results.side_effect = mock_results
+
+    aggregate_stats([1, 2], graph_type='pca')
+    results = json.loads(mock_save_results.call_args[0][0])
+    assert set(results.keys()) == {'layout', 'data'}