text:
#!/usr/bin/env python
import numpy as np
import pandas as pd
from sklearn.metrics import confusion_matrix
from os import path

_README="""
A script to generate confusion matrix and evaluate accuracy of rfmix.

-Matthew Aguirre (magu[at]stanford[dot]edu)
"""

_TODO="""
1. modify output directory (currently the output text file is stored in the directory that this script is called from).
"""

# quick check if we're on galangal
import platform
if platform.uname()[1]=='galangal.stanford.edu':
    print('error: script must be modified to be run on galangal')
else:
    # assume we're on sherlock -- load modules and check versions
    print('Assuming we are on sherlock')

# define functions
def load_data(args_dict):
    print('Loading in data')
    ## check if output file already exists
    if path.exists(args_dict['output-filename']):
        print('Error: output file already exists. Aborting script.')
        return '', '', True
    ## load y and yhat_raw
    data_dir='/scratch/users/magu/deepmix/data/simulated_chr20/'
    yhat_raw=pd.read_table(data_dir+'vcf/rf_out/'+args_dict['rfmix-result-filename'], skiprows=1)
    y=np.load(data_dir+'label/'+args_dict['gt-filename'])

    return y, yhat_raw, False

def expand_rfmix_windows(y, yhat_raw, S):
    print('Expanding rfmix windows')
    V_pos=y['V'][:,1].astype(int)
    yhat=pd.DataFrame(index=['_'.join(s) for s in y['V']], columns=S)

    for ix in range(yhat_raw.shape[0]):
        ids=(yhat_raw.iloc[ix,1] <= V_pos) & (V_pos <= yhat_raw.iloc[ix,2])
        yhat.iloc[ids,:]=np.vstack([yhat_raw.iloc[ix,6:] for _ in range(sum(ids))]).astype(int)+1

    return yhat

def evaluate_model(y, yhat, args_dict):
    print('Evaluating model and creating text file')
    ## create df of confusion matrices and evaluate accuracy
    # confusion
    cm=confusion_matrix(y['L'].flatten(), yhat.T.values.flatten().astype(int))
    # accuracy
    acc=np.sum(np.diag(cm))/np.sum(cm)
    anc_label=['AFR', 'EAS', 'EUR', 'NAT', 'SAS']
    row_normalized_df = pd.DataFrame(cm, index=anc_label, columns=anc_label).divide(cm.sum(axis=1), axis=0)
    col_normalized_df = pd.DataFrame(cm, index=anc_label, columns=anc_label).divide(cm.sum(axis=0), axis=1)
    bp_df = pd.DataFrame(cm, index=anc_label, columns=anc_label)

    ## write df and accuracy to text file
    output_file_handle = open(args_dict['output-filename'],"w")
    output_file_handle.writelines('Row-normalized confusion matrix:\n\n')
    output_file_handle.writelines(row_normalized_df.to_string())
    output_file_handle.writelines('\n\n\n\nCol-normalized confusion matrix:\n\n')
    output_file_handle.writelines(col_normalized_df.to_string())
    output_file_handle.writelines('\n\n\n\nBP confusion matrix:\n\n')
    output_file_handle.writelines(bp_df.to_string())
    output_file_handle.writelines('\n\n\n\nmodel accuracy = '+str(acc))
    output_file_handle.close()

def get_args():
    import argparse
    parser=argparse.ArgumentParser(description=_README)
    parser.add_argument('rfmix-result-filename', metavar='rfmix_result', type=str,
                         help='filename of rfmix output (.msp.tsv file)')
    parser.add_argument('gt-filename', metavar='groundtruth', type=str,
                         help='filename of ground truth labels (.npz file)')
    parser.add_argument('output-filename', metavar='output', type=str,
                         help='output filename (.txt file)')
    args=parser.parse_args()
    return args

def main():
    args=get_args()
#    print(args)
    y, yhat_raw, abort = load_data(vars(args))
    if abort: return
    S = np.array([s.replace('_S1','.0').replace('_S2','.1') for s in y['S']]) # match samples
    yhat = expand_rfmix_windows(y, yhat_raw, S) # expand rfmix windows
    evaluate_model(y, yhat, vars(args)) # evaluate model and save accuracy & confusion matrices to text file
    return

if __name__=='__main__':
    main()

reason: import numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# """
# - Author: steve simmert
# - Copyright: 2017
# """
# Model function of IRM intensity profile of a tilted MT and TIRF intensity profile
from scipy import sin, cos, arcsin, pi, sinc

def height_fun(x, x0, A, xt):
    """
    Model the height of a cantilevered microubule.

    The MT is fixed at x0 and pulled up. A is the shape giving/ scaling factor,
    equal to F/EI, where F ist the forc acting on the MT at xt, and EI is the
    flexural regidity of the mictotubule.
    """
    def heaviside(x, x0):
        return 1 * (x >= x0)

    return heaviside(x, x0) * A * ((xt - x0) / 2 * (x - x0)**2 - ((x - x0)**3) / 6)

def irm_intensity_vs_height(h, bg, d, INA, wl, n1):
    """
    Intensity of an IRM signal, given the background and difference
    signal and the height of the scatting object.

    All length units should be in meters!

    bg could be, for example, 0.0 for a perfect gray value.
    Rädler & Sackmann, Eq. 4 in Limozion et al. 2009
    """
    theta = arcsin(INA / n1)
    k = 2 * pi * n1 / wl
    # following division by pi is done because of the scipy.sinc
    # definition: sinc(x) = sin(pi*x)/(pi*x)
    y =  2 * k * h * (sin(theta / 2))**2 / pi

    i = bg - d * sinc(y) * cos(2 * k * (h * (1 - (sin(theta/2))**2)))

    return i

def irm_signal_bent_mt(x, x0, bg, d, A, xt, INA=1.33, wl=0.448, n1=1.33):
    """
    Intensity profile of a bent and cantilevered microtubule.

    Here the position and wavelength values should be in µm!
    """
    # conversion to meters
    h = height_fun(x, x0, A, xt) * 1e-6  # conversion to meters

    return irm_intensity_vs_height(h, bg, d, INA, wl * 1e-6, n1)

ext: py

reason: from scipy