<commit_before>import h5py as h5
import numpy as np
from scipy.sparse import vstack, coo_matrix
import sys
from math import sqrt
from sklearn.cluster import DBSCAN
import argparse
import multiprocessing
from functools import partial

from .ananke_database import TimeSeriesData

#Input:
#  - Tab-separated table file with sequence hash as rownames, time points as column names, sequence abundances as entries
#  - Absolute sequence abundance minimum threshold
#  - Number of threads
#  - Possible pre-computed distance matrix

#Output:
#  - Cluster label file with clustering parameter as column names, hash as row names, temporal cluster number as entries
#  - Distance matrix (for fast reclustering)


#  Calculate the slopes, mx
def calculate_slopes(matrix, time_points, mask):
    time_points = np.array(time_points)
    border = []
    for i in range(len(mask)-1):
        border.append(mask[i] == mask[i+1])
    border = np.array(border)
    time_difference = time_points[1:] - time_points[0:len(time_points)-1]
    time_difference = time_difference[border]
    if (min(time_difference) <= 0):
        raise ValueError("Minimum time difference is less than or equal to zero (may be caused by two consecutive samples with identical time points)")
    slope_matrix = matrix[:, 1:]-matrix[:, 0:matrix.shape[1]-1]
    slope_matrix = slope_matrix[:,border]
    slope_matrix = slope_matrix / time_difference
    return slope_matrix

#  Handles multithreading of STS distance matrix calculation
def generate_STS_distance_matrix(slope_matrix, nthreads=4):
    sts_dist_matrix = np.zeros(shape=(slope_matrix.shape[0],slope_matrix.shape[0]),dtype='float64')
    nrows = slope_matrix.shape[0]
    p = multiprocessing.Pool(nthreads)
    partial_sts_matrix_generator = partial(sts_matrix_generator, slope_matrix = slope_matrix)
    print("Beginning parallel calculations on %d threads"%nthreads)
    count = 1
    for result in p.imap_unordered(partial_sts_matrix_generator, range(0, nrows-1), 1000):
        ind = result[0]
        dists = result[1].flatten()
        sts_dist_matrix[ind,ind:] = dists
        sts_dist_matrix[ind:,ind] = dists
        count += 1
        sys.stdout.write("\r%d/%d" % (count, nrows))
        sys.stdout.flush()
    p.close()
    p.join()
    return sts_dist_matrix

#  Calculates short time-series distance
#  somewhat efficiently
def sts_matrix_generator(ind, slope_matrix):
    mx = slope_matrix[ind, :]
    mv = slope_matrix[ind:, :]
    mx_rep = np.vstack((mx,)*mv.shape[0])
    diff = mx_rep - mv
    diff = np.square(diff)
    sts_squared = diff.sum(axis=1)
    dists = np.sqrt(sts_squared)
    return (ind, dists)

#  DBSCAN from scikit learn     
def cluster_dbscan(dist_matrix, eps=1, distance_measure="sts"):
    if (distance_measure == "sts"):
        dbs = DBSCAN(eps=eps, metric='precomputed', min_samples=2)
    else:
        dbs = DBSCAN(eps=eps, metric=distance_measure, min_samples=2)
    cluster_labels = dbs.fit_predict(dist_matrix)
    return cluster_labels

#  Main method
def run_cluster(timeseriesdata_path, num_cores, param_min=0.01, param_max=1000, param_step=0.01, distance_measure="sts"):
    print("Loading time-series database file")
    timeseriesdb = TimeSeriesData(timeseriesdata_path)
    print("Importing time-series matrix")
    matrix = timeseriesdb.get_sparse_matrix()
    time_points = timeseriesdb.get_time_points()
    mask = timeseriesdb.get_mask()
    matrix = matrix.todense()
    nrows = matrix.shape[0]
    if nrows <= 1:
        raise ValueError("Time-series matrix contains no information. Was all of your data filtered out?")
    #Normalize the matrix for sequence depth then into Z-scores
    print("Normalizing matrix")
    matrix = matrix/matrix.sum(0)
    #Normalizing is an issue if you have a column that was completely filtered
    #Set these columns back to 0
    matrix[np.invert(np.isfinite(matrix))] = 0
    #Standardize to Z-scores
    norm_matrix = matrix
    for mask_val in np.unique(mask):
        y = norm_matrix[:,np.where(mask==mask_val)[0]]
        y = np.apply_along_axis(lambda x: (x-np.mean(x))/np.std(x),1,y)
        norm_matrix[:,np.where(mask==mask_val)[0]] = y
        del y
    if (distance_measure == "sts"):
        print("Calculating slopes")
        slope_matrix = calculate_slopes(norm_matrix, time_points, mask)
        print("Generating STS distance matrix")<commit_msg>Break out z-score into a proper function<commit_after>import h5py as h5
import numpy as np
from scipy.sparse import vstack, coo_matrix
import sys
from math import sqrt
from sklearn.cluster import DBSCAN
import argparse
import multiprocessing
from functools import partial

