<commit_before>#!/usr/bin/env python

import sys, os
from rasmus import util, env, treelib, spidirlib, stats, tablelib

import rpy


options = [
  ["s:", "stree=", "stree", "<species tree>",
    {"single": True}],
  ["p:", "params=", "params", "<sindir params file>",
    {"single": True}],
  ["l:", "lens=", "lens", "<sindir branch lengths file>",
    {"single": True}],
  ["o:", "outdir=", "outdir", "<output directory>",
    {"single": True,
     "default": "."}],
  ["b:", "nbins=", "nbins", "<number of bins>",
    {"single": True,
     "parser": int,
     "default": 20}]
]


# parse options
conf = util.parseOptions(sys.argv, options, quit=True)


# read data
env.addEnvPaths("DATAPATH")
stree = treelib.readTree(env.findFile(conf["stree"]))
params = spidirlib.readParams(env.findFile(conf["params"]))
lens = spidirlib.readTreeDistrib(env.findFile(conf["lens"]))
totals = map(sum, zip(* lens.values()))
rlens = util.mapdict(lens, valfunc=lambda x: util.vidiv(x, totals))

fitting = tablelib.Table(headers=["name", "type", "param1", "param2", "pval"])


def plotAbsLens(name, lens, low, high, step):
    p = util.Gnuplot()
    p.enableOutput(False)
    
    lens = filter(util.withinfunc(low, high), lens)
    
    lens2 = filter(util.gtfunc(.001), lens)
    mu = stats.mean(lens2)
    sigma2 = stats.variance(lens2)
    prms = [mu*mu/sigma2, mu/sigma2]
    
    prms, resid = stats.fitDistrib(stats.gammaPdf, prms, 
                                   lens2, low, high, step)
    
    util.plotdistrib(lens, low=low, width=step, plot=p)
    p.plotfunc(lambda x: stats.gammaPdf(x, prms), low, high, step/4.)
    
    
    bins, obs = util.hist(lens, low=low, width=step)
    tot = len(lens)
    ext = map(lambda x: tot * (stats.gammaCdf(x+step, prms) - 
                               stats.gammaCdf(x, prms)), bins)
    
    ind = util.find(lambda x: x>5, ext)
    ind = ind[1:]
    obs = util.mget(obs, ind)
    bins = util.mget(bins, ind)
    ext = util.mget(ext, ind)
    
    try:
        chisq = rpy.r.chisq_test(obs, p=util.oneNorm(ext))
        pval = chisq["p.value"]
    except rpy.RException, e:
        print e
        pval = 0.0
    
    
    util.logger("%s\t%.3e\t%s" % (name, pval, str(pval > .05)))
    
    p.set(xmin=low, xmax=high, 
          xlab="sub/site",
          main="%s absolute branch lengths (a=%f, b=%f, Pval=%e)" % 
          (str(name), prms[0], prms[1], pval))
    
    return p, util.Bundle(pval=pval, prms=prms)


def plotRelLens(name, params, lens, low, high, step):
    p = util.Gnuplot()
    p.enableOutput(False)
    
    lens = filter(util.withinfunc(low, high), lens)
    lens2 = filter(util.withinfunc(.00001, high), lens)
    
    prms, resid = stats.fitDistrib(stats.normalPdf, params, 
                                   lens2, low, high, step)    
    util.plotdistrib(lens, low=low, width=step, plot=p)
    p.plotfunc(lambda x: stats.normalPdf(x, prms), low, high, step/4.)
    
    bins, obs = util.hist(lens, low=low, width=step)
    tot = len(lens)
    ext = map(lambda x: tot * (stats.normalCdf(x+step, prms) - 
                               stats.normalCdf(x, prms)), bins)
    
    ind = util.find(lambda x: x>5, ext)
    ind = ind[1:]
    obs = util.mget(obs, ind)
    bins = util.mget(bins, ind)
    ext = util.mget(ext, ind)
    
    
    try:
        chisq = rpy.r.chisq_test(obs, p=util.oneNorm(ext))
        pval = chisq["p.value"]
    except rpy.RException, e:
        print e
        pval = 0.0
    
    #pval = rpy.r.shapiro_test(lens2)["p.value"]
    
    util.logger("%s\t%.3e\t%s" % (name, pval, str(pval > .05)))<commit_msg>*** empty log message ***<commit_after>#!/usr/bin/env python

import sys, os
from rasmus import util, env, treelib, spidirlib, stats, tablelib

import rpy