from .ananke_database import TimeSeriesData

#Input:
#  - Tab-separated table file with sequence hash as rownames, time points as column names, sequence abundances as entries
#  - Absolute sequence abundance minimum threshold
#  - Number of threads
#  - Possible pre-computed distance matrix

#Output:
#  - Cluster label file with clustering parameter as column names, hash as row names, temporal cluster number as entries
#  - Distance matrix (for fast reclustering)


#  Calculate the slopes, mx
def calculate_slopes(matrix, time_points, mask):
    time_points = np.array(time_points)
    border = []
    for i in range(len(mask)-1):
        border.append(mask[i] == mask[i+1])
    border = np.array(border)
    time_difference = time_points[1:] - time_points[0:len(time_points)-1]
    time_difference = time_difference[border]
    if (min(time_difference) <= 0):
        raise ValueError("Minimum time difference is less than or equal to zero (may be caused by two consecutive samples with identical time points)")
    slope_matrix = matrix[:, 1:]-matrix[:, 0:matrix.shape[1]-1]
    slope_matrix = slope_matrix[:,border]
    slope_matrix = slope_matrix / time_difference
    return slope_matrix

#  Handles multithreading of STS distance matrix calculation
def generate_STS_distance_matrix(slope_matrix, nthreads=4):
    sts_dist_matrix = np.zeros(shape=(slope_matrix.shape[0],slope_matrix.shape[0]),dtype='float64')
    nrows = slope_matrix.shape[0]
    p = multiprocessing.Pool(nthreads)
    partial_sts_matrix_generator = partial(sts_matrix_generator, slope_matrix = slope_matrix)
    print("Beginning parallel calculations on %d threads"%nthreads)
    count = 1
    for result in p.imap_unordered(partial_sts_matrix_generator, range(0, nrows-1), 1000):
        ind = result[0]
        dists = result[1].flatten()
        sts_dist_matrix[ind,ind:] = dists
        sts_dist_matrix[ind:,ind] = dists
        count += 1
        sys.stdout.write("\r%d/%d" % (count, nrows))
        sys.stdout.flush()
    p.close()
    p.join()
    return sts_dist_matrix

#  Calculates short time-series distance
#  somewhat efficiently
def sts_matrix_generator(ind, slope_matrix):
    mx = slope_matrix[ind, :]
    mv = slope_matrix[ind:, :]
    mx_rep = np.vstack((mx,)*mv.shape[0])
    diff = mx_rep - mv
    diff = np.square(diff)
    sts_squared = diff.sum(axis=1)
    dists = np.sqrt(sts_squared)
    return (ind, dists)

#  DBSCAN from scikit learn     
def cluster_dbscan(dist_matrix, eps=1, distance_measure="sts"):
    if (distance_measure == "sts"):
        dbs = DBSCAN(eps=eps, metric='precomputed', min_samples=2)
    else:
        dbs = DBSCAN(eps=eps, metric=distance_measure, min_samples=2)
    cluster_labels = dbs.fit_predict(dist_matrix)
    return cluster_labels

def zscore(x):
    mean = np.mean(x)
    sd = np.std(x)
    if sd == 0:
        z = np.zeros_like(x)
    else:
        z = (x-mean)/sd
    return z

#  Main method
def run_cluster(timeseriesdata_path, num_cores, param_min=0.01, param_max=1000, param_step=0.01, distance_measure="sts"):
    print("Loading time-series database file")
    timeseriesdb = TimeSeriesData(timeseriesdata_path)
    print("Importing time-series matrix")
    matrix = timeseriesdb.get_sparse_matrix()
    time_points = timeseriesdb.get_time_points()
    mask = timeseriesdb.get_mask()
    matrix = matrix.todense()
    nrows = matrix.shape[0]
    if nrows <= 1:
        raise ValueError("Time-series matrix contains no information. Was all of your data filtered out?")
    #Normalize the matrix for sequence depth then into Z-scores
    print("Normalizing matrix")
    matrix = matrix/matrix.sum(0)
    #Normalizing is an issue if you have a column that was completely filtered
    #Set these columns back to 0
    matrix[np.invert(np.isfinite(matrix))] = 0
    #Standardize to Z-scores
    norm_matrix = matrix
    for mask_val in np.unique(mask):
        y = norm_matrix[:,np.where(mask==mask_val)[0]]
        y = np.apply_along_axis(zscore,1,y)
        norm_matrix[:,np.where(mask==mask_val)[0]] = y
        del y
    if (distance_measure == "sts"):
        print("Calculating slopes")
        slope_matrix = calculate_slopes(norm_matrix, time_points, mask)
        print("Generating STS distance matrix")
{'fraction_non_alphanumeric': 0.10285200988097912, 'fraction_numerical': 0.011565236918931058, 'mean_word_length': 3.8276422764227642, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1968935', 'n_tokens_mistral': 2898, 'n_tokens_neox': 2823, 'n_words': 745}