options = [
  ["s:", "stree=", "stree", "<species tree>",
    {"single": True}],
  ["p:", "params=", "params", "<sindir params file>",
    {"single": True}],
  ["l:", "lens=", "lens", "<sindir branch lengths file>",
    {"single": True}],
  ["o:", "outdir=", "outdir", "<output directory>",
    {"single": True,
     "default": "."}],
  ["b:", "nbins=", "nbins", "<number of bins>",
    {"single": True,
     "parser": int,
     "default": 20}]
]


# parse options
conf = util.parseOptions(sys.argv, options, quit=True)


# read data
env.addEnvPaths("DATAPATH")
stree = treelib.readTree(env.findFile(conf["stree"]))
params = spidirlib.readParams(env.findFile(conf["params"]))
lens = spidirlib.readTreeDistrib(env.findFile(conf["lens"]))
totals = map(sum, zip(* lens.values()))
rlens = util.mapdict(lens, valfunc=lambda x: util.vidiv(x, totals))

fitting = tablelib.Table(headers=["name", "type", "param1", "param2", "pval"])


def plotAbsLens(name, lens, low, high, step):
    p = util.Gnuplot()
    p.enableOutput(False)
    
    lens = filter(util.withinfunc(low, high), lens)
    
    lens2 = filter(util.gtfunc(.001), lens)
    mu = stats.mean(lens2)
    sigma2 = stats.variance(lens2)
    prms = [mu*mu/sigma2, mu/sigma2]
    
    prms, resid = stats.fitDistrib(stats.gammaPdf, prms, 
                                   lens2, low, high, step)
    
    util.plotdistrib(lens, low=low, width=step, plot=p)
    p.plotfunc(lambda x: stats.gammaPdf(x, prms), low, high, step/4.)
    
    
    bins, obs = util.hist(lens, low=low, width=step)
    tot = len(lens)
    ext = map(lambda x: tot * (stats.gammaCdf(x+step, prms) - 
                               stats.gammaCdf(x, prms)), bins)
    
    ind = util.find(lambda x: x>5, ext)
    ind = ind[1:]
    obs = util.mget(obs, ind)
    bins = util.mget(bins, ind)
    ext = util.mget(ext, ind)
    
    try:
        chisq = rpy.r.chisq_test(obs, p=util.oneNorm(ext))
        pval = float(chisq["p.value"])
    except rpy.RException, e:
        print e
        pval = 0.0
    
    
    util.logger("%s\t%.3e\t%s" % (name, pval, str(pval > .05)))
    
    p.set(xmin=low, xmax=high, 
          xlab="sub/site",
          main="%s absolute branch lengths (a=%f, b=%f, Pval=%e)" % 
          (str(name), prms[0], prms[1], pval))
    
    return p, util.Bundle(pval=pval, prms=prms)


def plotRelLens(name, params, lens, low, high, step):
    p = util.Gnuplot()
    p.enableOutput(False)
    
    lens = filter(util.withinfunc(low, high), lens)
    lens2 = filter(util.withinfunc(.00001, high), lens)
    
    prms, resid = stats.fitDistrib(stats.normalPdf, params, 
                                   lens2, low, high, step)    
    util.plotdistrib(lens, low=low, width=step, plot=p)
    p.plotfunc(lambda x: stats.normalPdf(x, prms), low, high, step/4.)
    
    bins, obs = util.hist(lens, low=low, width=step)
    tot = len(lens)
    ext = map(lambda x: tot * (stats.normalCdf(x+step, prms) - 
                               stats.normalCdf(x, prms)), bins)
    
    ind = util.find(lambda x: x>5, ext)
    ind = ind[1:]
    obs = util.mget(obs, ind)
    bins = util.mget(bins, ind)
    ext = util.mget(ext, ind)
    
    
    try:
        chisq = rpy.r.chisq_test(obs, p=util.oneNorm(ext))
        pval = float(chisq["p.value"])
    except rpy.RException, e:
        print e
        pval = 0.0
    
    #pval = rpy.r.shapiro_test(lens2)["p.value"]
    
    util.logger("%s\t%.3e\t%s" % (name, pval, str(pval > .05)))
{'fraction_non_alphanumeric': 0.17751315425090003, 'fraction_numerical': 0.011077263915812794, 'mean_word_length': 2.613306653326663, 'pattern_counts': {'":': 16, '<': 13, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5410877', 'n_tokens_mistral': 2863, 'n_tokens_neox': 2742, 'n_words': 585}